Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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 #if 0
66 // Unused; commented out to avoid build warnings
67 
68 static bool verbose = false;
69 static bool SymmetricGallery = false;
70 static bool Solver = AZ_gmres;
71 #endif // 0
72 
73 //=============================================
74 // Test the BlockDiagMatrix
75 bool 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 //=============================================
169 void 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 //=============================================
207 void 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 //=============================================
246 void 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 //=============================================
315 void 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 //=============================================
335 double 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 
347  EpetraExt_PointToBlockDiagPermute Perm(MAT);
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 //=============================================
366 double 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 
375  EpetraExt_PointToBlockDiagPermute Perm(MAT);
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 //=============================================
393 bool 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 //=============================================
440 double 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 //=============================================
476 double 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 //=============================================
509 bool 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 //=============================================
590 int 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
599  Epetra_SerialComm Comm;
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 
636 int 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
643  Epetra_SerialComm Comm;
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
Epetra_MultiVector
Ifpack_ConfigDefs.h
Epetra_Map
EpetraExt_PointToBlockDiagPermute
Ifpack_Chebyshev.h
verbose
bool verbose
Definition: test/CompareWithAztecOO/cxx_main.cpp:66
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Epetra_CrsMatrix::InsertGlobalValues
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_CrsMatrix::Apply
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
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
Solver
static bool Solver
Definition: performance.cpp:70
Epetra_CrsMatrix::RowMap
const Epetra_Map & RowMap() const
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_BlockMap::GID64
long long GID64(int LID) const
Epetra_Comm
Ifpack_BlockRelaxation.h
Ifpack_Chebyshev
Ifpack_Chebyshev: class for preconditioning with Chebyshev polynomials in Ifpack. ...
Definition: Ifpack_Chebyshev.h:104
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
SymmetricGallery
static bool SymmetricGallery
Definition: performance.cpp:69
Epetra_CrsMatrix
Ifpack_PointRelaxation.h
Generated on Thu Apr 18 2024 09:19:03 for Ifpack Package Browser (Single Doxygen Collection) by   doxygen 1.8.5