Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
64
65#include "Ifpack_Amesos.h"
66#include "AztecOO.h"
67
68static bool verbose = false;
69static bool SymmetricGallery = false;
70static bool Solver = AZ_gmres;
71const int NumVectors = 3;
72
73// ======================================================================
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
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// ======================================================================
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");
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// ======================================================================
257int 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
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// ======================================================================
299int 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
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// ======================================================================
338int 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
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// ======================================================================
380int 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
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// ======================================================================
418int 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
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// ======================================================================
455bool 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
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// ======================================================================
553bool 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// ======================================================================
654bool 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// ======================================================================
711int main(int argc, char *argv[])
712{
713#ifdef HAVE_MPI
714 MPI_Init(&argc,&argv);
715 Epetra_MpiComm Comm( MPI_COMM_WORLD );
716#else
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;
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}
Copy
#define RHS(a)
Definition: MatGenFD.c:60
int GID(int LID) const
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_RowMatrix * GetMatrix() const
int NumProc() const
int MyPID() const
int PutScalar(double ScalarConstant)
virtual const Epetra_Comm & Comm() const=0
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual int Initialize()
Initializes the preconditioner.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the block Jacobi preconditioner to X, returns the result in Y.
virtual int Compute()
Computes the preconditioner.
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix's.
virtual int Compute()
Computes the preconditioners.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
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[])
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)
static bool Solver
bool KrylovTest(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, bool backward, bool reorder=false)
bool TestTriDiVariableBlocking(const Epetra_Comm &Comm)
const int NumVectors
static bool verbose
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