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