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