44 #if defined(HAVE_IFPACK_AZTECOO) && defined(HAVE_IFPACK_EPETRAEXT)
46 #include "Epetra_MpiComm.h"
48 #include "Epetra_SerialComm.h"
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"
60 #include "EpetraExt_BlockDiagMatrix.h"
61 #include "EpetraExt_PointToBlockDiagPermute.h"
64 using namespace Trilinos_Util;
70 static bool Solver = AZ_gmres;
79 const int NUM_BLOCKS=30;
81 long long my_blockgids[NUM_BLOCKS];
82 int my_blocksizes[NUM_BLOCKS];
84 for(
int i=0;i<NUM_BLOCKS;i++){
85 my_blockgids[i]=i+((
long long)NUM_BLOCKS)*Comm.
MyPID();
88 else if(i<2*NUM_BLOCKS/3)
96 Epetra_BlockMap BDMap(-1LL,NUM_BLOCKS,my_blockgids,my_blocksizes,0LL,Comm);
100 for(
int i=0;i<BMAT.NumMyBlocks();i++){
101 double*start=BMAT[i];
102 if(BMAT.BlockSize(i)==1)
104 else if(BMAT.BlockSize(i)==2){
110 else if(BMAT.BlockSize(i)==3){
138 List.
set(
"apply mode",
"invert");
139 BMAT.SetParameters(List);
140 X.SetSeed(24601); X.Random();
141 BMAT_forward.ApplyInverse(X,Y);
143 BMAT.ApplyInverse(Y,Z);
144 X.Update(1.0,Z,-1.0);
146 if(!Comm.
MyPID()) cout<<
"Forward/Invert Error = "<<norm2<<endl;
147 if(norm2 > 1e-12) TestPassed=
false;
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);
160 if(!Comm.
MyPID()) cout<<
"Forward/Factor Error = "<<norm2<<endl;
161 if(norm2 > 1e-12) TestPassed=
false;
169 void Build_Local_Contiguous_Size2_BlockMatrix(
const Epetra_Comm & Comm,
int NUM_ROWS,
int &NumBlocks,
171 double values[3]={-1.0,2.0,-1.0};
172 long long indices[3];
180 Blockstart_=
new int [NumBlocks+1];
187 int row_in_block=i%2;
189 Blockstart_[curr_block]=i;
190 indices[0]=i; indices[1]=i+1;
194 else if(row_in_block==1){
195 indices[0]=i-1; indices[1]=i;
207 void Build_Local_NonContiguous_Size2_BlockMatrix(
const Epetra_Comm & Comm,
int NUM_ROWS,
int &NumBlocks,
209 double values[3]={-1.0,2.0,-1.0};
210 long long indices[3];
218 Blockstart_=
new int [NumBlocks+1];
224 int row_in_block=(i%2)?1:0;
227 Blockstart_[curr_block]=i;
228 indices[0]=Map.GID64(idx); indices[1]=Map.GID64(idx+NumBlocks);
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);
246 void Build_NonLocal_BlockMatrix(
const Epetra_Comm & Comm,
int NUM_ROWS,
int &NumBlocks,
248 double values[3]={-1.0,2.0,-1.0};
249 long long indices[3];
255 int MyPID=Comm.
MyPID();
258 long long GID=Map.GID64(i);
260 if(i==0) values[1]=2.0;
261 else values[1]=NumProcs+1;
266 if(NumProcs>1 && MyPID==0 && i>0){
268 for(
int j=1;j<NumProcs;j++){
269 indices[1]+=NUM_ROWS;
273 else if(NumProcs>1 && MyPID!=0 && i>0){
274 indices[1]=GID%NUM_ROWS;
283 if(MyPID==0) NumBlocks=NUM_ROWS;
285 Blockstart_=
new int [NumBlocks+1];
287 int curr_idx,curr_block;
291 Blockids_=
new long long[1];
292 Blockstart_[0]=0; Blockstart_[1]=1;
293 Blockids_[0]=Map.GID64(0);
297 int nnz=NumProcs*NumBlocks;
298 Blockids_=
new long long[nnz+1];
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;
310 Blockstart_[curr_block]=curr_idx;
315 void Build_DiagonalStructure(
const Epetra_Map &Map,
int &NumBlocks,
int *&Blockstart_,
int *&Blockids_int_,
long long *&Blockids_LL_,
bool local_ids){
318 Blockids_int_=
new int[NumBlocks];
320 Blockids_LL_=
new long long[NumBlocks];
321 Blockstart_=
new int[NumBlocks+1];
322 for(
int i=0;i<NumBlocks;i++){
324 if(local_ids) Blockids_int_[i]=i;
325 else Blockids_LL_[i]=Map.
GID64(i);
327 Blockstart_[NumBlocks]=NumBlocks;
335 double Test_PTBDP(
const Epetra_CrsMatrix& MAT,
int NumBlocks,
int* Blockstart_,
int* Blockids_int_,
long long* Blockids_LL_,
bool is_lid){
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_);
343 Sublist.
set(
"apply mode",
"invert");
345 List.
set(
"blockdiagmatrix: list",Sublist);
348 Perm.SetParameters(List);
354 X.SetSeed(24601); X.Random();
357 Perm.ApplyInverse(X,Y);
359 X.Update(1.0,Z,-1.0);
369 List.
set(
"contiguous block size",BlockSize);
371 Sublist.
set(
"apply mode",
"invert");
373 List.
set(
"blockdiagmatrix: list",Sublist);
376 Perm.SetParameters(List);
382 X.SetSeed(24601); X.Random();
385 Perm.ApplyInverse(X,Y);
387 X.Update(1.0,Z,-1.0);
393 bool TestPointToBlockDiagPermute(
const Epetra_Comm & Comm){
397 const int NUM_ROWS=64;
399 bool TestPassed=
true;
401 int NumBlocks, *Blockstart_;
403 long long *Blockids_LL_;
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_;
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_;
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_;
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_;
440 double Test_Cheby(
const Epetra_CrsMatrix& MAT,
int NumBlocks,
int* Blockstart_,
int* Blockids_int_,
long long* Blockids_LL_,
int maxits,
bool is_lid){
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_);
449 Sublist.
set(
"apply mode",
"invert");
450 List.
set(
"blockdiagmatrix: list",Sublist);
452 ChebyList.
set(
"chebyshev: use block mode",
true);
453 ChebyList.
set(
"chebyshev: block list",List);
454 ChebyList.
set(
"chebyshev: eigenvalue autocompute ratio",30.0);
455 ChebyList.
set(
"chebyshev: degree",maxits);
459 Cheby.SetParameters(ChebyList);
465 X.SetSeed(24601); X.Random();
469 Cheby.ApplyInverse(Y,Z);
470 X.Update(1.0,Z,-1.0);
472 return norm2 / norm0;
480 List.
set(
"contiguous block size",BlockSize);
481 Sublist.
set(
"apply mode",
"invert");
482 List.
set(
"blockdiagmatrix: list",Sublist);
484 ChebyList.
set(
"chebyshev: use block mode",
true);
485 ChebyList.
set(
"chebyshev: block list",List);
486 ChebyList.
set(
"chebyshev: eigenvalue autocompute ratio",30.0);
487 ChebyList.
set(
"chebyshev: degree",maxits);
491 Cheby.SetParameters(ChebyList);
497 X.SetSeed(24601); X.Random();
501 Cheby.ApplyInverse(Y,Z);
502 X.Update(1.0,Z,-1.0);
504 return norm2 / norm0;
513 const int NUM_ROWS=100;
515 bool TestPassed=
true;
517 int NumBlocks, *Blockstart_;
519 long long *Blockids_LL_;
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_;
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_;
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_;
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_;
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_;
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_;
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_;
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_;
590 int main(
int argc,
char *argv[])
596 MPI_Init(&argc,&argv);
601 bool TestPassed=
true;
604 TestPassed=TestPassed && TestBlockDiagMatrix(Comm);
607 TestPassed=TestPassed && TestPointToBlockDiagPermute(Comm);
610 TestPassed=TestPassed && TestBlockChebyshev(Comm);
617 if(!Comm.
MyPID()) cout <<
"Test `BlockCheby.exe' FAILED!" << endl;
624 if(!Comm.
MyPID()) cout <<
"Test `BlockCheby.exe' passed!" << endl;
631 #include "Epetra_MpiComm.h"
633 #include "Epetra_SerialComm.h"
636 int main(
int argc,
char *argv[])
640 MPI_Init(&argc,&argv);
646 puts(
"please configure IFPACK with --eanble-aztecoo --enable-teuchos --enable-epetraext");
647 puts(
"to run this test");
652 return(EXIT_SUCCESS);
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
const Epetra_Map & RowMap() const
int NumMyElements() const
long long GID64(int LID) const
Ifpack_Chebyshev: class for preconditioning with Chebyshev polynomials in Ifpack. ...
int main(int argc, char *argv[])
virtual int NumProc() const =0