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;
78 const int NUM_BLOCKS=30;
80 int my_blockgids[NUM_BLOCKS];
81 int my_blocksizes[NUM_BLOCKS];
83 for(
int i=0;i<NUM_BLOCKS;i++){
84 my_blockgids[i]=i+NUM_BLOCKS*Comm.
MyPID();
87 else if(i<2*NUM_BLOCKS/3)
99 for(
int i=0;i<BMAT.NumMyBlocks();i++){
100 double*start=BMAT[i];
101 if(BMAT.BlockSize(i)==1)
103 else if(BMAT.BlockSize(i)==2){
109 else if(BMAT.BlockSize(i)==3){
137 List.
set(
"apply mode",
"invert");
138 BMAT.SetParameters(List);
139 X.SetSeed(24601); X.Random();
140 BMAT_forward.ApplyInverse(X,Y);
142 BMAT.ApplyInverse(Y,Z);
143 X.Update(1.0,Z,-1.0);
145 if(!Comm.
MyPID()) cout<<
"Forward/Invert Error = "<<norm2<<endl;
146 if(norm2 > 1e-12) TestPassed=
false;
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);
159 if(!Comm.
MyPID()) cout<<
"Forward/Factor Error = "<<norm2<<endl;
160 if(norm2 > 1e-12) TestPassed=
false;
168 void Build_Local_Contiguous_Size2_BlockMatrix(
const Epetra_Comm & Comm,
int NUM_ROWS,
int &NumBlocks,
170 double values[3]={-1.0,2.0,-1.0};
179 Blockstart_=
new int [NumBlocks+1];
186 int row_in_block=i%2;
188 Blockstart_[curr_block]=i;
189 indices[0]=i; indices[1]=i+1;
193 else if(row_in_block==1){
194 indices[0]=i-1; indices[1]=i;
206 void Build_Local_NonContiguous_Size2_BlockMatrix(
const Epetra_Comm & Comm,
int NUM_ROWS,
int &NumBlocks,
208 double values[3]={-1.0,2.0,-1.0};
217 Blockstart_=
new int [NumBlocks+1];
223 int row_in_block=(i%2)?1:0;
226 Blockstart_[curr_block]=i;
227 indices[0]=Map.GID(idx); indices[1]=Map.GID(idx+NumBlocks);
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);
245 void Build_NonLocal_BlockMatrix(
const Epetra_Comm & Comm,
int NUM_ROWS,
int &NumBlocks,
247 double values[3]={-1.0,2.0,-1.0};
254 int MyPID=Comm.
MyPID();
259 if(i==0) values[1]=2.0;
260 else values[1]=NumProcs+1;
265 if(NumProcs>1 && MyPID==0 && i>0){
267 for(
int j=1;j<NumProcs;j++){
268 indices[1]+=NUM_ROWS;
272 else if(NumProcs>1 && MyPID!=0 && i>0){
273 indices[1]=GID%NUM_ROWS;
282 if(MyPID==0) NumBlocks=NUM_ROWS;
284 Blockstart_=
new int [NumBlocks+1];
286 int curr_idx,curr_block;
290 Blockids_=
new int[1];
291 Blockstart_[0]=0; Blockstart_[1]=1;
292 Blockids_[0]=Map.GID(0);
296 int nnz=NumProcs*NumBlocks;
297 Blockids_=
new int[nnz+1];
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;
309 Blockstart_[curr_block]=curr_idx;
314 void Build_DiagonalStructure(
const Epetra_Map &Map,
int &NumBlocks,
int *&Blockstart_,
int *&Blockids_,
bool local_ids){
316 Blockstart_=
new int[NumBlocks+1];
317 Blockids_=
new int[NumBlocks];
318 for(
int i=0;i<NumBlocks;i++){
320 if(local_ids) Blockids_[i]=i;
321 else Blockids_[i]=Map.
GID(i);
323 Blockstart_[NumBlocks]=NumBlocks;
331 double Test_PTBDP(
const Epetra_CrsMatrix& MAT,
int NumBlocks,
int* Blockstart_,
int* Blockids_,
bool is_lid){
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_);
339 Sublist.
set(
"apply mode",
"invert");
341 List.
set(
"blockdiagmatrix: list",Sublist);
344 Perm.SetParameters(List);
350 X.SetSeed(24601); X.Random();
353 Perm.ApplyInverse(X,Y);
355 X.Update(1.0,Z,-1.0);
365 List.
set(
"contiguous block size",BlockSize);
367 Sublist.
set(
"apply mode",
"invert");
369 List.
set(
"blockdiagmatrix: list",Sublist);
372 Perm.SetParameters(List);
378 X.SetSeed(24601); X.Random();
381 Perm.ApplyInverse(X,Y);
383 X.Update(1.0,Z,-1.0);
389 bool TestPointToBlockDiagPermute(
const Epetra_Comm & Comm){
393 const int NUM_ROWS=64;
395 bool TestPassed=
true;
397 int NumBlocks, *Blockstart_,*Blockids_;
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_;
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_;
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_;
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_;
434 double Test_Cheby(
const Epetra_CrsMatrix& MAT,
int NumBlocks,
int* Blockstart_,
int* Blockids_,
int maxits,
bool is_lid){
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_);
443 Sublist.
set(
"apply mode",
"invert");
444 List.
set(
"blockdiagmatrix: list",Sublist);
446 ChebyList.
set(
"chebyshev: use block mode",
true);
447 ChebyList.
set(
"chebyshev: block list",List);
448 ChebyList.
set(
"chebyshev: eigenvalue autocompute ratio",30.0);
449 ChebyList.
set(
"chebyshev: degree",maxits);
453 Cheby.SetParameters(ChebyList);
459 X.SetSeed(24601); X.Random();
463 Cheby.ApplyInverse(Y,Z);
464 X.Update(1.0,Z,-1.0);
466 return norm2 / norm0;
474 List.
set(
"contiguous block size",BlockSize);
475 Sublist.
set(
"apply mode",
"invert");
476 List.
set(
"blockdiagmatrix: list",Sublist);
478 ChebyList.
set(
"chebyshev: use block mode",
true);
479 ChebyList.
set(
"chebyshev: block list",List);
480 ChebyList.
set(
"chebyshev: eigenvalue autocompute ratio",30.0);
481 ChebyList.
set(
"chebyshev: degree",maxits);
485 Cheby.SetParameters(ChebyList);
491 X.SetSeed(24601); X.Random();
495 Cheby.ApplyInverse(Y,Z);
496 X.Update(1.0,Z,-1.0);
498 return norm2 / norm0;
507 const int NUM_ROWS=100;
509 bool TestPassed=
true;
511 int NumBlocks, *Blockstart_,*Blockids_;
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_;
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_;
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_;
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_;
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_;
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_;
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_;
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_;
582 int main(
int argc,
char *argv[])
588 MPI_Init(&argc,&argv);
593 bool TestPassed=
true;
596 TestPassed=TestPassed && TestBlockDiagMatrix(Comm);
599 TestPassed=TestPassed && TestPointToBlockDiagPermute(Comm);
602 TestPassed=TestPassed && TestBlockChebyshev(Comm);
609 if(!Comm.
MyPID()) cout <<
"Test `BlockCheby.exe' FAILED!" << endl;
616 if(!Comm.
MyPID()) cout <<
"Test `BlockCheby.exe' passed!" << endl;
623 #include "Epetra_MpiComm.h"
625 #include "Epetra_SerialComm.h"
628 int main(
int argc,
char *argv[])
632 MPI_Init(&argc,&argv);
638 puts(
"please configure IFPACK with --eanble-aztecoo --enable-teuchos --enable-epetraext");
639 puts(
"to run this test");
644 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
Ifpack_Chebyshev: class for preconditioning with Chebyshev polynomials in Ifpack. ...
int main(int argc, char *argv[])
virtual int NumProc() const =0