Ifpack Package Browser (Single Doxygen Collection) Development
test/Relaxation/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 
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 "Epetra_Map.h"
00040 #include "Galeri_Maps.h"
00041 #include "Galeri_CrsMatrices.h"
00042 #include "Galeri_Utils.h"
00043 #include "Teuchos_ParameterList.hpp"
00044 #include "Teuchos_RefCountPtr.hpp"
00045 #include "Ifpack_PointRelaxation.h"
00046 #include "Ifpack_BlockRelaxation.h"
00047 #include "Ifpack_SparseContainer.h"
00048 #include "Ifpack_Amesos.h"
00049 #include "AztecOO.h"
00050 
00051 static bool verbose = false;
00052 static bool SymmetricGallery = false;
00053 static bool Solver = AZ_gmres;
00054 const int NumVectors = 3;
00055 
00056 // ====================================================================== 
00057 int CompareBlockOverlap(const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int Overlap)
00058 {
00059   Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
00060   Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
00061   LHS.PutScalar(0.0); RHS.Random();
00062 
00063   Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
00064 
00065   Teuchos::ParameterList List;
00066   List.set("relaxation: damping factor", 1.0);
00067   List.set("relaxation: type", "Jacobi");
00068   List.set("relaxation: sweeps",1);
00069   List.set("partitioner: overlap", Overlap);
00070   List.set("partitioner: type", "linear");
00071   List.set("partitioner: local parts", 16);
00072 
00073   RHS.PutScalar(1.0);
00074   LHS.PutScalar(0.0);
00075 
00076   Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
00077   Prec.SetParameters(List);
00078   Prec.Compute();
00079 
00080   // set AztecOO solver object
00081   AztecOO AztecOOSolver(Problem);
00082   AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00083   if (verbose)
00084     AztecOOSolver.SetAztecOption(AZ_output,32);
00085   else
00086     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00087   AztecOOSolver.SetPrecOperator(&Prec);
00088 
00089   AztecOOSolver.Iterate(1550,1e-5);
00090 
00091   return(AztecOOSolver.NumIters());
00092 }
00093 
00094 // ====================================================================== 
00095 int CompareBlockSizes(string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int NumParts)
00096 {
00097   Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
00098   Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
00099   LHS.PutScalar(0.0); RHS.Random();
00100 
00101   Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
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", "linear");
00108   List.set("partitioner: local parts", NumParts);
00109 
00110   RHS.PutScalar(1.0);
00111   LHS.PutScalar(0.0);
00112 
00113   Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > 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-5);
00127 
00128   return(AztecOOSolver.NumIters());
00129 }
00130 
00131 // ====================================================================== 
00132 bool ComparePointAndBlock(string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int sweeps)
00133 {
00134   Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
00135   Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
00136   LHS.PutScalar(0.0); RHS.Random();
00137 
00138   Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
00139 
00140   // Set up the list
00141   Teuchos::ParameterList List;
00142   List.set("relaxation: damping factor", 1.0);
00143   List.set("relaxation: type", PrecType);
00144   List.set("relaxation: sweeps",sweeps);
00145   List.set("partitioner: type", "linear");
00146   List.set("partitioner: local parts", A->NumMyRows());
00147 
00148   int ItersPoint, ItersBlock;
00149 
00150   // ================================================== //
00151   // get the number of iterations with point relaxation //
00152   // ================================================== //
00153   {
00154     RHS.PutScalar(1.0);
00155     LHS.PutScalar(0.0);
00156 
00157     Ifpack_PointRelaxation Point(&*A);
00158     Point.SetParameters(List);
00159     Point.Compute();
00160 
00161     // set AztecOO solver object
00162     AztecOO AztecOOSolver(Problem);
00163     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00164     if (verbose)
00165       AztecOOSolver.SetAztecOption(AZ_output,32);
00166     else
00167       AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00168     AztecOOSolver.SetPrecOperator(&Point);
00169 
00170     AztecOOSolver.Iterate(1550,1e-2);
00171 
00172     double TrueResidual = AztecOOSolver.TrueResidual();
00173     ItersPoint = AztecOOSolver.NumIters();
00174     // some output
00175     if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
00176       cout << "Iterations  = " << ItersPoint << endl;
00177       cout << "Norm of the true residual = " << TrueResidual << endl;
00178     }
00179   }
00180 
00181   // ================================================== //
00182   // get the number of iterations with block relaxation //
00183   // ================================================== //
00184   {
00185 
00186     RHS.PutScalar(1.0);
00187     LHS.PutScalar(0.0);
00188 
00189     Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Block(&*A);
00190     Block.SetParameters(List);
00191     Block.Compute();
00192 
00193     // set AztecOO solver object
00194     AztecOO AztecOOSolver(Problem);
00195     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00196     if (verbose)
00197       AztecOOSolver.SetAztecOption(AZ_output,32);
00198     else
00199       AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00200     AztecOOSolver.SetPrecOperator(&Block);
00201 
00202     AztecOOSolver.Iterate(1550,1e-2);
00203 
00204     double TrueResidual = AztecOOSolver.TrueResidual();
00205     ItersBlock = AztecOOSolver.NumIters();
00206     // some output
00207     if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
00208       cout << "Iterations " << ItersBlock << endl;
00209       cout << "Norm of the true residual = " << TrueResidual << endl;
00210     }
00211   }
00212 
00213   int diff = ItersPoint - ItersBlock;
00214   if (diff < 0) diff = -diff;
00215     
00216   if (diff > 10)
00217   {
00218     if (verbose)
00219       cout << "TEST FAILED!" << endl;
00220     return(false);
00221   }
00222   else {
00223     if (verbose)
00224       cout << "TEST PASSED" << endl;
00225     return(true);
00226   }
00227 }
00228 
00229 // ====================================================================== 
00230 bool KrylovTest(string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, bool backward)
00231 {
00232   Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
00233   Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
00234   LHS.PutScalar(0.0); RHS.Random();
00235 
00236   Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
00237 
00238   // Set up the list
00239   Teuchos::ParameterList List;
00240   List.set("relaxation: damping factor", 1.0);
00241   List.set("relaxation: type", PrecType);
00242   if(backward) List.set("relaxation: backward mode",backward);  
00243 
00244   int Iters1, Iters10;
00245 
00246   if (verbose) {
00247     cout << "Krylov test: Using " << PrecType 
00248          << " with AztecOO" << endl;
00249   }
00250 
00251   // ============================================== //
00252   // get the number of iterations with 1 sweep only //
00253   // ============================================== //
00254   {
00255 
00256     List.set("relaxation: sweeps",1);
00257     Ifpack_PointRelaxation Point(&*A);
00258     Point.SetParameters(List);
00259     Point.Compute();
00260 
00261     // set AztecOO solver object
00262     AztecOO AztecOOSolver(Problem);
00263     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00264     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00265     AztecOOSolver.SetPrecOperator(&Point);
00266 
00267     AztecOOSolver.Iterate(1550,1e-5);
00268 
00269     double TrueResidual = AztecOOSolver.TrueResidual();
00270     // some output
00271     if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
00272       cout << "Norm of the true residual = " << TrueResidual << endl;
00273     }
00274     Iters1 = AztecOOSolver.NumIters();
00275   }
00276 
00277   // ======================================================== //
00278   // now re-run with 10 sweeps, solver should converge faster
00279   // ======================================================== //
00280   {
00281     List.set("relaxation: sweeps",10);
00282     Ifpack_PointRelaxation Point(&*A);
00283     Point.SetParameters(List);
00284     Point.Compute();
00285     LHS.PutScalar(0.0);
00286 
00287     // set AztecOO solver object
00288     AztecOO AztecOOSolver(Problem);
00289     AztecOOSolver.SetAztecOption(AZ_solver,Solver);
00290     AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
00291     AztecOOSolver.SetPrecOperator(&Point);
00292     AztecOOSolver.Iterate(1550,1e-5);
00293 
00294     double TrueResidual = AztecOOSolver.TrueResidual();
00295     // some output
00296     if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
00297       cout << "Norm of the true residual = " << TrueResidual << endl;
00298     }
00299     Iters10 = AztecOOSolver.NumIters();
00300   }
00301 
00302   if (verbose) {
00303     cout << "Iters_1 = " << Iters1 << ", Iters_10 = " << Iters10 << endl;
00304     cout << "(second number should be smaller than first one)" << endl;
00305   }
00306 
00307   if (Iters10 > Iters1) {
00308     if (verbose)
00309       cout << "TEST FAILED!" << endl;
00310     return(false);
00311   }
00312   else {
00313     if (verbose)
00314       cout << "TEST PASSED" << endl;
00315     return(true);
00316   }
00317 }
00318 
00319 // ====================================================================== 
00320 bool BasicTest(string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,bool backward)
00321 {
00322   Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
00323   Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
00324   LHS.PutScalar(0.0); RHS.Random();
00325 
00326   double starting_residual = Galeri::ComputeNorm(&*A, &LHS, &RHS);
00327   Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
00328 
00329   // Set up the list
00330   Teuchos::ParameterList List;
00331   List.set("relaxation: damping factor", 1.0);
00332   List.set("relaxation: sweeps",1550);
00333   List.set("relaxation: type", PrecType);
00334   if(backward) List.set("relaxation: backward mode",backward);
00335   
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 = Galeri::ComputeNorm(&*A, &LHS, &RHS);
00347   
00348   if (A->Comm().MyPID() == 0 && verbose)
00349     cout << "||A * x - b||_2 (scaled) = " << residual / starting_residual << endl;
00350   
00351   // Jacobi is very slow to converge here
00352   if (residual / starting_residual < 1e-2) {
00353     if (verbose)
00354       cout << "Test passed" << endl;
00355     return(true);
00356   }
00357   else {
00358     if (verbose)
00359       cout << "Test failed!" << endl;
00360     return(false);
00361   }
00362 }
00363 
00364 // ====================================================================== 
00365 int main(int argc, char *argv[])
00366 {
00367 #ifdef HAVE_MPI
00368   MPI_Init(&argc,&argv);
00369   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00370 #else
00371   Epetra_SerialComm Comm;
00372 #endif
00373 
00374   verbose = (Comm.MyPID() == 0);
00375 
00376   for (int i = 1 ; i < argc ; ++i) {
00377     if (strcmp(argv[i],"-s") == 0) {
00378       SymmetricGallery = true;
00379       Solver = AZ_cg;
00380     }
00381   }
00382 
00383   // size of the global matrix. 
00384   Teuchos::ParameterList GaleriList;
00385   int nx = 30; 
00386   GaleriList.set("nx", nx);
00387   GaleriList.set("ny", nx * Comm.NumProc());
00388   GaleriList.set("mx", 1);
00389   GaleriList.set("my", Comm.NumProc());
00390   Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
00391   Teuchos::RefCountPtr<Epetra_CrsMatrix> A;
00392   if (SymmetricGallery)
00393     A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
00394   else
00395     A = Teuchos::rcp( Galeri::CreateCrsMatrix("Recirc2D", &*Map, GaleriList) );
00396 
00397   // test the preconditioner
00398   int TestPassed = true;
00399 
00400   // ======================================== //
00401   // first verify that we can get convergence //
00402   // with all point relaxation methods        //
00403   // ======================================== //
00404 
00405   if(!BasicTest("Jacobi",A,false))
00406     TestPassed = false;
00407 
00408   if(!BasicTest("symmetric Gauss-Seidel",A,false))
00409     TestPassed = false;
00410 
00411   if (!SymmetricGallery) {
00412     if(!BasicTest("Gauss-Seidel",A,false))
00413       TestPassed = false;
00414     if(!BasicTest("Gauss-Seidel",A,true))
00415       TestPassed = false;  
00416 
00417   }
00418 
00419   // ============================= //
00420   // check uses as preconditioners //
00421   // ============================= //
00422   
00423   if(!KrylovTest("symmetric Gauss-Seidel",A,false))
00424     TestPassed = false;
00425 
00426   if (!SymmetricGallery) {
00427     if(!KrylovTest("Gauss-Seidel",A,false))
00428       TestPassed = false;
00429     if(!KrylovTest("Gauss-Seidel",A,true))
00430       TestPassed = false;
00431 
00432   }
00433 
00434   // ================================== //
00435   // compare point and block relaxation //
00436   // ================================== //
00437 
00438   //TestPassed = TestPassed && 
00439    // ComparePointAndBlock("Jacobi",A,1);
00440 
00441   TestPassed = TestPassed && 
00442     ComparePointAndBlock("Jacobi",A,10);
00443 
00444   //TestPassed = TestPassed && 
00445     //ComparePointAndBlock("symmetric Gauss-Seidel",A,1);
00446 
00447   TestPassed = TestPassed && 
00448     ComparePointAndBlock("symmetric Gauss-Seidel",A,10);
00449 
00450   if (!SymmetricGallery) {
00451     //TestPassed = TestPassed && 
00452       //ComparePointAndBlock("Gauss-Seidel",A,1);
00453 
00454     TestPassed = TestPassed && 
00455       ComparePointAndBlock("Gauss-Seidel",A,10);
00456   }
00457 
00458   // ============================ //
00459   // verify effect of # of blocks //
00460   // ============================ //
00461   
00462   {
00463     int Iters4, Iters8, Iters16;
00464     Iters4 = CompareBlockSizes("Jacobi",A,4);
00465     Iters8 = CompareBlockSizes("Jacobi",A,8);
00466     Iters16 = CompareBlockSizes("Jacobi",A,16);
00467 
00468     if ((Iters16 > Iters8) && (Iters8 > Iters4)) {
00469       if (verbose)
00470         cout << "Test passed" << endl;
00471     }
00472     else {
00473       if (verbose) 
00474         cout << "TEST FAILED!" << endl;
00475       TestPassed = TestPassed && false;
00476     }
00477   }
00478 
00479   // ================================== //
00480   // verify effect of overlap in Jacobi //
00481   // ================================== //
00482 
00483   {
00484     int Iters0, Iters2, Iters4;
00485     Iters0 = CompareBlockOverlap(A,0);
00486     Iters2 = CompareBlockOverlap(A,2);
00487     Iters4 = CompareBlockOverlap(A,4);
00488     if ((Iters4 < Iters2) && (Iters2 < Iters0)) {
00489       if (verbose)
00490         cout << "Test passed" << endl;
00491     }
00492     else {
00493       if (verbose) 
00494         cout << "TEST FAILED!" << endl;
00495       TestPassed = TestPassed && false;
00496     }
00497   }
00498 
00499   // ============ //
00500   // final output //
00501   // ============ //
00502 
00503   if (!TestPassed) {
00504     cout << "Test `TestRelaxation.exe' failed!" << endl;
00505     exit(EXIT_FAILURE);
00506   }
00507   
00508 #ifdef HAVE_MPI
00509   MPI_Finalize(); 
00510 #endif
00511   
00512   cout << endl;
00513   cout << "Test `TestRelaxation.exe' passed!" << endl;
00514   cout << endl;
00515   return(EXIT_SUCCESS);
00516 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines