example/ReducedLinearProblem/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright (2001) 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 "Epetra_Comm.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Time.h"
00032 #include "Epetra_BlockMap.h"
00033 #include "Epetra_MultiVector.h"
00034 #include "Epetra_Vector.h"
00035 #include "Epetra_Export.h"
00036 #include "Epetra_LinearProblem.h"
00037 #include "Epetra_CrsSingletonFilter.h"
00038 
00039 #include "Epetra_VbrMatrix.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #ifdef EPETRA_MPI
00042 #include "mpi.h"
00043 #include "Epetra_MpiComm.h"
00044 #else
00045 #include "Epetra_SerialComm.h"
00046 #endif
00047 #include "Trilinos_Util.h"
00048 #include "Ifpack_IlukGraph.h"
00049 #include "Ifpack_CrsRiluk.h"
00050 
00051 void BiCGSTAB(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, 
00052         Ifpack_CrsRiluk *M, 
00053         int Maxiter, double Tolerance, 
00054         double *residual, bool verbose);
00055 int Statistics(const Epetra_CrsSingletonFilter & filter);
00056 int main(int argc, char *argv[]) {
00057 
00058 #ifdef EPETRA_MPI
00059   MPI_Init(&argc,&argv);
00060   Epetra_MpiComm Comm (MPI_COMM_WORLD);
00061 #else
00062   Epetra_SerialComm Comm;
00063 #endif
00064 
00065   cout << Comm << endl;
00066 
00067   int MyPID = Comm.MyPID();
00068 
00069   bool verbose = false;
00070   bool verbose1 = true;
00071   if (MyPID==0) verbose = true;
00072 
00073   if(argc < 2 && verbose) {
00074     cerr << "Usage: " << argv[0] 
00075    << " HPC_filename [level_fill [level_overlap [absolute_threshold [ relative_threshold]]]]" << endl
00076    << "where:" << endl
00077    << "HB_filename        - filename and path of a Harwell-Boeing data set" << endl
00078    << "level_fill         - The amount of fill to use for ILU(k) preconditioner (default 0)" << endl
00079    << "level_overlap      - The amount of overlap used for overlapping Schwarz subdomains (default 0)" << endl
00080    << "absolute_threshold - The minimum value to place on the diagonal prior to factorization (default 0.0)" << endl
00081    << "relative_threshold - The relative amount to perturb the diagonal prior to factorization (default 1.0)" << endl << endl
00082    << "To specify a non-default value for one of these parameters, you must specify all" << endl
00083    << " preceding values but not any subsequent parameters. Example:" << endl
00084    << "ifpackHpcSerialMsr.exe mymatrix.hpc 1  - loads mymatrix.hpc, uses level fill of one, all other values are defaults" << endl
00085    << endl;
00086     return(1);
00087 
00088   }
00089 
00090   // Uncomment the next three lines to debug in mpi mode
00091   int tmp;
00092   if (MyPID==0) cin >> tmp;
00093   Comm.Barrier();
00094 
00095   Epetra_Map * map;
00096   Epetra_CrsMatrix * A; 
00097   Epetra_Vector * x; 
00098   Epetra_Vector * b;
00099   Epetra_Vector * xexact;
00100 
00101   Trilinos_Util_ReadHpc2Epetra(argv[1], Comm, map, A, x, b, xexact);
00102 
00103   bool smallProblem = false;
00104   if (A->RowMap().NumGlobalElements()<100) smallProblem = true;
00105 
00106   if (smallProblem)
00107     cout << "Original Matrix = " << endl << *A   << endl;
00108 
00109   x->PutScalar(0.0);
00110 
00111   Epetra_LinearProblem FullProblem(A, x, b);
00112   double normb, norma;
00113   b->NormInf(&normb);
00114   norma = A->NormInf();
00115   if (verbose)
00116     cout << "Inf norm of Original Matrix = " << A->NormInf() << endl
00117    << "Inf norm of Original RHS    = " << normb << endl;
00118   
00119   Epetra_Time ReductionTimer(Comm);
00120   Epetra_CrsSingletonFilter SingletonFilter;
00121   Comm.Barrier();
00122   double reduceInitTime = ReductionTimer.ElapsedTime();
00123   SingletonFilter.Analyze(A);
00124   Comm.Barrier();
00125   double reduceAnalyzeTime = ReductionTimer.ElapsedTime() - reduceInitTime;
00126 
00127   if (SingletonFilter.SingletonsDetected())
00128     cout << "Singletons found" << endl;
00129   else {
00130     cout << "Singletons not found" << endl;
00131     exit(1);
00132   }
00133   SingletonFilter.ConstructReducedProblem(&FullProblem);
00134   Comm.Barrier();
00135   double reduceConstructTime = ReductionTimer.ElapsedTime() - reduceInitTime;
00136 
00137   double totalReduceTime = ReductionTimer.ElapsedTime();
00138 
00139   if (verbose)
00140     cout << "\n\n****************************************************" << endl
00141    << "    Reduction init  time (sec)           = " << reduceInitTime<< endl
00142    << "    Reduction Analyze time (sec)         = " << reduceAnalyzeTime << endl
00143    << "    Construct Reduced Problem time (sec) = " << reduceConstructTime << endl
00144    << "    Reduction Total time (sec)           = " << totalReduceTime << endl<< endl;
00145 
00146   Statistics(SingletonFilter);
00147 
00148   Epetra_LinearProblem * ReducedProblem = SingletonFilter.ReducedProblem();
00149 
00150   Epetra_CrsMatrix * Ap = dynamic_cast<Epetra_CrsMatrix *>(ReducedProblem->GetMatrix());
00151   Epetra_Vector * bp = (*ReducedProblem->GetRHS())(0);
00152   Epetra_Vector * xp = (*ReducedProblem->GetLHS())(0);
00153   
00154 
00155   if (smallProblem)
00156     cout << " Reduced Matrix = " << endl << *Ap << endl
00157    << " LHS before sol = " << endl << *xp << endl
00158    << " RHS            = " << endl << *bp << endl;
00159 
00160   // Construct ILU preconditioner
00161 
00162   double elapsed_time, total_flops, MFLOPs;
00163   Epetra_Time timer(Comm);
00164 
00165   int LevelFill = 0;
00166   if (argc > 2)  LevelFill = atoi(argv[2]);
00167   if (verbose) cout << "Using Level Fill = " << LevelFill << endl;
00168   int Overlap = 0;
00169   if (argc > 3) Overlap = atoi(argv[3]);
00170   if (verbose) cout << "Using Level Overlap = " << Overlap << endl;
00171   double Athresh = 0.0;
00172   if (argc > 4) Athresh = atof(argv[4]);
00173   if (verbose) cout << "Using Absolute Threshold Value of = " << Athresh << endl;
00174 
00175   double Rthresh = 1.0;
00176   if (argc > 5) Rthresh = atof(argv[5]);
00177   if (verbose) cout << "Using Relative Threshold Value of = " << Rthresh << endl;
00178 
00179   Ifpack_IlukGraph * IlukGraph = 0;
00180   Ifpack_CrsRiluk * ILUK = 0;
00181 
00182   if (LevelFill>-1) {
00183     elapsed_time = timer.ElapsedTime();
00184     IlukGraph = new Ifpack_IlukGraph(Ap->Graph(), LevelFill, Overlap);
00185     assert(IlukGraph->ConstructFilledGraph()==0);
00186     elapsed_time = timer.ElapsedTime() - elapsed_time;
00187     if (verbose) cout << "Time to construct ILUK graph = " << elapsed_time << endl;
00188 
00189 
00190     Epetra_Flops fact_counter;
00191   
00192     elapsed_time = timer.ElapsedTime();
00193     ILUK = new Ifpack_CrsRiluk(*IlukGraph);
00194     ILUK->SetFlopCounter(fact_counter);
00195     ILUK->SetAbsoluteThreshold(Athresh);
00196     ILUK->SetRelativeThreshold(Rthresh);
00197     //assert(ILUK->InitValues()==0);
00198     int initerr = ILUK->InitValues(*Ap);
00199     if (initerr!=0) {
00200       cout << endl << Comm << endl << "  InitValues error = " << initerr;
00201       if (initerr==1) cout << "  Zero diagonal found, warning error only";
00202       cout << endl << endl;
00203     }
00204     assert(ILUK->Factor()==0);
00205     elapsed_time = timer.ElapsedTime() - elapsed_time;
00206     total_flops = ILUK->Flops();
00207     MFLOPs = total_flops/elapsed_time/1000000.0;
00208     if (verbose) cout << "Time to compute preconditioner values = " 
00209         << elapsed_time << endl
00210         << "MFLOPS for Factorization = " << MFLOPs << endl;
00211     //cout << *ILUK << endl;
00212   double Condest;
00213   ILUK->Condest(false, Condest);
00214 
00215   if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
00216   }
00217   int Maxiter = 100;
00218   double Tolerance = 1.0E-8;
00219 
00220   Epetra_Flops counter;
00221   Ap->SetFlopCounter(counter);
00222   xp->SetFlopCounter(*Ap);
00223   bp->SetFlopCounter(*Ap);
00224   if (ILUK!=0) ILUK->SetFlopCounter(*Ap);
00225 
00226   elapsed_time = timer.ElapsedTime();
00227 
00228   double normreducedb, normreduceda;
00229   bp->NormInf(&normreducedb);
00230   normreduceda = Ap->NormInf();
00231   if (verbose) 
00232     cout << "Inf norm of Reduced Matrix = " << normreduceda << endl
00233    << "Inf norm of Reduced RHS    = " << normreducedb << endl;
00234 
00235   double residual;
00236   BiCGSTAB(*Ap, *xp, *bp, ILUK, Maxiter, Tolerance, &residual, verbose);
00237 
00238   elapsed_time = timer.ElapsedTime() - elapsed_time;
00239   total_flops = counter.Flops();
00240   MFLOPs = total_flops/elapsed_time/1000000.0;
00241   if (verbose) cout << "Time to compute solution = " 
00242         << elapsed_time << endl
00243         << "Number of operations in solve = " << total_flops << endl
00244         << "MFLOPS for Solve = " << MFLOPs<< endl << endl;
00245 
00246   SingletonFilter.ComputeFullSolution();
00247 
00248   if (smallProblem)
00249   cout << " Reduced LHS after sol = " << endl << *xp << endl
00250        << " Full    LHS after sol = " << endl << *x << endl
00251        << " Full  Exact LHS         = " << endl << *xexact << endl;
00252 
00253   Epetra_Vector resid(*x);
00254 
00255   resid.Update(1.0, *x, -1.0, *xexact, 0.0); // resid = xcomp - xexact
00256 
00257   resid.Norm2(&residual);
00258   double normx, normxexact;
00259   x->Norm2(&normx);
00260   xexact->Norm2(&normxexact);
00261 
00262   if (verbose) 
00263     cout << "2-norm of computed solution                               = " << normx << endl
00264    << "2-norm of exact solution                                  = " << normxexact << endl
00265    << "2-norm of difference between computed and exact solution  = " << residual << endl;
00266     
00267   if (verbose1 && residual>1.0e-5) {
00268     if (verbose)
00269       cout << "Difference between computed and exact solution appears large..." << endl
00270      << "Computing norm of A times this difference.  If this norm is small, then matrix is singular"
00271      << endl;
00272     Epetra_Vector bdiff(*b);
00273     assert(A->Multiply(false, resid, bdiff)==0);
00274     assert(bdiff.Norm2(&residual)==0);
00275     if (verbose) 
00276       cout << "2-norm of A times difference between computed and exact solution  = " << residual << endl;
00277     
00278   }
00279   if (verbose) 
00280     cout << "********************************************************" << endl
00281    << "              Solving again with 2*Ax=2*b" << endl
00282    << "********************************************************" << endl;
00283 
00284   A->Scale(1.0); // A = 2*A
00285   b->Scale(1.0); // b = 2*b
00286   x->PutScalar(0.0);
00287   b->NormInf(&normb);
00288   norma = A->NormInf();
00289   if (verbose)
00290     cout << "Inf norm of Original Matrix = " << norma << endl
00291    << "Inf norm of Original RHS    = " << normb << endl;
00292   double updateReducedProblemTime = ReductionTimer.ElapsedTime();
00293   SingletonFilter.UpdateReducedProblem(&FullProblem);
00294   Comm.Barrier();
00295   updateReducedProblemTime = ReductionTimer.ElapsedTime() - updateReducedProblemTime;
00296   if (verbose)
00297     cout << "\n\n****************************************************" << endl
00298    << "    Update Reduced Problem time (sec)           = " << updateReducedProblemTime<< endl
00299    << "****************************************************" << endl;
00300   Statistics(SingletonFilter);
00301 
00302   if (LevelFill>-1) {
00303 
00304     Epetra_Flops fact_counter;
00305   
00306     elapsed_time = timer.ElapsedTime();
00307 
00308     int initerr = ILUK->InitValues(*Ap);
00309     if (initerr!=0) {
00310       cout << endl << Comm << endl << "  InitValues error = " << initerr;
00311       if (initerr==1) cout << "  Zero diagonal found, warning error only";
00312       cout << endl << endl;
00313     }
00314     assert(ILUK->Factor()==0);
00315     elapsed_time = timer.ElapsedTime() - elapsed_time;
00316     total_flops = ILUK->Flops();
00317     MFLOPs = total_flops/elapsed_time/1000000.0;
00318     if (verbose) cout << "Time to compute preconditioner values = " 
00319         << elapsed_time << endl
00320         << "MFLOPS for Factorization = " << MFLOPs << endl;
00321     double Condest;
00322     ILUK->Condest(false, Condest);
00323     
00324     if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
00325   }
00326   bp->NormInf(&normreducedb);
00327   normreduceda = Ap->NormInf();
00328   if (verbose) 
00329     cout << "Inf norm of Reduced Matrix = " << normreduceda << endl
00330    << "Inf norm of Reduced RHS    = " << normreducedb << endl;
00331 
00332   BiCGSTAB(*Ap, *xp, *bp, ILUK, Maxiter, Tolerance, &residual, verbose);
00333 
00334   elapsed_time = timer.ElapsedTime() - elapsed_time;
00335   total_flops = counter.Flops();
00336   MFLOPs = total_flops/elapsed_time/1000000.0;
00337   if (verbose) cout << "Time to compute solution = " 
00338         << elapsed_time << endl
00339         << "Number of operations in solve = " << total_flops << endl
00340         << "MFLOPS for Solve = " << MFLOPs<< endl << endl;
00341 
00342   SingletonFilter.ComputeFullSolution();
00343 
00344   if (smallProblem)
00345   cout << " Reduced LHS after sol = " << endl << *xp << endl
00346        << " Full    LHS after sol = " << endl << *x << endl
00347        << " Full  Exact LHS         = " << endl << *xexact << endl;
00348 
00349   resid.Update(1.0, *x, -1.0, *xexact, 0.0); // resid = xcomp - xexact
00350 
00351   resid.Norm2(&residual);
00352   x->Norm2(&normx);
00353   xexact->Norm2(&normxexact);
00354 
00355   if (verbose) 
00356     cout << "2-norm of computed solution                               = " << normx << endl
00357    << "2-norm of exact solution                                  = " << normxexact << endl
00358    << "2-norm of difference between computed and exact solution  = " << residual << endl;
00359     
00360   if (verbose1 && residual>1.0e-5) {
00361     if (verbose)
00362       cout << "Difference between computed and exact solution appears large..." << endl
00363      << "Computing norm of A times this difference.  If this norm is small, then matrix is singular"
00364      << endl;
00365     Epetra_Vector bdiff(*b);
00366     assert(A->Multiply(false, resid, bdiff)==0);
00367     assert(bdiff.Norm2(&residual)==0);
00368     if (verbose) 
00369       cout << "2-norm of A times difference between computed and exact solution  = " << residual << endl;
00370     
00371   }
00372  
00373 
00374 
00375   if (ILUK!=0) delete ILUK;
00376   if (IlukGraph!=0) delete IlukGraph;
00377                
00378 #ifdef EPETRA_MPI
00379   MPI_Finalize() ;
00380 #endif
00381 
00382 return 0 ;
00383 }
00384 void BiCGSTAB(Epetra_CrsMatrix &A, 
00385         Epetra_Vector &x, 
00386         Epetra_Vector &b, 
00387         Ifpack_CrsRiluk *M, 
00388         int Maxiter, 
00389         double Tolerance, 
00390         double *residual, bool verbose) {
00391 
00392   // Allocate vectors needed for iterations
00393   Epetra_Vector phat(x.Map()); phat.SetFlopCounter(x);
00394   Epetra_Vector p(x.Map()); p.SetFlopCounter(x);
00395   Epetra_Vector shat(x.Map()); shat.SetFlopCounter(x);
00396   Epetra_Vector s(x.Map()); s.SetFlopCounter(x);
00397   Epetra_Vector r(b.Map()); r.SetFlopCounter(x);
00398   Epetra_Vector rtilde(x.Map()); rtilde.PutScalar(1.0); rtilde.SetFlopCounter(x);
00399   Epetra_Vector v(x.Map()); v.SetFlopCounter(x);
00400   
00401 
00402   A.Multiply(false, x, r); // r = A*x
00403 
00404   r.Update(1.0, b, -1.0); // r = b - A*x
00405 
00406   double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
00407   double alpha = 1.0, omega = 1.0, sigma;
00408   double omega_num, omega_den;
00409   r.Norm2(&r_norm);
00410   b.Norm2(&b_norm);
00411   scaled_r_norm = r_norm/b_norm;
00412   r.Dot(rtilde,&rhon);
00413 
00414   if (verbose) cout << "Initial residual = " << r_norm
00415         << " Scaled residual = " << scaled_r_norm << endl;
00416 
00417 
00418   for (int i=0; i<Maxiter; i++) { // Main iteration loop   
00419 
00420     double beta = (rhon/rhonm1) * (alpha/omega);
00421     rhonm1 = rhon;
00422 
00423     /* p    = r + beta*(p - omega*v)       */
00424     /* phat = M^-1 p                       */
00425     /* v    = A phat                       */
00426 
00427     double dtemp = - beta*omega;
00428 
00429     p.Update(1.0, r, dtemp, v, beta);
00430     if (M==0) 
00431       phat.Scale(1.0, p);
00432     else
00433       M->Solve(false, p, phat);
00434     A.Multiply(false, phat, v);
00435 
00436     
00437     rtilde.Dot(v,&sigma);
00438     alpha = rhon/sigma;    
00439 
00440     /* s = r - alpha*v                     */
00441     /* shat = M^-1 s                       */
00442     /* r = A shat (r is a tmp here for t ) */
00443 
00444     s.Update(-alpha, v, 1.0, r, 0.0);
00445     if (M==0) 
00446       shat.Scale(1.0, s);
00447     else
00448       M->Solve(false, s, shat);
00449     A.Multiply(false, shat, r);
00450 
00451     r.Dot(s, &omega_num);
00452     r.Dot(r, &omega_den);
00453     omega = omega_num / omega_den;
00454 
00455     /* x = x + alpha*phat + omega*shat */
00456     /* r = s - omega*r */
00457 
00458     x.Update(alpha, phat, omega, shat, 1.0);
00459     r.Update(1.0, s, -omega); 
00460     
00461     r.Norm2(&r_norm);
00462     scaled_r_norm = r_norm/b_norm;
00463     r.Dot(rtilde,&rhon);
00464 
00465     if (verbose) cout << "Iter "<< i << " residual = " << r_norm
00466           << " Scaled residual = " << scaled_r_norm << endl;
00467 
00468     if (r_norm < Tolerance) break;
00469   }
00470   return;
00471 }
00472 //==============================================================================
00473 int Statistics(const Epetra_CrsSingletonFilter & filter) {
00474 
00475   // Create local variables of some of the stats we need because we don't know if the
00476   // method calls are collective or not for the particular implementation of the Row Matrix
00477   int fmaxentries;
00478   int maxentries = filter.FullMatrix()->MaxNumEntries();
00479   filter.FullMatrix()->RowMatrixRowMap().Comm().SumAll(&maxentries, &fmaxentries, 1);
00480   int fnrows = filter.FullMatrix()->NumGlobalRows();
00481   int fnnz = filter.FullMatrix()->NumGlobalNonzeros();
00482   if (filter.FullMatrix()->RowMatrixRowMap().Comm().MyPID()!=0) return(0);
00483 
00484   cout << "Full System characteristics:" << endl << endl
00485        << "  Dimension                             = " << fnrows << endl
00486        << "  Number of nonzeros                    = " <<fnnz << endl
00487        << "  Maximum Number of Row Entries         = " << fmaxentries << endl << endl
00488        << "Reduced System characteristics:" << endl << endl
00489        << "  Dimension                             = " << filter.ReducedMatrix()->NumGlobalRows() << endl
00490        << "  Number of nonzeros                    = " << filter.ReducedMatrix()->NumGlobalNonzeros() << endl
00491        << "  Maximum Number of Row Entries         = " << filter.ReducedMatrix()->GlobalMaxNumEntries() << endl << endl
00492        << "Singleton information: " << endl
00493        << "  Total number of singletons            = " << filter.NumSingletons() << endl
00494        << "  Number of rows with single entries    = " << filter.NumRowSingletons() << endl
00495        << "  Number of columns with single entries " << endl
00496        << "    (that were not already counted as " << endl
00497        << "     row singletons)                    = " << filter.NumColSingletons() << endl << endl
00498        << "Ratios: " << endl
00499        << "  Percent reduction in dimension        = " << filter.RatioOfDimensions()*100.0 << endl
00500        << "  Percent reduction in nonzero count    = " << filter.RatioOfNonzeros()*100.0 << endl << endl;
00501 
00502   return(0);
00503 
00504 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 09:58:40 2011 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.6.3