Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/Relaxation/cxx_main.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 
45 #ifdef HAVE_MPI
46 #include "Epetra_MpiComm.h"
47 #else
48 #include "Epetra_SerialComm.h"
49 #endif
50 #include "Epetra_CrsMatrix.h"
51 #include "Epetra_MultiVector.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_LinearProblem.h"
54 #include "Epetra_Map.h"
55 #include "Galeri_Maps.h"
56 #include "Galeri_CrsMatrices.h"
57 #include "Galeri_Utils.h"
58 #include "Teuchos_ParameterList.hpp"
59 #include "Teuchos_RefCountPtr.hpp"
60 #include "Ifpack_PointRelaxation.h"
61 #include "Ifpack_BlockRelaxation.h"
62 #include "Ifpack_SparseContainer.h"
63 #include "Ifpack_TriDiContainer.h"
64 
65 #include "Ifpack_Amesos.h"
66 #include "AztecOO.h"
67 
68 static bool verbose = false;
69 static bool SymmetricGallery = false;
70 static bool Solver = AZ_gmres;
71 const int NumVectors = 3;
72 
73 // ======================================================================
74 bool TestTriDiVariableBlocking(const Epetra_Comm & Comm) {
75  // Basically each processor gets this 5x5 block lower-triangular matrix:
76  //
77  // [ 2 -1 0 0 0 ;...
78  // [-1 2 0 0 0 ;...
79  // [ 0 -1 3 -1 0 ;...
80  // [ 0 0 -1 3 -1 ;...
81  // [ 0 0 0 -1 2 ];
82  //
83 
84  Epetra_Map RowMap(-1,5,0,Comm); // 5 rows per proc
85 
86  Epetra_CrsMatrix A(Copy,RowMap,0);
87 
88  int num_entries;
89  int indices[5];
90  double values[5];
91  int rb = RowMap.GID(0);
92 
93  /*** Fill RHS / LHS ***/
94  Epetra_Vector rhs(RowMap), lhs(RowMap), exact_soln(RowMap);
95  rhs.PutScalar(2.0);
96  lhs.PutScalar(0.0);
97  exact_soln.PutScalar(2.0);
98 
99  /*** Fill Matrix ****/
100  // Row 0
101  num_entries=2;
102  indices[0]=rb; indices[1]=rb+1;
103  values[0] =2; values[1] =-1;
104  A.InsertGlobalValues(rb,num_entries,&values[0],&indices[0]);
105 
106  // Row 1
107  num_entries=2;
108  indices[0]=rb; indices[1]=rb+1;
109  values[0] =-1; values[1] =2;
110  A.InsertGlobalValues(rb+1,num_entries,&values[0],&indices[0]);
111 
112  // Row 2
113  num_entries=3;
114  indices[0]=rb+1; indices[1]=rb+2; indices[2]=rb+3;
115  values[0] =-1; values[1] = 3; values[2] =-1;
116  A.InsertGlobalValues(rb+2,num_entries,&values[0],&indices[0]);
117 
118  // Row 3
119  num_entries=3;
120  indices[0]=rb+2; indices[1]=rb+3; indices[2]=rb+4;
121  values[0] =-1; values[1] = 3; values[2] =-1;
122  A.InsertGlobalValues(rb+3,num_entries,&values[0],&indices[0]);
123 
124  // Row 4
125  num_entries=2;
126  indices[0]=rb+3; indices[1]=rb+4;
127  values[0] =-1; values[1] = 2;
128  A.InsertGlobalValues(rb+4,num_entries,&values[0],&indices[0]);
129  A.FillComplete();
130 
131  /* Setup Block Relaxation */
132  int PartMap[5]={0,0,1,1,1};
133 
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");
140 
141  Ifpack_BlockRelaxation<Ifpack_TriDiContainer> TDRelax(&A);
142 
143  TDRelax.SetParameters(ilist);
144  TDRelax.Initialize();
145  TDRelax.Compute();
146  TDRelax.ApplyInverse(rhs,lhs);
147 
148  double norm;
149  lhs.Update(1.0,exact_soln,-1.0);
150  lhs.Norm2(&norm);
151 
152  if(verbose) cout<<"Variable Block Partitioning Test"<<endl;
153 
154  if(norm < 1e-14) {
155  if(verbose) cout << "Test passed" << endl;
156  return true;
157  }
158  else {
159  if(verbose) cout << "Test failed" << endl;
160  return false;
161  }
162 }
163 // ======================================================================
164 bool TestVariableBlocking(const Epetra_Comm & Comm) {
165  // Basically each processor gets this 5x5 block lower-triangular matrix:
166  //
167  // [ 2 -1 0 0 0 ;...
168  // [-1 2 0 0 0 ;...
169  // [-1 -1 5 -1 -1 ;...
170  // [-1 -1 -1 5 -1 ;...
171  // [-1 -1 -1 -1 5 ];
172  //
173  // The nice thing about this matrix is that if the RHS is a constant,the solution is the same constant...
174 
175  Epetra_Map RowMap(-1,5,0,Comm); // 5 rows per proc
176 
177  Epetra_CrsMatrix A(Copy,RowMap,0);
178 
179  int num_entries;
180  int indices[5];
181  double values[5];
182  int rb = RowMap.GID(0);
183 
184  /*** Fill RHS / LHS ***/
185  Epetra_Vector rhs(RowMap), lhs(RowMap), exact_soln(RowMap);
186  rhs.PutScalar(2.0);
187  lhs.PutScalar(0.0);
188  exact_soln.PutScalar(2.0);
189 
190  /*** Fill Matrix ****/
191  // Row 0
192  num_entries=2;
193  indices[0]=rb; indices[1]=rb+1;
194  values[0] =2; values[1] =-1;
195  A.InsertGlobalValues(rb,num_entries,&values[0],&indices[0]);
196 
197  // Row 1
198  num_entries=2;
199  indices[0]=rb; indices[1]=rb+1;
200  values[0] =-1; values[1] =2;
201  A.InsertGlobalValues(rb+1,num_entries,&values[0],&indices[0]);
202 
203  // Row 2
204  num_entries=5;
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;
207  A.InsertGlobalValues(rb+2,num_entries,&values[0],&indices[0]);
208 
209  // Row 3
210  num_entries=5;
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;
213  A.InsertGlobalValues(rb+3,num_entries,&values[0],&indices[0]);
214 
215  // Row 4
216  num_entries=5;
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;
219  A.InsertGlobalValues(rb+4,num_entries,&values[0],&indices[0]);
220  A.FillComplete();
221 
222 
223  /* Setup Block Relaxation */
224  int PartMap[5]={0,0,1,1,1};
225 
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);
234  Relax.Initialize();
235  Relax.Compute();
236 
237  Relax.ApplyInverse(rhs,lhs);
238 
239 
240  double norm;
241  lhs.Update(1.0,exact_soln,-1.0);
242  lhs.Norm2(&norm);
243 
244  if(verbose) cout<<"Variable Block Partitioning Test"<<endl;
245 
246  if(norm < 1e-14) {
247  if(verbose) cout << "Test passed" << endl;
248  return true;
249  }
250  else {
251  if(verbose) cout << "Test failed" << endl;
252  return false;
253  }
254 }
255 
256 // ======================================================================
257 int CompareLineSmoother(const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, Teuchos::RCP<Epetra_MultiVector> coord)
258 {
259  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
260  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
261  LHS.PutScalar(0.0); RHS.Random();
262 
263  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
264 
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);
275 
276  RHS.PutScalar(1.0);
277  LHS.PutScalar(0.0);
278 
279  Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
280  Prec.SetParameters(List);
281  Prec.Compute();
282 
283  // set AztecOO solver object
284  AztecOO AztecOOSolver(Problem);
285  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
286  if (verbose)
287  AztecOOSolver.SetAztecOption(AZ_output,32);
288  else
289  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
290  AztecOOSolver.SetPrecOperator(&Prec);
291 
292  AztecOOSolver.Iterate(2550,1e-5);
293 
294  return(AztecOOSolver.NumIters());
295 }
296 
297 
298 // ======================================================================
299 int CompareLineSmootherEntries(const Teuchos::RefCountPtr<Epetra_RowMatrix>& A)
300 {
301  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
302  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
303  LHS.PutScalar(0.0); RHS.Random();
304 
305  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
306 
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);
315 
316  RHS.PutScalar(1.0);
317  LHS.PutScalar(0.0);
318 
319  Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
320  Prec.SetParameters(List);
321  Prec.Compute();
322 
323  // set AztecOO solver object
324  AztecOO AztecOOSolver(Problem);
325  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
326  if (verbose)
327  AztecOOSolver.SetAztecOption(AZ_output,32);
328  else
329  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
330  AztecOOSolver.SetPrecOperator(&Prec);
331 
332  AztecOOSolver.Iterate(2550,1e-5);
333 
334  return(AztecOOSolver.NumIters());
335 }
336 
337 // ======================================================================
338 int AllSingle(const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, Teuchos::RCP<Epetra_MultiVector> coord)
339 {
340  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
341  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
342  LHS.PutScalar(0.0); RHS.Random();
343 
344  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
345 
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);
356 
357  RHS.PutScalar(1.0);
358  LHS.PutScalar(0.0);
359 
360  Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
361  Prec.SetParameters(List);
362  Prec.Compute();
363 
364  // set AztecOO solver object
365  AztecOO AztecOOSolver(Problem);
366  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
367  if (verbose)
368  AztecOOSolver.SetAztecOption(AZ_output,32);
369  else
370  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
371  AztecOOSolver.SetPrecOperator(&Prec);
372 
373  AztecOOSolver.Iterate(2550,1e-5);
374 
375  printf(" AllSingle iters %d \n",AztecOOSolver.NumIters());
376  return(AztecOOSolver.NumIters());
377 }
378 
379 // ======================================================================
380 int CompareBlockOverlap(const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int Overlap)
381 {
382  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
383  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
384  LHS.PutScalar(0.0); RHS.Random();
385 
386  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
387 
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);
395 
396  RHS.PutScalar(1.0);
397  LHS.PutScalar(0.0);
398 
399  Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
400  Prec.SetParameters(List);
401  Prec.Compute();
402 
403  // set AztecOO solver object
404  AztecOO AztecOOSolver(Problem);
405  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
406  if (verbose)
407  AztecOOSolver.SetAztecOption(AZ_output,32);
408  else
409  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
410  AztecOOSolver.SetPrecOperator(&Prec);
411 
412  AztecOOSolver.Iterate(2550,1e-5);
413 
414  return(AztecOOSolver.NumIters());
415 }
416 
417 // ======================================================================
418 int CompareBlockSizes(std::string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int NumParts)
419 {
420  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
421  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
422  LHS.PutScalar(0.0); RHS.Random();
423 
424  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
425 
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);
432 
433  RHS.PutScalar(1.0);
434  LHS.PutScalar(0.0);
435 
436  Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Prec(&*A);
437  Prec.SetParameters(List);
438  Prec.Compute();
439 
440  // set AztecOO solver object
441  AztecOO AztecOOSolver(Problem);
442  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
443  if (verbose)
444  AztecOOSolver.SetAztecOption(AZ_output,32);
445  else
446  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
447  AztecOOSolver.SetPrecOperator(&Prec);
448 
449  AztecOOSolver.Iterate(2550,1e-5);
450 
451  return(AztecOOSolver.NumIters());
452 }
453 
454 // ======================================================================
455 bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int sweeps)
456 {
457  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
458  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
459  LHS.PutScalar(0.0); RHS.Random();
460 
461  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
462 
463  // Set up the list
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());
470 
471  int ItersPoint, ItersBlock;
472 
473  // ================================================== //
474  // get the number of iterations with point relaxation //
475  // ================================================== //
476  {
477  RHS.PutScalar(1.0);
478  LHS.PutScalar(0.0);
479 
480  Ifpack_PointRelaxation Point(&*A);
481  Point.SetParameters(List);
482  Point.Compute();
483 
484  // set AztecOO solver object
485  AztecOO AztecOOSolver(Problem);
486  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
487  if (verbose)
488  AztecOOSolver.SetAztecOption(AZ_output,32);
489  else
490  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
491  AztecOOSolver.SetPrecOperator(&Point);
492 
493  AztecOOSolver.Iterate(2550,1e-2);
494 
495  double TrueResidual = AztecOOSolver.TrueResidual();
496  ItersPoint = AztecOOSolver.NumIters();
497  // some output
498  if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
499  cout << "Iterations = " << ItersPoint << endl;
500  cout << "Norm of the true residual = " << TrueResidual << endl;
501  }
502  }
503 
504  // ================================================== //
505  // get the number of iterations with block relaxation //
506  // ================================================== //
507  {
508 
509  RHS.PutScalar(1.0);
510  LHS.PutScalar(0.0);
511 
512  Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > Block(&*A);
513  Block.SetParameters(List);
514  Block.Compute();
515 
516  // set AztecOO solver object
517  AztecOO AztecOOSolver(Problem);
518  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
519  if (verbose)
520  AztecOOSolver.SetAztecOption(AZ_output,32);
521  else
522  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
523  AztecOOSolver.SetPrecOperator(&Block);
524 
525  AztecOOSolver.Iterate(2550,1e-2);
526 
527  double TrueResidual = AztecOOSolver.TrueResidual();
528  ItersBlock = AztecOOSolver.NumIters();
529  // some output
530  if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
531  cout << "Iterations " << ItersBlock << endl;
532  cout << "Norm of the true residual = " << TrueResidual << endl;
533  }
534  }
535 
536  int diff = ItersPoint - ItersBlock;
537  if (diff < 0) diff = -diff;
538 
539  if (diff > 10)
540  {
541  if (verbose)
542  cout << "ComparePointandBlock TEST FAILED!" << endl;
543  return(false);
544  }
545  else {
546  if (verbose)
547  cout << "ComparePointandBlock TEST PASSED" << endl;
548  return(true);
549  }
550 }
551 
552 // ======================================================================
553 bool KrylovTest(std::string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, bool backward, bool reorder=false)
554 {
555  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
556  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
557  LHS.PutScalar(0.0); RHS.Random();
558 
559  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
560 
561  // Set up the list
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);
566 
567  // Reordering if needed
568  int NumRows=A->NumMyRows();
569  std::vector<int> RowList(NumRows);
570  if(reorder) {
571  for(int i=0; i<NumRows; i++)
572  RowList[i]=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);
575  }
576 
577 
578  int Iters1, Iters10;
579 
580  if (verbose) {
581  cout << "Krylov test: Using " << PrecType
582  << " with AztecOO" << endl;
583  }
584 
585  // ============================================== //
586  // get the number of iterations with 1 sweep only //
587  // ============================================== //
588  {
589 
590  List.set("relaxation: sweeps",1);
591  Ifpack_PointRelaxation Point(&*A);
592  Point.SetParameters(List);
593  Point.Compute();
594 
595  // set AztecOO solver object
596  AztecOO AztecOOSolver(Problem);
597  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
598  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
599  AztecOOSolver.SetPrecOperator(&Point);
600 
601  AztecOOSolver.Iterate(2550,1e-5);
602 
603  double TrueResidual = AztecOOSolver.TrueResidual();
604  // some output
605  if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
606  cout << "Norm of the true residual = " << TrueResidual << endl;
607  }
608  Iters1 = AztecOOSolver.NumIters();
609  }
610 
611  // ======================================================== //
612  // now re-run with 10 sweeps, solver should converge faster
613  // ======================================================== //
614  {
615  List.set("relaxation: sweeps",10);
616  Ifpack_PointRelaxation Point(&*A);
617  Point.SetParameters(List);
618  Point.Compute();
619  LHS.PutScalar(0.0);
620 
621  // set AztecOO solver object
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);
627 
628  double TrueResidual = AztecOOSolver.TrueResidual();
629  // some output
630  if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
631  cout << "Norm of the true residual = " << TrueResidual << endl;
632  }
633  Iters10 = AztecOOSolver.NumIters();
634  }
635 
636  if (verbose) {
637  cout << "Iters_1 = " << Iters1 << ", Iters_10 = " << Iters10 << endl;
638  cout << "(second number should be smaller than first one)" << endl;
639  }
640 
641  if (Iters10 > Iters1) {
642  if (verbose)
643  cout << "KrylovTest TEST FAILED!" << endl;
644  return(false);
645  }
646  else {
647  if (verbose)
648  cout << "KrylovTest TEST PASSED" << endl;
649  return(true);
650  }
651 }
652 
653 // ======================================================================
654 bool BasicTest(std::string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A,bool backward, bool reorder=false)
655 {
656  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
657  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
658  LHS.PutScalar(0.0); RHS.Random();
659 
660  double starting_residual = Galeri::ComputeNorm(&*A, &LHS, &RHS);
661  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
662 
663  // Set up the list
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);
669 
670  // Reordering if needed
671  int NumRows=A->NumMyRows();
672  std::vector<int> RowList(NumRows);
673  if(reorder) {
674  for(int i=0; i<NumRows; i++)
675  RowList[i]=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);
678  }
679 
680  Ifpack_PointRelaxation Point(&*A);
681 
682  Point.SetParameters(List);
683  Point.Compute();
684  // use the preconditioner as solver, with 1550 iterations
685  Point.ApplyInverse(RHS,LHS);
686 
687  // compute the real residual
688 
689  double residual = Galeri::ComputeNorm(&*A, &LHS, &RHS);
690 
691  if (A->Comm().MyPID() == 0 && verbose)
692  cout << "||A * x - b||_2 (scaled) = " << residual / starting_residual << endl;
693 
694  // Jacobi is very slow to converge here
695  if (residual / starting_residual < 1e-2) {
696  if (verbose)
697  cout << "BasicTest Test passed" << endl;
698  return(true);
699  }
700  else {
701  if (verbose)
702  cout << "BasicTest Test failed!" << endl;
703  return(false);
704  }
705 }
706 
707 
708 
709 
710 // ======================================================================
711 int main(int argc, char *argv[])
712 {
713 #ifdef HAVE_MPI
714  MPI_Init(&argc,&argv);
715  Epetra_MpiComm Comm( MPI_COMM_WORLD );
716 #else
717  Epetra_SerialComm Comm;
718 #endif
719 
720  verbose = (Comm.MyPID() == 0);
721 
722  for (int i = 1 ; i < argc ; ++i) {
723  if (strcmp(argv[i],"-s") == 0) {
724  SymmetricGallery = true;
725  Solver = AZ_cg;
726  }
727  }
728 
729  // size of the global matrix.
730  Teuchos::ParameterList GaleriList;
731  int nx = 30;
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;
738  if (SymmetricGallery)
739  A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
740  else
741  A = Teuchos::rcp( Galeri::CreateCrsMatrix("Recirc2D", &*Map, GaleriList) );
742 
743  // coordinates
744  Teuchos::RCP<Epetra_MultiVector> coord = Teuchos::rcp( Galeri::CreateCartesianCoordinates("2D",&*Map,GaleriList));
745 
746  // test the preconditioner
747  int TestPassed = true;
748 
749  // ======================================== //
750  // first verify that we can get convergence //
751  // with all point relaxation methods //
752  // ======================================== //
753 
754  if(!BasicTest("Jacobi",A,false))
755  TestPassed = false;
756 
757  if(!BasicTest("symmetric Gauss-Seidel",A,false))
758  TestPassed = false;
759 
760  if(!BasicTest("symmetric Gauss-Seidel",A,false,true))
761  TestPassed = false;
762 
763  if (!SymmetricGallery) {
764  if(!BasicTest("Gauss-Seidel",A,false))
765  TestPassed = false;
766  if(!BasicTest("Gauss-Seidel",A,true))
767  TestPassed = false;
768 
769  if(!BasicTest("Gauss-Seidel",A,false,true))
770  TestPassed = false;
771  if(!BasicTest("Gauss-Seidel",A,true,true))
772  TestPassed = false;
773 
774  }
775 
776  // ============================= //
777  // check uses as preconditioners //
778  // ============================= //
779 
780  if(!KrylovTest("symmetric Gauss-Seidel",A,false))
781  TestPassed = false;
782 
783  if(!KrylovTest("symmetric Gauss-Seidel",A,false,true))
784  TestPassed = false;
785 
786 
787  if (!SymmetricGallery) {
788  if(!KrylovTest("Gauss-Seidel",A,false))
789  TestPassed = false;
790  if(!KrylovTest("Gauss-Seidel",A,true))
791  TestPassed = false;
792 
793  if(!KrylovTest("Gauss-Seidel",A,false,true))
794  TestPassed = false;
795  if(!KrylovTest("Gauss-Seidel",A,true,true))
796  TestPassed = false;
797 
798  }
799 
800  // ================================== //
801  // compare point and block relaxation //
802  // ================================== //
803 
804  //TestPassed = TestPassed &&
805  // ComparePointAndBlock("Jacobi",A,1);
806 
807  TestPassed = TestPassed &&
808  ComparePointAndBlock("Jacobi",A,10);
809 
810  //TestPassed = TestPassed &&
811  //ComparePointAndBlock("symmetric Gauss-Seidel",A,1);
812 
813  TestPassed = TestPassed &&
814  ComparePointAndBlock("symmetric Gauss-Seidel",A,10);
815 
816  if (!SymmetricGallery) {
817  //TestPassed = TestPassed &&
818  //ComparePointAndBlock("Gauss-Seidel",A,1);
819 
820  TestPassed = TestPassed &&
821  ComparePointAndBlock("Gauss-Seidel",A,10);
822  }
823 
824  // ============================ //
825  // verify effect of # of blocks //
826  // ============================ //
827 
828  {
829  int Iters4, Iters8, Iters16;
830  Iters4 = CompareBlockSizes("Jacobi",A,4);
831  Iters8 = CompareBlockSizes("Jacobi",A,8);
832  Iters16 = CompareBlockSizes("Jacobi",A,16);
833 
834  if ((Iters16 > Iters8) && (Iters8 > Iters4)) {
835  if (verbose)
836  cout << "CompareBlockSizes Test passed" << endl;
837  }
838  else {
839  if (verbose)
840  cout << "CompareBlockSizes TEST FAILED!" << endl;
841  TestPassed = TestPassed && false;
842  }
843  }
844 
845  // ================================== //
846  // verify effect of overlap in Jacobi //
847  // ================================== //
848 
849  {
850  int Iters0, Iters2, Iters4;
851  Iters0 = CompareBlockOverlap(A,0);
852  Iters2 = CompareBlockOverlap(A,2);
853  Iters4 = CompareBlockOverlap(A,4);
854  if ((Iters4 < Iters2) && (Iters2 < Iters0)) {
855  if (verbose)
856  cout << "CompareBlockOverlap Test passed" << endl;
857  }
858  else {
859  if (verbose)
860  cout << "CompareBlockOverlap TEST FAILED!" << endl;
861  TestPassed = TestPassed && false;
862  }
863  }
864 
865  // ================================== //
866  // check if (coordinate) line smoothing works //
867  // ================================== //
868  {
869  int Iters1= CompareLineSmoother(A,coord);
870  printf(" comparelinesmoother (coordinate) iters %d \n",Iters1);
871  }
872 
873 
874  // ================================= //
875  // check if (entry) line smoothing works //
876  // ================================== //
877  {
878  int Iters1= CompareLineSmootherEntries(A);
879  printf(" comparelinesmoother (entries) iters %d \n",Iters1);
880  }
881 
882  // ================================== //
883  // check if All singleton version of CompareLineSmoother //
884  // ================================== //
885  {
886 
887  AllSingle(A,coord);
888 
889  }
890 
891  // ================================== //
892  // test variable blocking //
893  // ================================== //
894  {
895  TestPassed = TestPassed && TestVariableBlocking(A->Comm());
896  }
897 
898  // ================================== //
899  // test variable blocking //
900  // ================================== //
901  {
902  TestPassed = TestPassed && TestTriDiVariableBlocking(A->Comm());
903  }
904 
905  // ============ //
906  // final output //
907  // ============ //
908 
909  if (!TestPassed) {
910  cout << "Test `TestRelaxation.exe' failed!" << endl;
911  exit(EXIT_FAILURE);
912  }
913 
914 #ifdef HAVE_MPI
915  MPI_Finalize();
916 #endif
917 
918  cout << endl;
919  cout << "Test `TestRelaxation.exe' passed!" << endl;
920  cout << endl;
921  return(EXIT_SUCCESS);
922 }
Epetra_MultiVector
Ifpack_ConfigDefs.h
Epetra_Map
CompareLineSmoother
int CompareLineSmoother(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, Teuchos::RCP< Epetra_MultiVector > coord)
Definition: test/Relaxation/cxx_main.cpp:257
Ifpack_BlockRelaxation
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix&#39;s.
Definition: Ifpack_BlockRelaxation.h:139
verbose
bool verbose
Definition: test/CompareWithAztecOO/cxx_main.cpp:66
SymmetricGallery
static bool SymmetricGallery
Definition: test/Relaxation/cxx_main.cpp:69
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Epetra_CrsMatrix::InsertGlobalValues
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
BasicTest
bool BasicTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward, bool reorder=false)
Definition: test/Relaxation/cxx_main.cpp:654
KrylovTest
bool KrylovTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward, bool reorder=false)
Definition: test/Relaxation/cxx_main.cpp:553
Ifpack_BlockRelaxation::ApplyInverse
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the block Jacobi preconditioner to X, returns the result in Y.
Definition: Ifpack_BlockRelaxation.h:614
NumVectors
const int NumVectors
Definition: performance.cpp:71
Epetra_Vector
Epetra_SerialComm::MyPID
int MyPID() const
Epetra_MpiComm
Epetra_Comm::MyPID
virtual int MyPID() const =0
Epetra_CrsMatrix::FillComplete
int FillComplete(bool OptimizeDataStorage=true)
Ifpack_PointRelaxation::SetParameters
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Definition: Ifpack_PointRelaxation.cpp:97
TestVariableBlocking
bool TestVariableBlocking(const Epetra_Comm &Comm)
Definition: test/Relaxation/cxx_main.cpp:164
Solver
static bool Solver
Definition: performance.cpp:70
Epetra_Operator::Comm
virtual const Epetra_Comm & Comm() const =0
Ifpack_BlockRelaxation::Compute
virtual int Compute()
Computes the preconditioner.
Definition: Ifpack_BlockRelaxation.h:581
Epetra_Comm
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ifpack_BlockRelaxation.h
AllSingle
int AllSingle(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, Teuchos::RCP< Epetra_MultiVector > coord)
Definition: test/Relaxation/cxx_main.cpp:338
Epetra_BlockMap::GID
int GID(int LID) const
Ifpack_PointRelaxation
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix&#39;s...
Definition: Ifpack_PointRelaxation.h:130
Ifpack_BlockRelaxation::SetParameters
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Definition: Ifpack_BlockRelaxation.h:1166
main
int main(int argc, char *argv[])
Definition: Ifpack_ex_Amesos.cpp:61
Ifpack_PointRelaxation::Compute
virtual int Compute()
Computes the preconditioners.
Definition: Ifpack_PointRelaxation.cpp:201
Epetra_SerialComm::NumProc
int NumProc() const
Epetra_SerialComm
CompareLineSmootherEntries
int CompareLineSmootherEntries(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A)
Definition: test/Relaxation/cxx_main.cpp:299
Ifpack_Amesos.h
ComparePointAndBlock
bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int sweeps)
Definition: performance.cpp:74
Epetra_LinearProblem::GetMatrix
Epetra_RowMatrix * GetMatrix() const
Ifpack_SparseContainer.h
Ifpack_BlockRelaxation::Initialize
virtual int Initialize()
Initializes the preconditioner.
Definition: Ifpack_BlockRelaxation.h:1238
Epetra_CrsMatrix
CompareBlockOverlap
int CompareBlockOverlap(const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int Overlap)
Definition: test/Relaxation/cxx_main.cpp:380
Epetra_LinearProblem
Ifpack_TriDiContainer.h
TestTriDiVariableBlocking
bool TestTriDiVariableBlocking(const Epetra_Comm &Comm)
Definition: test/Relaxation/cxx_main.cpp:74
Ifpack_PointRelaxation::ApplyInverse
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
Definition: Ifpack_PointRelaxation.cpp:414
Ifpack_PointRelaxation.h
CompareBlockSizes
int CompareBlockSizes(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int NumParts)
Definition: test/Relaxation/cxx_main.cpp:418
RHS
#define RHS(a)
Definition: MatGenFD.c:60
Generated by   doxygen 1.8.5