Ifpack_IC_Utils.cpp

Go to the documentation of this file.
00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_IC_Utils.h"
00003 
00004 #define SYMSTR 1 /* structurally symmetric version */
00005 #include <stdio.h>
00006 
00007 #define MIN(a,b) ((a)<=(b) ? (a) : (b))
00008 #define MAX(a,b) ((a)>=(b) ? (a) : (b))
00009 #define ABS(a) ((a)>=0 ? (a) : -(a))
00010 
00011 #define SHORTCUT(p, a, ja, ia) \
00012         (a)  = (p)->val; \
00013         (ja) = (p)->col; \
00014         (ia) = (p)->ptr;
00015 
00016 #define MATNULL(p) \
00017         (p).val = NULL; \
00018         (p).col = NULL; \
00019         (p).ptr = NULL;
00020 
00021 void Ifpack_AIJMatrix_alloc(Ifpack_AIJMatrix *a, int n, int nnz)
00022 {
00023     a->val = new double[nnz];
00024     a->col = new int[nnz];
00025     a->ptr = new int[n+1];
00026 }
00027  
00028 void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
00029 {
00030     delete [] (a->val);
00031     delete [] (a->col);
00032     delete [] (a->ptr);
00033     a->val = 0;
00034     a->col = 0;
00035     a->ptr = 0;
00036 }
00037 
00038 static void qsplit(double *a, int *ind, int n, int ncut)
00039 {
00040     double tmp, abskey;
00041     int itmp, first, last, mid;
00042     int j;
00043  
00044     ncut--;
00045     first = 0;
00046     last = n-1;
00047     if (ncut < first || ncut > last) 
00048         return;
00049  
00050     /* outer loop while mid != ncut */
00051     while (1)
00052     {
00053         mid = first;
00054         abskey = ABS(a[mid]);
00055         for (j=first+1; j<=last; j++)
00056         {
00057             if (ABS(a[j]) > abskey)
00058             {
00059                 mid = mid+1;
00060                 /* interchange */
00061                 tmp = a[mid];
00062                 itmp = ind[mid];
00063                 a[mid] = a[j];
00064                 ind[mid] = ind[j];
00065                 a[j]  = tmp;
00066                 ind[j] = itmp;
00067             }
00068         }
00069        
00070         /* interchange */
00071         tmp = a[mid];
00072         a[mid] = a[first];
00073         a[first]  = tmp;
00074         itmp = ind[mid];
00075         ind[mid] = ind[first];
00076         ind[first] = itmp;
00077        
00078         /* test for while loop */
00079         if (mid == ncut) 
00080             return;
00081 
00082         if (mid > ncut)
00083             last = mid-1;
00084         else
00085             first = mid+1;
00086     }
00087 }
00088 
00089 /* update column k using previous columns */
00090 /* assumes that column of A resides in the work vector */
00091 /* this function can also be used to update rows */
00092  
00093 static void update_column(
00094     int k,
00095     const int *ia, const int *ja, const double *a,
00096     const int *ifirst,
00097     const int *ifirst2,
00098     const int *list2,
00099     const double *multipliers, /* the val array of the other factor */
00100     const double *d, /* diagonal of factorization */
00101     int *marker,
00102     double *ta,
00103     int *itcol,
00104     int *ptalen)
00105 {
00106     int j, i, isj, iej, row;
00107     int talen, pos;
00108     double mult;
00109  
00110     talen = *ptalen;
00111 
00112     j = list2[k];
00113     while (j != -1)
00114     {
00115         /* update column k using column j */
00116  
00117         isj = ifirst[j];
00118 
00119         /* skip first nonzero in column, since it does not contribute */
00120         /* if symmetric structure */
00121         /* isj++; */
00122  
00123         /* now do the actual update */
00124        if (isj != -1)
00125        {
00126         /* multiplier */
00127         mult = multipliers[ifirst2[j]] * d[j];
00128 
00129         /* section of column used for update */
00130         iej = ia[j+1]-1;
00131 
00132         for (i=isj; i<=iej; i++)
00133         {
00134             row = ja[i];
00135 
00136             /* if nonsymmetric structure */
00137             if (row <= k)
00138                 continue;
00139  
00140             if ((pos = marker[row]) != -1)
00141             {
00142                 /* already in pattern of working row */
00143                 ta[pos] -= mult*a[i];
00144             }
00145             else
00146             {
00147                 /* not yet in pattern of working row */
00148                 itcol[talen] = row;
00149                 ta[talen] = -mult*a[i];
00150                 marker[row] = talen++;
00151             }
00152         }
00153        }
00154  
00155         j = list2[j];
00156     }
00157 
00158     *ptalen = talen;
00159 }
00160 
00161 /* update ifirst and list */
00162  
00163 static void update_lists(
00164     int k,
00165     const int *ia,
00166     const int *ja,
00167     int *ifirst,
00168     int *list)
00169 {
00170     int j, isj, iej, iptr;
00171 
00172     j = list[k];
00173     while (j != -1)
00174     {
00175         isj = ifirst[j];
00176         iej = ia[j+1]-1;
00177  
00178         isj++;
00179  
00180         if (isj != 0 && isj <= iej) /* nonsymmetric structure */
00181         {
00182             /* increment ifirst for column j */
00183             ifirst[j] = isj;
00184 
00185             /* add j to head of list for list[ja[isj]] */
00186             iptr = j;
00187             j = list[j];
00188             list[iptr] = list[ja[isj]];
00189             list[ja[isj]] = iptr;
00190         }
00191         else
00192         {
00193             j = list[j];
00194         }
00195     }
00196 }
00197 
00198 static void update_lists_newcol(
00199     int k,
00200     int isk,
00201     int iptr,
00202     int *ifirst,
00203     int *list)
00204 {
00205     ifirst[k] = isk;
00206     list[k] = list[iptr];
00207     list[iptr] = k;
00208 }
00209 
00210 /*
00211  * crout_ict - Crout version of ICT - Incomplete Cholesky by Threshold
00212  *
00213  * The incomplete factorization L D L^T is computed,
00214  * where L is unit triangular, and D is diagonal
00215  *
00216  * INPUTS
00217  * n = number of rows or columns of the matrices
00218  * AL = unit lower triangular part of A stored by columns
00219  *            the unit diagonal is implied (not stored)
00220  * Adiag = diagonal part of A
00221  * droptol = drop tolerance
00222  * lfil  = max nonzeros per col in L factor or per row in U factor
00223  *
00224  * OUTPUTS
00225  * L     = lower triangular factor stored by columns
00226  * pdiag = diagonal factor stored in an array
00227  *
00228  * NOTE: calling function must free the memory allocated by this
00229  * function, i.e., L, pdiag.
00230  */
00231 
00232 void crout_ict(
00233     int n,
00234 #ifdef IFPACK
00235     void * A,
00236     int maxentries;
00237     int (*getcol)( void * A, int col, int ** nentries, double * val, int * ind),
00238     int (*getdiag)( void *A, double * diag),
00239 #else
00240     const Ifpack_AIJMatrix *AL,
00241     const double *Adiag,
00242 #endif
00243     double droptol,
00244     int lfil,
00245     Ifpack_AIJMatrix *L,
00246     double **pdiag)
00247 {
00248     int k, j, i, index;
00249     int count_l;
00250     double norm_l;
00251 
00252     /* work arrays; work_l is list of nonzero values */
00253     double *work_l = new double[n];
00254     int    *ind_l  = new int[n];
00255     int     len_l;
00256 
00257     /* list and ifirst data structures */
00258     int *list_l    = new int[n];
00259     int *ifirst_l  = new int[n];
00260     
00261     /* aliases */
00262     int *marker_l  = ifirst_l;
00263 
00264     /* matrix arrays */
00265     double *al; int *jal, *ial;
00266     double *l;  int *jl,  *il;
00267 
00268     int kl = 0;
00269 
00270     double *diag = new double[n];
00271     *pdiag = diag;
00272 
00273     Ifpack_AIJMatrix_alloc(L,  n, lfil*n);
00274 
00275     SHORTCUT(AL, al, jal, ial);
00276     SHORTCUT(L,  l,  jl,  il);
00277 
00278     /* initialize work arrays */
00279     for (i=0; i<n; i++)
00280     {
00281         list_l[i]    = -1;
00282         ifirst_l[i]  = -1;
00283         marker_l[i]  = -1;
00284     }
00285 
00286     /* copy the diagonal */
00287 #ifdef IFPACK
00288     getdiag(A, diag);
00289 #else
00290     for (i=0; i<n; i++)
00291         diag[i] = Adiag[i];
00292 #endif
00293 
00294     /* start off the matrices right */
00295     il[0]  = 0;
00296 
00297     /*
00298      * Main loop over columns and rows
00299      */
00300 
00301     for (k=0; k<n; k++)
00302     {
00303         /*
00304          * lower triangular factor update by columns
00305          * (need ifirst for L and lists for U)
00306          */
00307  
00308         /* copy column of A into work vector */
00309         norm_l = 0.;
00310 #ifdef IFPACK
00311       getcol(A, k, len_l, work_l, ind_l);
00312       for (j=0; j<len_l; j++)
00313   {
00314     norm_l += ABS(work_l[j]);
00315     marker_l[ind_l[j]] = j;
00316   }
00317 #else
00318         len_l = 0;
00319         for (j=ial[k]; j<ial[k+1]; j++)
00320         {
00321             index = jal[j];
00322             work_l[len_l] = al[j];
00323             norm_l += ABS(al[j]);
00324             ind_l[len_l] = index;
00325             marker_l[index] = len_l++; /* points into work array */
00326         }
00327 #endif
00328         norm_l = (len_l == 0) ? 0.0 : norm_l/((double) len_l);
00329  
00330         /* update and scale */
00331 
00332         update_column(k, il, jl, l, ifirst_l, ifirst_l, list_l, l,
00333             diag, marker_l, work_l, ind_l, &len_l);
00334 
00335         for (j=0; j<len_l; j++)
00336             work_l[j] /= diag[k];
00337 
00338   i = 0;
00339         for (j=0; j<len_l; j++)
00340   {
00341       if (ABS(work_l[j]) < droptol * norm_l)
00342       {
00343           /* zero out marker array for this */
00344     marker_l[ind_l[j]] = -1;
00345       }
00346       else
00347       {
00348           work_l[i] = work_l[j];
00349     ind_l[i]  = ind_l[j];
00350     i++;
00351       }
00352   }
00353   len_l = i;
00354 
00355   /*
00356    * dropping:  for each work vector, perform qsplit, and then
00357    * sort each part by increasing index; then copy each sorted
00358    * part into the factors
00359    */
00360 
00361   count_l = MIN(len_l, lfil);
00362   qsplit(work_l, ind_l, len_l, count_l);
00363   quicksort(ind_l, work_l, count_l);
00364 
00365   for (j=0; j<count_l; j++)
00366   {
00367       l[kl] = work_l[j];
00368       jl[kl++] = ind_l[j];
00369   }
00370   il[k+1] = kl;
00371 
00372   /*
00373    * update lists
00374    */
00375   
00376         update_lists(k, il, jl, ifirst_l, list_l);
00377 
00378   if (kl - il[k] > SYMSTR)
00379       update_lists_newcol(k, il[k], jl[il[k]], ifirst_l, list_l);
00380 
00381         /* zero out the marker arrays */
00382         for (j=0; j<len_l; j++)
00383             marker_l[ind_l[j]] = -1;
00384 
00385         /*
00386          * update the diagonal (after dropping)
00387          */
00388 
00389         for (j=0; j<count_l; j++)
00390         {
00391             index = ind_l[j];
00392             diag[index] -= work_l[j] * work_l[j] * diag[k];
00393         }
00394     }
00395 
00396     delete [] work_l;
00397     delete [] ind_l;
00398     delete [] list_l;
00399     delete [] ifirst_l;
00400 }

Generated on Thu Sep 18 12:37:21 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1