ifp_biluk.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 "ifp_LocalMat.h"
00031 #include "ifp_BlockVec.h"
00032 #include "ifp_BlockMat.h"
00033 #include "ifp_biluk.h"
00034 #include "ifp_ifpack.h"
00035 #include "ifp_SparseUtil.h"
00036 #include "ifp_DenseMat.h"
00037 #include <iostream>
00038 using namespace std;
00039 /*int ifp_biluk::growth = 10; */
00040 
00041 ifp_biluk::ifp_biluk()
00042 {
00043     diag = (ifp_LocalMat **) NULL;
00044      al = NULL;
00045     jal = NULL;
00046     ial = NULL;
00047      au = NULL;
00048     jau = NULL;
00049     iau = NULL;
00050     NumEntries_ = 0;
00051     NumNonzeros_ = 0;
00052 }
00053 
00054 ifp_biluk::~ifp_biluk()
00055 {
00056     int i;
00057 
00058     if (diag != NULL)
00059         for (i=0; i<Ap->numrow(); i++)
00060             delete diag[i];
00061     delete [] diag;
00062 
00063     if (al != NULL)
00064         for (i=0; i<ial[Ap->numrow()]; i++)
00065             delete al[i];
00066     delete [] al;
00067 
00068     if (au != NULL)
00069         for (i=0; i<iau[Ap->numrow()]; i++)
00070             delete au[i];
00071     delete [] au;
00072 
00073     delete [] jal;
00074     delete [] ial;
00075     delete [] jau;
00076     delete [] iau;
00077 }
00078 
00079 // if levfill < 0, then use existing factorization for the pattern
00080 
00081 void ifp_biluk::setup(const ifp_BlockMat& A, int levfill)
00082 {
00083     Ap = &A;
00084 
00085     if (levfill >= 0)
00086     {
00087     int nzl, nzu, ierr;
00088 
00089         delete jal;
00090         delete ial;
00091         delete jau;
00092         delete iau;
00093     int growthl = 2;
00094     int growthu = 2;
00095     ierr = 1;
00096 
00097     while (ierr !=0)
00098       {
00099         allocate_ilu(levfill, Ap->numrow(), &nzl, &nzu, 
00100              &A.row_ptr(0), &A.col_ind(0), &ial, &jal, &iau, &jau, 
00101              growthl, growthu);
00102 
00103         ierr = symbolic_ilu(levfill, Ap->numrow(), &nzl, &nzu, 
00104                 &A.row_ptr(0), &A.col_ind(0), ial, jal, iau, jau);
00105         if (ierr !=0 )
00106           {
00107         cout << "Doubling storage and trying again\n" << endl;
00108         delete [] ial;
00109         delete [] jal;
00110         delete [] iau;
00111         delete [] jau;
00112         if (ierr == -1) // Not enough for L
00113           {
00114             growthl *=2;
00115             growthu *=2; // Assume current U size is too small
00116           }
00117         else // ierr = -2
00118           growthu *=2;
00119           }
00120       }
00121             
00122 
00123         // allocate for values
00124     delete [] al;
00125     delete [] au;
00126         al = new ifp_LocalMatp[nzl];
00127         au = new ifp_LocalMatp[nzu];
00128     }
00129 
00130     int i, j, k, kk, id, idd;
00131 
00132     int *marker = new int[Ap->numrow()];
00133     for (i=0; i<Ap->numrow(); i++)
00134         marker[i] = -1; // negative flag
00135 
00136     // full length work array of matrices
00137     ifp_LocalMat **row = new ifp_LocalMatp[A.numcol()];  // array of pointers to matrices
00138 
00139     diag = new ifp_LocalMatp[A.numrow()];
00140 
00141     NumEntries_ = ial[A.numrow()] - ial[0] + iau[A.numrow()] - iau[0];
00142     NumNonzeros_ = 0;
00143 
00144     for (i=0; i<A.numrow(); i++)
00145     {
00146     int neqr = A.kvst_row(i+1) - A.kvst_row(i);
00147 
00148         // scatter data structure of L and U
00149     j = 0;
00150         for (k=ial[i]; k<ial[i+1]; k++)
00151     {
00152             marker[jal[k]] = j;
00153         row[j] = A.val(0).CreateEmpty();
00154         const int nvars = A.kvst_col(jal[k]+1)-A.kvst_col(jal[k]);
00155         NumNonzeros_ += neqr*nvars;
00156         row[j++]->SetToZero(neqr, nvars);
00157     }
00158         for (k=iau[i]; k<iau[i+1]; k++)
00159     {
00160             marker[jau[k]] = j;
00161         row[j] = A.val(0).CreateEmpty();
00162         const int nvars = A.kvst_col(jau[k]+1)-A.kvst_col(jau[k]);
00163         NumNonzeros_ += neqr*nvars;
00164         row[j++]->SetToZero(neqr, nvars);
00165     }
00166         
00167     //  cout << "============== Block Row "<< i << "==================" << endl;
00168     //  cout << "Number of rows in block row = "<< neqr << endl;
00169         // scatter row of A
00170         for (k=A.row_ptr(i); k<A.row_ptr(i+1); k++)
00171       {
00172             row[marker[A.col_ind(k)]]->MatCopy(A.val(k));
00173         /*
00174         int nrow = ((ifp_DenseMat *)&(A.val(k)))->nrow;
00175         int ncol = ((ifp_DenseMat *)&(A.val(k)))->ncol;
00176           cout << "Block Col "<< A.col_ind(k) << ",  Number of cols = " <<ncol << endl;
00177         register double *qxx = ((ifp_DenseMat *)&(A.val(k)))->a;
00178         for (int ixx=0; ixx<nrow; ixx++)
00179           {
00180         for (int jxx=0; jxx<ncol; jxx++)
00181           cout <<setiosflags(ios::scientific) <<*(qxx+jxx*nrow+ixx)<<",  ";
00182         cout << endl;
00183           }
00184         */
00185       }
00186         ifp_LocalMat *mult = A.val(0).CreateEmpty();
00187 
00188         // eliminate the elements in L in order
00189         for (k=ial[i]; k<ial[i+1]; k++)
00190         {
00191             id = jal[k];
00192 
00193             // mult = row[id] / au[idiag[id]];
00194             row[marker[id]]->Mat_Mat_Mult(diag[id], mult, 1.0, 0.0);
00195 
00196             row[marker[id]]->MatCopy(*mult);
00197 
00198             for (kk=iau[id]+1; kk<iau[id+1]; kk++)
00199             {
00200                 idd = jau[kk];
00201                 if (marker[idd] >= 0) 
00202         {
00203                     // row[idd] = row[idd] - mult * au[kk];
00204                     mult->Mat_Mat_Mult(au[kk], row[marker[idd]], -1.0, 1.0);
00205         }
00206             }
00207         }
00208 
00209     delete mult;
00210 
00211         // gather resulting rows in L and U
00212         for (k=ial[i]; k<ial[i+1]; k++)
00213         {
00214             al[k] = row[marker[jal[k]]];
00215             marker[jal[k]] = -1;
00216         }
00217         for (k=iau[i]; k<iau[i+1]; k++)
00218         {
00219             au[k] = row[marker[jau[k]]];
00220             marker[jau[k]] = -1;
00221         }
00222 
00223         diag[i] = au[iau[i]]->CreateInv(local_precon);  // inverse of diagonal
00224     }
00225 
00226     delete [] row;
00227     delete [] marker;
00228 }
00229 
00230 void ifp_biluk::apply(int nr, int nc, const double *u, int ldu, 
00231     double *v, int ldv)
00232 {
00233     applyl(nr, nc, u, ldu, v, ldv);
00234     applyr(nr, nc, v, ldv, v, ldv);
00235 }
00236 
00237 void ifp_biluk::applyl(int nr, int nc, const double *u, int ldu, 
00238     double *v, int ldv)
00239 {
00240     int i, j;
00241 
00242     ifp_BlockVec  V(nr, nc, v, ldv, &Ap->kvst_col(0));
00243     ifp_BlockVec V2(nr, nc, v, ldv, &Ap->kvst_col(0));
00244     ifp_BlockVec  U(nr, nc, u, ldu, &Ap->kvst_col(0));
00245     V.VecCopy(U);
00246 
00247     // forward solve with lower triang factor (identity on diagonal assumed)
00248 
00249     for (i=0; i<Ap->numrow(); i++)
00250     {
00251         for (j=ial[i]; j<ial[i+1]; j++)
00252         {
00253             // V(i) = V(i) - al[j] * V(jal[j])
00254             al[j]->Mat_Vec_Mult(V(jal[j]), V2(i), -1.0, 1.0);
00255         }
00256     }
00257 
00258 }
00259 
00260 void ifp_biluk::applyr(int nr, int nc, const double *u, int ldu, 
00261     double *v, int ldv)
00262 {
00263     int i, j;
00264 
00265     ifp_BlockVec  V(nr, nc, v, ldv, &Ap->kvst_col(0));
00266     ifp_BlockVec V2(nr, nc, v, ldv, &Ap->kvst_col(0));
00267     ifp_BlockVec  U(nr, nc, u, ldu, &Ap->kvst_col(0));
00268     V.VecCopy(U);
00269 
00270     // backward solve with upper triang factor
00271 
00272     for (i=Ap->numrow()-1; i>=0; i--)
00273     {
00274         for (j=iau[i]+1; j<iau[i+1]; j++)
00275         {
00276             // V(i) = V(i) - au[j] * V(jau[j])
00277             au[j]->Mat_Vec_Mult(V(jau[j]), V2(i), -1.0, 1.0);
00278         }
00279         diag[i]->Mat_Vec_Solve(V(i), V(i));
00280     }
00281 }
00282 
00283 void ifp_biluk::multiply(int nr, int nc, const double *u, int ldu, 
00284     double *v, int ldv)
00285 {
00286     int i, j;
00287 
00288     ifp_BlockVec U(nr, nc, u, ldu, &Ap->kvst_col(0));
00289     ifp_BlockVec V(nr, nc, v, ldv, &Ap->kvst_col(0));
00290     ifp_BlockVec T(U);
00291     V.VecSetToZero();
00292 
00293     // multiply with upper triang factor
00294 
00295     for (i=0; i<Ap->numrow(); i++)
00296     {
00297         for (j=iau[i]; j<iau[i+1]; j++)
00298         {
00299             // V(i) = V(i) + au[j] * U(jau[j])
00300             au[j]->Mat_Vec_Mult(U(jau[j]), V(i), 1.0, 1.0);
00301         }
00302     }
00303 
00304     // multiply with lower triang factor (unit diagonal assumed)
00305     T.VecCopy(V);
00306 
00307     for (i=0; i<Ap->numrow(); i++)
00308     {
00309         for (j=ial[i]; j<ial[i+1]; j++)
00310         {
00311             // V(i) = V(i) + al[j] * temp(jal[j])
00312             al[j]->Mat_Vec_Mult(T(jal[j]), V(i), 1.0, 1.0);
00313         }
00314     }
00315 }
00316 double ifp_biluk::condest()
00317 {
00318   // This routine computes a  bound of the infinity-norm condition number 
00319   // of the preconditioner.
00320   // It is equal to infinity-norm of M^-1 e, where e is the vector of all
00321   // ones.
00322 
00323   int i;
00324 
00325   int m = Ap->dimrow(); // Object is a matrix, so it has row/col dimensions.
00326   int n = Ap->dimcol();
00327   double *u = new double[n];
00328   double *v = new double[m];
00329 
00330     for (i=0; i<n; i++) u[i] = 1.0;
00331 
00332     apply(n, 1, u, n, v, n);
00333 
00334     double cond_number = 0.0;
00335     for (i=0; i<m; i++)
00336         if (ABS(v[i]) > cond_number)
00337            cond_number = ABS(v[i]);
00338 
00339     delete [] u;
00340     delete [] v;
00341 
00342     return(cond_number);
00343 }

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