Ifpack Package Browser (Single Doxygen Collection) Development
icrout_cholesky_mex.c
Go to the documentation of this file.
00001 /*BHEADER**********************************************************************
00002  * (c) 2002   The Regents of the University of California
00003  *
00004  * See the file COPYRIGHT_and_DISCLAIMER for a complete copyright
00005  * notice, contact person, and disclaimer.
00006  *
00007  * $Revision$
00008 *********************************************************************EHEADER*/
00009 
00010 /*
00011  * Crout-form incomplete factorization
00012  *
00013  * To create a Matlab interface for icrout_cholesky_mex, define the symbol
00014  * MATLAB and compile this file as a Matlab mex program.
00015  * 
00016  * Contact Info for this code
00017  * --------------------------
00018  * Edmond Chow, Lawrence Livermore National Laboratory, echow@llnl.gov
00019  *
00020  * Revision History
00021  * ----------------
00022  * 06/04/02 Cleaned up for Boeing
00023  * 11/06/02 Added icrout_cholesky_mex
00024  */
00025 
00026 #undef IFPACK
00027 
00028 #define SYMSTR 1 /* structurally symmetric version */
00029 
00030 #ifdef MATLAB
00031 
00032 #include "matrix.h"
00033 #include "mex.h"
00034 #define malloc mxMalloc
00035 #define calloc mxCalloc
00036 #define realloc mxRealloc
00037 #define free mxFree
00038 #define printf mexPrintf
00039 
00040 #else
00041 
00042 #include <stdio.h>
00043 
00044 #endif
00045 
00046 #include "Ifpack_config.h"
00047 
00048 #include <stdlib.h>
00049 #include <math.h>
00050 #include <assert.h>
00051 
00052 #define MIN(a,b) ((a)<=(b) ? (a) : (b))
00053 #define MAX(a,b) ((a)>=(b) ? (a) : (b))
00054 #define ABS(a) ((a)>=0 ? (a) : -(a))
00055 
00056 /*
00057  * Data structure for sparse matrices is CSR, 0-based indexing.
00058  */
00059 typedef struct {
00060     double *val;  /* also known as A  */
00061     int    *col;  /* also known as JA; first column is column 0 */
00062     int    *ptr;  /* also known as IA; with ptr[0] = 0 */
00063 } Matrix;
00064 
00065 IFPACK_DEPRECATED void ifpack_quicksort (int *const pbase, double *const daux, size_t total_elems);
00066 
00067 #define SHORTCUT(p, a, ja, ia) \
00068         (a)  = (p)->val; \
00069         (ja) = (p)->col; \
00070         (ia) = (p)->ptr;
00071 
00072 #define MATNULL(p) \
00073         (p).val = NULL; \
00074         (p).col = NULL; \
00075         (p).ptr = NULL;
00076 
00077 IFPACK_DEPRECATED void Matrix_alloc(Matrix *a, int n, int nnz)
00078 {
00079     a->val = (double *) malloc(nnz * sizeof(double));
00080     a->col = (int *) malloc(nnz * sizeof(int));
00081     a->ptr = (int *) malloc((n+1) * sizeof(int));
00082 }
00083  
00084 IFPACK_DEPRECATED void Matrix_dealloc(Matrix *a)
00085 {
00086     free(a->val);
00087     free(a->col);
00088     free(a->ptr);
00089     a->val = NULL;
00090     a->col = NULL;
00091     a->ptr = NULL;
00092 }
00093 
00094 static void qsplit(double *a, int *ind, int n, int ncut)
00095 {
00096     double tmp, abskey;
00097     int itmp, first, last, mid;
00098     int j;
00099  
00100     ncut--;
00101     first = 0;
00102     last = n-1;
00103     if (ncut < first || ncut > last) 
00104         return;
00105  
00106     /* outer loop while mid != ncut */
00107     while (1)
00108     {
00109         mid = first;
00110         abskey = ABS(a[mid]);
00111         for (j=first+1; j<=last; j++)
00112         {
00113             if (ABS(a[j]) > abskey)
00114             {
00115                 mid = mid+1;
00116                 /* interchange */
00117                 tmp = a[mid];
00118                 itmp = ind[mid];
00119                 a[mid] = a[j];
00120                 ind[mid] = ind[j];
00121                 a[j]  = tmp;
00122                 ind[j] = itmp;
00123             }
00124         }
00125        
00126         /* interchange */
00127         tmp = a[mid];
00128         a[mid] = a[first];
00129         a[first]  = tmp;
00130         itmp = ind[mid];
00131         ind[mid] = ind[first];
00132         ind[first] = itmp;
00133        
00134         /* test for while loop */
00135         if (mid == ncut) 
00136             return;
00137 
00138         if (mid > ncut)
00139             last = mid-1;
00140         else
00141             first = mid+1;
00142     }
00143 }
00144 
00145 /* update column k using previous columns */
00146 /* assumes that column of A resides in the work vector */
00147 /* this function can also be used to update rows */
00148  
00149 static void update_column(
00150     int k,
00151     const int *ia, const int *ja, const double *a,
00152     const int *ifirst,
00153     const int *ifirst2,
00154     const int *list2,
00155     const double *multipliers, /* the val array of the other factor */
00156     const double *d, /* diagonal of factorization */
00157     int *marker,
00158     double *ta,
00159     int *itcol,
00160     int *ptalen)
00161 {
00162     int j, i, isj, iej, row;
00163     int talen, pos;
00164     double mult;
00165  
00166     talen = *ptalen;
00167 
00168     j = list2[k];
00169     while (j != -1)
00170     {
00171         /* update column k using column j */
00172  
00173         isj = ifirst[j];
00174 
00175         /* skip first nonzero in column, since it does not contribute */
00176         /* if symmetric structure */
00177         /* isj++; */
00178  
00179         /* now do the actual update */
00180        if (isj != -1)
00181        {
00182         /* multiplier */
00183         mult = multipliers[ifirst2[j]] * d[j];
00184 
00185         /* section of column used for update */
00186         iej = ia[j+1]-1;
00187 
00188         for (i=isj; i<=iej; i++)
00189         {
00190             row = ja[i];
00191 
00192             /* if nonsymmetric structure */
00193             if (row <= k)
00194                 continue;
00195  
00196             if ((pos = marker[row]) != -1)
00197             {
00198                 /* already in pattern of working row */
00199                 ta[pos] -= mult*a[i];
00200             }
00201             else
00202             {
00203                 /* not yet in pattern of working row */
00204                 itcol[talen] = row;
00205                 ta[talen] = -mult*a[i];
00206                 marker[row] = talen++;
00207             }
00208         }
00209        }
00210  
00211         j = list2[j];
00212     }
00213 
00214     *ptalen = talen;
00215 }
00216 
00217 /* update ifirst and list */
00218  
00219 static void update_lists(
00220     int k,
00221     const int *ia,
00222     const int *ja,
00223     int *ifirst,
00224     int *list)
00225 {
00226     int j, isj, iej, iptr;
00227 
00228     j = list[k];
00229     while (j != -1)
00230     {
00231         isj = ifirst[j];
00232         iej = ia[j+1]-1;
00233  
00234         isj++;
00235  
00236         if (isj != 0 && isj <= iej) /* nonsymmetric structure */
00237         {
00238             /* increment ifirst for column j */
00239             ifirst[j] = isj;
00240 
00241             /* add j to head of list for list[ja[isj]] */
00242             iptr = j;
00243             j = list[j];
00244             list[iptr] = list[ja[isj]];
00245             list[ja[isj]] = iptr;
00246         }
00247         else
00248         {
00249             j = list[j];
00250         }
00251     }
00252 }
00253 
00254 static void update_lists_newcol(
00255     int k,
00256     int isk,
00257     int iptr,
00258     int *ifirst,
00259     int *list)
00260 {
00261     ifirst[k] = isk;
00262     list[k] = list[iptr];
00263     list[iptr] = k;
00264 }
00265 
00266 /*
00267  * crout_ict - Crout version of ICT - Incomplete Cholesky by Threshold
00268  *
00269  * The incomplete factorization L D L^T is computed,
00270  * where L is unit triangular, and D is diagonal
00271  *
00272  * INPUTS
00273  * n = number of rows or columns of the matrices
00274  * AL = unit lower triangular part of A stored by columns
00275  *            the unit diagonal is implied (not stored)
00276  * Adiag = diagonal part of A
00277  * droptol = drop tolerance
00278  * lfil  = max nonzeros per col in L factor or per row in U factor
00279  *
00280  * OUTPUTS
00281  * L     = lower triangular factor stored by columns
00282  * pdiag = diagonal factor stored in an array
00283  *
00284  * NOTE: calling function must free the memory allocated by this
00285  * function, i.e., L, pdiag.
00286  */
00287 
00288 IFPACK_DEPRECATED void crout_ict(
00289     int n,
00290 #ifdef IFPACK
00291     void * A,
00292     int maxentries;
00293     int (*getcol)( void * A, int col, int ** nentries, double * val, int * ind),
00294     int (*getdiag)( void *A, double * diag),
00295 #else
00296     const Matrix *AL,
00297     const double *Adiag,
00298 #endif
00299     double droptol,
00300     int lfil,
00301     Matrix *L,
00302     double **pdiag)
00303 {
00304     int k, j, i, index;
00305     int count_l;
00306     double norm_l;
00307 
00308     /* work arrays; work_l is list of nonzero values */
00309     double *work_l = (double *) malloc(n * sizeof(double));
00310     int    *ind_l  = (int *)    malloc(n * sizeof(int));
00311     int     len_l;
00312 
00313     /* list and ifirst data structures */
00314     int *list_l    = (int *) malloc(n * sizeof(int));
00315     int *ifirst_l  = (int *) malloc(n * sizeof(int));
00316     
00317     /* aliases */
00318     int *marker_l  = ifirst_l;
00319 
00320     /* matrix arrays */
00321     double *al; int *jal, *ial;
00322     double *l;  int *jl,  *il;
00323 
00324     int kl = 0;
00325 
00326     double *diag = (double *) malloc(n * sizeof(double));
00327     *pdiag = diag;
00328 
00329     Matrix_alloc(L,  n, lfil*n);
00330 
00331     SHORTCUT(AL, al, jal, ial);
00332     SHORTCUT(L,  l,  jl,  il);
00333 
00334     /* initialize work arrays */
00335     for (i=0; i<n; i++)
00336     {
00337         list_l[i]    = -1;
00338         ifirst_l[i]  = -1;
00339         marker_l[i]  = -1;
00340     }
00341 
00342     /* copy the diagonal */
00343 #ifdef IFPACK
00344     getdiag(A, diag);
00345 #else
00346     for (i=0; i<n; i++)
00347         diag[i] = Adiag[i];
00348 #endif
00349 
00350     /* start off the matrices right */
00351     il[0]  = 0;
00352 
00353     /*
00354      * Main loop over columns and rows
00355      */
00356 
00357     for (k=0; k<n; k++)
00358     {
00359         /*
00360          * lower triangular factor update by columns
00361          * (need ifirst for L and lists for U)
00362          */
00363  
00364         /* copy column of A into work vector */
00365         norm_l = 0.;
00366 #ifdef IFPACK
00367       getcol(A, k, len_l, work_l, ind_l);
00368       for (j=0; j<len_l; j++)
00369   {
00370     norm_l += ABS(work_l[j]);
00371     marker_l[ind_l[j]] = j;
00372   }
00373 #else
00374         len_l = 0;
00375         for (j=ial[k]; j<ial[k+1]; j++)
00376         {
00377             index = jal[j];
00378             work_l[len_l] = al[j];
00379             norm_l += ABS(al[j]);
00380             ind_l[len_l] = index;
00381             marker_l[index] = len_l++; /* points into work array */
00382         }
00383 #endif
00384         norm_l = (len_l == 0) ? 0.0 : norm_l/((double) len_l);
00385  
00386         /* update and scale */
00387 
00388         update_column(k, il, jl, l, ifirst_l, ifirst_l, list_l, l,
00389             diag, marker_l, work_l, ind_l, &len_l);
00390 
00391         for (j=0; j<len_l; j++)
00392             work_l[j] /= diag[k];
00393 
00394   i = 0;
00395         for (j=0; j<len_l; j++)
00396   {
00397       if (ABS(work_l[j]) < droptol * norm_l)
00398       {
00399           /* zero out marker array for this */
00400     marker_l[ind_l[j]] = -1;
00401       }
00402       else
00403       {
00404           work_l[i] = work_l[j];
00405     ind_l[i]  = ind_l[j];
00406     i++;
00407       }
00408   }
00409   len_l = i;
00410 
00411   /*
00412    * dropping:  for each work vector, perform qsplit, and then
00413    * sort each part by increasing index; then copy each sorted
00414    * part into the factors
00415    */
00416 
00417   count_l = MIN(len_l, lfil);
00418   qsplit(work_l, ind_l, len_l, count_l);
00419   ifpack_quicksort(ind_l, work_l, count_l);
00420 
00421   for (j=0; j<count_l; j++)
00422   {
00423       l[kl] = work_l[j];
00424       jl[kl++] = ind_l[j];
00425   }
00426   il[k+1] = kl;
00427 
00428   /*
00429    * update lists
00430    */
00431   
00432         update_lists(k, il, jl, ifirst_l, list_l);
00433 
00434   if (kl - il[k] > SYMSTR)
00435       update_lists_newcol(k, il[k], jl[il[k]], ifirst_l, list_l);
00436 
00437         /* zero out the marker arrays */
00438         for (j=0; j<len_l; j++)
00439             marker_l[ind_l[j]] = -1;
00440 
00441         /*
00442          * update the diagonal (after dropping)
00443          */
00444 
00445         for (j=0; j<count_l; j++)
00446         {
00447             index = ind_l[j];
00448             diag[index] -= work_l[j] * work_l[j] * diag[k];
00449         }
00450     }
00451 
00452     free(work_l);
00453     free(ind_l);
00454     free(list_l);
00455     free(ifirst_l);
00456 }
00457 
00458 #ifdef MATLAB
00459 /*
00460  * [l, d] = icrout_cholesky_mex(AL, Adiag, droptol, lfil)
00461  * 
00462  * where AL is a sparse matrix (strictly lower triangular part)
00463  * Adiag is a vector
00464  * droptol is a scalar real
00465  * lfil is a scalar integer
00466  *
00467  * d must be converted to a matrix (from a vector) by icrout_cholesky.m
00468  */
00469 IFPACK_DEPRECATED void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
00470 {
00471     int n;
00472     Matrix AL;
00473     double *Adiag;
00474     double droptol;
00475     int lfil;
00476     Matrix L;
00477     double *D;
00478 
00479     if (nrhs != 4)
00480         mexErrMsgTxt("mex function called with bad number of arguments");
00481 
00482     n = mxGetN(prhs[0]);
00483     AL.ptr = mxGetJc(prhs[0]);
00484     AL.col = mxGetIr(prhs[0]);
00485     AL.val = mxGetPr(prhs[0]);
00486     Adiag  = mxGetPr(prhs[1]);
00487 
00488     droptol = (double) *mxGetPr(prhs[2]);
00489     lfil   = (int) *mxGetPr(prhs[3]);
00490 
00491     crout_ict(n, &AL, Adiag, droptol, lfil, &L, &D);
00492 
00493     /* create output matrices */
00494 
00495     /* L */
00496     plhs[0] = mxCreateSparse(n, n, 1, mxREAL);
00497     mxFree(mxGetJc(plhs[0]));
00498     mxFree(mxGetIr(plhs[0]));
00499     mxFree(mxGetPr(plhs[0]));
00500     mxSetJc(plhs[0], L.ptr);
00501     mxSetIr(plhs[0], L.col);
00502     mxSetPr(plhs[0], L.val);
00503     mxSetNzmax(plhs[0], n*lfil); /* must agree */
00504 
00505     /* D */
00506     plhs[1] = mxCreateDoubleMatrix(n, 1, mxREAL);
00507     mxFree(mxGetPr(plhs[1]));
00508     mxSetPr(plhs[1], D);
00509 }
00510 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines