ifp_brelax.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_LocalMat.h"
00030 #include "ifp_BlockVec.h"
00031 #include "ifp_BlockMat.h"
00032 #include "ifp_brelax.h"
00033 #include "ifp_ifpack.h"
00034 
00035 void ifp_None::apply(int nr, int nc, const double *u, int ldu, 
00036     double *v, int ldv)
00037 {
00038     int i, j;
00039     double *p;
00040     const double *q;
00041 
00042     for (i=0; i<nc; i++)
00043     {
00044         p = v + i*ldv;
00045         q = u + i*ldu;
00046         for (j=0; j<nr; j++)
00047             *p++ = *q++;
00048     }
00049 }
00050 
00051 ifp_BJacobi::ifp_BJacobi()
00052 {
00053     Ap = (ifp_BlockMat *) NULL;
00054     diag = (ifp_LocalMat **) NULL;
00055 }
00056 
00057 ifp_BJacobi::~ifp_BJacobi()
00058 {
00059     if (Ap != NULL)
00060         for (int i=0; i<Ap->numrow(); i++)
00061             delete diag[i];
00062     delete [] diag;
00063 }
00064 
00065 void ifp_BJacobi::setup(const ifp_BlockMat& A)
00066 {
00067     int i, j;
00068     int got_diag;
00069 
00070     Ap = &A;
00071     diag = new ifp_LocalMatp[A.numrow()];
00072 
00073     // search for diagonal blocks
00074     for (i=0; i<A.numrow(); i++)
00075     {
00076         got_diag = FALSE;
00077         for (j=A.row_ptr(i); j<A.row_ptr(i+1); j++)
00078         {
00079             if (A.col_ind(j) == i)
00080             {
00081                 diag[i] = A.val(j).CreateInv(local_precon);
00082                 got_diag = TRUE;
00083             }
00084         }
00085         if (!got_diag)
00086             ifp_error("ifp_brelax: matrix does not have diagonal block", i);
00087     }
00088 }
00089 
00090 void ifp_BJacobi::apply(int nr, int nc, const double *u, int ldu,
00091     double *v, int ldv)
00092 {
00093     ifp_BlockVec U(nr, nc, u, ldu, &Ap->kvst_row(0));
00094     ifp_BlockVec V(nr, nc, v, ldv, &Ap->kvst_row(0));
00095 
00096     for (int i=0; i<Ap->numrow(); i++)
00097     {
00098         diag[i]->Mat_Vec_Solve(U(i), V(i));
00099     }
00100 }
00101 
00102 double ifp_BJacobi::condest()
00103 {
00104   // This routine computes a  bound of the infinity-norm condition number 
00105   // of the preconditioner.
00106   // It is equal to infinity-norm of M^-1 e, where e is the vector of all
00107   // ones.
00108 
00109   int i;
00110 
00111   int m = Ap->dimrow(); // Object is a matrix, so it has row/col dimensions.
00112   int n = Ap->dimcol();
00113   double *u = new double[n];
00114   double *v = new double[m];
00115 
00116     for (i=0; i<n; i++) u[i] = 1.0;
00117 
00118     apply(n, 1, u, n, v, n);
00119 
00120     double cond_number = 0.0;
00121     for (i=0; i<m; i++)
00122         if (ABS(v[i]) > cond_number)
00123            cond_number = ABS(v[i]);
00124 
00125     delete [] u;
00126     delete [] v;
00127 
00128     return(cond_number);
00129 }
00130 
00131 // note:
00132 // assumes each matrix row is stored such that lower triangular elements,
00133 // diagonal elements, upper triangular elements are stored in order
00134 
00135 ifp_BSOR_Base::ifp_BSOR_Base()
00136 {
00137     Ap = (ifp_BlockMat *) NULL;
00138     diag = (ifp_LocalMat **) NULL;
00139     idiag = (int *) NULL;
00140 }
00141 
00142 ifp_BSOR_Base::~ifp_BSOR_Base()
00143 {
00144     if (Ap != NULL)
00145         for (int i=0; i<Ap->numrow(); i++)
00146             delete diag[i];
00147     delete [] diag;
00148     delete [] idiag;
00149 }
00150 
00151 void ifp_BSOR_Base::setup(const ifp_BlockMat& A, double omega, int iterations)
00152 {
00153     int i, j, nrow;
00154     int got_diag;
00155 
00156     omega_ = omega;
00157     iterations_ = iterations;
00158     Ap = &A;
00159     nrow = A.numrow();
00160     diag = new ifp_LocalMatp[nrow];
00161     idiag = new int[nrow];
00162 
00163     // search for diagonal blocks
00164     for (i=0; i<nrow; i++)
00165     {
00166         got_diag = FALSE;
00167         for (j=A.row_ptr(i); j<A.row_ptr(i+1); j++)
00168         {
00169             if (A.col_ind(j) == i)
00170             {
00171                 diag[i] = A.val(j).CreateInv(local_precon);
00172                 idiag[i] = j;
00173                 got_diag = TRUE;
00174             }
00175         }
00176         if (!got_diag)
00177             ifp_error("ifp_brelax: matrix does not have diagonal block", i);
00178     }
00179 }
00180 
00181 // c = alpha * a + beta * b
00182 // works on Blocks within a ifp_BlockVec
00183 
00184 static void gaxpy(const double& alpha, const ifp_BlockVec& a,
00185     const double& beta, const ifp_BlockVec& b, ifp_BlockVec& c)
00186 {
00187     double *ap, *bp, *cp;
00188     for (int j=0; j<a.dim1; j++)
00189     {
00190         ap = a.v + j * a.ld;
00191         bp = b.v + j * b.ld;
00192         cp = c.v + j * c.ld;
00193 
00194         for (int i=0; i<a.size0; i++) 
00195         *cp++ = alpha * *ap++ + beta * *bp++;
00196             // c(i,j) = alpha*a(i,j) + beta*b(i,j);
00197     }
00198 }
00199 
00200 void ifp_BSOR::apply(int nr, int nc, const double *u, int ldu,
00201     double *v, int ldv)
00202 {
00203     const int *ia = &Ap->row_ptr(0);
00204     const int *ja = &Ap->col_ind(0);
00205     int it, i, j;
00206 
00207 #if 0
00208     ifp_BlockVec  V(nr, nc, v, ldv, &Ap->kvst_col(0));
00209     ifp_BlockVec V2(nr, nc, v, ldv, &Ap->kvst_col(0));
00210     ifp_BlockVec  U(nr, nc, u, ldu, &Ap->kvst_col(0));
00211     V.VecCopy(U);
00212 
00213     // Specialized code for first step.
00214 
00215     for (i=0; i<Ap->numrow(); i++)
00216     {
00217         for (j=ia[i]; j<idiag[i]; j++)
00218         {
00219             // V(i) = V(i) - omega_ a[j] * V(ja[j])
00220             Ap->val(j).Mat_Vec_Mult(V(ja[j]), V2(i), -omega_, 1.0);
00221         }
00222         diag[i]->Mat_Vec_Solve(V(i), V(i));
00223     }
00224 
00225     // After first step....
00226 #endif
00227     ifp_BlockVec  V(nr, nc, v, ldv, &Ap->kvst_col(0));
00228     ifp_BlockVec  U(nr, nc, u, ldu, &Ap->kvst_col(0));
00229     V.VecSetToZero();
00230 
00231     for (it=1; it<=iterations_; it++)
00232     {
00233         for (i=0; i<Ap->numrow(); i++)
00234         {
00235             ifp_BlockVec temp(U,i);
00236 
00237             for (j=ia[i]; j<idiag[i]; j++)
00238             {
00239                 // temp = temp - a[j] * v[ja[j]];
00240                 Ap->val(j).Mat_Vec_Mult(V(ja[j]), temp, -1.0, 1.0);
00241             }
00242             for (j=idiag[i]+1; j<ia[i+1]; j++)
00243             {
00244                 // temp = temp - a[j] * v[ja[j]];
00245                 Ap->val(j).Mat_Vec_Mult(V(ja[j]), temp, -1.0, 1.0);
00246             }
00247 
00248             diag[i]->Mat_Vec_Solve(temp, temp);
00249 
00250             // v[i] = (1.0-omega_) * v[i] + omega_ * temp
00251             gaxpy(1.0-omega_, V(i), omega_, temp, V(i));
00252         }
00253     }
00254 }
00255 
00256 void ifp_BSSOR::apply(int nr, int nc, const double *u, int ldu,
00257     double *v, int ldv)
00258 {
00259     const int *ia = &Ap->row_ptr(0);
00260     const int *ja = &Ap->col_ind(0);
00261     int it, i, j;
00262 
00263 #if 0
00264     ifp_BlockVec  V(nr, nc, v, ldv, &Ap->kvst_col(0));
00265     ifp_BlockVec V2(nr, nc, v, ldv, &Ap->kvst_col(0));
00266     ifp_BlockVec  U(nr, nc, u, ldu, &Ap->kvst_col(0));
00267     V.VecCopy(U);
00268 
00269     // Specialized code for first step.
00270 
00271     // lower sweep
00272     for (i=0; i<Ap->numrow(); i++)
00273     {
00274         for (j=ia[i]; j<idiag[i]; j++)
00275         {
00276             // V(i) = V(i) - omega_ a[j] * V(ja[j])
00277             Ap->val(j).Mat_Vec_Mult(V(ja[j]), V2(i), -omega_, 1.0);
00278         }
00279         diag[i]->Mat_Vec_Solve(V(i), V(i));
00280     }
00281 
00282     // multiply by diagonal blocks
00283     for (i=0; i<Ap->numrow(); i++)
00284     {
00285         // V(i) = diag[i] * V(i)
00286         // this cannot be done in place
00287         ifp_BlockVec y(V,i); // make a copy of V(i)
00288 
00289         Ap->val(idiag[i]).Mat_Vec_Mult(y, V(i));
00290     }
00291 
00292     // upper sweep
00293     for (i=Ap->numrow()-1; i>=0; i--)
00294     {
00295         for (j=idiag[i]+1; j<ia[i+1]; j++)
00296         {
00297             // V(i) = V(i) - omega_ a[j] * V(ja[j])
00298             Ap->val(j).Mat_Vec_Mult(V(ja[j]), V2(i), -omega_, 1.0);
00299         }
00300         diag[i]->Mat_Vec_Solve(V(i), V(i));
00301     }
00302 
00303     // After first step....
00304 #endif
00305 
00306     ifp_BlockVec  V(nr, nc, v, ldv, &Ap->kvst_col(0));
00307     ifp_BlockVec  U(nr, nc, u, ldu, &Ap->kvst_col(0));
00308     V.VecSetToZero();
00309 
00310     for (it=1; it<=iterations_; it++)
00311     {
00312         for (i=0; i<Ap->numrow(); i++)
00313         {
00314             ifp_BlockVec temp(U,i);
00315 
00316             for (j=ia[i]; j<idiag[i]; j++)
00317             {
00318                 // temp = temp - a[j] * v[ja[j]];
00319                 Ap->val(j).Mat_Vec_Mult(V(ja[j]), temp, -1.0, 1.0);
00320             }
00321             for (j=idiag[i]+1; j<ia[i+1]; j++)
00322             {
00323                 // temp = temp - a[j] * v[ja[j]];
00324                 Ap->val(j).Mat_Vec_Mult(V(ja[j]), temp, -1.0, 1.0);
00325             }
00326 
00327             diag[i]->Mat_Vec_Solve(temp, temp);
00328 
00329             // v[i] = (1.0-omega_) * v[i] + omega_ * temp
00330             gaxpy(1.0-omega_, V(i), omega_, temp, V(i));
00331         }
00332 
00333         for (i=Ap->numrow()-1; i>=0; i--)
00334         {
00335             ifp_BlockVec temp(U,i);
00336 
00337             for (j=ia[i]; j<idiag[i]; j++)
00338             {
00339                 // temp = temp - a[j] * v[ja[j]];
00340                 Ap->val(j).Mat_Vec_Mult(V(ja[j]), temp, -1.0, 1.0);
00341             }
00342             for (j=idiag[i]+1; j<ia[i+1]; j++)
00343             {
00344                 // temp = temp - a[j] * v[ja[j]];
00345                 Ap->val(j).Mat_Vec_Mult(V(ja[j]), temp, -1.0, 1.0);
00346             }
00347 
00348             diag[i]->Mat_Vec_Solve(temp, temp);
00349 
00350             // v[i] = (1.0-omega_) * v[i] + omega_ * temp
00351             gaxpy(1.0-omega_, V(i), omega_, temp, V(i));
00352         }
00353     }
00354 }

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1