ifp_SparseUtil.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 #include "ifp_ifpack.h"
00030 #include "ifp_SparseUtil.h"
00031 #include "stdio.h" // kludge
00032 using namespace std;
00033 /* #define DEBUG */
00034 // shell sort
00035 // stable, so it is fast if already sorted
00036 // sorts x[0:n-1] in place, ascending order.
00037 
00038 void shell_sort(
00039   const int n,
00040   int x[])
00041 {
00042     int m, max, j, k, itemp;
00043     
00044     m = n/2;
00045 
00046     while (m > 0) {
00047         max = n - m;
00048         for (j=0; j<max; j++)
00049         {
00050             for (k=j; k>=0; k-=m)
00051             {
00052                 if (x[k+m] >= x[k])
00053                     break;
00054                 itemp = x[k+m];
00055                 x[k+m] = x[k];
00056                 x[k] = itemp;
00057             }
00058         }   
00059         m = m/2;
00060     }
00061 }
00062 
00063 // allocate space for integer data for level ILU
00064 
00065 void allocate_ilu(
00066   const int levfill,                 // level of fill
00067   const int n,                       // order of matrix
00068   int *nzl, int *nzu,                // space allocated
00069   const int ia[], const int ja[],    // input
00070   int *ial[], int *jal[],            // output lower factor structure
00071   int *iau[], int *jau[],            // output upper factor structure
00072   int growthl, int growthu)              // storage parameters
00073 {
00074     int i, j, nzla, nzua;
00075     int maxl, maxu;
00076 
00077     maxu = n*(n+1)/2;
00078     maxl = n*n - maxu;
00079 
00080     // count number of entries in lower and upper triangular parts
00081 
00082     nzla = 0;
00083     nzua = 0;
00084     for (i=0; i<n; i++)
00085     {
00086         for (j=ia[i]; j<ia[i+1]; j++)
00087         {
00088             if (ja[j] < i)
00089                nzla++;
00090             else
00091                nzua++;
00092         }
00093     }
00094 
00095     *ial = new int[n+1];
00096     *iau = new int[n+1];
00097 
00098     if (levfill == 0) // ILU(0)
00099     {
00100         *nzl = nzla;
00101         *nzu = nzua + n;
00102     }
00103     else if (levfill == n) // full factorization
00104     {
00105         *nzl = maxl;
00106         *nzu = maxu;
00107     }
00108     else
00109     {
00110         *nzl = MIN((levfill+growthl)*nzla, maxl);
00111         *nzu = MIN((levfill+growthu)*nzua + n, maxu);
00112     }
00113 
00114     *jal = new int[*nzl];
00115     *jau = new int[*nzu];
00116 
00117     if (levfill != 0)
00118     {
00119       // cerr << "nnz in unfactored mat: " << nzla + nzua << endl;
00120       // cerr << "nnz allocated for ILU: " << *nzl + *nzu << endl;
00121     }
00122 }
00123 
00124 // symbolic level ILU
00125 // fortran-style data structures and algorithms
00126 // factors into separate upper and lower parts
00127 // assumes rows are sorted  *** symbolic ILU
00128 // assumes no zero rows
00129 // Returns an integer:
00130 //  0 - No error
00131 // -1 - Not enough memory for L
00132 // -2 - Not enough memory for U
00133 
00134 int symbolic_ilu(
00135   const int levfill,                 // level of fill
00136   const int n,                       // order of matrix
00137   int *nzl,                          // input-output
00138   int *nzu,                          // input-output
00139   const int ia[], const int ja[],    // input
00140   int ial[], int jal[],              // output lower factor structure
00141   int iau[], int jau[])              // output upper factor structure
00142 {
00143 #if 0
00144     // copy if levfill is 0
00145     if (levfill == 0)
00146     {
00147         int kl = 0;
00148         int ku = 0;
00149         ial[0] = 0;
00150         iau[0] = 0;
00151         for (int i=0; i<n; i++)
00152         {
00153             for (int j=ia[i]; j<ia[i+1]; j++)
00154             {
00155                 if (ja[j] < i)
00156                    jal[kl++] = ja[j];
00157                 else
00158                    jau[ku++] = ja[j];
00159             }
00160             ial[i+1] = kl;
00161             iau[i+1] = ku;
00162             shell_sort(ial[i+1]-ial[i], &jal[ial[i]]);
00163             shell_sort(iau[i+1]-iau[i], &jau[iau[i]]);
00164         }
00165         assert (kl <= *nzl); 
00166         assert (ku <= *nzu);
00167         *nzl = kl;
00168         *nzu = ku;
00169     }
00170     else
00171 #endif
00172     {
00173 
00174     int *lnklst = new int[n];
00175     int *curlev = new int[n];
00176     int *levels = new int[*nzu];
00177     int *iwork = new int[n];
00178 
00179     int knzl = 0;
00180     int knzu = 0;
00181 
00182     ial[0] = 0;
00183     iau[0] = 0;
00184 
00185     for (int i=0; i<n; i++)
00186     {
00187         int first, next, j;
00188 
00189         // copy column indices of row into workspace and sort them
00190 
00191         int len = ia[i+1] - ia[i];
00192     if (len == 0 ) cerr << " Row " << i << " is zero " << endl;
00193         next = 0;
00194         for (j=ia[i]; j<ia[i+1]; j++)
00195             iwork[next++] = ja[j];
00196         shell_sort(len, iwork);
00197 
00198         // construct implied linked list for row
00199 
00200         first = iwork[0];
00201         curlev[first] = 0;
00202 
00203         for (j=0; j<=len-2; j++)
00204         {
00205             lnklst[iwork[j]] = iwork[j+1];
00206             curlev[iwork[j]] = 0;
00207         }
00208 
00209         lnklst[iwork[len-1]] = n;
00210         curlev[iwork[len-1]] = 0;
00211 
00212         // merge with rows in U
00213 
00214         next = first;
00215         while (next < i)
00216         {
00217             int oldlst = next;
00218             int nxtlst = lnklst[next];
00219             int row = next;
00220             int ii;
00221 
00222             // scan row
00223 
00224             for (ii=iau[row]+1; ii<iau[row+1]; /*nop*/)
00225             {
00226                 if (jau[ii] < nxtlst)
00227                 {
00228                     // new fill-in
00229                     int newlev = curlev[row] + levels[ii] + 1;
00230                     if (newlev <= levfill)
00231                     {
00232                         lnklst[oldlst]  = jau[ii];
00233                         lnklst[jau[ii]] = nxtlst;
00234                         oldlst = jau[ii];
00235                         curlev[jau[ii]] = newlev;
00236                     }
00237                     ii++;
00238                 }
00239                 else if (jau[ii] == nxtlst)
00240                 {
00241                     oldlst = nxtlst;
00242                     nxtlst = lnklst[oldlst];
00243                     int newlev = curlev[row] + levels[ii] + 1;
00244                     curlev[jau[ii]] = MIN(curlev[jau[ii]], newlev);
00245                     ii++;
00246                 }
00247                 else // (jau[ii] > nxtlst)
00248                 {
00249                     oldlst = nxtlst;
00250                     nxtlst = lnklst[oldlst];
00251                 }
00252             }
00253             next = lnklst[next];
00254         }
00255         
00256         // gather the pattern into L and U
00257 
00258         next = first;
00259         while (next < i)
00260         {
00261            if (knzl >= *nzl)
00262          {
00263            cerr << "For row "<< i << endl;
00264            cerr << " nzl = "<< *nzl <<" knzl = "<< knzl << endl;
00265            printf("Not enough space allocated for lower symbolic factor\n");
00266            return(-1);
00267          }
00268 
00269             jal[knzl++] = next;
00270             next = lnklst[next];
00271         }
00272         ial[i+1] = knzl;
00273 
00274         if (next != i)
00275         {
00276             cerr << i << "  U has zero on diag, forcing nonzero" << endl;
00277             if (knzu >= *nzu)
00278          {
00279            cerr << " nzl = "<< *nzu <<" knzl = "<< knzu << endl;
00280            printf("Not enough space allocated for symbolic factor\n");
00281            return(-2);
00282          }
00283             levels[knzu] = 2*n; // infinity
00284             jau[knzu++] = i;
00285         }
00286 
00287         while (next < n)
00288         {
00289             if (knzu >= *nzu)
00290          {
00291            cerr << "For row "<< i << endl;
00292            cerr << " nzu = "<< *nzu <<" knzu = "<< knzu << endl;
00293            printf("Not enough space allocated for upper symbolic factor\n");
00294            return(-2);
00295     
00296          }
00297             levels[knzu] = curlev[next];
00298             jau[knzu++] = next;
00299             next = lnklst[next];
00300         }
00301         iau[i+1] = knzu;
00302     }
00303 
00304     delete [] lnklst;
00305     delete [] curlev;
00306     delete [] levels;
00307     delete [] iwork;
00308 
00309     *nzl = knzl;
00310     *nzu = knzu;
00311 
00312 #ifdef DEBUG
00313     cerr << "Actual nnz for ILU: " << *nzl + *nzu << endl;
00314 #endif
00315     }
00316     return(0);
00317 }

Generated on Wed May 12 21:30:17 2010 for IFPACK by  doxygen 1.4.7