91 int rb = RowMap.
GID(0);
97 exact_soln.PutScalar(2.0);
102 indices[0]=rb; indices[1]=rb+1;
103 values[0] =2; values[1] =-1;
108 indices[0]=rb; indices[1]=rb+1;
109 values[0] =-1; values[1] =2;
114 indices[0]=rb+1; indices[1]=rb+2; indices[2]=rb+3;
115 values[0] =-1; values[1] = 3; values[2] =-1;
120 indices[0]=rb+2; indices[1]=rb+3; indices[2]=rb+4;
121 values[0] =-1; values[1] = 3; values[2] =-1;
126 indices[0]=rb+3; indices[1]=rb+4;
127 values[0] =-1; values[1] = 2;
132 int PartMap[5]={0,0,1,1,1};
134 Teuchos::ParameterList ilist;
135 ilist.set(
"partitioner: type",
"user");
136 ilist.set(
"partitioner: map",&PartMap[0]);
137 ilist.set(
"partitioner: local parts",2);
138 ilist.set(
"relaxation: sweeps",1);
139 ilist.set(
"relaxation: type",
"Gauss-Seidel");
141 Ifpack_BlockRelaxation<Ifpack_TriDiContainer> TDRelax(&A);
143 TDRelax.SetParameters(ilist);
144 TDRelax.Initialize();
146 TDRelax.ApplyInverse(rhs,lhs);
149 lhs.Update(1.0,exact_soln,-1.0);
152 if(
verbose) cout<<
"Variable Block Partitioning Test"<<endl;
155 if(
verbose) cout <<
"Test passed" << endl;
159 if(
verbose) cout <<
"Test failed" << endl;
182 int rb = RowMap.
GID(0);
188 exact_soln.PutScalar(2.0);
193 indices[0]=rb; indices[1]=rb+1;
194 values[0] =2; values[1] =-1;
199 indices[0]=rb; indices[1]=rb+1;
200 values[0] =-1; values[1] =2;
205 indices[0]=rb; indices[1]=rb+1; indices[2]=rb+2; indices[3]=rb+3; indices[4]=rb+4;
206 values[0] =-1; values[1] =-1; values[2] = 5; values[3] =-1; values[4] =-1;
211 indices[0]=rb; indices[1]=rb+1; indices[2]=rb+2; indices[3]=rb+3; indices[4]=rb+4;
212 values[0] =-1; values[1] =-1; values[2] =-1; values[3] = 5; values[4] =-1;
217 indices[0]=rb; indices[1]=rb+1; indices[2]=rb+2; indices[3]=rb+3; indices[4]=rb+4;
218 values[0] =-1; values[1] =-1; values[2] =-1; values[3] =-1; values[4] = 5;
224 int PartMap[5]={0,0,1,1,1};
226 Teuchos::ParameterList ilist;
227 ilist.set(
"partitioner: type",
"user");
228 ilist.set(
"partitioner: map",&PartMap[0]);
229 ilist.set(
"partitioner: local parts",2);
230 ilist.set(
"relaxation: sweeps",1);
231 ilist.set(
"relaxation: type",
"Gauss-Seidel");
232 Ifpack_BlockRelaxation<Ifpack_DenseContainer> Relax(&A);
233 Relax.SetParameters(ilist);
237 Relax.ApplyInverse(rhs,lhs);
241 lhs.Update(1.0,exact_soln,-1.0);
244 if(
verbose) cout<<
"Variable Block Partitioning Test"<<endl;
247 if(
verbose) cout <<
"Test passed" << endl;
251 if(
verbose) cout <<
"Test failed" << endl;
257int CompareLineSmoother(
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, Teuchos::RCP<Epetra_MultiVector> coord)
265 Teuchos::ParameterList List;
266 List.set(
"relaxation: damping factor", 1.0);
267 List.set(
"relaxation: type",
"symmetric Gauss-Seidel");
268 List.set(
"relaxation: sweeps",1);
269 List.set(
"partitioner: overlap",0);
270 List.set(
"partitioner: type",
"line");
271 List.set(
"partitioner: line detection threshold",0.1);
272 List.set(
"partitioner: x-coordinates",&(*coord)[0][0]);
273 List.set(
"partitioner: y-coordinates",&(*coord)[1][0]);
274 List.set(
"partitioner: z-coordinates",(
double*) 0);
279 Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
280 Prec.SetParameters(List);
284 AztecOO AztecOOSolver(Problem);
285 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
287 AztecOOSolver.SetAztecOption(AZ_output,32);
289 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
290 AztecOOSolver.SetPrecOperator(&Prec);
292 AztecOOSolver.Iterate(2550,1e-5);
294 return(AztecOOSolver.NumIters());
307 Teuchos::ParameterList List;
308 List.set(
"relaxation: damping factor", 1.0);
309 List.set(
"relaxation: type",
"symmetric Gauss-Seidel");
310 List.set(
"relaxation: sweeps",1);
311 List.set(
"partitioner: overlap",0);
312 List.set(
"partitioner: type",
"line");
313 List.set(
"partitioner: line mode",
"matrix entries");
314 List.set(
"partitioner: line detection threshold",10.0);
319 Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
320 Prec.SetParameters(List);
324 AztecOO AztecOOSolver(Problem);
325 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
327 AztecOOSolver.SetAztecOption(AZ_output,32);
329 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
330 AztecOOSolver.SetPrecOperator(&Prec);
332 AztecOOSolver.Iterate(2550,1e-5);
334 return(AztecOOSolver.NumIters());
338int AllSingle(
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, Teuchos::RCP<Epetra_MultiVector> coord)
346 Teuchos::ParameterList List;
347 List.set(
"relaxation: damping factor", 1.0);
348 List.set(
"relaxation: type",
"symmetric Gauss-Seidel");
349 List.set(
"relaxation: sweeps",1);
350 List.set(
"partitioner: overlap",0);
351 List.set(
"partitioner: type",
"line");
352 List.set(
"partitioner: line detection threshold",1.0);
353 List.set(
"partitioner: x-coordinates",&(*coord)[0][0]);
354 List.set(
"partitioner: y-coordinates",&(*coord)[1][0]);
355 List.set(
"partitioner: z-coordinates",(
double*) 0);
360 Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
361 Prec.SetParameters(List);
365 AztecOO AztecOOSolver(Problem);
366 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
368 AztecOOSolver.SetAztecOption(AZ_output,32);
370 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
371 AztecOOSolver.SetPrecOperator(&Prec);
373 AztecOOSolver.Iterate(2550,1e-5);
375 printf(
" AllSingle iters %d \n",AztecOOSolver.NumIters());
376 return(AztecOOSolver.NumIters());
388 Teuchos::ParameterList List;
389 List.set(
"relaxation: damping factor", 1.0);
390 List.set(
"relaxation: type",
"Jacobi");
391 List.set(
"relaxation: sweeps",1);
392 List.set(
"partitioner: overlap", Overlap);
393 List.set(
"partitioner: type",
"linear");
394 List.set(
"partitioner: local parts", 16);
399 Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
400 Prec.SetParameters(List);
404 AztecOO AztecOOSolver(Problem);
405 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
407 AztecOOSolver.SetAztecOption(AZ_output,32);
409 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
410 AztecOOSolver.SetPrecOperator(&Prec);
412 AztecOOSolver.Iterate(2550,1e-5);
414 return(AztecOOSolver.NumIters());
418int CompareBlockSizes(std::string PrecType,
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,
int NumParts)
426 Teuchos::ParameterList List;
427 List.set(
"relaxation: damping factor", 1.0);
428 List.set(
"relaxation: type", PrecType);
429 List.set(
"relaxation: sweeps",1);
430 List.set(
"partitioner: type",
"linear");
431 List.set(
"partitioner: local parts", NumParts);
436 Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
437 Prec.SetParameters(List);
441 AztecOO AztecOOSolver(Problem);
442 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
444 AztecOOSolver.SetAztecOption(AZ_output,32);
446 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
447 AztecOOSolver.SetPrecOperator(&Prec);
449 AztecOOSolver.Iterate(2550,1e-5);
451 return(AztecOOSolver.NumIters());
464 Teuchos::ParameterList List;
465 List.set(
"relaxation: damping factor", 1.0);
466 List.set(
"relaxation: type", PrecType);
467 List.set(
"relaxation: sweeps",sweeps);
468 List.set(
"partitioner: type",
"linear");
469 List.set(
"partitioner: local parts", A->NumMyRows());
471 int ItersPoint, ItersBlock;
481 Point.SetParameters(List);
485 AztecOO AztecOOSolver(Problem);
486 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
488 AztecOOSolver.SetAztecOption(AZ_output,32);
490 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
491 AztecOOSolver.SetPrecOperator(&Point);
493 AztecOOSolver.Iterate(2550,1e-2);
495 double TrueResidual = AztecOOSolver.TrueResidual();
496 ItersPoint = AztecOOSolver.NumIters();
498 if (
verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
499 cout <<
"Iterations = " << ItersPoint << endl;
500 cout <<
"Norm of the true residual = " << TrueResidual << endl;
512 Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Block(&*A);
513 Block.SetParameters(List);
517 AztecOO AztecOOSolver(Problem);
518 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
520 AztecOOSolver.SetAztecOption(AZ_output,32);
522 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
523 AztecOOSolver.SetPrecOperator(&Block);
525 AztecOOSolver.Iterate(2550,1e-2);
527 double TrueResidual = AztecOOSolver.TrueResidual();
528 ItersBlock = AztecOOSolver.NumIters();
530 if (
verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
531 cout <<
"Iterations " << ItersBlock << endl;
532 cout <<
"Norm of the true residual = " << TrueResidual << endl;
536 int diff = ItersPoint - ItersBlock;
537 if (diff < 0) diff = -diff;
542 cout <<
"ComparePointandBlock TEST FAILED!" << endl;
547 cout <<
"ComparePointandBlock TEST PASSED" << endl;
553bool KrylovTest(std::string PrecType,
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,
bool backward,
bool reorder=
false)
562 Teuchos::ParameterList List;
563 List.set(
"relaxation: damping factor", 1.0);
564 List.set(
"relaxation: type", PrecType);
565 if(backward) List.set(
"relaxation: backward mode",backward);
568 int NumRows=A->NumMyRows();
569 std::vector<int> RowList(NumRows);
571 for(
int i=0; i<NumRows; i++)
573 List.set(
"relaxation: number of local smoothing indices",NumRows);
574 List.set(
"relaxation: local smoothing indices",RowList.size()>0? &RowList[0] : (
int*)0);
581 cout <<
"Krylov test: Using " << PrecType
582 <<
" with AztecOO" << endl;
590 List.set(
"relaxation: sweeps",1);
592 Point.SetParameters(List);
596 AztecOO AztecOOSolver(Problem);
597 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
598 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
599 AztecOOSolver.SetPrecOperator(&Point);
601 AztecOOSolver.Iterate(2550,1e-5);
603 double TrueResidual = AztecOOSolver.TrueResidual();
605 if (
verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
606 cout <<
"Norm of the true residual = " << TrueResidual << endl;
608 Iters1 = AztecOOSolver.NumIters();
615 List.set(
"relaxation: sweeps",10);
617 Point.SetParameters(List);
622 AztecOO AztecOOSolver(Problem);
623 AztecOOSolver.SetAztecOption(AZ_solver,
Solver);
624 AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
625 AztecOOSolver.SetPrecOperator(&Point);
626 AztecOOSolver.Iterate(2550,1e-5);
628 double TrueResidual = AztecOOSolver.TrueResidual();
630 if (
verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
631 cout <<
"Norm of the true residual = " << TrueResidual << endl;
633 Iters10 = AztecOOSolver.NumIters();
637 cout <<
"Iters_1 = " << Iters1 <<
", Iters_10 = " << Iters10 << endl;
638 cout <<
"(second number should be smaller than first one)" << endl;
641 if (Iters10 > Iters1) {
643 cout <<
"KrylovTest TEST FAILED!" << endl;
648 cout <<
"KrylovTest TEST PASSED" << endl;
654bool BasicTest(std::string PrecType,
const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,
bool backward,
bool reorder=
false)
660 double starting_residual = Galeri::ComputeNorm(&*A, &LHS, &
RHS);
664 Teuchos::ParameterList List;
665 List.set(
"relaxation: damping factor", 1.0);
666 List.set(
"relaxation: sweeps",2550);
667 List.set(
"relaxation: type", PrecType);
668 if(backward) List.set(
"relaxation: backward mode",backward);
671 int NumRows=A->NumMyRows();
672 std::vector<int> RowList(NumRows);
674 for(
int i=0; i<NumRows; i++)
676 List.set(
"relaxation: number of local smoothing indices",NumRows);
677 List.set(
"relaxation: local smoothing indices",RowList.size()>0? &RowList[0] : (
int*)0);
682 Point.SetParameters(List);
685 Point.ApplyInverse(
RHS,LHS);
689 double residual = Galeri::ComputeNorm(&*A, &LHS, &
RHS);
691 if (A->Comm().MyPID() == 0 &&
verbose)
692 cout <<
"||A * x - b||_2 (scaled) = " << residual / starting_residual << endl;
695 if (residual / starting_residual < 1e-2) {
697 cout <<
"BasicTest Test passed" << endl;
702 cout <<
"BasicTest Test failed!" << endl;
711int main(
int argc,
char *argv[])
714 MPI_Init(&argc,&argv);
722 for (
int i = 1 ; i < argc ; ++i) {
723 if (strcmp(argv[i],
"-s") == 0) {
730 Teuchos::ParameterList GaleriList;
732 GaleriList.set(
"nx", nx);
733 GaleriList.set(
"ny", nx * Comm.
NumProc());
734 GaleriList.set(
"mx", 1);
735 GaleriList.set(
"my", Comm.
NumProc());
736 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap(
"Cartesian2D", Comm, GaleriList) );
737 Teuchos::RefCountPtr<Epetra_CrsMatrix> A;
739 A = Teuchos::rcp( Galeri::CreateCrsMatrix(
"Laplace2D", &*Map, GaleriList) );
741 A = Teuchos::rcp( Galeri::CreateCrsMatrix(
"Recirc2D", &*Map, GaleriList) );
744 Teuchos::RCP<Epetra_MultiVector> coord = Teuchos::rcp( Galeri::CreateCartesianCoordinates(
"2D",&*Map,GaleriList));
747 int TestPassed =
true;
757 if(!
BasicTest(
"symmetric Gauss-Seidel",A,
false))
760 if(!
BasicTest(
"symmetric Gauss-Seidel",A,
false,
true))
769 if(!
BasicTest(
"Gauss-Seidel",A,
false,
true))
771 if(!
BasicTest(
"Gauss-Seidel",A,
true,
true))
780 if(!
KrylovTest(
"symmetric Gauss-Seidel",A,
false))
783 if(!
KrylovTest(
"symmetric Gauss-Seidel",A,
false,
true))
807 TestPassed = TestPassed &&
813 TestPassed = TestPassed &&
820 TestPassed = TestPassed &&
829 int Iters4, Iters8, Iters16;
834 if ((Iters16 > Iters8) && (Iters8 > Iters4)) {
836 cout <<
"CompareBlockSizes Test passed" << endl;
840 cout <<
"CompareBlockSizes TEST FAILED!" << endl;
841 TestPassed = TestPassed &&
false;
850 int Iters0, Iters2, Iters4;
854 if ((Iters4 < Iters2) && (Iters2 < Iters0)) {
856 cout <<
"CompareBlockOverlap Test passed" << endl;
860 cout <<
"CompareBlockOverlap TEST FAILED!" << endl;
861 TestPassed = TestPassed &&
false;
870 printf(
" comparelinesmoother (coordinate) iters %d \n",Iters1);
879 printf(
" comparelinesmoother (entries) iters %d \n",Iters1);
910 cout <<
"Test `TestRelaxation.exe' failed!" << endl;
919 cout <<
"Test `TestRelaxation.exe' passed!" << endl;
921 return(EXIT_SUCCESS);
int main(int argc, char *argv[])
int CompareBlockOverlap(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int Overlap)
int CompareLineSmoother(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, Teuchos::RCP< Epetra_MultiVector > coord)
int AllSingle(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, Teuchos::RCP< Epetra_MultiVector > coord)
int CompareBlockSizes(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int NumParts)
int CompareLineSmootherEntries(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A)
bool TestVariableBlocking(const Epetra_Comm &Comm)
bool KrylovTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward, bool reorder=false)
bool TestTriDiVariableBlocking(const Epetra_Comm &Comm)
bool BasicTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward, bool reorder=false)
bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int sweeps)
static bool SymmetricGallery