71int main(
int argc,
char *argv[])
73 int ierr = 0, forierr = 0;
80 MPI_Init(&argc,&argv);
83 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
96 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
98 int verbose_int = verbose ? 1 : 0;
100 verbose = verbose_int==1 ? true :
false;
109 int MyPID = Comm.
MyPID();
112 if(verbose && MyPID==0)
115 if (verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
116 <<
" is alive."<<endl;
118 bool verbose1 = verbose;
121 if(verbose && rank!=0)
124 int NumMyEquations = 10000;
125 int NumGlobalEquations = (NumMyEquations * NumProc) +
EPETRA_MIN(NumProc,3);
131 Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
140 int* NumNz =
new int[NumMyEquations];
145 for (
int i = 0; i < NumMyEquations; i++)
146 if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
162 double* Values =
new double[2];
165 int* Indices =
new int[2];
170 for (
int i = 0; i < NumMyEquations; i++) {
171 if(MyGlobalElements[i] == 0) {
175 else if (MyGlobalElements[i] == NumGlobalEquations-1) {
176 Indices[0] = NumGlobalEquations-2;
180 Indices[0] = MyGlobalElements[i]-1;
181 Indices[1] = MyGlobalElements[i]+1;
184 forierr += !(A.
InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
185 forierr += !(A.
InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i)>0);
189 int * indexOffsetTmp;
207 int NumMyNonzeros = 3 * NumMyEquations;
210 if(A.
LRID(NumGlobalEquations-1) >= 0)
212 EPETRA_TEST_ERR(
check(A, NumMyEquations, NumGlobalEquations, NumMyNonzeros, 3*NumGlobalEquations-2,
213 MyGlobalElements, verbose),ierr);
215 for (
int i = 0; i < NumMyEquations; i++)
219 for (
int i = 0; i < NumMyEquations; i++)
223 if (verbose) cout <<
"\n\nNumEntries function check OK" << std::endl<< std::endl;
237 double tolerance = 1.0e-1;
254 double MFLOPs = total_flops/elapsed_time/1000000.0;
256 if (verbose) cout <<
"\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
262 if (verbose) cout <<
"\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
271 MFLOPs = total_flops/elapsed_time/1000000.0;
273 if (verbose) cout <<
"\n\nTotal MFLOPs for transpose solve = " << MFLOPs << std::endl<< endl;
279 if (verbose) cout <<
"\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
285 double * Rowvals =
new double [numvals];
286 int * Rowinds =
new int [numvals];
289 for (
int i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
302 MFLOPs = total_flops/elapsed_time/1000000.0;
304 if (verbose) cout <<
"\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
310 if (verbose) cout <<
"\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
320 MFLOPs = total_flops/elapsed_time/1000000.0;
323 if (verbose) cout <<
"\n\nTotal MFLOPs for tranpose of second solve = " << MFLOPs << endl<< endl;
325 if (verbose) cout <<
"\n\n*****Testing constant entry constructor" << endl<< endl;
331 double dble_one = 1.0;
332 for (
int i=0; i< NumMyEquations; i++) AA.
InsertGlobalValues(MyGlobalElements[i], 1, &dble_one, MyGlobalElements+i);
348 EPETRA_TEST_ERR(
check(AA, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
349 MyGlobalElements, verbose),ierr);
354 for (
int i=0; i<NumMyEquations; i++) forierr += !(AA.
NumGlobalEntries(MyGlobalElements[i])==1);
357 if (verbose) cout <<
"\n\nNumEntries function check OK" << endl<< endl;
361 if (verbose) cout <<
"\n\n*****Testing copy constructor" << endl<< endl;
364 EPETRA_TEST_ERR(
check(B, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
365 MyGlobalElements, verbose),ierr);
368 for (
int i=0; i<NumMyEquations; i++) forierr += !(B.
NumGlobalEntries(MyGlobalElements[i])==1);
371 if (verbose) cout <<
"\n\nNumEntries function check OK" << endl<< endl;
375 if (verbose) cout <<
"\n\n*****Testing local view constructor" << endl<< endl;
382 for (
int i = 0; i < NumMyEquations; i++) {
387 EPETRA_TEST_ERR(
check(BV, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
388 MyGlobalElements, verbose),ierr);
391 for (
int i=0; i<NumMyEquations; i++) forierr += !(BV.
NumGlobalEntries(MyGlobalElements[i])==1);
394 if (verbose) cout <<
"\n\nNumEntries function check OK" << endl<< endl;
397 if (verbose) cout <<
"\n\n*****Testing post construction modifications" << endl<< endl;
406 delete [] MyGlobalElements;
413 int NumMyElements1 = 2;
414 int NumMyEquations1 = NumMyElements1;
415 int NumGlobalEquations1 = NumMyEquations1*NumProc;
426 int * NumNz1 =
new int[NumMyEquations1];
431 for (
int i=0; i<NumMyEquations1; i++)
432 if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
446 double *Values1 =
new double[2];
447 Values1[0] = -1.0; Values1[1] = -1.0;
448 int *Indices1 =
new int[2];
453 for (
int i=0; i<NumMyEquations1; i++)
455 if (MyGlobalElements1[i]==0)
460 else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
462 Indices1[0] = NumGlobalEquations1-2;
467 Indices1[0] = MyGlobalElements1[i]-1;
468 Indices1[1] = MyGlobalElements1[i]+1;
471 forierr += !(A1.
InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1)==0);
472 forierr += !(A1.
InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i)>0);
487 for (
int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==two1);
493 for (
int i=0; i<NumMyEquations1; i++) checkDiag[i]=two1*two1;
502 for (
int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
505 if (verbose) cout <<
"\n\nDiagonal extraction and replacement OK.\n\n" << endl;
507 double orignorm = A1.
NormOne();
511 if (verbose) cout <<
"\n\nMatrix scale OK.\n\n" << endl;
513 if (verbose) cout <<
"\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
519 delete [] MyGlobalElements1;
523 if (verbose) cout <<
"\n\n*****Testing LeftScale and RightScale" << endl << endl;
525 int NumMyElements2 = 7;
529 Epetra_Map ColMap(NumMyElements2,NumMyElements2,0,Comm);
533 int NumMyRangeElements2 = 0;
536 NumMyRangeElements2 = NumMyRows2*2;
537 if (NumProc % 2 == 1 && MyPID == NumProc-1)
538 NumMyRangeElements2 = NumMyRows2;
540 Epetra_Map RangeMap(-1,NumMyRangeElements2,0,Comm);
542 double * Values2 =
new double[NumMyElements2];
543 int * Indices2 =
new int[NumMyElements2];
545 for (
int i=0; i<NumMyElements2; i++) {
546 Values2[i] = i+MyPID;
554 double * RowLeftScaleValues =
new double[NumMyRows2];
555 double * ColRightScaleValues =
new double[NumMyElements2];
557 for (
int i=0; i<RowLoopLength; i++)
558 RowLeftScaleValues[i] = (i + RowMap.
MinMyGID() ) % 2 + 1;
560 for (
int i=0; i<NumMyElements2;i++)
561 ColRightScaleValues[i] = i % 2 + 1;
564 double * RangeLeftScaleValues =
new double[RangeLoopLength];
566 double * DomainRightScaleValues =
new double[DomainLoopLength];
567 for (
int i=0; i<RangeLoopLength; i++)
568 RangeLeftScaleValues[i] = 1.0/((i + RangeMap.
MinMyGID() ) % 2 + 1);
569 for (
int i=0; i<DomainLoopLength;i++)
570 DomainRightScaleValues[i] = 1.0/((i + DomainMap.
MinMyGID() ) % 2 + 1);
577 double A2infNorm = A2.
NormInf();
578 double A2oneNorm = A2.
NormOne();
580 if (verbose1) cout << A2;
582 double A2infNorm1 = A2.
NormInf();
583 double A2oneNorm1 = A2.
NormOne();
584 bool ScalingBroke =
false;
585 if (A2infNorm1>2*A2infNorm||A2infNorm1<A2infNorm) {
589 if (A2oneNorm1>2*A2oneNorm||A2oneNorm1<A2oneNorm) {
594 if (verbose1) cout << A2;
596 double A2infNorm2 = A2.
NormInf();
597 double A2oneNorm2 = A2.
NormOne();
598 if (A2infNorm2>=2*A2infNorm1||A2infNorm2<=A2infNorm1) {
602 if (A2oneNorm2>2*A2oneNorm1||A2oneNorm2<=A2oneNorm1) {
606 if (verbose1) cout << A2;
608 double A2infNorm3 = A2.
NormInf();
609 double A2oneNorm3 = A2.
NormOne();
611 if (A2infNorm3!=A2infNorm1) {
615 if (A2oneNorm3!=A2oneNorm1) {
619 if (verbose1) cout << A2;
621 double A2infNorm4 = A2.
NormInf();
622 double A2oneNorm4 = A2.
NormOne();
624 if (A2infNorm4!=A2infNorm) {
628 if (A2oneNorm4!=A2oneNorm) {
643 for (
int i=0 ; i< num_my_rows; i++ ) {
645 for (
int j = 0 ; j <num_entries; j++ ) {
654 double A2infNorm5 = A2.
NormInf();
655 double A2oneNorm5 = A2.
NormOne();
657 if (A2infNorm5!=2.0 * A2infNorm4) {
661 if (A2oneNorm5!= 2.0 * A2oneNorm4) {
669 for (
int i=0 ; i< num_my_rows; i++ ) {
671 for (
int j = 0 ; j <num_entries; j++ ) {
676 if (verbose1) cout << A2;
679 if (verbose) cout << endl <<
"LeftScale and RightScale tests FAILED" << endl << endl;
682 if (verbose) cout << endl <<
"LeftScale and RightScale tests PASSED" << endl << endl;
687 if (verbose) cout <<
"\n\n*****Testing InvRowMaxs and InvColMaxs" << endl << endl;
689 if (verbose1) cout << A2 << endl;
692 if (verbose1) cout << xRow << endl << xRange << endl;
694 if (verbose) cout <<
"\n\n*****Testing InvRowSums and InvColSums" << endl << endl;
695 bool InvSumsBroke =
false;
698 if (verbose1) cout << xRow;
700 float A2infNormFloat = A2.
NormInf();
701 if (verbose1) cout << A2 << endl;
702 if (fabs(1.0-A2infNormFloat) > 1.e-5) {
708 int expectedcode = 1;
709 if (Comm.
NumProc()>1) expectedcode = 0;
711 if (verbose1) cout << xDomain << endl;
713 float A2oneNormFloat2 = A2.
NormOne();
714 if (verbose1) cout << A2;
715 if (fabs(1.0-A2oneNormFloat2)>1.e-5) {
723 if (verbose1) cout << xRange;
725 float A2infNormFloat2 = A2.
NormInf();
727 if (verbose1) cout << A2;
728 if (fabs(1.0-A2infNormFloat2)>1.e-5) {
729 cout <<
"InfNorm should be = 1, but InfNorm = " << A2infNormFloat2 << endl;
745 delete [] ColRightScaleValues;
746 delete [] DomainRightScaleValues;
747 if (verbose) cout <<
"Begin partial sum testing." << endl;
751 int * myGlobalElements =
new int[NumMyRows3];
752 for (
int i=0; i<NumMyRows3; i++) myGlobalElements[i] = MyPID+i;
753 Epetra_Map RowMap3(NumProc*2, NumMyRows3, myGlobalElements, 0, Comm);
754 int NumMyElements3 = 5;
756 double * Values3 =
new double[NumMyElements3];
757 int * Indices3 =
new int[NumMyElements3];
758 for (
int i=0; i < NumMyElements3; i++) {
759 Values3[i] = (int) (MyPID + (i+1));
762 for (
int i=0; i<NumMyRows3; i++) {
766 Epetra_Map DomainMap3(NumMyElements3, 0, Comm);
768 if (verbose1) cout << A3;
774 if (verbose1) cout << xRange3;
776 float A3infNormFloat = A3.
NormInf();
777 if (verbose1) cout << A3;
778 if (1.0!=A3infNormFloat) {
779 cout <<
"InfNorm should be = 1, but InfNorm = " << A3infNormFloat <<endl;
784 int NumMyColumns3 = NumMyRows3;
789 double *Values3cm =
new double[NumMyColumns3];
790 int * Indices3cm =
new int[NumMyColumns3];
791 for (
int i=0; i<NumMyColumns3; i++) {
792 Values3cm[i] = MyPID + i + 1;
793 Indices3cm[i]= i + MyPID;
795 for (
int ii=0; ii<NumMyElements3; ii++) {
806 if (verbose1) cout << A3cm << endl;
814 if (verbose1) cout << xDomain3cm << endl;
817 float A3cmOneNormFloat = A3cm.
NormOne();
818 if (verbose1) cout << A3cm << endl;
819 if (1.0!=A3cmOneNormFloat) {
820 cout <<
"OneNorm should be = 1, but OneNorm = " << A3cmOneNormFloat << endl;
825 if (verbose) cout <<
"End partial sum testing" << endl;
826 if (verbose) cout <<
"Begin replicated testing" << endl;
833 for (
int ii=0; ii < NumMyElements3; ii++) {
834 Values3[ii] = (int)((ii*.6)+1.0);
836 for (
int ii=0; ii<NumMyRows3; ii++) {
840 if (verbose1) cout << A4 << endl;
844 if (verbose1) cout << xRow3 << xRange3;
847 if (verbose1) cout << xRow3;
849 float A4infNormFloat = A4.
NormInf();
850 if (verbose1) cout << A4;
851 if (2.0!=A4infNormFloat && NumProc != 1) {
852 if (verbose1) cout <<
"InfNorm should be = 2 (because one column is replicated on two processors and NormOne() does not handle replication), but InfNorm = " << A4infNormFloat <<endl;
856 else if (1.0!=A4infNormFloat && NumProc == 1) {
857 if (verbose1) cout <<
"InfNorm should be = 1, but InfNorm = " << A4infNormFloat <<endl;
865 for (
int ii=0; ii<NumMyElements3; ii++) {
869 if (verbose1) cout << A4cm << endl;
873 if (verbose1) cout << xCol3cm << xDomain3cm;
877 if (verbose1) cout << xCol3cm << endl;
879 float A4cmOneNormFloat = A4cm.
NormOne();
880 if (verbose1) cout << A4cm << endl;
881 if (2.0!=A4cmOneNormFloat && NumProc != 1) {
882 if (verbose1) cout <<
"OneNorm should be = 2 (because one column is replicated on two processors and NormOne() does not handle replication), but OneNorm = " << A4cmOneNormFloat << endl;
886 else if (1.0!=A4cmOneNormFloat && NumProc == 1) {
887 if (verbose1) cout <<
"OneNorm should be = 1, but OneNorm = " << A4infNormFloat <<endl;
892 if (verbose) cout <<
"End replicated testing" << endl;
895 if (verbose) cout << endl <<
"InvRowSums tests FAILED" << endl << endl;
898 if (verbose) cout << endl <<
"InvRowSums tests PASSED" << endl << endl;
902 double check_frobnorm = sqrt(nnz_A3cm*4.0);
905 bool frobnorm_test_failed =
false;
906 if (fabs(check_frobnorm-frobnorm) > 5.e-5) {
907 frobnorm_test_failed =
true;
910 if (frobnorm_test_failed) {
911 if (verbose) std::cout <<
"Frobenius-norm test FAILED."<<std::endl;
918 int NumMyRows = (MyPID==0) ? NumGlobalEquations : 0;
923 for(
int i=0; i<NumMyRows; i++) {
924 int GID = Map1.
GID(i);
932 if (verbose) std::cout <<
"Subcommunicator test FAILED."<<std::endl;
942 delete [] myGlobalElements;
946 delete [] Indices3cm;
947 delete [] RangeLeftScaleValues;
948 delete [] RowLeftScaleValues;