68 #include "../epetra_test_err.h" 69 #include "../src/Epetra_matrix_data.h" 74 int checkValues(
double x,
double y,
string message =
"",
bool verbose =
false) {
75 if (fabs((x-y)/x) > 0.01) {
77 if (verbose) cout <<
"********** " << message <<
" check failed.********** " << endl;
80 if (verbose) cout << message <<
" check OK." << endl;
89 int globalbadvalue = 0;
90 for (
int j=0; j<numVectors; j++)
91 for (
int i=0; i< length; i++)
96 if (globalbadvalue==0) cout << message <<
" check OK." << endl;
97 else cout <<
"********* " << message <<
" check failed.********** " << endl;
99 return(globalbadvalue);
106 int CompareValues(
double *
A,
int LDA,
int NumRowsA,
int NumColsA,
107 double *
B,
int LDB,
int NumRowsB,
int NumColsB);
110 int NumMyRows1,
int NumGlobalRows1,
int NumMyNonzeros1,
int NumGlobalNonzeros1,
111 int NumMyBlockRows1,
int NumGlobalBlockRows1,
int NumMyBlockNonzeros1,
int NumGlobalBlockNonzeros1,
112 int * MyGlobalElements,
bool verbose);
118 double * lambda,
int niters,
double tolerance,
143 vector<int> MyGlobalElements( NumMyElements );
147 vector<int> MyGlobalColumns( NumMyColumns );
152 vector<int> LocalColumnIndices(MaxNumIndices);
153 vector<int> GlobalColumnIndices(MaxNumIndices);
154 vector<double> MatrixValues(MaxNumIndices);
156 for(
int LocalRow=0; LocalRow<NumMyElements; ++LocalRow ) {
162 &LocalColumnIndices[0] );
164 for (
int j = 0 ; j < NumIndices ; j++ )
166 GlobalColumnIndices[j] = MyGlobalColumns[ LocalColumnIndices[j] ] ;
173 &GlobalColumnIndices[0] )!=0)abort();
178 &LocalColumnIndices[0] )!= 0) abort();
201 assert( (vecX && vecY) || (!vecX && !vecY) ) ;
203 if ( vecX && vecY ) {
204 double normY, NormError;
209 A.Multiply1( transpose, *vecX, Y ) ;
211 Error.
Update( 1.0, Y, -1.0, *vecY, 0.0 ) ;
214 if ( NormError / normY > 1e-13 ) {
226 A.Multiply1( transpose, Z, Z ) ;
227 Error.
Update( 1.0, Z, -1.0, *vecY, 0.0 ) ;
231 if ( NormError / normY > 1e-13 ) numerrors++;
238 A.Multiply( transpose, X, Y ) ;
241 vector<double> NormError(NumVecs);
242 vector<double> Normy(NumVecs);
246 Error.
Update( 1.0, Y, -1.0, Check_Y, 0.0 ) ;
247 Error.
NormInf( &NormError[0] ) ;
249 bool LoopError = false ;
250 for (
int ii = 0 ; ii < NumVecs ; ii++ ) {
251 if( NormError[ii] / Normy[ii] > 1e-13 ) {
296 int NumMyElements,
int MinSize,
int MaxSize,
297 ConsType ConstructorType,
bool ExtraBlocks,
304 cout <<
"MinSize = " << MinSize << endl
305 <<
"MaxSize = " << MaxSize << endl
306 <<
"ConstructorType = " << ConstructorType << endl
307 <<
"ExtraBlocks = " << ExtraBlocks << endl
308 <<
"insertlocal = " << insertlocal << endl
309 <<
"symmetric = " << symmetric << endl
310 <<
"PreviousA = " << PreviousA << endl;
312 int ierr = 0, forierr = 0;
313 int MyPID = Comm.
MyPID();
314 if (MyPID < 3) NumMyElements++;
315 if (NumMyElements<2) NumMyElements = 2;
318 Epetra_Map randmap(-1, NumMyElements, 0, Comm);
321 int * ElementSizeList =
new int[NumMyElements];
322 int SizeRange = MaxSize - MinSize + 1;
323 double DSizeRange = SizeRange;
328 if ( *PreviousA != 0 ) {
330 rowmap = &((*PreviousA)->RowMap());
331 colmap = &((*PreviousA)->ColMap());
334 ElementSizeList[0] = MaxSize;
335 for (
int i=1; i<NumMyElements-1; i++) {
336 int curSize = MinSize + (int) (DSizeRange * fabs(randvec[i]) + .99);
339 ElementSizeList[NumMyElements-1] = MaxSize;
354 Epetra_BlockMap Map (-1, NumMyElements, randMyGlobalElements, ElementSizeList, 0, Comm);
363 int * NumNz =
new int[NumMyElements];
367 for (
int i=0; i<NumMyElements; i++)
368 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
374 bool FixedNumEntries = false ;
376 bool HaveColMap = false ;
377 bool FixedBlockSize = ( MinSize == MaxSize ) ;
378 bool HaveGraph = false ;
383 switch( ConstructorType ) {
389 FixedNumEntries =
true;
400 FixedNumEntries =
true;
421 if ( insertlocal ) assert( HaveColMap );
422 if ( ExtraBlocks ) assert( FixedBlockSize );
423 if ( ExtraBlocks ) assert( ! HaveColMap );
424 if ( FixedNumEntries) assert( FixedBlockSize ) ;
425 if ( insertlocal && HaveGraph ) assert( ! FixedNumEntries ) ;
426 if ( insertlocal && HaveGraph ) assert( ! ExtraBlocks ) ;
440 for (
int kr=0; kr<SizeRange; kr++) {
442 int RowDim = MinSize+kr;
443 for (
int kc = 0; kc<SizeRange; kc++) {
444 int ColDim = MinSize+kc;
446 curmat->
Shape(RowDim,ColDim);
447 for (
int j=0; j < ColDim; j++)
448 for (
int i=0; i < RowDim; i++) {
449 BlockEntries[kr][kc][j][i] = -1.0;
450 if (i==j && kr==kc) BlockEntries[kr][kc][j][i] = 9.0;
451 else BlockEntries[kr][kc][j][i] = -1.0;
453 if ( ! symmetric ) BlockEntries[kr][kc][j][i] += ((double) j)/10000.0;
460 int *Indices =
new int[3];
461 int *MyIndices =
new int[3];
462 int *ColDims =
new int[3];
464 int NumMyNonzeros = 0, NumMyEquations = 0;
468 for (
int i=0; i<NumMyElements; i++) {
469 int CurRow = MyGlobalElements[i];
471 assert ( i == rowmap->
LID( CurRow ) ) ;
472 assert ( CurRow == rowmap->
GID( i ) ) ;
474 int RowDim = ElementSizeList[i]-MinSize;
475 NumMyEquations += BlockEntries[RowDim][RowDim].M();
480 Indices[1] = CurRow+1;
481 if ( FixedNumEntries ) {
482 Indices[2] = NumGlobalElements-1;
484 assert( ElementSizeList[i] == MinSize );
489 ColDims[0] = ElementSizeList[i] - MinSize;
490 ColDims[1] = ElementSizeList[i+1] - MinSize;
492 else if (CurRow == NumGlobalElements-1)
494 Indices[0] = CurRow-1;
496 if ( FixedNumEntries ) {
499 assert( ElementSizeList[i] == MinSize );
504 ColDims[0] = ElementSizeList[i-1] - MinSize;
505 ColDims[1] = ElementSizeList[i] - MinSize;
508 Indices[0] = CurRow-1;
510 Indices[2] = CurRow+1;
512 if (i==0) ColDims[0] = MaxSize - MinSize;
513 else ColDims[0] = ElementSizeList[i-1] - MinSize;
514 ColDims[1] = ElementSizeList[i] - MinSize;
516 if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
517 else ColDims[2] = ElementSizeList[i+1] - MinSize;
520 for (
int ii=0; ii < NumEntries; ii++ )
521 MyIndices[ii] = colmap->
LID( Indices[ii] ) ;
525 EPETRA_TEST_ERR(!(
A->BeginReplaceMyValues(rowmap->
LID(CurRow), NumEntries, MyIndices)==0),ierr);
527 EPETRA_TEST_ERR(!(
A->BeginInsertMyValues(rowmap->
LID(CurRow), NumEntries, MyIndices)==0),ierr);
532 EPETRA_TEST_ERR(!(
A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
541 EPETRA_TEST_ERR(!(
A->BeginInsertGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
545 for (
int j=0; j < NumEntries; j++) {
547 NumMyNonzeros += AD->
M() * AD->
N();
548 forierr += !(
A->SubmitBlockEntry(AD->
A(), AD->
LDA(), AD->
M(), AD->
N())==0);
552 A->EndSubmitEntries();
555 int NumMyBlockEntries = 3*NumMyElements;
556 if ( ! FixedNumEntries ) {
557 if (
A->LRID(0)>=0) NumMyBlockEntries--;
558 if (
A->LRID(NumGlobalElements-1)>=0) NumMyBlockEntries--;
574 const int NumTries = 100;
581 BlockIindex.
Scale( 1.0 * NumMyElements );
582 BlockJindex.
Scale( 1.0 *
A->NumGlobalBlockCols() );
583 bool OffDiagonalBlockAdded = false ;
584 for (
int ii=0; ii < NumTries && ! OffDiagonalBlockAdded; ii++ ) {
585 int i = (int) BlockIindex[0][ii];
586 int j = (int) BlockJindex[0][ii];
587 if ( i < 0 ) i = - i ;
588 if ( j < 0 ) j = - j ;
591 assert( i < NumMyElements ) ;
592 assert( j < A->NumGlobalBlockCols() ) ;
593 int CurRow = MyGlobalElements[i];
597 EPETRA_TEST_ERR(!(
A->BeginInsertGlobalValues( CurRow, 1, IndicesL)==0),ierr);
599 if ( CurRow < j-1 || CurRow > j+1 ) {
600 OffDiagonalBlockAdded = true ;
601 NumMyNonzeros += AD->
M() * AD->
N();
602 NumMyBlockEntries++ ;
607 A->EndSubmitEntries();
612 if ( ! HaveGraph && ! insertlocal )
619 A->OptimizeStorage();
620 if ( FixedBlockSize ) {
634 int NumGlobalBlockEntries ;
635 int NumGlobalNonzeros, NumGlobalEquations;
636 Comm.
SumAll(&NumMyBlockEntries, &NumGlobalBlockEntries, 1);
637 Comm.
SumAll(&NumMyNonzeros, &NumGlobalNonzeros, 1);
639 Comm.
SumAll(&NumMyEquations, &NumGlobalEquations, 1);
641 if (! ExtraBlocks ) {
642 if ( FixedNumEntries ) {
643 EPETRA_TEST_ERR( !( NumGlobalBlockEntries == (3*NumGlobalElements)), ierr );
645 EPETRA_TEST_ERR( !( NumGlobalBlockEntries == (3*NumGlobalElements-2)), ierr );
651 NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries,
652 MyGlobalElements, verbose)==0),ierr);
655 if ( ! ExtraBlocks ) {
656 if ( FixedNumEntries )
657 for (
int i=0; i<NumMyElements; i++) forierr += !(
A->NumGlobalBlockEntries(MyGlobalElements[i])==3);
659 for (
int i=0; i<NumMyElements; i++) forierr += !(
A->NumGlobalBlockEntries(MyGlobalElements[i])==NumNz[i]);
663 if ( ! ExtraBlocks ) {
664 if ( FixedNumEntries )
665 for (
int i=0; i<NumMyElements; i++) forierr += !(
A->NumMyBlockEntries(i)==3);
667 for (
int i=0; i<NumMyElements; i++) forierr += !(
A->NumMyBlockEntries(i)==NumNz[i]);
671 if (verbose) cout <<
"\n\nNumEntries function check OK" << endl<< endl;
691 double tolerance = 1.0e-3;
697 A->SetFlopCounter(flopcounter);
703 int ierr1 =
power_method(
false, *
A, q, z, resid, &lambda, niters, tolerance, verbose);
705 double total_flops = flopcounter.
Flops();
706 double MFLOPs = total_flops/elapsed_time/1000000.0;
708 if (verbose) cout <<
"\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
709 if (verbose && ierr1==1) cout <<
"***** Power Method did not converge. *****" << endl << endl;
715 if (verbose) cout <<
"\n\nUsing transpose of matrix and solving again (should give same result).\n\n" 722 ierr1 =
power_method(
true, *
A, q, z, resid, &lambda, niters, tolerance, verbose);
724 total_flops = flopcounter.
Flops();
725 MFLOPs = total_flops/elapsed_time/1000000.0;
727 if (verbose) cout <<
"\n\nTotal MFLOPs for transpose solve = " << MFLOPs << endl<< endl;
728 if (verbose && ierr1==1) cout <<
"***** Power Method did not converge. *****" << endl << endl;
737 if (verbose) cout <<
"\n\nIncreasing the magnitude of first diagonal term and solving again\n\n" 740 double AnormInf = -13 ;
741 double AnormOne = -13 ;
742 AnormInf =
A->NormInf( ) ;
743 AnormOne =
A->NormOne( ) ;
748 if (
A->MyGlobalBlockRow(0)) {
749 int numvals =
A->NumGlobalBlockEntries(0);
751 int* Rowinds =
new int[numvals];
753 A->ExtractGlobalBlockRowPointers(0, numvals, RowDim, numvals, Rowinds,
756 for (
int i=0; i<numvals; i++) {
757 if (Rowinds[i] == 0) {
759 Rowvals[i]->
A()[0] += 1000.0;
787 ierr1 =
power_method(
false, *
A, q, z, resid, &lambda, niters, tolerance, verbose);
789 total_flops = flopcounter.
Flops();
790 MFLOPs = total_flops/elapsed_time/1000000.0;
792 if (verbose) cout <<
"\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
793 if (verbose && ierr1==1) cout <<
"***** Power Method did not converge. *****" << endl << endl;
803 if (verbose) cout <<
"\n\nUsing transpose of matrix and solving again (should give same result).\n\n" 811 ierr1 =
power_method(
true, *
A, q, z, resid, &lambda, niters, tolerance, verbose);
813 total_flops = flopcounter.
Flops();
814 MFLOPs = total_flops/elapsed_time/1000000.0;
816 if (verbose) cout <<
"\n\nTotal MFLOPs for tranpose of second solve = " << MFLOPs << endl<< endl;
817 if (verbose && ierr1==1) cout <<
"***** Power Method did not converge. *****" << endl << endl;
822 if (verbose) cout <<
"\n\n*****Comparing against CrsMatrix " << endl<< endl;
841 CrsA->
Multiply(
false, CrsX, CrsY ) ;
842 OrigCrsA->
Multiply(
false, CrsX, OrigCrsY ) ;
846 orig_check_y = OrigCrsY ;
847 CrsA->
Multiply(
true, CrsX, CrsY ) ;
848 check_ytranspose = CrsY ;
859 if (verbose) cout <<
"\n\n*****Try the Apply method " << endl<< endl;
862 A->Apply( CrsX, Y_Apply ) ;
863 Apply_check_y = Y_Apply ;
866 if (verbose) cout <<
"\n\n*****Multiply multivectors " << endl<< endl;
868 const int NumVecs = 4 ;
882 CrsA->
Multiply(
false, CrsMX, CrsMY ) ;
885 CrsA->
Multiply(
false, CrsMY, CrsMY ) ;
889 CrsA->
Multiply(
true, CrsMY, CrsMY ) ;
890 check_mytranspose = CrsMY ;
900 if (verbose) cout <<
"\n\n*****Testing copy constructor" << endl<< endl;
909 NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries,
910 MyGlobalElements, verbose)==0),ierr);
920 NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries,
921 MyGlobalElements, verbose)==0),ierr);
927 AnormInf =
A->NormInf( );
928 AnormOne =
A->NormOne( );
936 if (verbose) cout <<
"\n\n*****Testing PutScalar, LeftScale, RightScale, and ReplaceDiagonalValues" << endl<< endl;
940 B->PutScalar( 1.0 ) ;
959 double B_norm_frob =
B->NormFrobenius();
966 if (fabs(B_norm_frob-CrsB_norm_frob) > 5.e-5) {
967 std::cout <<
"fabs(B_norm_frob-CrsB_norm_frob): " 968 << fabs(B_norm_frob-CrsB_norm_frob) << std::endl;
969 std::cout <<
"VbrMatrix::NormFrobenius test FAILED."<<std::endl;
972 if (verbose) std::cout <<
"\n\nVbrMatrix::NormFrobenius OK"<<std::endl<<std::endl;
976 if (verbose) cout <<
"\n\n*****Testing post construction modifications" << endl<< endl;
977 if (verbose) cout <<
"\n\n*****Replace methods should be OK" << endl<< endl;
989 NumMyNonzeros = NumMyEquations = 0;
991 for (
int i=0; i<NumMyElements; i++) {
992 int CurRow = MyGlobalElements[i];
993 int RowDim = ElementSizeList[i]-MinSize;
994 NumMyEquations += BlockEntries[RowDim][RowDim].M();
999 Indices[1] = CurRow+1;
1001 ColDims[0] = ElementSizeList[i] - MinSize;
1002 ColDims[1] = ElementSizeList[i+1] - MinSize;
1004 else if (CurRow == NumGlobalElements-1)
1006 Indices[0] = CurRow-1;
1007 Indices[1] = CurRow;
1009 ColDims[0] = ElementSizeList[i-1] - MinSize;
1010 ColDims[1] = ElementSizeList[i] - MinSize;
1013 Indices[0] = CurRow-1;
1014 Indices[1] = CurRow;
1015 Indices[2] = CurRow+1;
1017 if (i==0) ColDims[0] = MaxSize - MinSize;
1018 else ColDims[0] = ElementSizeList[i-1] - MinSize;
1019 ColDims[1] = ElementSizeList[i] - MinSize;
1021 if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
1022 else ColDims[2] = ElementSizeList[i+1] - MinSize;
1026 EPETRA_TEST_ERR(!(
A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
1028 for (
int j=0; j < NumEntries; j++) {
1030 NumMyNonzeros += AD->
M() * AD->
N();
1031 forierr += !(
A->SubmitBlockEntry(AD->
A(), AD->
LDA(), AD->
M(), AD->
N())==0);
1035 A->EndSubmitEntries();
1043 if ( ! ExtraBlocks ) {
1046 NumMyNonzeros = NumMyEquations = 0;
1048 for (
int i=0; i<NumMyElements; i++) {
1049 int CurRow = MyGlobalElements[i];
1050 int RowDim = ElementSizeList[i]-MinSize;
1051 NumMyEquations += BlockEntries[RowDim][RowDim].M();
1055 Indices[0] = CurRow;
1056 Indices[1] = CurRow+1;
1058 ColDims[0] = ElementSizeList[i] - MinSize;
1059 ColDims[1] = ElementSizeList[i+1] - MinSize;
1061 else if (CurRow == NumGlobalElements-1)
1063 Indices[0] = CurRow-1;
1064 Indices[1] = CurRow;
1066 ColDims[0] = ElementSizeList[i-1] - MinSize;
1067 ColDims[1] = ElementSizeList[i] - MinSize;
1070 Indices[0] = CurRow-1;
1071 Indices[1] = CurRow;
1072 Indices[2] = CurRow+1;
1074 if (i==0) ColDims[0] = MaxSize - MinSize;
1075 else ColDims[0] = ElementSizeList[i-1] - MinSize;
1076 ColDims[1] = ElementSizeList[i] - MinSize;
1078 if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
1079 else ColDims[2] = ElementSizeList[i+1] - MinSize;
1081 if ( insertlocal ) {
1082 for (
int ii=0; ii < NumEntries; ii++ )
1083 MyIndices[ii] = colmap->
LID( Indices[ii] ) ;
1084 EPETRA_TEST_ERR(!(
A->BeginSumIntoMyValues( rowmap->
LID(CurRow), NumEntries, MyIndices)==0),ierr);
1086 EPETRA_TEST_ERR(!(
A->BeginSumIntoGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
1089 for (
int j=0; j < NumEntries; j++) {
1091 NumMyNonzeros += AD->
M() * AD->
N();
1093 if ( insertlocal ) {
1094 forierr += !(
A->SubmitBlockEntry( *AD )==0);
1096 forierr += !(
A->SubmitBlockEntry(AD->
A(), AD->
LDA(), AD->
M(), AD->
N())==0);
1101 A->EndSubmitEntries();
1104 orig_check_y.
Scale(2.0) ;
1107 if ( ! FixedNumEntries )
1112 {
for (
int kr=0; kr<SizeRange; kr++)
delete [] BlockEntries[kr];}
1113 delete [] BlockEntries;
1115 delete [] MyIndices ;
1117 delete [] ElementSizeList;
1122 if (verbose) cout <<
"\n\n*****Insert methods should not be accepted" << endl<< endl;
1125 if (
B->MyGRID(0))
EPETRA_TEST_ERR(!(
B->BeginInsertGlobalValues(0, 1, &One)==-2),ierr);
1129 int NumMyEquations1 =
B->NumMyRows();
1135 for (
int i=0; i<NumMyEquations1; i++) checkDiag[i]=two1;
1144 for (
int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
1147 if (verbose) cout <<
"\n\nDiagonal extraction and replacement OK.\n\n" << endl;
1149 double originfnorm =
B->NormInf();
1150 double origonenorm =
B->NormOne();
1156 if (verbose) cout <<
"\n\nMatrix scale OK.\n\n" << endl;
1194 MPI_Init(&argc,&argv);
1200 int MyPID = Comm.
MyPID();
1202 bool verbose =
false;
1205 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
1216 if (verbose && MyPID==0)
1219 if (verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
1220 <<
" is alive."<<endl;
1223 if (verbose && Comm.
MyPID()!=0) verbose =
false;
1226 int NumMyElements = 400;
1230 bool NoExtraBlocks =
false;
1231 bool symmetric =
true;
1232 bool NonSymmetric =
false;
1233 bool NoInsertLocal = false ;
1234 bool InsertLocal = true ;
1238 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 1, 1,
VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1244 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize,
VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1246 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1248 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize,
FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1250 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1252 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
RowMapColMap_VEPR, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1254 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
RowMapColMap_NEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1256 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
RowMapColMap_VEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1258 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1260 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
WithGraph, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1261 assert ( PreviousA == 0 ) ;
1267 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1269 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1271 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize,
WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1272 assert ( PreviousA == 0 ) ;
1274 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1276 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4,
RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1278 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4,
WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1279 assert ( PreviousA == 0 ) ;
1283 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize,
FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1285 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize,
RowMapColMap_FEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1289 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1291 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 3, 3,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1293 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1295 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 5, 5,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1297 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 6, 6,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1299 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 7, 7,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1301 ierr +=
TestMatrix( Comm, verbose, debug, NumMyElements, 8, 8,
NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1399 if (ierr==0) cout <<
"All VbrMatrix tests OK" << endl;
1400 else cout << ierr <<
" errors encountered." << endl;
1416 double * lambda,
int niters,
double tolerance,
1420 double normz, residual;
1424 for (
int iter = 0; iter < niters; iter++)
1427 q.
Scale(1.0/normz, z);
1428 A.Multiply(TransA, q, z);
1430 if (iter%100==0 || iter+1==niters)
1432 resid.
Update(1.0, z, -(*lambda), q, 0.0);
1433 resid.
Norm2(&residual);
1434 if (verbose) cout <<
"Iter = " << iter <<
" Lambda = " << *lambda
1435 <<
" Residual of A*q - lambda*q = " << residual << endl;
1437 if (residual < tolerance) {
1445 int NumMyRows1,
int NumGlobalRows1,
int NumMyNonzeros1,
int NumGlobalNonzeros1,
1446 int NumMyBlockRows1,
int NumGlobalBlockRows1,
int NumMyBlockNonzeros1,
int NumGlobalBlockNonzeros1,
1447 int * MyGlobalElements,
bool verbose) {
1449 int ierr = 0, forierr = 0;
1452 int NumMyRows =
A.NumMyRows();
1453 if (verbose) cout <<
"\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
1456 if (verbose) cout <<
"\n\nNumber of local Rows is = " << NumMyRows << endl<< endl;
1457 if (verbose) cout <<
"\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
1461 int NumMyNonzeros =
A.NumMyNonzeros();
1462 if (verbose) cout <<
"\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
1466 if ( NumMyNonzeros != NumMyNonzeros1 ) {
1467 cout <<
" MyPID = " <<
A.Comm().MyPID()
1468 <<
" NumMyNonzeros = " << NumMyNonzeros
1469 <<
" NumMyNonzeros1 = " << NumMyNonzeros1 << endl ;
1474 int NumGlobalRows =
A.NumGlobalRows();
1475 if (verbose) cout <<
"\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
1479 int NumGlobalNonzeros =
A.NumGlobalNonzeros();
1480 if (verbose) cout <<
"\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
1482 if ( NumGlobalNonzeros != NumGlobalNonzeros1 ) {
1483 cout <<
" MyPID = " <<
A.Comm().MyPID()
1484 <<
" NumGlobalNonzeros = " << NumGlobalNonzeros
1485 <<
" NumGlobalNonzeros1 = " << NumGlobalNonzeros1 << endl ;
1489 int NumMyBlockRows =
A.NumMyBlockRows();
1490 if (verbose) cout <<
"\n\nNumber of local Block Rows = " << NumMyBlockRows << endl<< endl;
1494 int NumMyBlockNonzeros =
A.NumMyBlockEntries();
1496 if (verbose) cout <<
"\n\nNumber of local Nonzero Block entries = " << NumMyBlockNonzeros << endl<< endl;
1497 if (verbose) cout <<
"\n\nNumber of local Nonzero Block entries 1 = " << NumMyBlockNonzeros1 << endl<< endl;
1501 int NumGlobalBlockRows =
A.NumGlobalBlockRows();
1502 if (verbose) cout <<
"\n\nNumber of global Block Rows = " << NumGlobalBlockRows << endl<< endl;
1506 int NumGlobalBlockNonzeros =
A.NumGlobalBlockEntries();
1507 if (verbose) cout <<
"\n\nNumber of global Nonzero Block entries = " << NumGlobalBlockNonzeros << endl<< endl;
1509 EPETRA_TEST_ERR(!(NumGlobalBlockNonzeros==NumGlobalBlockNonzeros1),ierr);
1513 int RowDim, NumBlockEntries, * BlockIndices;
1516 A.ExtractMyBlockRowView(NumMyBlockRows-1, RowDim, NumBlockEntries,
1517 BlockIndices, Values);
1518 int NumMyEntries1 = 0;
1519 {
for (
int i=0; i < NumBlockEntries; i++) NumMyEntries1 += Values[i]->N();}
1521 A.NumMyRowEntries(NumMyRows-1, NumMyEntries);
1523 cout <<
"\n\nNumber of nonzero values in last row = " 1524 << NumMyEntries << endl<< endl;
1543 int MaxNumBlockEntries =
A.MaxNumBlockEntries();
1548 int MyPointersRowDim, MyPointersNumBlockEntries;
1549 int * MyPointersBlockIndices =
new int[MaxNumBlockEntries];
1552 int GlobalPointersRowDim, GlobalPointersNumBlockEntries;
1553 int * GlobalPointersBlockIndices =
new int[MaxNumBlockEntries];
1559 int MyCopyRowDim, MyCopyNumBlockEntries;
1560 int * MyCopyBlockIndices =
new int[MaxNumBlockEntries];
1561 int * MyCopyColDims =
new int[MaxNumBlockEntries];
1562 int * MyCopyLDAs =
new int[MaxNumBlockEntries];
1563 int MaxRowDim =
A.MaxRowDim();
1564 int MaxColDim =
A.MaxColDim();
1565 int MyCopySizeOfValues = MaxRowDim*MaxColDim;
1566 double ** MyCopyValuesPointers =
new double*[MaxNumBlockEntries];
1567 for (
int i=0; i<MaxNumBlockEntries; i++)
1568 MyCopyValuesPointers[i] =
new double[MaxRowDim*MaxColDim];
1571 int GlobalCopyRowDim, GlobalCopyNumBlockEntries;
1572 int * GlobalCopyBlockIndices =
new int[MaxNumBlockEntries];
1573 int * GlobalCopyColDims =
new int[MaxNumBlockEntries];
1574 int * GlobalCopyLDAs =
new int[MaxNumBlockEntries];
1576 int GlobalMaxRowDim =
A.GlobalMaxRowDim();
1577 int GlobalMaxColDim =
A.GlobalMaxColDim();
1578 int GlobalCopySizeOfValues = GlobalMaxRowDim*GlobalMaxColDim;
1579 double ** GlobalCopyValuesPointers =
new double*[MaxNumBlockEntries];
1580 for (
int i=0; i<MaxNumBlockEntries; i++)
1581 GlobalCopyValuesPointers[i] =
new double[GlobalMaxRowDim*GlobalMaxColDim];
1586 int MyView1RowDim, MyView1NumBlockEntries;
1587 int * MyView1BlockIndices;
1591 int MyView2RowDim, MyView2NumBlockEntries;
1592 int * MyView2BlockIndices;
1598 for (
int i=0; i<NumMyBlockRows; i++) {
1600 int GlobalRow =
A.GRID(i);
1602 EPETRA_TEST_ERR(
A.ExtractMyBlockRowPointers(MyRow, MaxNumBlockEntries, MyPointersRowDim,
1603 MyPointersNumBlockEntries, MyPointersBlockIndices,
1604 MyPointersValuesPointers), ierr );
1606 EPETRA_TEST_ERR(
A.ExtractGlobalBlockRowPointers(GlobalRow, MaxNumBlockEntries, GlobalPointersRowDim,
1607 GlobalPointersNumBlockEntries, GlobalPointersBlockIndices,
1608 GlobalPointersValuesPointers), ierr ) ;
1611 EPETRA_TEST_ERR(
A.BeginExtractMyBlockRowCopy(MyRow, MaxNumBlockEntries, MyCopyRowDim,
1612 MyCopyNumBlockEntries, MyCopyBlockIndices,
1613 MyCopyColDims), ierr);
1615 for (
int j=0; j<MyCopyNumBlockEntries; j++) {
1616 EPETRA_TEST_ERR(
A.ExtractEntryCopy(MyCopySizeOfValues, MyCopyValuesPointers[j], MaxRowDim,
false), ierr) ;
1617 MyCopyLDAs[j] = MaxRowDim;
1621 EPETRA_TEST_ERR(
A.BeginExtractGlobalBlockRowCopy(GlobalRow, MaxNumBlockEntries, GlobalCopyRowDim,
1622 GlobalCopyNumBlockEntries, GlobalCopyBlockIndices,
1623 GlobalCopyColDims), ierr ) ;
1625 for (
int j=0; j<GlobalCopyNumBlockEntries; j++) {
1626 EPETRA_TEST_ERR(
A.ExtractEntryCopy(GlobalCopySizeOfValues, GlobalCopyValuesPointers[j], GlobalMaxRowDim,
false), ierr );
1627 GlobalCopyLDAs[j] = GlobalMaxRowDim;
1632 MyView1NumBlockEntries, MyView1BlockIndices), ierr ) ;
1634 for (
int j=0; j<MyView1NumBlockEntries; j++)
1640 MyView2NumBlockEntries, MyView2BlockIndices,
1641 MyView2ValuesPointers), ierr );
1643 forierr += !(MyPointersNumBlockEntries==GlobalPointersNumBlockEntries);
1644 forierr += !(MyPointersNumBlockEntries==MyCopyNumBlockEntries);
1645 forierr += !(MyPointersNumBlockEntries==GlobalCopyNumBlockEntries);
1646 forierr += !(MyPointersNumBlockEntries==MyView1NumBlockEntries);
1647 forierr += !(MyPointersNumBlockEntries==MyView2NumBlockEntries);
1648 for (
int j=1; j<MyPointersNumBlockEntries; j++) {
1649 forierr += !(MyCopyBlockIndices[j-1]<MyCopyBlockIndices[j]);
1650 forierr += !(MyView1BlockIndices[j-1]<MyView1BlockIndices[j]);
1651 forierr += !(MyView2BlockIndices[j-1]<MyView2BlockIndices[j]);
1653 forierr += !(GlobalPointersBlockIndices[j]==
A.GCID(MyPointersBlockIndices[j]));
1654 forierr += !(
A.LCID(GlobalPointersBlockIndices[j])==MyPointersBlockIndices[j]);
1655 forierr += !(GlobalPointersBlockIndices[j]==GlobalCopyBlockIndices[j]);
1664 MyPointersRowDim, my->
N(),
1665 global->
A(), global->
LDA(),
1666 GlobalPointersRowDim, global->
N())==0);
1668 MyPointersRowDim, my->
N(),
1669 MyCopyValuesPointers[j], MyCopyLDAs[j],
1670 MyCopyRowDim, MyCopyColDims[j])==0);
1672 MyPointersRowDim, my->
N(),
1673 GlobalCopyValuesPointers[j], GlobalCopyLDAs[j],
1674 GlobalCopyRowDim, GlobalCopyColDims[j])==0);
1676 MyPointersRowDim, my->
N(),
1677 myview1->
A(), myview1->
LDA(),
1678 MyView1RowDim, myview1->
N())==0);
1680 MyPointersRowDim, my->
N(),
1681 myview2->
A(), myview2->
LDA(),
1682 MyView2RowDim, myview2->
N())==0);
1689 MyView1NumBlockEntries,
1690 MyView1BlockIndices)==-2),ierr);
1694 MyView2NumBlockEntries, MyView2BlockIndices,
1695 MyView2ValuesPointers)==-2),ierr);
1697 delete [] MyPointersBlockIndices;
1698 delete [] GlobalPointersBlockIndices;
1699 delete [] MyCopyBlockIndices;
1700 delete [] MyCopyColDims;
1701 delete [] MyCopyLDAs;
1702 for (
int i=0; i<MaxNumBlockEntries; i++)
delete [] MyCopyValuesPointers[i];
1703 delete [] MyCopyValuesPointers;
1704 delete [] GlobalCopyBlockIndices;
1705 delete [] GlobalCopyColDims;
1706 delete [] GlobalCopyLDAs;
1707 for (
int i=0; i<MaxNumBlockEntries; i++)
delete [] GlobalCopyValuesPointers[i];
1708 delete [] GlobalCopyValuesPointers;
1709 delete [] MyView1ValuesPointers;
1710 if (verbose) cout <<
"\n\nRows sorted check OK" << endl<< endl;
1717 double *
B,
int LDB,
int NumRowsB,
int NumColsB) {
1719 int ierr = 0, forierr = 0;
1728 for (
int j=0; j<NumColsA; j++) {
1731 for (
int i=0; i<NumRowsA; i++) forierr += (*ptr1++ != *ptr2++);
1739 int numProcs = comm.
NumProc();
1740 int localProc = comm.
MyPID();
1742 int myFirstRow = localProc*3;
1743 int myLastRow = myFirstRow+2;
1744 int numMyRows = myLastRow - myFirstRow + 1;
1745 int numGlobalRows = numProcs*numMyRows;
1752 int numCols = 2*numMyRows;
1753 int* myCols =
new int[numCols];
1755 int col = myFirstRow;
1756 for(
int i=0; i<numCols; ++i) {
1758 if (col > myLastRow) col = myFirstRow;
1765 elemSize, indexBase, comm);
1775 double* coef =
new double[elemSize*elemSize];
1776 for(
int i=0; i<elemSize*elemSize; ++i) {
1784 for(
int i=myFirstRow; i<=myLastRow; ++i) {
1787 for(
int j=0; j<numCols; ++j) {
1789 elemSize, elemSize), ierr);
1799 if (verbose) cout <<
"Multiply x"<<endl;
1810 int numBlockEntries = 0;
1812 int** BlockIndices =
new int*[numMyRows];
1816 for(
int i=myFirstRow; i<=myLastRow; ++i) {
1817 BlockIndices[i-myFirstRow] =
new int[numCols];
1819 RowDim, numBlockEntries,
1820 BlockIndices[i-myFirstRow],
1824 BlockIndices[i-myFirstRow]), ierr);
1826 if (numMyRows != numBlockEntries)
return(-1);
1827 if (RowDim != elemSize)
return(-2);
1828 for(
int j=0; j<numBlockEntries; ++j) {
1829 if (Values[j]->
A()[0] != 1.0) {
1830 cout <<
"Row " << i <<
" Values["<<j<<
"][0]: "<< Values[j][0]
1831 <<
" should be 1.0" << endl;
1838 Values[j]->
N()), ierr);
1849 for(
int i=myFirstRow; i<=myLastRow; ++i) {
1851 RowDim, numBlockEntries,
1852 BlockIndices[i-myFirstRow],
1855 if (numMyRows != numBlockEntries)
return(-1);
1856 if (RowDim != elemSize)
return(-2);
1857 for(
int j=0; j<numBlockEntries; ++j) {
1858 if (Values[j]->
A()[0] != 1.0) {
1859 cout <<
"Aview: Row " << i <<
" Values["<<j<<
"][0]: "<< Values[j][0]
1860 <<
" should be 1.0" << endl;
1865 delete [] BlockIndices[i-myFirstRow];
1869 if (verbose&&localProc==0) {
1870 cout <<
"checkMergeRedundantEntries, A:" << endl;
1874 delete [] BlockIndices;
1882 int numProcs = comm.
NumProc();
1883 int localProc = comm.
MyPID();
1885 int myFirstRow = localProc*3;
1886 int myLastRow = myFirstRow+2;
1887 int numMyRows = myLastRow - myFirstRow + 1;
1888 int numGlobalRows = numProcs*numMyRows;
1891 int numCols = numMyRows;
1892 int* myCols =
new int[numCols];
1894 int col = myFirstRow;
1895 for(
int i=0; i<numCols; ++i) {
1897 if (col > myLastRow) col = myFirstRow;
1904 elemSize, indexBase, comm);
1908 double* coef =
new double[elemSize*elemSize];
1910 for(
int i=myFirstRow; i<=myLastRow; ++i) {
1911 int myPointRow = i*elemSize;
1915 for(
int ii=0; ii<elemSize; ++ii) {
1916 for(
int jj=0; jj<elemSize; ++jj) {
1917 double val = (myPointRow+ii)*1.0;
1918 coef[ii+elemSize*jj] = val;
1924 for(
int j=0; j<numCols; ++j) {
1926 elemSize, elemSize), ierr);
1938 int len = elemSize*numCols, checkLen;
1939 double* values =
new double[len];
1940 int* indices =
new int[len];
1941 int RowDim, numBlockEntries;
1943 for(
int i=myFirstRow; i<=myLastRow; ++i) {
1945 RowDim, numBlockEntries,
1947 blockEntries), ierr);
1948 if (numMyRows != numBlockEntries)
return(-1);
1949 if (RowDim != elemSize)
return(-2);
1951 int myPointRow = i*elemSize - myFirstRow*elemSize;
1952 for(
int ii=0; ii<elemSize; ++ii) {
1954 checkLen, values, indices), ierr);
1955 if (len != checkLen)
return(-3);
1957 double val = (i*elemSize+ii)*1.0;
1958 double blockvalue = blockEntries[0]->
A()[ii];
1960 for(
int jj=0; jj<len; ++jj) {
1961 if (values[jj] != val)
return(-4);
1962 if (values[jj] != blockvalue)
return(-5);
1975 int numProcs = comm.
NumProc();
1976 int localProc = comm.
MyPID();
1978 int myFirstRow = localProc*3;
1979 int myLastRow = myFirstRow+2;
1980 int numMyRows = myLastRow - myFirstRow + 1;
1981 int numGlobalRows = numProcs*numMyRows;
1985 int num_off_diagonals = 1;
1988 num_off_diagonals, elemSize);
1997 for(
int i=myFirstRow; i<=myLastRow; ++i) {
2000 colindices[i]), ierr);
2002 for(
int j=0; j<rowlengths[i]; ++j) {
2004 elemSize, elemSize), ierr);
2016 A.Multiply(
false, x, y);
2017 A.Multiply(
false, x, x);
2019 double* xptr = x.Values();
2020 double* yptr = y.Values();
2022 for(
int i=0; i<numMyRows; ++i) {
2023 if (xptr[i] != yptr[i]) {
2033 int localProc = comm.
MyPID();
2034 int numProcs = comm.
NumProc();
2035 int myFirstRow = localProc*3;
2036 int myLastRow = myFirstRow+2;
2037 int numMyRows = myLastRow - myFirstRow + 1;
2038 int numGlobalRows = numProcs*numMyRows;
2042 int num_off_diagonals = 1;
2045 num_off_diagonals, elemSize);
2054 for(
int i=myFirstRow; i<=myLastRow; ++i) {
2057 colindices[i]), ierr);
2059 for(
int j=0; j<rowlengths[i]; ++j) {
2061 elemSize, elemSize, elemSize), ierr);
2069 int errcode =
A->BeginReplaceMyValues(0, 0, 0);
2078 for(
int i=myFirstRow; i<=myLastRow; ++i) {
2081 colindices[i]), ierr);
2083 for(
int j=0; j<rowlengths[i]; ++j) {
2085 elemSize, elemSize, elemSize), ierr);
2102 const int node = comm.
MyPID();
2103 const int nodes = comm.
NumProc();
2120 else if (nodes == 2)
2128 for (
int j = 0; j < 4; ++j)
2130 for (
int i = 0; i < 2; ++i)
2132 l2g[i + 2 * j] = (i + i_off) + (j + j_off) * Gi;
2137 else if (nodes == 4)
2155 for (
int j = 0; j < 2; ++j)
2157 for (
int i = 0; i < 2; ++i)
2159 l2g[i + 2 * j] = (i + i_off) + (j + j_off) * Gi;
2173 for (
int j = 0; j < Nj; ++j)
2175 for (
int i = 0; i < Ni; ++i)
2180 indices[ctr++] = gi + gj * Gi;
2182 indices[ctr++] = (gi - 1) + gj * Gi;
2184 indices[ctr++] = (gi + 1) + gj * Gi;
2186 indices[ctr++] = gi + (gj - 1) * Gi;
2188 indices[ctr++] = gi + (gj + 1) * Gi;
2217 std::fill(
C.A(),
C.A()+9, -4.0);
2218 std::fill(L.
A(), L.
A()+9, 2.0);
2219 std::fill(R.
A(), R.
A()+9, 2.0);
2223 for (
int j = 0; j < Nj; ++j)
2225 for (
int i = 0; i < Ni; ++i)
2231 int local = i + j * Ni;
2232 int global = gi + gj * Gi;
2234 int left = (gi - 1) + gj * Gi;
2235 int right = (gi + 1) + gj * Gi;
2236 int bottom = gi + (gj - 1) * Gi;
2237 int top = gi + (gj + 1) * Gi;
2241 indices[ctr++] = local;
2243 indices[ctr++] = left;
2245 indices[ctr++] = right;;
2247 indices[ctr++] = bottom;
2249 indices[ctr++] = top;
2251 matrix->BeginReplaceMyValues(local, ctr, &indices[0]);
2252 matrix->SubmitBlockEntry(
C);
2253 if (gi > first) matrix->SubmitBlockEntry(L);
2254 if (gi < last) matrix->SubmitBlockEntry(R);
2255 if (gj > first) matrix->SubmitBlockEntry(L);
2256 if (gj < last) matrix->SubmitBlockEntry(R);
2257 matrix->EndSubmitEntries();
2261 matrix->FillComplete();
2278 if (verbose) cout <<
"Checking VbrRowMatrix Adapter..." << endl;
2294 for (
int i=0; i<
A.NumMyRows(); i++) {
2296 A.NumMyRowEntries(i,nA);
B.NumMyRowEntries(i,nB);
2316 bool transA =
false;
2317 A.SetUseTranspose(transA);
2318 B.SetUseTranspose(transA);
2320 A.Multiply(transA, X, YA2);
2324 B.Multiply(transA, X, YB2);
2337 A.SetUseTranspose(transA);
2338 B.SetUseTranspose(transA);
2340 A.Multiply(transA, X, YA2);
2344 B.Multiply(transA, X,YB2);
2382 vector<double> valuesA(
A.MaxNumEntries());
2383 vector<int> indicesA(
A.MaxNumEntries());
2384 vector<double> valuesB(
B.MaxNumEntries());
2385 vector<int> indicesB(
B.MaxNumEntries());
2387 for (
int i=0; i<
A.NumMyRows(); i++) {
2389 EPETRA_TEST_ERR(
A.ExtractMyRowCopy(i,
A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr);
2390 EPETRA_TEST_ERR(
B.ExtractMyRowCopy(i,
B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
2392 for (
int j=0; j<nA; j++) {
2393 double curVal = valuesA[j];
2394 int curIndex = indicesA[j];
2395 bool notfound =
true;
2397 while (notfound && jj< nB) {
2398 if (!
checkValues(curVal, valuesB[jj])) notfound =
false;
2402 vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex);
2410 cout <<
"RowMatrix Methods check OK" << endl;
2412 cout <<
"ierr = " << ierr <<
". RowMatrix Methods check failed." << endl;
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int Random()
Set multi-vector values to random numbers.
int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
int N() const
Returns column dimension of system.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int M() const
Returns row dimension of system.
int Random()
Set matrix values to random numbers.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
#define EPETRA_TEST_ERR(a, b)
const Epetra_Map & RowMatrixColMap() const
Returns the Epetra_Map object associated with columns of this matrix.
int CompareValues(double *A, int LDA, int NumRowsA, int NumColsA, double *B, int LDB, int NumRowsB, int NumColsB)
double NormOne() const
Returns the one norm of the global matrix.
int checkMatvecSameVectors(Epetra_Comm &comm, bool verbose)
int checkEarlyDelete(Epetra_Comm &comm, bool verbose)
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int checkVbrMatrixOptimizedGraph(Epetra_Comm &comm, bool verbose)
Epetra_VbrMatrix: A class for the construction and use of real-valued double-precision variable block...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])
Epetra_MpiComm: The Epetra MPI Communication Class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
std::string Epetra_Version()
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int MyPID() const =0
Return my process ID.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
Epetra_Time: The Epetra Timing Class.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
int NumMyRows() const
Returns the number of matrix rows on this processor.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
double Flops() const
Returns the number of floating point operations with this object and resets the count.
void ConvertVbrToCrs(Epetra_VbrMatrix *VbrIn, Epetra_CrsMatrix *&CrsOut)
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int NumMyBlockRows() const
Returns the number of block matrix rows on this processor.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
Epetra_VbrRowMatrix: A class for using an existing Epetra_VbrMatrix object as an Epetra_RowMatrix obj...
int checkMergeRedundantEntries(Epetra_Comm &comm, bool verbose)
Epetra_Comm: The Epetra Communication Abstract Base Class.
double * A() const
Returns pointer to the this matrix.
int NumVectors() const
Returns the number of vectors in the multi-vector.
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
int NumMyElements() const
Number of elements on the calling processor.
int TestMatrix(Epetra_Comm &Comm, bool verbose, bool debug, int NumMyElements, int MinSize, int MaxSize, ConsType ConstructorType, bool ExtraBlocks, bool insertlocal, bool symmetric, Epetra_VbrMatrix **PreviousA)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false...
int ExtractGlobalBlockRowPointers(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Copy the block indices into user-provided array, set pointers for rest of data for specified global b...
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the left with a Epetra_Vector x.
int checkmultiply(bool transpose, Epetra_VbrMatrix &A, Epetra_MultiVector &X, Epetra_MultiVector &Check_Y)
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
Epetra_Flops: The Epetra Floating Point Operations Class.
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
matrix_data is a very simple data source to be used for filling test matrices.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int checkValues(double x, double y, string message="", bool verbose=false)
virtual int NumProc() const =0
Returns total number of processes.
int MyPID() const
Return my process ID.
int checkMultiVectors(Epetra_MultiVector &X, Epetra_MultiVector &Y, string message="", bool verbose=false)
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int power_method(bool TransA, Epetra_VbrMatrix &A, Epetra_MultiVector &q, Epetra_MultiVector &z, Epetra_MultiVector &resid, double *lambda, int niters, double tolerance, bool verbose)
int NumGlobalElements() const
Number of elements across all processors.
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int FillComplete()
Signal that data entry is complete, perform transformations to local index space. ...
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given global row of the matrix, values are inserted via...
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
Submit a block entry to the indicated block row and column specified in the Begin routine...
int check(Epetra_VbrMatrix &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1, int *MyGlobalElements, bool verbose)
int LDA() const
Returns the leading dimension of the this matrix.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int InsertMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given local row of the matrix.
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the right with a Epetra_Vector x.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
const Epetra_Map & RowMatrixRowMap() const
Returns the EpetraMap object associated with the rows of this matrix.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
int checkVbrRowMatrix(Epetra_RowMatrix &A, Epetra_RowMatrix &B, bool verbose)
double ElapsedTime(void) const
Epetra_Time elapsed time function.
int checkExtractMyRowCopy(Epetra_Comm &comm, bool verbose)