Ifpack_IC_Utils.cpp

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

Generated on Tue Jul 13 09:27:13 2010 for IFPACK by  doxygen 1.4.7