test/PointPreconditioner/cxx_main.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                IFPACK
00005 //                 Copyright (2004) 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 "Ifpack_ConfigDefs.h"
00030 #ifdef HAVE_IFPACK_AZTECOO 
00031 #ifdef HAVE_MPI
00032 #include "Epetra_MpiComm.h"
00033 #else
00034 #include "Epetra_SerialComm.h"
00035 #endif
00036 #include "Epetra_CrsMatrix.h"
00037 #include "Epetra_Vector.h"
00038 #include "Epetra_LinearProblem.h"
00039 #include "Trilinos_Util_CrsMatrixGallery.h"
00040 #include "Teuchos_ParameterList.hpp"
00041 #include "Ifpack_PointRelaxation.h"
00042 #include "Ifpack_BlockRelaxation.h"
00043 #include "Ifpack_DenseContainer.h"
00044 #include "Ifpack_AdditiveSchwarz.h"
00045 #include "AztecOO.h"
00046 
00047 using namespace Trilinos_Util;
00048 static bool verbose = false;
00049 static bool SymmetricGallery = false;
00050 static bool Solver = AZ_gmres;
00051 
00052 
00053 // ====================================================================== 
00054 int CompareBlockOverlap(CrsMatrixGallery& Gallery, int Overlap)
00055 {
00056 
00057   Epetra_RowMatrix* A = Gallery.GetMatrix();
00058   Epetra_LinearProblem* Problem = Gallery.GetLinearProblem();
00059 
00060   Epetra_MultiVector& RHS = *(Problem->GetRHS());
00061   Epetra_MultiVector& LHS = *(Problem->GetLHS());
00062 
00063   Teuchos::ParameterList List;
00064   List.set("relaxation: damping factor", 1.0);
00065   List.set("relaxation: type", "Jacobi");
00066   List.set("relaxation: sweeps",1);
00067   List.set("partitioner: overlap", Overlap);
00068   List.set("partitioner: type", "greedy");
00069   List.set("partitioner: local parts", 16);
00070 
00071   RHS.PutScalar(1.0);
00072   LHS.PutScalar(0.0);
00073 
00074   Ifpack_AdditiveSchwarz<Ifpack_BlockRelaxation<Ifpack_DenseContainer> > Prec(A);
00075   Prec.SetParameters(List);
00076   Prec.Compute();
00077 
00078   // set AztecOO solver object
00079   AztecOO AztecOOSolver(*Problem);
00080   AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00081   if (verbose)
00082     AztecOOSolver.SetAztecOption(AZ_output,32);
00083   else
00084     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00085   AztecOOSolver.SetPrecOperator(&Prec);
00086 
00087   AztecOOSolver.Iterate(1550,1e-8);
00088 
00089   return(AztecOOSolver.NumIters());
00090 }
00091 
00092 // ====================================================================== 
00093 int CompareBlockSizes(string PrecType,
00094                       CrsMatrixGallery& Gallery, int NumParts)
00095 {
00096 
00097   Epetra_RowMatrix* A = Gallery.GetMatrix();
00098   Epetra_LinearProblem* Problem = Gallery.GetLinearProblem();
00099 
00100   Epetra_MultiVector& RHS = *(Problem->GetRHS());
00101   Epetra_MultiVector& LHS = *(Problem->GetLHS());
00102 
00103   Teuchos::ParameterList List;
00104   List.set("relaxation: damping factor", 1.0);
00105   List.set("relaxation: type", PrecType);
00106   List.set("relaxation: sweeps",1);
00107   List.set("partitioner: type", "greedy");
00108   List.set("partitioner: local parts", NumParts);
00109 
00110   RHS.PutScalar(1.0);
00111   LHS.PutScalar(0.0);
00112 
00113   Ifpack_AdditiveSchwarz<Ifpack_BlockRelaxation<Ifpack_DenseContainer> > Prec(A);
00114   Prec.SetParameters(List);
00115   Prec.Compute();
00116 
00117   // set AztecOO solver object
00118   AztecOO AztecOOSolver(*Problem);
00119   AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00120   if (verbose)
00121     AztecOOSolver.SetAztecOption(AZ_output,32);
00122   else
00123     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00124   AztecOOSolver.SetPrecOperator(&Prec);
00125 
00126   AztecOOSolver.Iterate(1550,1e-8);
00127 
00128   return(AztecOOSolver.NumIters());
00129 }
00130 
00131 // ====================================================================== 
00132 bool ComparePointAndBlock(string PrecType,
00133                           CrsMatrixGallery& Gallery, int sweeps)
00134 {
00135 
00136   Epetra_RowMatrix* A = Gallery.GetMatrix();
00137   Epetra_LinearProblem* Problem = Gallery.GetLinearProblem();
00138 
00139   Epetra_MultiVector& RHS = *(Problem->GetRHS());
00140   Epetra_MultiVector& LHS = *(Problem->GetLHS());
00141 
00142   // Set up the list
00143   Teuchos::ParameterList List;
00144   List.set("relaxation: damping factor", 1.0);
00145   List.set("relaxation: type", PrecType);
00146   List.set("relaxation: sweeps",sweeps);
00147   List.set("partitioner: type", "linear");
00148   List.set("partitioner: local parts", A->NumMyRows());
00149 
00150   int ItersPoint, ItersBlock;
00151 
00152   // ================================================== //
00153   // get the number of iterations with point relaxation //
00154   // ================================================== //
00155   {
00156 
00157     RHS.PutScalar(1.0);
00158     LHS.PutScalar(0.0);
00159 
00160     Ifpack_AdditiveSchwarz<Ifpack_PointRelaxation> Point(A);
00161     Point.SetParameters(List);
00162     Point.Compute();
00163 
00164     // set AztecOO solver object
00165     AztecOO AztecOOSolver(*Problem);
00166     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00167     if (verbose)
00168       AztecOOSolver.SetAztecOption(AZ_output,32);
00169     else
00170       AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00171     AztecOOSolver.SetPrecOperator(&Point);
00172 
00173     AztecOOSolver.Iterate(1550,1e-8);
00174 
00175     double TrueResidual = AztecOOSolver.TrueResidual();
00176     ItersPoint = AztecOOSolver.NumIters();
00177     // some output
00178     if (verbose && Problem->GetMatrix()->Comm().MyPID() == 0) {
00179       cout << "Iterations  = " << ItersPoint << endl;
00180       cout << "Norm of the true residual = " << TrueResidual << endl;
00181     }
00182   }
00183 
00184   // ================================================== //
00185   // get the number of iterations with block relaxation //
00186   // ================================================== //
00187   {
00188 
00189     RHS.PutScalar(1.0);
00190     LHS.PutScalar(0.0);
00191 
00192     Ifpack_AdditiveSchwarz<Ifpack_BlockRelaxation<Ifpack_DenseContainer> > Block(A);
00193     Block.SetParameters(List);
00194     Block.Compute();
00195 
00196     // set AztecOO solver object
00197     AztecOO AztecOOSolver(*Problem);
00198     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00199     if (verbose)
00200       AztecOOSolver.SetAztecOption(AZ_output,32);
00201     else
00202       AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00203     AztecOOSolver.SetPrecOperator(&Block);
00204 
00205     AztecOOSolver.Iterate(1550,1e-8);
00206 
00207     double TrueResidual = AztecOOSolver.TrueResidual();
00208     ItersBlock = AztecOOSolver.NumIters();
00209     // some output
00210     if (verbose && Problem->GetMatrix()->Comm().MyPID() == 0) {
00211       cout << "Iterations " << ItersBlock << endl;
00212       cout << "Norm of the true residual = " << TrueResidual << endl;
00213     }
00214   }
00215 
00216   if (ItersPoint != ItersBlock) {
00217     if (verbose)
00218       cout << "TEST FAILED!" << endl;
00219     return(false);
00220   }
00221   else {
00222     if (verbose)
00223       cout << "TEST PASSED" << endl;
00224     return(true);
00225   }
00226 }
00227 
00228 // ====================================================================== 
00229 bool KrylovTest(string PrecType,
00230                 CrsMatrixGallery& Gallery)
00231 {
00232 
00233   // The following methods of CrsMatrixGallery are used to get pointers
00234   // to internally stored Epetra_RowMatrix and Epetra_LinearProblem.
00235   Epetra_RowMatrix* A = Gallery.GetMatrix();
00236   Epetra_LinearProblem* Problem = Gallery.GetLinearProblem();
00237 
00238   Epetra_MultiVector& LHS = *(Problem->GetLHS());
00239   LHS.PutScalar(0.0);
00240 
00241   // Set up the list
00242   Teuchos::ParameterList List;
00243   List.set("relaxation: damping factor", 1.0);
00244   List.set("relaxation: type", PrecType);
00245 
00246   int Iters1, Iters10;
00247 
00248   // ============================================== //
00249   // get the number of iterations with 1 sweep only //
00250   // ============================================== //
00251   {
00252 
00253     List.set("relaxation: sweeps",1);
00254     Ifpack_AdditiveSchwarz<Ifpack_PointRelaxation> Point(A);
00255     Point.SetParameters(List);
00256     Point.Compute();
00257 
00258     // set AztecOO solver object
00259     AztecOO AztecOOSolver(*Problem);
00260     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00261     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00262     AztecOOSolver.SetPrecOperator(&Point);
00263 
00264     AztecOOSolver.Iterate(1550,1e-8);
00265 
00266     double TrueResidual = AztecOOSolver.TrueResidual();
00267     // some output
00268     if (verbose && Problem->GetMatrix()->Comm().MyPID() == 0) {
00269       cout << "Norm of the true residual = " << TrueResidual << endl;
00270     }
00271     Iters1 = AztecOOSolver.NumIters();
00272   }
00273 
00274   // ======================================================== //
00275   // now re-run with 10 sweeps, solver should converge faster
00276   // ======================================================== //
00277   {
00278     List.set("relaxation: sweeps",10);
00279     Ifpack_AdditiveSchwarz<Ifpack_PointRelaxation> Point(A);
00280     Point.SetParameters(List);
00281     Point.Compute();
00282     LHS.PutScalar(0.0);
00283 
00284     // set AztecOO solver object
00285     AztecOO AztecOOSolver(*Problem);
00286     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00287     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00288     AztecOOSolver.SetPrecOperator(&Point);
00289     AztecOOSolver.Iterate(1550,1e-8);
00290 
00291     double TrueResidual = AztecOOSolver.TrueResidual();
00292     // some output
00293     if (verbose && Problem->GetMatrix()->Comm().MyPID() == 0) {
00294       cout << "Norm of the true residual = " << TrueResidual << endl;
00295     }
00296     Iters10 = AztecOOSolver.NumIters();
00297   }
00298 
00299   if (verbose) {
00300     cout << "Iters_1 = " << Iters1 << ", Iters_10 = " << Iters10 << endl;
00301     cout << "(second number should be smaller than first one)" << endl;
00302   }
00303 
00304   if (Iters10 > Iters1) {
00305     if (verbose)
00306       cout << "TEST FAILED!" << endl;
00307     return(false);
00308   }
00309   else {
00310     if (verbose)
00311       cout << "TEST PASSED" << endl;
00312     return(true);
00313   }
00314 }
00315 
00316 // ====================================================================== 
00317 bool BasicTest(string PrecType,
00318                CrsMatrixGallery& Gallery)
00319 {
00320 
00321   // The following methods of CrsMatrixGallery are used to get pointers
00322   // to internally stored Epetra_RowMatrix and Epetra_LinearProblem.
00323   Epetra_RowMatrix* A = Gallery.GetMatrix();
00324   Epetra_LinearProblem* Problem = Gallery.GetLinearProblem();
00325 
00326   Epetra_MultiVector& RHS = *(Problem->GetRHS());
00327   Epetra_MultiVector& LHS = *(Problem->GetLHS());
00328 
00329   LHS.PutScalar(0.0);
00330 
00331   // Set up the list
00332   Teuchos::ParameterList List;
00333   List.set("relaxation: damping factor", 1.0);
00334   List.set("relaxation: sweeps",1550);
00335   List.set("relaxation: type", PrecType);
00336 
00337   Ifpack_PointRelaxation Point(A);
00338 
00339   Point.SetParameters(List);
00340   Point.Compute();
00341   // use the preconditioner as solver, with 1550 iterations
00342   Point.ApplyInverse(RHS,LHS);
00343 
00344   // compute the real residual
00345 
00346   double residual, diff;
00347   Gallery.ComputeResidual(&residual);
00348   Gallery.ComputeDiffBetweenStartingAndExactSolutions(&diff);
00349 
00350   if (verbose && A->Comm().MyPID()==0) {
00351     cout << "||b-Ax||_2 = " << residual << endl;
00352     cout << "||x_exact - x||_2 = " << diff << endl;
00353   }
00354   
00355   // Jacobi is very slow to converge here
00356   if (residual < 1e-2) {
00357     if (verbose)
00358       cout << "Test passed" << endl;
00359     return(true);
00360   }
00361   else {
00362     if (verbose)
00363       cout << "Test failed!" << endl;
00364     return(false);
00365   }
00366 }
00367 
00368 // ====================================================================== 
00369 int main(int argc, char *argv[])
00370 {
00371 
00372 #ifdef HAVE_MPI
00373   MPI_Init(&argc,&argv);
00374   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00375 #else
00376   Epetra_SerialComm Comm;
00377 #endif
00378 
00379   bool verbose = (Comm.MyPID() == 0);
00380 
00381   for (int i = 1 ; i < argc ; ++i) {
00382     if (strcmp(argv[i],"-s") == 0) {
00383       SymmetricGallery = true;
00384       Solver = AZ_cg;
00385     }
00386   }
00387 
00388   // size of the global matrix. 
00389   const int NumPoints = 900;
00390 
00391   CrsMatrixGallery Gallery("", Comm);
00392   Gallery.Set("problem_size", NumPoints);
00393   Gallery.Set("map_type", "linear");
00394   if (SymmetricGallery)
00395     Gallery.Set("problem_type", "laplace_2d");
00396   else
00397     Gallery.Set("problem_type", "recirc_2d");
00398 
00399   // test the preconditioner
00400   int TestPassed = true;
00401 
00402   // ======================================== //
00403   // first verify that we can get convergence //
00404   // with all point relaxation methods        //
00405   // ======================================== //
00406   if(!BasicTest("Jacobi",Gallery))
00407     TestPassed = false;
00408 
00409   if(!BasicTest("symmetric Gauss-Seidel",Gallery))
00410     TestPassed = false;
00411 
00412   if (!SymmetricGallery) {
00413     if(!BasicTest("Gauss-Seidel",Gallery))
00414       TestPassed = false;
00415   }
00416 
00417   // ============================= //
00418   // check uses as preconditioners //
00419   // ============================= //
00420   
00421   if(!KrylovTest("symmetric Gauss-Seidel",Gallery))
00422     TestPassed = false;
00423 
00424   if (!SymmetricGallery) {
00425     if(!KrylovTest("Gauss-Seidel",Gallery))
00426       TestPassed = false;
00427   }
00428 
00429   // ================================== //
00430   // compare point and block relaxation //
00431   // ================================== //
00432 
00433   TestPassed = TestPassed && 
00434     ComparePointAndBlock("Jacobi",Gallery,1);
00435 
00436   TestPassed = TestPassed && 
00437     ComparePointAndBlock("Jacobi",Gallery,10);
00438 
00439   TestPassed = TestPassed && 
00440     ComparePointAndBlock("symmetric Gauss-Seidel",Gallery,1);
00441 
00442   TestPassed = TestPassed && 
00443     ComparePointAndBlock("symmetric Gauss-Seidel",Gallery,10);
00444 
00445   if (!SymmetricGallery) {
00446     TestPassed = TestPassed && 
00447       ComparePointAndBlock("Gauss-Seidel",Gallery,1);
00448 
00449     TestPassed = TestPassed && 
00450       ComparePointAndBlock("Gauss-Seidel",Gallery,10);
00451   }
00452 
00453   // ============================ //
00454   // verify effect of # of blocks //
00455   // ============================ //
00456   
00457   {
00458     int Iters4, Iters8, Iters16;
00459     Iters4 = CompareBlockSizes("Jacobi",Gallery,4);
00460     Iters8 = CompareBlockSizes("Jacobi",Gallery,8);
00461     Iters16 = CompareBlockSizes("Jacobi",Gallery,16);
00462     if ((Iters16 > Iters8) && (Iters8 > Iters4)) {
00463       if (verbose)
00464         cout << "Test passed" << endl;
00465     }
00466     else {
00467       if (verbose) 
00468         cout << "TEST FAILED!" << endl;
00469       TestPassed = TestPassed && false;
00470     }
00471   }
00472 
00473   // ================================== //
00474   // verify effect of overlap in Jacobi //
00475   // ================================== //
00476   
00477   {
00478     int Iters0, Iters2, Iters4;
00479     Iters0 = CompareBlockOverlap(Gallery,0);
00480     Iters2 = CompareBlockOverlap(Gallery,2);
00481     Iters4 = CompareBlockOverlap(Gallery,4);
00482     if ((Iters4 < Iters2) && (Iters2 < Iters0)) {
00483       if (verbose)
00484         cout << "Test passed" << endl;
00485     }
00486     else {
00487       if (verbose) 
00488         cout << "TEST FAILED!" << endl;
00489       TestPassed = TestPassed && false;
00490     }
00491   }
00492 
00493   // ============ //
00494   // final output //
00495   // ============ //
00496 
00497   if (!TestPassed) {
00498     cout << "Test `PointPreconditioner.exe' FAILED!" << endl;
00499     exit(EXIT_FAILURE);
00500   }
00501   
00502 #ifdef HAVE_MPI
00503   MPI_Finalize(); 
00504 #endif
00505 
00506   cout << "Test `PointPreconditioner.exe' passed!" << endl;
00507   exit(EXIT_SUCCESS);
00508 }
00509 
00510 #else
00511 
00512 #ifdef HAVE_MPI
00513 #include "Epetra_MpiComm.h"
00514 #else
00515 #include "Epetra_SerialComm.h"
00516 #endif
00517 
00518 int main(int argc, char *argv[])
00519 {
00520 
00521 #ifdef HAVE_MPI
00522   MPI_Init(&argc,&argv);
00523   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00524 #else
00525   Epetra_SerialComm Comm;
00526 #endif
00527 
00528   puts("please configure IFPACK with --eanble-aztecoo --enable-teuchos");
00529   puts("to run this test");
00530 
00531 #ifdef HAVE_MPI
00532   MPI_Finalize() ;
00533 #endif
00534   return(EXIT_SUCCESS);
00535 }
00536 
00537 #endif
 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