Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/BlockCheby/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#if defined(HAVE_IFPACK_AZTECOO) && defined(HAVE_IFPACK_EPETRAEXT)
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_Vector.h"
53#include "Trilinos_Util_CrsMatrixGallery.h"
54#include "Teuchos_ParameterList.hpp"
59#include "AztecOO.h"
62#include "Ifpack_Chebyshev.h"
63
64using namespace Trilinos_Util;
65/*
66(unused; commented out to avoid build warnings)
67static bool verbose = false;
68static bool SymmetricGallery = false;
69static bool Solver = AZ_gmres;
70*/
71
72//=============================================
73// Test the BlockDiagMatrix
74bool TestBlockDiagMatrix(const Epetra_Comm& Comm){
75 using std::cout;
76 using std::endl;
77
78 const int NUM_BLOCKS=30;
79 bool TestPassed=true;
80 int my_blockgids[NUM_BLOCKS];
81 int my_blocksizes[NUM_BLOCKS];
82
83 for(int i=0;i<NUM_BLOCKS;i++){
84 my_blockgids[i]=i+NUM_BLOCKS*Comm.MyPID();
85 if(i<NUM_BLOCKS/3)
86 my_blocksizes[i]=1;
87 else if(i<2*NUM_BLOCKS/3)
88 my_blocksizes[i]=2;
89 else
90 my_blocksizes[i]=3;
91 }
92
93
94 // Build me a map and a DBM to go with it...
95 Epetra_BlockMap BDMap(-1,NUM_BLOCKS,my_blockgids,my_blocksizes,0,Comm);
96 EpetraExt_BlockDiagMatrix BMAT(BDMap,true);
97
98 // Fill the matrix - This tests all three block size cases in the code, 1x1, 2x2 and larger.
99 for(int i=0;i<BMAT.NumMyBlocks();i++){
100 double*start=BMAT[i];
101 if(BMAT.BlockSize(i)==1)
102 *start=2.0;
103 else if(BMAT.BlockSize(i)==2){
104 start[0]=2.0;
105 start[1]=-1.0;
106 start[2]=-1.0;
107 start[3]=2.0;
108 }
109 else if(BMAT.BlockSize(i)==3){
110 start[0]=4.0;
111 start[1]=-1.0;
112 start[2]=-1.0;
113 start[3]=-1.0;
114 start[4]=4.0;
115 start[5]=-1.0;
116 start[2]=-1.0;
117 start[7]=-1.0;
118 start[8]=4.0;
119 }
120 else
121 exit(1);
122 }
123
124
125 // Allocations for tests
126 double norm2;
127 Epetra_MultiVector X(BDMap,1);
128 Epetra_MultiVector Y(BDMap,1);
129 Epetra_MultiVector Z(BDMap,1);
130 EpetraExt_BlockDiagMatrix BMAT_forward(BMAT);
131 EpetraExt_BlockDiagMatrix BMAT_factor(BMAT);
132 Teuchos::ParameterList List;
133
134 //***************************
135 // Test the Forward/Invert
136 //***************************
137 List.set("apply mode","invert");
138 BMAT.SetParameters(List);
139 X.SetSeed(24601); X.Random();
140 BMAT_forward.ApplyInverse(X,Y);
141 BMAT.Compute();
142 BMAT.ApplyInverse(Y,Z);
143 X.Update(1.0,Z,-1.0);
144 X.Norm2(&norm2);
145 if(!Comm.MyPID()) cout<<"Forward/Invert Error = "<<norm2<<endl;
146 if(norm2 > 1e-12) TestPassed=false;
147
148 //***************************
149 // Test the Forward/Factor
150 //***************************
151 List.set("apply mode","factor");
152 BMAT_factor.SetParameters(List);
153 X.SetSeed(24601); X.Random();
154 BMAT_forward.ApplyInverse(X,Y);
155 BMAT_factor.Compute();
156 BMAT_factor.ApplyInverse(Y,Z);
157 X.Update(1.0,Z,-1.0);
158 X.Norm2(&norm2);
159 if(!Comm.MyPID()) cout<<"Forward/Factor Error = "<<norm2<<endl;
160 if(norm2 > 1e-12) TestPassed=false;
161
162 return TestPassed;
163}
164
165//=============================================
166//=============================================
167//=============================================
168void Build_Local_Contiguous_Size2_BlockMatrix(const Epetra_Comm & Comm, int NUM_ROWS,int &NumBlocks,
169 int *&Blockstart_, int *&Blockids_,Epetra_CrsMatrix *&MAT){
170 double values[3]={-1.0,2.0,-1.0};
171 int indices[3];
172
173 // Build me a CrsMatrix
174 Epetra_Map Map(-1,NUM_ROWS,0,Comm);
175 MAT=new Epetra_CrsMatrix(Copy,Map,2);
176 assert(MAT->NumMyRows()%2==0);
177
178 NumBlocks=MAT->NumMyRows()/2;
179 Blockstart_=new int [NumBlocks+1];
180 Blockids_=new int [MAT->NumMyRows()];
181 Blockstart_[0]=0;
182 int curr_block=0;
183
184 for(int i=0;i<MAT->NumMyRows();i++){
185 // local contiguous blocks of constant size 2
186 int row_in_block=i%2;
187 if(row_in_block==0){
188 Blockstart_[curr_block]=i;
189 indices[0]=i; indices[1]=i+1;
190 Blockids_[i]=i;
191 MAT->InsertGlobalValues(Map.GID(i),2,&values[1],&indices[0]);
192 }
193 else if(row_in_block==1){
194 indices[0]=i-1; indices[1]=i;
195 MAT->InsertGlobalValues(Map.GID(i),2,&values[0],&indices[0]);
196 Blockids_[i]=i;
197 curr_block++;
198 }
199 }
200 Blockstart_[NumBlocks]=MAT->NumMyRows();
201
202 MAT->FillComplete();
203}
204
205//=============================================
206void Build_Local_NonContiguous_Size2_BlockMatrix(const Epetra_Comm & Comm, int NUM_ROWS,int &NumBlocks,
207 int *&Blockstart_, int *&Blockids_,Epetra_CrsMatrix *&MAT){
208 double values[3]={-1.0,2.0,-1.0};
209 int indices[3];
210
211 // Build me a CrsMatrix
212 Epetra_Map Map(-1,NUM_ROWS,0,Comm);
213 MAT=new Epetra_CrsMatrix(Copy,Map,2);
214 assert(MAT->NumMyRows()%2==0);
215
216 NumBlocks=MAT->NumMyRows()/2;
217 Blockstart_=new int [NumBlocks+1];
218 Blockids_=new int [MAT->NumMyRows()];
219 Blockstart_[0]=0;
220 int curr_block=0;
221
222 for(int i=0;i<MAT->NumMyRows();i++){
223 int row_in_block=(i%2)?1:0;
224 if(row_in_block==0){
225 int idx=i/2;
226 Blockstart_[curr_block]=i;
227 indices[0]=Map.GID(idx); indices[1]=Map.GID(idx+NumBlocks);
228 Blockids_[i]=idx;
229 MAT->InsertGlobalValues(Map.GID(idx),2,&values[1],&indices[0]);
230
231 }
232 else if(row_in_block==1){
233 int idx=(i-1)/2+NumBlocks;
234 indices[0]=Map.GID(idx-NumBlocks); indices[1]=Map.GID(idx);
235 MAT->InsertGlobalValues(Map.GID(idx),2,&values[0],&indices[0]);
236 Blockids_[i]=idx;
237 curr_block++;
238 }
239 }
240 Blockstart_[NumBlocks]=MAT->NumMyRows();
241
242 MAT->FillComplete();
243}
244//=============================================
245void Build_NonLocal_BlockMatrix(const Epetra_Comm & Comm, int NUM_ROWS,int &NumBlocks,
246 int *&Blockstart_, int *&Blockids_,Epetra_CrsMatrix *&MAT){
247 double values[3]={-1.0,2.0,-1.0};
248 int indices[3];
249 int NumProcs=Comm.NumProc();
250
251 // Build me a CrsMatrix
252 Epetra_Map Map(-1,NUM_ROWS,0,Comm);
253 MAT=new Epetra_CrsMatrix(Copy,Map,2);
254 int MyPID=Comm.MyPID();
255
256 for(int i=0;i<MAT->NumMyRows();i++){
257 int GID=Map.GID(i);
258 indices[0]=GID;
259 if(i==0) values[1]=2.0;
260 else values[1]=NumProcs+1;
261
262 MAT->InsertGlobalValues(GID,1,&values[1],&indices[0]);
263
264 // A little off-diagonal for good measure
265 if(NumProcs>1 && MyPID==0 && i>0){
266 indices[1]=GID;
267 for(int j=1;j<NumProcs;j++){
268 indices[1]+=NUM_ROWS;//HAQ
269 MAT->InsertGlobalValues(GID,1,&values[0],&indices[1]);
270 }
271 }
272 else if(NumProcs>1 && MyPID!=0 && i>0){
273 indices[1]=GID%NUM_ROWS;//HAQ
274 MAT->InsertGlobalValues(GID,1,&values[0],&indices[1]);
275 }
276 }
277 MAT->FillComplete();
278
279 // Build me some block structure
280 // The first row on each proc is a solo block. All others get blocked to
281 // PID=0's second block
282 if(MyPID==0) NumBlocks=NUM_ROWS;
283 else NumBlocks=1;
284 Blockstart_=new int [NumBlocks+1];
285 Blockstart_[0]=0;
286 int curr_idx,curr_block;
287
288 if(MyPID){
289 // PID > 0
290 Blockids_=new int[1];
291 Blockstart_[0]=0; Blockstart_[1]=1;
292 Blockids_[0]=Map.GID(0);
293 }
294 else{
295 // PID 0
296 int nnz=NumProcs*NumBlocks;
297 Blockids_=new int[nnz+1];
298 Blockstart_[0]=0;
299 Blockids_[0]=Map.GID(0);
300 curr_idx=1; curr_block=1;
301 for(int j=1;j<NUM_ROWS;j++){
302 Blockstart_[curr_block]=curr_idx;
303 for(int i=0;i<NumProcs;i++){
304 Blockids_[curr_idx]=NUM_ROWS*i+j;//FIX: THIS IS A HACK
305 curr_idx++;
306 }
307 curr_block++;
308 }
309 Blockstart_[curr_block]=curr_idx;
310 }
311}
312
313//=============================================
314void Build_DiagonalStructure(const Epetra_Map &Map,int &NumBlocks,int *&Blockstart_, int *&Blockids_,bool local_ids){
315 NumBlocks=Map.NumMyElements();
316 Blockstart_=new int[NumBlocks+1];
317 Blockids_=new int[NumBlocks];
318 for(int i=0;i<NumBlocks;i++){
319 Blockstart_[i]=i;
320 if(local_ids) Blockids_[i]=i;
321 else Blockids_[i]=Map.GID(i);
322 }
323 Blockstart_[NumBlocks]=NumBlocks;
324}
325
326
327
328//=============================================
329//=============================================
330//=============================================
331double Test_PTBDP(const Epetra_CrsMatrix& MAT, int NumBlocks,int* Blockstart_,int* Blockids_,bool is_lid){
332 // Build the block lists
333 Teuchos::ParameterList List,Sublist;
334 List.set("number of local blocks",NumBlocks);
335 List.set("block start index",Blockstart_);
336 if(is_lid) List.set("block entry lids",Blockids_);
337 else List.set("block entry gids",Blockids_);
338
339 Sublist.set("apply mode","invert");
340 //Sublist.set("apply mode","multiply");
341 List.set("blockdiagmatrix: list",Sublist);
342
344 Perm.SetParameters(List);
345
346 Perm.Compute();
347 Epetra_MultiVector X(MAT.RowMap(),1);
348 Epetra_MultiVector Y(MAT.RowMap(),1);
349 Epetra_MultiVector Z(MAT.RowMap(),1);
350 X.SetSeed(24601); X.Random();
351
352 double norm2;
353 Perm.ApplyInverse(X,Y);
354 MAT.Apply(Y,Z);
355 X.Update(1.0,Z,-1.0);
356 X.Norm2(&norm2);
357 return norm2;
358}
359
360
361//=============================================
362double Test_PTBDP_C(const Epetra_CrsMatrix& MAT,int BlockSize){
363 // Build the block lists
364 Teuchos::ParameterList List,Sublist;
365 List.set("contiguous block size",BlockSize);
366
367 Sublist.set("apply mode","invert");
368 //Sublist.set("apply mode","multiply");
369 List.set("blockdiagmatrix: list",Sublist);
370
372 Perm.SetParameters(List);
373
374 Perm.Compute();
375 Epetra_MultiVector X(MAT.RowMap(),1);
376 Epetra_MultiVector Y(MAT.RowMap(),1);
377 Epetra_MultiVector Z(MAT.RowMap(),1);
378 X.SetSeed(24601); X.Random();
379
380 double norm2;
381 Perm.ApplyInverse(X,Y);
382 MAT.Apply(Y,Z);
383 X.Update(1.0,Z,-1.0);
384 X.Norm2(&norm2);
385 return norm2;
386}
387
388//=============================================
389bool TestPointToBlockDiagPermute(const Epetra_Comm & Comm){
390 using std::cout;
391 using std::endl;
392
393 const int NUM_ROWS=64;
394
395 bool TestPassed=true;
396 Epetra_CrsMatrix *MAT;
397 int NumBlocks, *Blockstart_,*Blockids_;
398 double norm2;
399
400 // TEST #1 - Local, Contiguous
401 Build_Local_Contiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
402 norm2=Test_PTBDP(*MAT,NumBlocks,Blockstart_,Blockids_,true);
403 if(norm2 > 1e-12) TestPassed=false;
404 if(!Comm.MyPID()) cout<<"P2BDP LCMAT Error = "<<norm2<<endl;
405 delete MAT; delete [] Blockstart_; delete [] Blockids_;
406
407 // TEST #2 - Local, Non-Contiguous
408 Build_Local_NonContiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
409 norm2=Test_PTBDP(*MAT,NumBlocks,Blockstart_,Blockids_,true);
410 if(norm2 > 1e-12) TestPassed=false;
411 if(!Comm.MyPID()) cout<<"P2BDP LNCMat Error = "<<norm2<<endl;
412 delete MAT; delete [] Blockstart_; delete [] Blockids_;
413
414 // TEST #3 - Non-Local
415 Build_NonLocal_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
416 norm2=Test_PTBDP(*MAT,NumBlocks,Blockstart_,Blockids_,false);
417 if(norm2 > 1e-12) TestPassed=false;
418 if(!Comm.MyPID()) cout<<"P2BDP NLMat Error = "<<norm2<<endl;
419 delete MAT; delete [] Blockstart_; delete [] Blockids_;
420
421 // TEST #4 - Local, Contiguous in ContiguousMode
422 Build_Local_Contiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
423 norm2=Test_PTBDP_C(*MAT,2);
424 if(norm2 > 1e-12) TestPassed=false;
425 if(!Comm.MyPID()) cout<<"P2BDP LCMAT-C Error = "<<norm2<<endl;
426 delete MAT; delete [] Blockstart_; delete [] Blockids_;
427
428
429 return TestPassed;
430}
431
432
433//=============================================
434double Test_Cheby(const Epetra_CrsMatrix& MAT, int NumBlocks,int* Blockstart_,int* Blockids_,int maxits,bool is_lid){
435 double norm2,norm0;
436 // Build the block lists
437 Teuchos::ParameterList ChebyList,List,Sublist;
438 List.set("number of local blocks",NumBlocks);
439 List.set("block start index",Blockstart_);
440 if(is_lid) List.set("block entry lids",Blockids_);
441 else List.set("block entry gids",Blockids_);
442
443 Sublist.set("apply mode","invert");
444 List.set("blockdiagmatrix: list",Sublist);
445
446 ChebyList.set("chebyshev: use block mode",true);
447 ChebyList.set("chebyshev: block list",List);
448 ChebyList.set("chebyshev: eigenvalue autocompute ratio",30.0);//HAQ
449 ChebyList.set("chebyshev: degree",maxits);
450
451 // Build a Chebyshev
452 Ifpack_Chebyshev Cheby(&MAT);
453 Cheby.SetParameters(ChebyList);
454 Cheby.Compute();
455
456 Epetra_MultiVector X(MAT.RowMap(),1);
457 Epetra_MultiVector Y(MAT.RowMap(),1);
458 Epetra_MultiVector Z(MAT.RowMap(),1);
459 X.SetSeed(24601); X.Random();
460 MAT.Apply(X,Y);
461 Y.Norm2(&norm0);
462
463 Cheby.ApplyInverse(Y,Z);
464 X.Update(1.0,Z,-1.0);
465 X.Norm2(&norm2);
466 return norm2 / norm0;
467}
468
469//=============================================
470double Test_Cheby_C(const Epetra_CrsMatrix& MAT, int BlockSize,int maxits){
471 double norm2,norm0;
472 // Build the block lists
473 Teuchos::ParameterList ChebyList,List,Sublist;
474 List.set("contiguous block size",BlockSize);
475 Sublist.set("apply mode","invert");
476 List.set("blockdiagmatrix: list",Sublist);
477
478 ChebyList.set("chebyshev: use block mode",true);
479 ChebyList.set("chebyshev: block list",List);
480 ChebyList.set("chebyshev: eigenvalue autocompute ratio",30.0);//HAQ
481 ChebyList.set("chebyshev: degree",maxits);
482
483 // Build a Chebyshev
484 Ifpack_Chebyshev Cheby(&MAT);
485 Cheby.SetParameters(ChebyList);
486 Cheby.Compute();
487
488 Epetra_MultiVector X(MAT.RowMap(),1);
489 Epetra_MultiVector Y(MAT.RowMap(),1);
490 Epetra_MultiVector Z(MAT.RowMap(),1);
491 X.SetSeed(24601); X.Random();
492 MAT.Apply(X,Y);
493 Y.Norm2(&norm0);
494
495 Cheby.ApplyInverse(Y,Z);
496 X.Update(1.0,Z,-1.0);
497 X.Norm2(&norm2);
498 return norm2 / norm0;
499}
500
501
502//=============================================
503bool TestBlockChebyshev(const Epetra_Comm & Comm){
504 using std::cout;
505 using std::endl;
506
507 const int NUM_ROWS=100;
508
509 bool TestPassed=true;
510 Epetra_CrsMatrix *MAT;
511 int NumBlocks, *Blockstart_,*Blockids_;
512 double norm2;
513
514 // Test #1 - Local, Contiguous matrix w/ diagonal precond
515 Build_Local_Contiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
516 delete [] Blockstart_; delete [] Blockids_;
517 Build_DiagonalStructure(MAT->RowMap(),NumBlocks,Blockstart_,Blockids_,true);
518 norm2=Test_Cheby(*MAT,NumBlocks,Blockstart_,Blockids_,100,true);
519 if(norm2 > 1e-12) TestPassed=false;
520 if(!Comm.MyPID()) cout<<"Cheby LC-D nrm-red = "<<norm2<<endl;
521 delete MAT; delete [] Blockstart_; delete [] Blockids_;
522
523 // Test #2 - Local, Non-Contiguous matrix w/ diagonal precond
524 Build_Local_NonContiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
525 delete [] Blockstart_; delete [] Blockids_;
526 Build_DiagonalStructure(MAT->RowMap(),NumBlocks,Blockstart_,Blockids_,true);
527 norm2=Test_Cheby(*MAT,NumBlocks,Blockstart_,Blockids_,100,true);
528 if(norm2 > 1e-12) TestPassed=false;
529 if(!Comm.MyPID()) cout<<"Cheby LNC-D nrm-red = "<<norm2<<endl;
530 delete MAT; delete [] Blockstart_; delete [] Blockids_;
531
532 // TEST #3 - Non-Local matrix w/ diagonal precond
533 Build_NonLocal_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
534 delete [] Blockstart_; delete [] Blockids_;
535 Build_DiagonalStructure(MAT->RowMap(),NumBlocks,Blockstart_,Blockids_,false);
536 norm2=Test_Cheby(*MAT,NumBlocks,Blockstart_,Blockids_,100,false);
537 if(norm2 > 1e-12) TestPassed=false;
538 if(!Comm.MyPID()) cout<<"Cheby NL-D nrm-red = "<<norm2<<endl;
539 delete MAT; delete [] Blockstart_; delete [] Blockids_;
540
541 // Test #4 - Local, Contiguous matrix w/ exact precond
542 Build_Local_Contiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
543 norm2=Test_Cheby(*MAT,NumBlocks,Blockstart_,Blockids_,1,true);
544 if(norm2 > 1e-12) TestPassed=false;
545 if(!Comm.MyPID()) cout<<"Cheby LC-E nrm-red = "<<norm2<<endl;
546 delete MAT; delete [] Blockstart_; delete [] Blockids_;
547
548 // Test #5 - Local, Non-Contiguous matrix w/ exact precond
549 Build_Local_NonContiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
550 norm2=Test_Cheby(*MAT,NumBlocks,Blockstart_,Blockids_,1,true);
551 if(norm2 > 1e-12) TestPassed=false;
552 if(!Comm.MyPID()) cout<<"Cheby LNC-E nrm-red = "<<norm2<<endl;
553 delete MAT; delete [] Blockstart_; delete [] Blockids_;
554
555 // TEST #6 - Non-Local matrix w/ exact precond
556 Build_NonLocal_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
557 norm2=Test_Cheby(*MAT,NumBlocks,Blockstart_,Blockids_,1,false);
558 if(norm2 > 1e-12) TestPassed=false;
559 if(!Comm.MyPID()) cout<<"Cheby NL-E nrm-red = "<<norm2<<endl;
560 delete MAT; delete [] Blockstart_; delete [] Blockids_;
561
562 // Test #7 - Local, Contiguous matrix w/ diagonal precond (contiguous mode)
563 Build_Local_Contiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
564 norm2=Test_Cheby_C(*MAT,1,100);
565 if(norm2 > 1e-12) TestPassed=false;
566 if(!Comm.MyPID()) cout<<"Cheby LC-Dc nrm-red = "<<norm2<<endl;
567 delete MAT; delete [] Blockstart_; delete [] Blockids_;
568
569 // Test #8 - Local, Contiguous matrix w/ exact precond (contiguous mode)
570 Build_Local_Contiguous_Size2_BlockMatrix(Comm,NUM_ROWS,NumBlocks,Blockstart_,Blockids_,MAT);
571 norm2=Test_Cheby_C(*MAT,2,1);
572 if(norm2 > 1e-12) TestPassed=false;
573 if(!Comm.MyPID()) cout<<"Cheby LC-Ec nrm-red = "<<norm2<<endl;
574 delete MAT; delete [] Blockstart_; delete [] Blockids_;
575
576 return TestPassed;
577}
578
579//=============================================
580//=============================================
581//=============================================
582int main(int argc, char *argv[])
583{
584 using std::cout;
585 using std::endl;
586
587#ifdef HAVE_MPI
588 MPI_Init(&argc,&argv);
589 Epetra_MpiComm Comm( MPI_COMM_WORLD );
590#else
592#endif
593 bool TestPassed=true;
594
595 // BlockDiagMatrix test
596 TestPassed=TestPassed && TestBlockDiagMatrix(Comm);
597
598 // PointToBlockDiagPermute tests
599 TestPassed=TestPassed && TestPointToBlockDiagPermute(Comm);
600
601 // Block Chebyshev Tests
602 TestPassed=TestPassed && TestBlockChebyshev(Comm);
603
604 // ============ //
605 // final output //
606 // ============ //
607
608 if (!TestPassed) {
609 if(!Comm.MyPID()) cout << "Test `BlockCheby.exe' FAILED!" << endl;
610 exit(EXIT_FAILURE);
611 }
612
613#ifdef HAVE_MPI
614 MPI_Finalize();
615#endif
616 if(!Comm.MyPID()) cout << "Test `BlockCheby.exe' passed!" << endl;
617 exit(EXIT_SUCCESS);
618}
619
620#else
621
622#ifdef HAVE_MPI
623#include "Epetra_MpiComm.h"
624#else
625#include "Epetra_SerialComm.h"
626#endif
627
628int main(int argc, char *argv[])
629{
630
631#ifdef HAVE_MPI
632 MPI_Init(&argc,&argv);
633 Epetra_MpiComm Comm( MPI_COMM_WORLD );
634#else
636#endif
637
638 puts("please configure IFPACK with --eanble-aztecoo --enable-teuchos --enable-epetraext");
639 puts("to run this test");
640
641#ifdef HAVE_MPI
642 MPI_Finalize() ;
643#endif
644 return(EXIT_SUCCESS);
645}
646
647#endif
int GID(int LID) const
int NumMyElements() const
virtual int NumProc() const=0
virtual int MyPID() const=0
int MyPID() const
Ifpack_Chebyshev: class for preconditioning with Chebyshev polynomials in Ifpack.
int main(int argc, char *argv[])