example/ifpack_hb/cxx_main.cpp

Go to the documentation of this file.
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 <stdio.h>
00032 #include <stdlib.h>
00033 #include <ctype.h>
00034 #include <assert.h>
00035 #include <string.h>
00036 #include <math.h>
00037 #ifdef EPETRA_MPI
00038 #include "mpi.h"
00039 #include "Epetra_MpiComm.h"
00040 #else
00041 #include "Epetra_SerialComm.h"
00042 #endif
00043 #include "Trilinos_Util.h"
00044 #include "Epetra_Comm.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Time.h"
00047 #include "Epetra_BlockMap.h"
00048 #include "Epetra_MultiVector.h"
00049 #include "Epetra_Vector.h"
00050 #include "Epetra_Export.h"
00051 
00052 #include "Epetra_VbrMatrix.h"
00053 #include "Epetra_CrsMatrix.h"
00054 #include "Ifpack_IlukGraph.h"
00055 #include "Ifpack_CrsRiluk.h"
00056 
00057 void BiCGSTAB(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, 
00058         Ifpack_CrsRiluk *M, 
00059         int Maxiter, double Tolerance, 
00060         double *residual, bool verbose);
00061 
00062 int main(int argc, char *argv[]) {
00063 
00064 #ifdef EPETRA_MPI
00065   MPI_Init(&argc,&argv);
00066   Epetra_MpiComm Comm (MPI_COMM_WORLD);
00067 #else
00068   Epetra_SerialComm Comm;
00069 #endif
00070 
00071   cout << Comm << endl;
00072 
00073   int MyPID = Comm.MyPID();
00074 
00075   bool verbose = false; 
00076   if (MyPID==0) verbose = true;
00077 
00078   if(argc < 2 && verbose) {
00079     cerr << "Usage: " << argv[0] 
00080    << " HB_filename [level_fill [level_overlap [absolute_threshold [ relative_threshold]]]]" << endl
00081    << "where:" << endl
00082    << "HB_filename        - filename and path of a Harwell-Boeing data set" << endl
00083    << "level_fill         - The amount of fill to use for ILU(k) preconditioner (default 0)" << endl
00084    << "level_overlap      - The amount of overlap used for overlapping Schwarz subdomains (default 0)" << endl
00085    << "absolute_threshold - The minimum value to place on the diagonal prior to factorization (default 0.0)" << endl
00086    << "relative_threshold - The relative amount to perturb the diagonal prior to factorization (default 1.0)" << endl << endl
00087    << "To specify a non-default value for one of these parameters, you must specify all" << endl
00088    << " preceding values but not any subsequent parameters. Example:" << endl
00089    << "ifpackHbSerialMsr.exe mymatrix.hb 1  - loads mymatrix.hb, uses level fill of one, all other values are defaults" << endl
00090    << endl;
00091     return(1);
00092 
00093   }
00094 
00095   // Uncomment the next three lines to debug in mpi mode
00096   //int tmp;
00097   //if (MyPID==0) cin >> tmp;
00098   //Comm.Barrier();
00099 
00100   Epetra_Map * readMap;
00101   Epetra_CrsMatrix * readA; 
00102   Epetra_Vector * readx; 
00103   Epetra_Vector * readb;
00104   Epetra_Vector * readxexact;
00105    
00106   // Call routine to read in HB problem
00107   Trilinos_Util_ReadHb2Epetra(argv[1], Comm, readMap, readA, readx, readb, readxexact);
00108 
00109   // Create uniform distributed map
00110   Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
00111 
00112   // Create Exporter to distribute read-in matrix and vectors
00113 
00114   Epetra_Export exporter(*readMap, map);
00115   Epetra_CrsMatrix A(Copy, map, 0);
00116   Epetra_Vector x(map);
00117   Epetra_Vector b(map);
00118   Epetra_Vector xexact(map);
00119 
00120   Epetra_Time FillTimer(Comm);
00121   x.Export(*readx, exporter, Add);
00122   b.Export(*readb, exporter, Add);
00123   xexact.Export(*readxexact, exporter, Add);
00124   Comm.Barrier();
00125   double vectorRedistributeTime = FillTimer.ElapsedTime();
00126   A.Export(*readA, exporter, Add);
00127   Comm.Barrier();
00128   double matrixRedistributeTime = FillTimer.ElapsedTime() - vectorRedistributeTime;
00129   assert(A.FillComplete()==0);    
00130   Comm.Barrier();
00131   double fillCompleteTime = FillTimer.ElapsedTime() - matrixRedistributeTime;
00132   if (Comm.MyPID()==0)  {
00133     cout << "\n\n****************************************************" << endl;
00134     cout << "\n Vector redistribute  time (sec) = " << vectorRedistributeTime<< endl;
00135     cout << "    Matrix redistribute time (sec) = " << matrixRedistributeTime << endl;
00136     cout << "    Transform to Local  time (sec) = " << fillCompleteTime << endl<< endl;
00137   }
00138   Epetra_Vector tmp1(*readMap);
00139   Epetra_Vector tmp2(map);
00140   readA->Multiply(false, *readxexact, tmp1);
00141 
00142   A.Multiply(false, xexact, tmp2);
00143   double residual;
00144   tmp1.Norm2(&residual);
00145   if (verbose) cout << "Norm of Ax from file            = " << residual << endl;
00146   tmp2.Norm2(&residual);
00147   if (verbose) cout << "Norm of Ax after redistribution = " << residual << endl << endl << endl;
00148 
00149   //cout << "A from file = " << *readA << endl << endl << endl;
00150 
00151   //cout << "A after dist = " << A << endl << endl << endl;
00152 
00153   delete readA;
00154   delete readx;
00155   delete readb;
00156   delete readxexact;
00157   delete readMap;
00158 
00159   Comm.Barrier();
00160 
00161   // Construct ILU preconditioner
00162 
00163   double elapsed_time, total_flops, MFLOPs;
00164   Epetra_Time timer(Comm);
00165 
00166   int LevelFill = 0;
00167   if (argc > 2)  LevelFill = atoi(argv[2]);
00168   if (verbose) cout << "Using Level Fill = " << LevelFill << endl;
00169   int Overlap = 0;
00170   if (argc > 3) Overlap = atoi(argv[3]);
00171   if (verbose) cout << "Using Level Overlap = " << Overlap << endl;
00172   double Athresh = 0.0;
00173   if (argc > 4) Athresh = atof(argv[4]);
00174   if (verbose) cout << "Using Absolute Threshold Value of = " << Athresh << endl;
00175 
00176   double Rthresh = 1.0;
00177   if (argc > 5) Rthresh = atof(argv[5]);
00178   if (verbose) cout << "Using Relative Threshold Value of = " << Rthresh << endl;
00179 
00180   Ifpack_IlukGraph * IlukGraph = 0;
00181   Ifpack_CrsRiluk * ILUK = 0;
00182 
00183   if (LevelFill>-1) {
00184     elapsed_time = timer.ElapsedTime();
00185     IlukGraph = new Ifpack_IlukGraph(A.Graph(), LevelFill, Overlap);
00186     assert(IlukGraph->ConstructFilledGraph()==0);
00187     elapsed_time = timer.ElapsedTime() - elapsed_time;
00188     if (verbose) cout << "Time to construct ILUK graph = " << elapsed_time << endl;
00189 
00190 
00191     Epetra_Flops fact_counter;
00192   
00193     elapsed_time = timer.ElapsedTime();
00194     ILUK = new Ifpack_CrsRiluk(*IlukGraph);
00195     ILUK->SetFlopCounter(fact_counter);
00196     ILUK->SetAbsoluteThreshold(Athresh);
00197     ILUK->SetRelativeThreshold(Rthresh);
00198     //assert(ILUK->InitValues()==0);
00199     int initerr = ILUK->InitValues(A);
00200     if (initerr!=0) cout << Comm << "InitValues error = " << initerr;
00201     assert(ILUK->Factor()==0);
00202     elapsed_time = timer.ElapsedTime() - elapsed_time;
00203     total_flops = ILUK->Flops();
00204     MFLOPs = total_flops/elapsed_time/1000000.0;
00205     if (verbose) cout << "Time to compute preconditioner values = " 
00206         << elapsed_time << endl
00207         << "MFLOPS for Factorization = " << MFLOPs << endl;
00208     //cout << *ILUK << endl;
00209   }
00210   double Condest;
00211   ILUK->Condest(false, Condest);
00212 
00213   if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
00214   int Maxiter = 500;
00215   double Tolerance = 1.0E-14;
00216 
00217   Epetra_Vector xcomp(map);
00218   Epetra_Vector resid(map);
00219 
00220   Epetra_Flops counter;
00221   A.SetFlopCounter(counter);
00222   xcomp.SetFlopCounter(A);
00223   b.SetFlopCounter(A);
00224   resid.SetFlopCounter(A);
00225   ILUK->SetFlopCounter(A);
00226 
00227   elapsed_time = timer.ElapsedTime();
00228 
00229   BiCGSTAB(A, xcomp, b, ILUK, Maxiter, Tolerance, &residual, verbose);
00230 
00231   elapsed_time = timer.ElapsedTime() - elapsed_time;
00232   total_flops = counter.Flops();
00233   MFLOPs = total_flops/elapsed_time/1000000.0;
00234   if (verbose) cout << "Time to compute solution = " 
00235         << elapsed_time << endl
00236         << "Number of operations in solve = " << total_flops << endl
00237         << "MFLOPS for Solve = " << MFLOPs<< endl << endl;
00238 
00239   resid.Update(1.0, xcomp, -1.0, xexact, 0.0); // resid = xcomp - xexact
00240 
00241   resid.Norm2(&residual);
00242 
00243   if (verbose) cout << "Norm of the difference between exact and computed solutions = " << residual << endl;
00244 
00245   
00246 
00247 
00248   if (ILUK!=0) delete ILUK;
00249   if (IlukGraph!=0) delete IlukGraph;
00250                
00251 #ifdef EPETRA_MPI
00252   MPI_Finalize() ;
00253 #endif
00254 
00255 return 0 ;
00256 }
00257 void BiCGSTAB(Epetra_CrsMatrix &A, 
00258         Epetra_Vector &x, 
00259         Epetra_Vector &b, 
00260         Ifpack_CrsRiluk *M, 
00261         int Maxiter, 
00262         double Tolerance, 
00263         double *residual, bool verbose) {
00264 
00265   // Allocate vectors needed for iterations
00266   Epetra_Vector phat(x.Map()); phat.SetFlopCounter(x);
00267   Epetra_Vector p(x.Map()); p.SetFlopCounter(x);
00268   Epetra_Vector shat(x.Map()); shat.SetFlopCounter(x);
00269   Epetra_Vector s(x.Map()); s.SetFlopCounter(x);
00270   Epetra_Vector r(x.Map()); r.SetFlopCounter(x);
00271   Epetra_Vector rtilde(x.Map()); rtilde.Random(); rtilde.SetFlopCounter(x);
00272   Epetra_Vector v(x.Map()); v.SetFlopCounter(x);
00273   
00274 
00275   A.Multiply(false, x, r); // r = A*x
00276 
00277   r.Update(1.0, b, -1.0); // r = b - A*x
00278 
00279   double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
00280   double alpha = 1.0, omega = 1.0, sigma;
00281   double omega_num, omega_den;
00282   r.Norm2(&r_norm);
00283   b.Norm2(&b_norm);
00284   scaled_r_norm = r_norm/b_norm;
00285   r.Dot(rtilde,&rhon);
00286 
00287   if (verbose) cout << "Initial residual = " << r_norm
00288         << " Scaled residual = " << scaled_r_norm << endl;
00289 
00290 
00291   for (int i=0; i<Maxiter; i++) { // Main iteration loop   
00292 
00293     double beta = (rhon/rhonm1) * (alpha/omega);
00294     rhonm1 = rhon;
00295 
00296     /* p    = r + beta*(p - omega*v)       */
00297     /* phat = M^-1 p                       */
00298     /* v    = A phat                       */
00299 
00300     double dtemp = - beta*omega;
00301 
00302     p.Update(1.0, r, dtemp, v, beta);
00303     if (M==0) 
00304       phat.Scale(1.0, p);
00305     else
00306       M->Solve(false, p, phat);
00307     A.Multiply(false, phat, v);
00308 
00309     
00310     rtilde.Dot(v,&sigma);
00311     alpha = rhon/sigma;    
00312 
00313     /* s = r - alpha*v                     */
00314     /* shat = M^-1 s                       */
00315     /* r = A shat (r is a tmp here for t ) */
00316 
00317     s.Update(-alpha, v, 1.0, r, 0.0);
00318     if (M==0) 
00319       shat.Scale(1.0, s);
00320     else
00321       M->Solve(false, s, shat);
00322     A.Multiply(false, shat, r);
00323 
00324     r.Dot(s, &omega_num);
00325     r.Dot(r, &omega_den);
00326     omega = omega_num / omega_den;
00327 
00328     /* x = x + alpha*phat + omega*shat */
00329     /* r = s - omega*r */
00330 
00331     x.Update(alpha, phat, omega, shat, 1.0);
00332     r.Update(1.0, s, -omega); 
00333     
00334     r.Norm2(&r_norm);
00335     scaled_r_norm = r_norm/b_norm;
00336     r.Dot(rtilde,&rhon);
00337 
00338     if (verbose) cout << "Iter "<< i << " residual = " << r_norm
00339           << " Scaled residual = " << scaled_r_norm << endl;
00340 
00341     if (scaled_r_norm < Tolerance) break;
00342   }
00343   return;
00344 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:33 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3