46 #include "Epetra_MpiComm.h" 48 #include "Epetra_SerialComm.h" 50 #include "Epetra_CrsMatrix.h" 51 #include "Epetra_Vector.h" 52 #include "Epetra_LinearProblem.h" 53 #include "Epetra_Map.h" 54 #include "Galeri_Maps.h" 55 #include "Galeri_CrsMatrices.h" 56 #include "Galeri_Utils.h" 57 #include "Teuchos_ParameterList.hpp" 58 #include "Teuchos_RefCountPtr.hpp" 75 LHS.PutScalar(0.0);
RHS.Random();
79 Teuchos::ParameterList List;
80 List.set(
"relaxation: damping factor", 1.0);
81 List.set(
"relaxation: type",
"Jacobi");
82 List.set(
"relaxation: sweeps",1);
83 List.set(
"partitioner: overlap", Overlap);
84 List.set(
"partitioner: type",
"linear");
85 List.set(
"partitioner: local parts", 16);
95 AztecOO AztecOOSolver(Problem);
96 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
98 AztecOOSolver.SetAztecOption(AZ_output,32);
100 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
101 AztecOOSolver.SetPrecOperator(&Prec);
103 AztecOOSolver.Iterate(1550,1e-5);
105 return(AztecOOSolver.NumIters());
109 int CompareBlockSizes(std::string PrecType,
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,
int NumParts)
117 Teuchos::ParameterList List;
118 List.set(
"relaxation: damping factor", 1.0);
119 List.set(
"relaxation: type", PrecType);
120 List.set(
"relaxation: sweeps",1);
121 List.set(
"partitioner: type",
"linear");
122 List.set(
"partitioner: local parts", NumParts);
132 AztecOO AztecOOSolver(Problem);
133 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
135 AztecOOSolver.SetAztecOption(AZ_output,32);
137 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
138 AztecOOSolver.SetPrecOperator(&Prec);
140 AztecOOSolver.Iterate(1550,1e-5);
142 return(AztecOOSolver.NumIters());
155 Teuchos::ParameterList List;
156 List.set(
"relaxation: damping factor", 1.0);
157 List.set(
"relaxation: type", PrecType);
158 List.set(
"relaxation: sweeps",sweeps);
159 List.set(
"partitioner: type",
"linear");
160 List.set(
"partitioner: local parts", A->NumMyRows());
162 int ItersPoint, ItersBlock;
176 AztecOO AztecOOSolver(Problem);
177 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
179 AztecOOSolver.SetAztecOption(AZ_output,32);
181 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
182 AztecOOSolver.SetPrecOperator(&Point);
184 AztecOOSolver.Iterate(1550,1e-2);
186 double TrueResidual = AztecOOSolver.TrueResidual();
187 ItersPoint = AztecOOSolver.NumIters();
190 cout <<
"Iterations = " << ItersPoint << endl;
191 cout <<
"Norm of the true residual = " << TrueResidual << endl;
208 AztecOO AztecOOSolver(Problem);
209 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
211 AztecOOSolver.SetAztecOption(AZ_output,32);
213 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
214 AztecOOSolver.SetPrecOperator(&Block);
216 AztecOOSolver.Iterate(1550,1e-2);
218 double TrueResidual = AztecOOSolver.TrueResidual();
219 ItersBlock = AztecOOSolver.NumIters();
222 cout <<
"Iterations " << ItersBlock << endl;
223 cout <<
"Norm of the true residual = " << TrueResidual << endl;
227 int diff = ItersPoint - ItersBlock;
228 if (diff < 0) diff = -diff;
233 cout <<
"TEST FAILED!" << endl;
238 cout <<
"TEST PASSED" << endl;
244 bool KrylovTest(std::string PrecType,
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,
bool backward)
248 LHS.PutScalar(0.0);
RHS.Random();
253 Teuchos::ParameterList List;
254 List.set(
"relaxation: damping factor", 1.0);
255 List.set(
"relaxation: type", PrecType);
256 if(backward) List.set(
"relaxation: backward mode",backward);
261 cout <<
"Krylov test: Using " << PrecType
262 <<
" with AztecOO" << endl;
270 List.set(
"relaxation: sweeps",1);
276 AztecOO AztecOOSolver(Problem);
277 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
278 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
279 AztecOOSolver.SetPrecOperator(&Point);
281 AztecOOSolver.Iterate(1550,1e-5);
283 double TrueResidual = AztecOOSolver.TrueResidual();
286 cout <<
"Norm of the true residual = " << TrueResidual << endl;
288 Iters1 = AztecOOSolver.NumIters();
295 List.set(
"relaxation: sweeps",10);
302 AztecOO AztecOOSolver(Problem);
303 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
304 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
305 AztecOOSolver.SetPrecOperator(&Point);
306 AztecOOSolver.Iterate(1550,1e-5);
308 double TrueResidual = AztecOOSolver.TrueResidual();
311 cout <<
"Norm of the true residual = " << TrueResidual << endl;
313 Iters10 = AztecOOSolver.NumIters();
317 cout <<
"Iters_1 = " << Iters1 <<
", Iters_10 = " << Iters10 << endl;
318 cout <<
"(second number should be smaller than first one)" << endl;
321 if (Iters10 > Iters1) {
323 cout <<
"TEST FAILED!" << endl;
328 cout <<
"TEST PASSED" << endl;
334 bool BasicTest(std::string PrecType,
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,
bool backward)
338 LHS.PutScalar(0.0);
RHS.Random();
340 double starting_residual = Galeri::ComputeNorm(&*A, &LHS, &
RHS);
344 Teuchos::ParameterList List;
345 List.set(
"relaxation: damping factor", 1.0);
346 List.set(
"relaxation: sweeps",1550);
347 List.set(
"relaxation: type", PrecType);
348 if(backward) List.set(
"relaxation: backward mode",backward);
360 double residual = Galeri::ComputeNorm(&*A, &LHS, &
RHS);
362 if (A->Comm().MyPID() == 0 &&
verbose)
363 cout <<
"||A * x - b||_2 (scaled) = " << residual / starting_residual << endl;
366 if (residual / starting_residual < 1e-2) {
368 cout <<
"Test passed" << endl;
373 cout <<
"Test failed!" << endl;
379 int main(
int argc,
char *argv[])
382 MPI_Init(&argc,&argv);
390 for (
int i = 1 ; i < argc ; ++i) {
391 if (strcmp(argv[i],
"-s") == 0) {
398 Teuchos::ParameterList GaleriList;
400 GaleriList.set(
"nx", nx);
401 GaleriList.set(
"ny", nx * Comm.
NumProc());
402 GaleriList.set(
"mx", 1);
403 GaleriList.set(
"my", Comm.
NumProc());
404 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap64(
"Cartesian2D", Comm, GaleriList) );
405 Teuchos::RefCountPtr<Epetra_CrsMatrix> A;
407 A = Teuchos::rcp( Galeri::CreateCrsMatrix(
"Laplace2D", &*Map, GaleriList) );
409 A = Teuchos::rcp( Galeri::CreateCrsMatrix(
"Recirc2D", &*Map, GaleriList) );
412 int TestPassed =
true;
422 if(!
BasicTest(
"symmetric Gauss-Seidel",A,
false))
437 if(!
KrylovTest(
"symmetric Gauss-Seidel",A,
false))
455 TestPassed = TestPassed &&
461 TestPassed = TestPassed &&
468 TestPassed = TestPassed &&
477 int Iters4, Iters8, Iters16;
482 if ((Iters16 > Iters8) && (Iters8 > Iters4)) {
484 cout <<
"Test passed" << endl;
488 cout <<
"TEST FAILED!" << endl;
489 TestPassed = TestPassed &&
false;
498 int Iters0, Iters2, Iters4;
502 if ((Iters4 < Iters2) && (Iters2 < Iters0)) {
504 cout <<
"Test passed" << endl;
508 cout <<
"TEST FAILED!" << endl;
509 TestPassed = TestPassed &&
false;
518 cout <<
"Test `TestRelaxation.exe' failed!" << endl;
527 cout <<
"Test `TestRelaxation.exe' passed!" << endl;
529 return(EXIT_SUCCESS);
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
int PutScalar(double ScalarConstant)
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
int main(int argc, char *argv[])
virtual int MyPID() const=0
bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int sweeps)
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
bool KrylovTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward)
virtual const Epetra_Comm & Comm() const=0
virtual int Compute()
Computes the preconditioner.
bool BasicTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward)
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix's...
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual int Compute()
Computes the preconditioners.
int CompareBlockOverlap(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int Overlap)
int CompareBlockSizes(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int NumParts)
static bool SymmetricGallery
Epetra_RowMatrix * GetMatrix() const