Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
52 #include "Epetra_LinearProblem.h"
53 #include "Trilinos_Util_CrsMatrixGallery.h"
54 #include "Teuchos_ParameterList.hpp"
55 #include "Ifpack_PointRelaxation.h"
56 #include "Ifpack_BlockRelaxation.h"
57 #include "Ifpack_DenseContainer.h"
58 #include "Ifpack_AdditiveSchwarz.h"
59 #include "AztecOO.h"
60 #include "EpetraExt_BlockDiagMatrix.h"
61 #include "EpetraExt_PointToBlockDiagPermute.h"
62 #include "Ifpack_Chebyshev.h"
63 
64 using namespace Trilinos_Util;
65 /*
66 (unused; commented out to avoid build warnings)
67 static bool verbose = false;
68 static bool SymmetricGallery = false;
69 static bool Solver = AZ_gmres;
70 */
71 
72 //=============================================
73 // Test the BlockDiagMatrix
74 bool 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 //=============================================
168 void 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 //=============================================
206 void 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 //=============================================
245 void 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 //=============================================
314 void 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 //=============================================
331 double 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 
343  EpetraExt_PointToBlockDiagPermute Perm(MAT);
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 //=============================================
362 double 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 
371  EpetraExt_PointToBlockDiagPermute Perm(MAT);
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 //=============================================
389 bool 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 //=============================================
434 double 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 //=============================================
470 double 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 //=============================================
503 bool 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 //=============================================
582 int 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
591  Epetra_SerialComm Comm;
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 
628 int 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
635  Epetra_SerialComm Comm;
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
Epetra_MultiVector
Ifpack_ConfigDefs.h
Epetra_Map
EpetraExt_PointToBlockDiagPermute
Ifpack_Chebyshev.h
Epetra_CrsMatrix::InsertGlobalValues
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_CrsMatrix::Apply
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Epetra_SerialComm::MyPID
int MyPID() const
Epetra_MpiComm
Epetra_Comm::MyPID
virtual int MyPID() const =0
Epetra_CrsMatrix::FillComplete
int FillComplete(bool OptimizeDataStorage=true)
EpetraExt_BlockDiagMatrix
Epetra_CrsMatrix::RowMap
const Epetra_Map & RowMap() const
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_Comm
Ifpack_BlockRelaxation.h
Ifpack_Chebyshev
Ifpack_Chebyshev: class for preconditioning with Chebyshev polynomials in Ifpack. ...
Definition: Ifpack_Chebyshev.h:110
Epetra_BlockMap::GID
int GID(int LID) const
Epetra_CrsMatrix::NumMyRows
int NumMyRows() const
Ifpack_AdditiveSchwarz.h
Epetra_BlockMap
main
int main(int argc, char *argv[])
Definition: Ifpack_ex_Amesos.cpp:61
Epetra_SerialComm
Epetra_Comm::NumProc
virtual int NumProc() const =0
Ifpack_DenseContainer.h
Epetra_CrsMatrix
Ifpack_PointRelaxation.h
Generated on Thu Nov 21 2024 09:19:18 for Ifpack Package Browser (Single Doxygen Collection) by   doxygen 1.8.5