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

Generated on Thu Sep 18 12:37:21 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1