test/VbrMatrix/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 
00030 //
00031 //
00032 //      valgrind usage:
00033 //         valgrind --suppressions=Suppressions --leak-check=yes --show-reachable=yes ./VbrMatrix_test.exe
00034 //
00035 //      mpirun -np 2 valgrind --suppressions=Suppressions --logfile=valg.out --leak-check=yes --show-reachable=yes ./VbrMatrix_test.exe
00036 //      The file Suppressions can be found in packages/epetra/test/VbrMatrix/Suppressions.in
00037 //
00038 //
00039 #include "Epetra_SerialDenseMatrix.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_Time.h"
00042 #include "Epetra_Vector.h"
00043 #include "Epetra_Flops.h"
00044 #include "Epetra_VbrMatrix.h"
00045 #include "Epetra_VbrRowMatrix.h"
00046 #include "Epetra_CrsMatrix.h"
00047 #include <vector>
00048 #ifdef EPETRA_MPI
00049 #include "Epetra_MpiComm.h"
00050 #include "mpi.h"
00051 #else
00052 #include "Epetra_SerialComm.h"
00053 #endif
00054 #include "../epetra_test_err.h"
00055 #include "../src/Epetra_matrix_data.h"
00056 #include "Epetra_Version.h"
00057 
00058 // prototypes
00059 
00060 int checkValues( double x, double y, string message = "", bool verbose = false) {
00061   if (fabs((x-y)/x) > 0.01) {
00062     return(1);
00063     if (verbose) cout << "********** " << message << " check failed.********** " << endl;
00064   }
00065   else {
00066     if (verbose) cout << message << " check OK." << endl;
00067     return(0);
00068   }
00069 }
00070 
00071 int checkMultiVectors( Epetra_MultiVector & X, Epetra_MultiVector & Y, string message = "", bool verbose = false) {
00072   int numVectors = X.NumVectors();
00073   int length = Y.MyLength();
00074   int badvalue = 0;
00075   int globalbadvalue = 0;
00076   for (int j=0; j<numVectors; j++)
00077     for (int i=0; i< length; i++)
00078       if (checkValues(X[j][i], Y[j][i])==1) badvalue = 1;
00079   X.Map().Comm().MaxAll(&badvalue, &globalbadvalue, 1);
00080 
00081   if (verbose) {
00082     if (globalbadvalue==0) cout << message << " check OK." << endl;
00083     else cout << "********* " << message << " check failed.********** " << endl;
00084   }
00085   return(globalbadvalue);
00086 }
00087 
00088 int checkVbrRowMatrix(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose);
00089 
00090 int CompareValues(double * A, int LDA, int NumRowsA, int NumColsA, 
00091       double * B, int LDB, int NumRowsB, int NumColsB);
00092 
00093 int check(Epetra_VbrMatrix& A, 
00094     int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, 
00095     int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1, 
00096     int * MyGlobalElements, bool verbose);
00097 
00098 int power_method(bool TransA, Epetra_VbrMatrix& A, 
00099      Epetra_MultiVector& q,
00100      Epetra_MultiVector& z, 
00101      Epetra_MultiVector& resid, 
00102      double * lambda, int niters, double tolerance,
00103      bool verbose);
00104 
00105 int checkMergeRedundantEntries(Epetra_Comm& comm, bool verbose);
00106 
00107 int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose);
00108 
00109 int checkMatvecSameVectors(Epetra_Comm& comm, bool verbose);
00110 
00111 int checkEarlyDelete(Epetra_Comm& comm, bool verbose);
00112 
00113 //
00114 //  ConvertVbrToCrs is a crude but effective way to convert a Vbr matrix to a Crs matrix
00115 //  Caller is responsible for deleting the CrsMatrix CrsOut
00116 //
00117 void ConvertVbrToCrs( Epetra_VbrMatrix* VbrIn, Epetra_CrsMatrix*& CrsOut ) {
00118 
00119   const Epetra_Map &E_Vbr_RowMap = VbrIn->RowMatrixRowMap() ; 
00120   const Epetra_Map &E_Vbr_ColMap = VbrIn->RowMatrixColMap() ; 
00121   //    const Epetra_Map &E_Vbr_RowMap = VbrIn->OperatorRangeMap() ; 
00122   //    const Epetra_Map &E_Vbr_ColMap = VbrIn->OperatorDomainMap() ; 
00123     
00124     CrsOut = new Epetra_CrsMatrix( Copy, E_Vbr_RowMap, E_Vbr_ColMap, 0 ); 
00125     //  CrsOut = new Epetra_CrsMatrix( Copy, E_Vbr_RowMap, 0 ); 
00126     int NumMyElements = VbrIn->RowMatrixRowMap().NumMyElements() ;
00127     vector<int> MyGlobalElements( NumMyElements );
00128     VbrIn->RowMatrixRowMap().MyGlobalElements( &MyGlobalElements[0] ) ;
00129 
00130     int NumMyColumns = VbrIn->RowMatrixColMap().NumMyElements() ;
00131     vector<int> MyGlobalColumns( NumMyColumns );
00132     VbrIn->RowMatrixColMap().MyGlobalElements( &MyGlobalColumns[0] ) ;
00133 
00134     int MaxNumIndices = VbrIn->MaxNumEntries();
00135     int NumIndices;
00136     vector<int> LocalColumnIndices(MaxNumIndices);
00137     vector<int> GlobalColumnIndices(MaxNumIndices);
00138     vector<double> MatrixValues(MaxNumIndices); 
00139 
00140     for( int LocalRow=0; LocalRow<NumMyElements; ++LocalRow ) {
00141 
00142       VbrIn->ExtractMyRowCopy( LocalRow, 
00143           MaxNumIndices,
00144           NumIndices, 
00145           &MatrixValues[0],
00146           &LocalColumnIndices[0] );
00147       
00148       for (int j = 0 ; j < NumIndices ; j++ ) 
00149   { 
00150     GlobalColumnIndices[j] = MyGlobalColumns[ LocalColumnIndices[j] ]  ;
00151   }
00152       
00153 #if 0
00154       if ( CrsOut->InsertGlobalValues( MyGlobalElements[LocalRow], 
00155               NumIndices, 
00156               &MatrixValues[0],
00157               &GlobalColumnIndices[0] )!=0)abort();
00158 #else
00159       if ( CrsOut->InsertMyValues( LocalRow, 
00160               NumIndices, 
00161               &MatrixValues[0],
00162               &LocalColumnIndices[0] )!= 0) abort();
00163 #endif
00164       
00165       
00166     }
00167     CrsOut->FillComplete();
00168 }
00169 
00170 //
00171 //  checkmultiply checks to make sure that AX=Y using both multiply 
00172 //  and both Crs Matrices, multiply1
00173 //
00174 
00175 int checkmultiply( bool transpose, Epetra_VbrMatrix& A, Epetra_MultiVector& X, Epetra_MultiVector& Check_Y ) {
00176 
00177   int numerrors = 0 ; 
00178 
00179 #if 1
00180   //
00181   //  If X and Y are Epetra_Vectors, we first test Multiply1
00182   //
00183   Epetra_Vector *vecY =  dynamic_cast<Epetra_Vector *>( &Check_Y );
00184   Epetra_Vector *vecX =  dynamic_cast<Epetra_Vector *>( &X );
00185   assert( vecX && vecY || (!vecX && !vecY) ) ;
00186 
00187   if ( vecX && vecY ) {
00188     double normY, NormError;
00189     Check_Y.NormInf( &normY ) ; 
00190     Epetra_Vector Y = *vecX ; 
00191     Y.PutScalar( -13.0 ) ; 
00192     Epetra_Vector Error = *vecX ; 
00193     A.Multiply1( transpose, *vecX, Y ) ;  
00194 
00195     Error.Update( 1.0, Y, -1.0, *vecY, 0.0 ) ; 
00196       
00197     Error.NormInf( &NormError ) ; 
00198     if ( NormError / normY > 1e-13 ) {
00199        numerrors++; 
00200        //cout << "Y = " << Y << endl;
00201        //cout << "vecY " << *vecY << endl;
00202        //cout << "Error " << Error << endl;
00203        //abort();
00204     }
00205     //
00206     //  Check x = Ax
00207     //
00208     Epetra_Vector Z = *vecX ; 
00209 
00210     A.Multiply1( transpose, Z, Z ) ;  
00211     Error.Update( 1.0, Z, -1.0, *vecY, 0.0 ) ; 
00212     //    Error.Update( 1.0, Y, -1.0, *vecY, 0.0 ) ; 
00213       
00214     Error.NormInf( &NormError ) ; 
00215     if ( NormError / normY > 1e-13 ) numerrors++; 
00216   }
00217   //
00218   //  Here we test Multiply 
00219   //
00220   Epetra_MultiVector Y = X ; 
00221   Epetra_MultiVector Error = X ; 
00222   A.Multiply( transpose, X, Y ) ; 
00223 
00224   int NumVecs = X.NumVectors() ; 
00225   vector<double> NormError(NumVecs);
00226   vector<double> Normy(NumVecs);
00227 
00228   Check_Y.NormInf( &Normy[0] ) ;
00229 
00230   Error.Update( 1.0, Y, -1.0, Check_Y, 0.0 ) ; 
00231   Error.NormInf( &NormError[0] ) ; 
00232 
00233   bool LoopError = false ; 
00234   for ( int ii = 0 ; ii < NumVecs ; ii++ ) {
00235     if( NormError[ii] / Normy[ii] > 1e-13 ) {
00236       LoopError = true ;
00237     }
00238   }
00239   if ( LoopError ) {
00240     numerrors++ ; 
00241   }
00242 
00243   //
00244   //  Check X = AX
00245   //
00246 
00247 #endif
00248   
00249   return numerrors;
00250 }
00251 //
00252 //  TestMatrix contains the bulk of the testing.  
00253 //    MinSize and MaxSize control the range of the block sizes - which are chosen randomly
00254 //    ConstructWithNumNz controls whether a Column Map is passed to the VbrMatrix contructor
00255 //      if false, the underlying graph will not be optimized
00256 //    ExtraBlocks, if true, causes extra blocks to be added to the matrix, further 
00257 //      guaranteeing that the matrix and graph will not have optimal storage.  
00258 //      ExtraBlocks is only supported for fixed block sizes (i.e. when minsize=maxsize)
00259 //
00260 //  If TestMatrix() is called with any of the ConsTypes that use PreviousA, the values used to 
00261 //  build A, i.e. MinSize, MaxSize must be the same as they were on the previous call to 
00262 //  TestMatrix (i.e. the call that built PreviousA).  Furthermore, the ConsType must be a 
00263 //  matching ConsType.  
00264 //  The ConsType values that cause TestMatrix to
00265 //  use PreviousA are:  
00266 //    RowMapColMap_VEPR, RowMapColMap_FEPR, RowMapColMap_NEPR, 
00267 //    WithGraph, CopyConstructor
00268 //  The matching ConsTypes are:  
00269 //               VariableEntriesPerRow, RowMapColMap_VEPR, NoEntriesPerRow, RowMapColMap_NEPR, WithGraph
00270 //               FixedEntriesPerRow, RowMapColMap_FEPR
00271 //
00272 //  TestMatrix() when called with ConsType=WithGraph, returns with PreviousA = 0 ;
00273 //              
00274 //
00275 enum ConsType { VariableEntriesPerRow, FixedEntriesPerRow, NoEntriesPerRow,
00276     RowMapColMap_VEPR, RowMapColMap_FEPR, RowMapColMap_NEPR, 
00277     WithGraph, CopyConstructor } ;
00278 
00279 int TestMatrix( Epetra_Comm& Comm, bool verbose, bool debug, 
00280     int NumMyElements, int MinSize, int MaxSize, 
00281     ConsType ConstructorType, bool ExtraBlocks, 
00282     bool insertlocal, 
00283     bool symmetric,
00284     Epetra_VbrMatrix** PreviousA ) {
00285 
00286 
00287   if (verbose) 
00288     cout << "MinSize         = " << MinSize << endl
00289    << "MaxSize         = " << MaxSize << endl
00290    << "ConstructorType = " << ConstructorType << endl
00291    << "ExtraBlocks     = " << ExtraBlocks << endl
00292    << "insertlocal     = " << insertlocal << endl
00293    << "symmetric       = " << symmetric << endl
00294    << "PreviousA       = " << PreviousA << endl;
00295 
00296   int ierr = 0, i, j, forierr = 0;
00297   int MyPID = Comm.MyPID();
00298   if (MyPID < 3) NumMyElements++;
00299   if (NumMyElements<2) NumMyElements = 2; // This value must be greater than one on each processor
00300 
00301   // Define pseudo-random block sizes using an Epetra Vector of random numbers
00302   Epetra_Map randmap(-1, NumMyElements, 0, Comm);
00303   Epetra_Vector randvec(randmap);
00304   randvec.Random(); // Fill with random numbers
00305   int * ElementSizeList = new int[NumMyElements];
00306   int SizeRange = MaxSize - MinSize + 1;
00307   double DSizeRange = SizeRange;
00308   
00309   const Epetra_BlockMap* rowmap = 0 ; 
00310   const Epetra_BlockMap* colmap = 0 ; 
00311   Epetra_CrsGraph* graph = 0 ; 
00312   if ( *PreviousA != 0 ) {
00313     
00314     rowmap = &((*PreviousA)->RowMap());
00315     colmap = &((*PreviousA)->ColMap());
00316   }
00317 
00318   ElementSizeList[0] = MaxSize;
00319   for (i=1; i<NumMyElements-1; i++) {
00320     int curSize = MinSize + (int) (DSizeRange * fabs(randvec[i]) + .99);
00321     ElementSizeList[i] = EPETRA_MAX(MinSize, EPETRA_MIN(MaxSize, curSize));
00322   }
00323   ElementSizeList[NumMyElements-1] = MaxSize;
00324 
00325   
00326 
00327   // Construct a Map
00328 
00329   int *randMyGlobalElements = randmap.MyGlobalElements();
00330 
00331   if ( ConstructorType == RowMapColMap_VEPR || 
00332        ConstructorType == RowMapColMap_FEPR ||
00333        ConstructorType == RowMapColMap_NEPR || 
00334        ConstructorType == WithGraph ) {
00335     rowmap->ElementSizeList( ElementSizeList ) ; 
00336   }
00337     
00338   Epetra_BlockMap Map (-1, NumMyElements, randMyGlobalElements, ElementSizeList, 0, Comm);
00339   
00340   // Get update list and number of local elements from newly created Map
00341   int NumGlobalElements = Map.NumGlobalElements();
00342   int * MyGlobalElements = Map.MyGlobalElements();
00343 
00344   // Create an integer vector NumNz that is used to build the Petra Matrix.
00345   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00346 
00347   int * NumNz = new int[NumMyElements];
00348 
00349   // We are building a block tridiagonal matrix
00350 
00351   for (i=0; i<NumMyElements; i++)
00352     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
00353       NumNz[i] = 2;
00354     else
00355       NumNz[i] = 3;
00356   // Create a Epetra_Matrix
00357 
00358   bool FixedNumEntries = false ;    // If FixedNumEntries == true, we add the upper left and lower right hand corners to 
00359                                     // our tri-diagonal matrix so that each row has exactly three entries
00360   bool HaveColMap = false ;         // Matrices constructed without a column map cannot use insertmy to create the matrix
00361   bool FixedBlockSize = ( MinSize == MaxSize ) ; 
00362   bool HaveGraph = false ; 
00363 
00364 
00365   Epetra_VbrMatrix* A ; 
00366 
00367   switch( ConstructorType ) { 
00368   case VariableEntriesPerRow:
00369     A = new Epetra_VbrMatrix( Copy, Map, NumNz ) ; 
00370     break; 
00371   case FixedEntriesPerRow:
00372     A = new Epetra_VbrMatrix( Copy, Map, 3 ) ; 
00373     FixedNumEntries = true; 
00374     break; 
00375   case NoEntriesPerRow:
00376     A = new Epetra_VbrMatrix( Copy, Map, 0 ) ; 
00377     break; 
00378   case RowMapColMap_VEPR: 
00379     A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, NumNz ) ; 
00380     HaveColMap = true; 
00381     break; 
00382   case RowMapColMap_FEPR: 
00383     A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, 3 ) ; 
00384     FixedNumEntries = true; 
00385     HaveColMap = true; 
00386     break; 
00387   case RowMapColMap_NEPR: 
00388     A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, 0 ) ; 
00389     HaveColMap = true; 
00390     break; 
00391   case WithGraph:
00392     graph =  new Epetra_CrsGraph( (*PreviousA)->Graph() );
00393     A = new Epetra_VbrMatrix( Copy, *graph );
00394     HaveGraph = true; 
00395     HaveColMap = true; 
00396     break; 
00397   case CopyConstructor:
00398     A = new Epetra_VbrMatrix( (**PreviousA) );
00399     HaveColMap = true; 
00400     break; 
00401   default:
00402     assert(false); 
00403  } 
00404 
00405   if ( insertlocal ) assert( HaveColMap );            // you can't insert local without a column map
00406   if ( ExtraBlocks ) assert( FixedBlockSize );        // ExtraBlocks is only supported for fixed block sizes
00407   if ( ExtraBlocks ) assert( ! HaveColMap );          // Matrices constructed with a column map cannot be extended
00408   if ( FixedNumEntries) assert( FixedBlockSize ) ;   // Can't handle a Fixed Number of Entries and a variable block size
00409   if ( insertlocal && HaveGraph ) assert( ! FixedNumEntries ) ;   // HaveGraph assumes the standard matrix shape
00410   if ( insertlocal && HaveGraph ) assert( ! ExtraBlocks ) ;   // HaveGraph assumes the standard matrix shape
00411 
00412 
00413   // WORK    Insert/Replace/Suminto  MY should fail here when there is no colmap 
00414 
00415 
00416   EPETRA_TEST_ERR(A->IndicesAreGlobal(),ierr);
00417   if ( ! HaveGraph ) EPETRA_TEST_ERR(A->IndicesAreLocal(),ierr);
00418   
00419   // Use an array of Epetra_SerialDenseMatrix objects to build VBR matrix
00420 
00421   Epetra_SerialDenseMatrix ** BlockEntries = new Epetra_SerialDenseMatrix*[SizeRange];
00422 
00423   // The array of dense matrices will increase in size from MinSize to MaxSize (defined above)
00424   for (int kr=0; kr<SizeRange; kr++) {
00425     BlockEntries[kr] = new Epetra_SerialDenseMatrix[SizeRange];
00426     int RowDim = MinSize+kr;
00427     for (int kc = 0; kc<SizeRange; kc++) {
00428       int ColDim = MinSize+kc;
00429       Epetra_SerialDenseMatrix * curmat = &(BlockEntries[kr][kc]);
00430       curmat->Shape(RowDim,ColDim);
00431       for (j=0; j < ColDim; j++)
00432   for (i=0; i < RowDim; i++) {
00433     BlockEntries[kr][kc][j][i] = -1.0;
00434     if (i==j && kr==kc) BlockEntries[kr][kc][j][i] = 9.0;
00435     else BlockEntries[kr][kc][j][i] = -1.0;
00436 
00437     if ( ! symmetric )  BlockEntries[kr][kc][j][i] += ((double) j)/10000.0;
00438   }
00439     }
00440   }
00441 
00442   // Add  rows one-at-a-time
00443 
00444   int *Indices = new int[3];
00445   int *MyIndices = new int[3];
00446   int *ColDims = new int[3];
00447   int NumEntries;
00448   int NumMyNonzeros = 0, NumMyEquations = 0;
00449 
00450 
00451 
00452   for (i=0; i<NumMyElements; i++) {
00453     int CurRow = MyGlobalElements[i];
00454     if ( HaveColMap ) { 
00455       assert ( i == rowmap->LID( CurRow ) ) ; 
00456       assert ( CurRow == rowmap->GID( i ) ) ; 
00457     }
00458     int RowDim = ElementSizeList[i]-MinSize;
00459     NumMyEquations += BlockEntries[RowDim][RowDim].M();
00460 
00461     if (CurRow==0)
00462       {
00463   Indices[0] = CurRow;
00464   Indices[1] = CurRow+1;
00465   if ( FixedNumEntries ) {
00466     Indices[2] = NumGlobalElements-1;
00467     ColDims[2] = 0 ; 
00468     assert( ElementSizeList[i] == MinSize );   // This is actually enforced above as well 
00469     NumEntries = 3;
00470   } else {
00471     NumEntries = 2;
00472   }
00473   ColDims[0] = ElementSizeList[i] - MinSize;
00474   ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
00475       }
00476     else if (CurRow == NumGlobalElements-1)
00477       {
00478   Indices[0] = CurRow-1;
00479   Indices[1] = CurRow;
00480   if ( FixedNumEntries ) {
00481     Indices[2] = 0;
00482     ColDims[2] = 0 ; 
00483     assert( ElementSizeList[i] == MinSize );   // This is actually enforced above as well 
00484     NumEntries = 3;
00485   } else {
00486     NumEntries = 2;
00487   }
00488   ColDims[0] = ElementSizeList[i-1] - MinSize;
00489   ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
00490       }
00491       else {
00492   Indices[0] = CurRow-1;
00493   Indices[1] = CurRow;
00494   Indices[2] = CurRow+1;
00495   NumEntries = 3;
00496   if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
00497   else ColDims[0] = ElementSizeList[i-1] - MinSize;
00498   ColDims[1] = ElementSizeList[i] - MinSize;
00499   // ElementSize on MyPID+1
00500   if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
00501   else ColDims[2] = ElementSizeList[i+1] - MinSize;
00502       }
00503     if ( insertlocal ) { 
00504       for ( int ii=0; ii < NumEntries; ii++ ) 
00505   MyIndices[ii] = colmap->LID( Indices[ii] ) ; 
00506       //      Epetra_MpiComm* MComm = dynamic_cast<Epetra_MpiComm*>( &Comm ) ;
00507       //  MComm->SetTracebackMode(1); // This should enable error traceback reporting
00508       if ( HaveGraph ) {
00509   EPETRA_TEST_ERR(!(A->BeginReplaceMyValues(rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
00510       } else { 
00511   EPETRA_TEST_ERR(!(A->BeginInsertMyValues(rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
00512       }
00513       //  MComm->SetTracebackMode(0); // This should shut down any error traceback reporting
00514     } else { 
00515       if ( HaveGraph ) { 
00516   EPETRA_TEST_ERR(!(A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
00517       } else { 
00518   //
00519   //  I, Ken Stanley, think the following call should return an error since it 
00520   //  makes no sense to insert a value with a local index in the absence of a 
00521   //  map indicating what that index means.  Instead, this call returns with an 
00522   //  error code of 0, but problems appear later.  
00523   //
00524   //  EPETRA_TEST_ERR((A->BeginInsertMyValues(CurRow, NumEntries, Indices)==0),ierr);  // Should fail
00525   EPETRA_TEST_ERR(!(A->BeginInsertGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
00526       }
00527     }
00528     forierr = 0;
00529     for (j=0; j < NumEntries; j++) {
00530       Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
00531       NumMyNonzeros += AD->M() * AD->N();   
00532       forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
00533     }
00534     EPETRA_TEST_ERR(forierr,ierr);
00535 
00536     A->EndSubmitEntries();
00537   }
00538 
00539   int NumMyBlockEntries = 3*NumMyElements;
00540   if ( ! FixedNumEntries ) { 
00541     if (A->LRID(0)>=0) NumMyBlockEntries--; // If I own first global row, then there is one less nonzero
00542     if (A->LRID(NumGlobalElements-1)>=0) NumMyBlockEntries--; // If I own last global row, then there is one less nonzero
00543   }
00544 
00545   if ( ExtraBlocks ) { 
00546 
00547 
00548     //
00549     //  Add a block to the matrix on each process.  
00550     //  The i index is chosen from among the block rows that this process owns (i.e. MyGlobalElements[0..NumMyElements-1])
00551     //  The j index is chosen from among all block columns 
00552     //  
00553     //  Bugs - does not support non-contiguous matrices
00554     //
00555     //  Adding more than one off diagonal block could have resulted in adding an off diagonal block
00556     //  twice, resulting in errors in NumMyNonzeros and NumMyBlockEntries
00557     //
00558     const int NumTries = 100; 
00559     Epetra_SerialDenseMatrix BlockIindex = Epetra_SerialDenseMatrix( NumTries, 1 );
00560     Epetra_SerialDenseMatrix BlockJindex = Epetra_SerialDenseMatrix( NumTries, 1 );
00561 
00562     BlockIindex.Random(); 
00563     BlockJindex.Random();
00564 
00565     BlockIindex.Scale( 1.0 * NumMyElements );
00566     BlockJindex.Scale( 1.0 * A->NumGlobalBlockCols() );
00567     bool OffDiagonalBlockAdded = false ; 
00568     for ( int ii=0; ii < NumTries && ! OffDiagonalBlockAdded; ii++ ) {
00569       int i = (int) BlockIindex[0][ii]; 
00570       int j = (int) BlockJindex[0][ii]; 
00571       if ( i < 0 ) i = - i ; 
00572       if ( j < 0 ) j = - j ; 
00573       assert( i >= 0 ) ;
00574       assert( j >= 0 ) ;
00575       assert( i < NumMyElements ) ;
00576       assert( j < A->NumGlobalBlockCols() ) ;
00577       int CurRow = MyGlobalElements[i];
00578       int Indices[1] ; 
00579       Indices[0] = j ; 
00580       Epetra_SerialDenseMatrix * AD = &(BlockEntries[0][0]);
00581       EPETRA_TEST_ERR(!(A->BeginInsertGlobalValues( CurRow, 1, Indices)==0),ierr);
00582 
00583       if ( CurRow < j-1 || CurRow > j+1 ) {
00584   OffDiagonalBlockAdded = true ; 
00585   NumMyNonzeros += AD->M() * AD->N();   
00586   NumMyBlockEntries++ ; 
00587       }
00588 
00589       //  EPETRA_TEST_ERR(!(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0), ierr);
00590       EPETRA_TEST_ERR(!(A->SubmitBlockEntry(*AD)==0), ierr);
00591       A->EndSubmitEntries();
00592     }
00593   }
00594 
00595   // Finish up
00596   if ( ! HaveGraph && ! insertlocal ) 
00597     EPETRA_TEST_ERR(!(A->IndicesAreGlobal()),ierr);
00598   EPETRA_TEST_ERR(!(A->FillComplete()==0),ierr);
00599   EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
00600   EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
00601   EPETRA_TEST_ERR(A->StorageOptimized(),ierr);
00602 
00603   A->OptimizeStorage();
00604   if ( FixedBlockSize ) {
00605     EPETRA_TEST_ERR(!(A->StorageOptimized()),ierr);
00606   } else { 
00607     //    EPETRA_TEST_ERR(A->StorageOptimized(),ierr);    //  Commented out until I figure out why it occasionally fails on one process
00608   }
00609   EPETRA_TEST_ERR(A->UpperTriangular(),ierr);
00610   EPETRA_TEST_ERR(A->LowerTriangular(),ierr);
00611 
00612 
00613 
00614 
00615 
00616 
00617 
00618   int NumGlobalBlockEntries ;
00619   int NumGlobalNonzeros, NumGlobalEquations;
00620   Comm.SumAll(&NumMyBlockEntries, &NumGlobalBlockEntries, 1);
00621   Comm.SumAll(&NumMyNonzeros, &NumGlobalNonzeros, 1);
00622 
00623   Comm.SumAll(&NumMyEquations, &NumGlobalEquations, 1);
00624 
00625   if (! ExtraBlocks ) {
00626     if ( FixedNumEntries ) {
00627       EPETRA_TEST_ERR( !( NumGlobalBlockEntries == (3*NumGlobalElements)), ierr );
00628     } else {
00629       EPETRA_TEST_ERR( !( NumGlobalBlockEntries == (3*NumGlobalElements-2)), ierr ); 
00630     }
00631   }
00632 
00633 
00634   EPETRA_TEST_ERR(!(check(*A, NumMyEquations, NumGlobalEquations, NumMyNonzeros, NumGlobalNonzeros, 
00635          NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries, 
00636          MyGlobalElements, verbose)==0),ierr);
00637 
00638   forierr = 0;
00639   if ( ! ExtraBlocks ) {
00640     if ( FixedNumEntries ) 
00641       for (i=0; i<NumMyElements; i++) forierr += !(A->NumGlobalBlockEntries(MyGlobalElements[i])==3);
00642     else
00643       for (i=0; i<NumMyElements; i++) forierr += !(A->NumGlobalBlockEntries(MyGlobalElements[i])==NumNz[i]);
00644   }
00645   EPETRA_TEST_ERR(forierr,ierr);
00646   forierr = 0;
00647   if ( ! ExtraBlocks ) 
00648     if ( FixedNumEntries ) 
00649       for (i=0; i<NumMyElements; i++) forierr += !(A->NumMyBlockEntries(i)==3);
00650     else
00651       for (i=0; i<NumMyElements; i++) forierr += !(A->NumMyBlockEntries(i)==NumNz[i]);
00652   EPETRA_TEST_ERR(forierr,ierr);
00653 
00654   if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
00655 
00656   delete [] NumNz;
00657 
00658 
00659   // Create vectors for Power method
00660 
00661   Epetra_Vector q(Map);
00662   Epetra_Vector z(Map);
00663   Epetra_Vector z_initial(Map);
00664   Epetra_Vector resid(Map);
00665 
00666   
00667   // Fill z with random Numbers 
00668   z_initial.Random();
00669 
00670   // variable needed for iteration
00671   double lambda = 0.0;
00672   int niters = 100;
00673   // int niters = 200;
00674   double tolerance = 1.0e-3;
00675 
00677 
00678   // Iterate
00679   Epetra_Flops flopcounter;
00680   A->SetFlopCounter(flopcounter);
00681   q.SetFlopCounter(*A);
00682   z.SetFlopCounter(*A);
00683   resid.SetFlopCounter(*A);
00684   z = z_initial;  // Start with common initial guess
00685   Epetra_Time timer(Comm);
00686   int ierr1 = power_method(false, *A, q, z, resid, &lambda, niters, tolerance, verbose);
00687   double elapsed_time = timer.ElapsedTime();
00688   double total_flops = flopcounter.Flops();
00689   double MFLOPs = total_flops/elapsed_time/1000000.0;
00690 
00691   if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
00692   if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
00693 
00695 
00696   // Solve transpose problem
00697 
00698   if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
00699         << endl;
00700   // Iterate
00701   lambda = 0.0;
00702   z = z_initial;  // Start with common initial guess
00703   flopcounter.ResetFlops();
00704   timer.ResetStartTime();
00705   ierr1 = power_method(true, *A, q, z, resid, &lambda, niters, tolerance, verbose);
00706   elapsed_time = timer.ElapsedTime();
00707   total_flops = flopcounter.Flops();
00708   MFLOPs = total_flops/elapsed_time/1000000.0;
00709 
00710   if (verbose) cout << "\n\nTotal MFLOPs for transpose solve = " << MFLOPs << endl<< endl;
00711   if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
00712 
00714 
00715   Epetra_CrsMatrix* OrigCrsA;
00716   ConvertVbrToCrs( A, OrigCrsA ) ; 
00717 
00718   // Increase diagonal dominance
00719 
00720   if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
00721         << endl;
00722 
00723   double AnormInf = -13 ;
00724   double AnormOne = -13 ;
00725   AnormInf = A->NormInf( ) ; 
00726   AnormOne = A->NormOne( ) ; 
00727 
00728   EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
00729   EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
00730 
00731   if (A->MyGlobalBlockRow(0)) {
00732     int numvals = A->NumGlobalBlockEntries(0);
00733     Epetra_SerialDenseMatrix ** Rowvals;
00734     int* Rowinds = new int[numvals];
00735     int  RowDim;
00736     A->ExtractGlobalBlockRowPointers(0, numvals, RowDim, numvals, Rowinds, 
00737             Rowvals); // Get A[0,:]
00738 
00739     for (i=0; i<numvals; i++) {
00740       if (Rowinds[i] == 0) {
00741   //  Rowvals[i]->A()[0] *= 10.0; // Multiply first diag value by 10.0
00742   Rowvals[i]->A()[0] += 1000.0; // Add 1000 to first diag value
00743       }
00744     }
00745     delete [] Rowinds;
00746   }
00747   //
00748   //  NormOne() and NormInf() will NOT return cached values
00749   //  See bug #1151
00750   //
00751 
00752   EPETRA_TEST_ERR( ! (AnormOne != A->NormOne( )), ierr ); 
00753   EPETRA_TEST_ERR( ! (AnormInf != A->NormInf( )), ierr );
00754   //
00755   //  On Process 0, let the class know that NormInf_ and NormOne_ are
00756   //  out of date.  
00757   // 
00758   if ( MyPID == 0 ) {
00759     EPETRA_TEST_ERR(!(A->BeginSumIntoGlobalValues( 0, 0, 0 )==0),ierr);
00760     EPETRA_TEST_ERR(  A->EndSubmitEntries(), ierr );
00761   }
00762   EPETRA_TEST_ERR( ! (AnormOne != A->NormOne( )), ierr ); 
00763   EPETRA_TEST_ERR( ! (AnormInf != A->NormInf( )), ierr ); 
00764   if ( MyPID == 0 ) 
00765   // Iterate (again)
00766   lambda = 0.0;
00767   z = z_initial;  // Start with common initial guess
00768   flopcounter.ResetFlops();
00769   timer.ResetStartTime();
00770   ierr1 = power_method(false, *A, q, z, resid, &lambda, niters, tolerance, verbose);
00771   elapsed_time = timer.ElapsedTime();
00772   total_flops = flopcounter.Flops();
00773   MFLOPs = total_flops/elapsed_time/1000000.0;
00774 
00775   if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
00776   if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
00777 
00778 
00779 
00781 
00782   EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
00783   EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
00784   // Solve transpose problem
00785 
00786   if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
00787         << endl;
00788 
00789   // Iterate (again)
00790   lambda = 0.0;
00791   z = z_initial;  // Start with common initial guess
00792   flopcounter.ResetFlops();
00793   timer.ResetStartTime();
00794   ierr1 = power_method(true, *A, q, z, resid, &lambda, niters, tolerance, verbose);
00795   elapsed_time = timer.ElapsedTime();
00796   total_flops = flopcounter.Flops();
00797   MFLOPs = total_flops/elapsed_time/1000000.0;
00798 
00799   if (verbose) cout << "\n\nTotal MFLOPs for tranpose of second solve = " << MFLOPs << endl<< endl;
00800   if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
00801 
00802 
00803   if (debug) Comm.Barrier();
00804 
00805   if (verbose) cout << "\n\n*****Comparing against CrsMatrix " << endl<< endl;
00806 
00807   Epetra_CrsMatrix* CrsA;
00808   ConvertVbrToCrs( A, CrsA ) ; 
00809 
00810   Epetra_Vector CrsX = Epetra_Vector( A->OperatorDomainMap(), false ) ; 
00811   Epetra_Vector CrsY = Epetra_Vector( A->OperatorRangeMap(), false ) ; 
00812   Epetra_Vector OrigCrsY = Epetra_Vector( A->OperatorRangeMap(), false ) ; 
00813   Epetra_Vector Y_Apply = Epetra_Vector( A->OperatorRangeMap(), false ) ; 
00814   Epetra_Vector x(Map);
00815   Epetra_Vector y(Map);
00816   Epetra_Vector orig_check_y(Map);
00817   Epetra_Vector Apply_check_y(Map);
00818   Epetra_Vector check_y(Map);
00819   Epetra_Vector check_ytranspose(Map);
00820 
00821   x.Random() ; 
00822   CrsX = x; 
00823   CrsY = CrsX;
00824   CrsA->Multiply( false, CrsX, CrsY ) ; 
00825   OrigCrsA->Multiply( false, CrsX, OrigCrsY ) ; 
00826 
00827 
00828   check_y = CrsY ; 
00829   orig_check_y = OrigCrsY ; 
00830   CrsA->Multiply( true, CrsX, CrsY ) ; 
00831   check_ytranspose = CrsY ; 
00832 
00833   EPETRA_TEST_ERR( checkmultiply( false, *A, x, check_y ), ierr ); 
00834 
00835   EPETRA_TEST_ERR( checkmultiply( true, *A, x, check_ytranspose ), ierr ); 
00836 
00837   if (! symmetric ) EPETRA_TEST_ERR( !checkmultiply( false, *A, x, check_ytranspose ), ierr );   // Just to confirm that A is not symmetric
00838 
00839   EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
00840   EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
00841 
00842   if (verbose) cout << "\n\n*****Try the Apply method " << endl<< endl;
00843 
00844 
00845   A->Apply( CrsX, Y_Apply ) ; 
00846   Apply_check_y = Y_Apply ; 
00847   EPETRA_TEST_ERR( checkmultiply( false, *A, x, Apply_check_y ), ierr ); 
00848 
00849   if (verbose) cout << "\n\n*****Multiply multivectors " << endl<< endl;
00850 
00851   const int NumVecs = 4 ; 
00852   
00853   Epetra_Map Amap = A->OperatorDomainMap() ; 
00854 
00855   //  Epetra_MultiVector CrsMX = Epetra_MultiVector( A->OperatorDomainMap(), false ) ; 
00856   Epetra_MultiVector CrsMX( Amap, NumVecs, false ) ; 
00857   Epetra_MultiVector CrsMY( A->OperatorRangeMap(), NumVecs, false ) ; 
00858   Epetra_MultiVector mx(Map, NumVecs);
00859   Epetra_MultiVector my(Map, NumVecs);
00860   Epetra_MultiVector check_my(Map, NumVecs);
00861   Epetra_MultiVector check_mytranspose(Map, NumVecs);
00862   mx.Random(); // Fill mx with random numbers
00863 #if 0 
00864   CrsMX = mx; 
00865   CrsA->Multiply( false, CrsMX, CrsMY ) ; 
00866 #else
00867   CrsMY = mx; 
00868   CrsA->Multiply( false, CrsMY, CrsMY ) ; 
00869 #endif
00870   check_my = CrsMY ; 
00871   CrsMY = mx; 
00872   CrsA->Multiply( true, CrsMY, CrsMY ) ; 
00873   check_mytranspose = CrsMY ; 
00874 
00875 
00876   EPETRA_TEST_ERR( checkmultiply( false, *A, mx, check_my ), ierr ); 
00877 
00878   EPETRA_TEST_ERR( checkmultiply( true, *A, mx, check_mytranspose ), ierr ); 
00879 
00880   
00881   EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
00882   EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
00883   if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
00884 
00885   //  B was changed to a pointer so that we could delete it before the
00886   //  underlying graph is deleted.  This was necessary before Bug #1116 was
00887   //  fixed, to avoid seg-faults. Bug #1116 is now fixed so this is no longer
00888   //  an issue.
00889   Epetra_VbrMatrix* B = new Epetra_VbrMatrix(*A);
00890 
00891   EPETRA_TEST_ERR(!(check(*B, NumMyEquations, NumGlobalEquations, NumMyNonzeros, NumGlobalNonzeros, 
00892          NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries, 
00893          MyGlobalElements, verbose)==0),ierr);
00894 
00895   EPETRA_TEST_ERR( ! ( A->StorageOptimized() == B->StorageOptimized() ), ierr ) ; 
00896 
00897   EPETRA_TEST_ERR( checkmultiply( false, *B, mx, check_my ), ierr ); 
00898 
00899 
00900   *B = *B;    // Should be harmless - check to make sure that it is
00901 
00902   EPETRA_TEST_ERR(!(check(*B, NumMyEquations, NumGlobalEquations, NumMyNonzeros, NumGlobalNonzeros, 
00903          NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries, 
00904          MyGlobalElements, verbose)==0),ierr);
00905 
00906   EPETRA_TEST_ERR( ! ( A->StorageOptimized() == B->StorageOptimized() ), ierr ) ; 
00907 
00908   EPETRA_TEST_ERR( checkmultiply( false, *B, mx, check_my ), ierr ); 
00909 
00910   AnormInf =  A->NormInf( );
00911   AnormOne =  A->NormOne( );
00912   EPETRA_TEST_ERR( ! (AnormOne == B->NormOne( )), ierr ); 
00913   EPETRA_TEST_ERR( ! (AnormInf == B->NormInf( )), ierr ); 
00914 
00915 
00916   Epetra_CrsMatrix* CrsB;
00917   ConvertVbrToCrs( B, CrsB ) ; 
00918 
00919   if (verbose) cout << "\n\n*****Testing PutScalar, LeftScale, RightScale, and ReplaceDiagonalValues" << endl<< endl;
00920   //
00921   //  Check PutScalar, 
00922   //
00923   B->PutScalar( 1.0 ) ; 
00924   CrsB->PutScalar( 1.0 ) ; 
00925   EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ; 
00926 
00927 
00928   check_y = CrsY ; 
00929   //
00930   EPETRA_TEST_ERR( B->ReplaceDiagonalValues( check_y ), ierr ) ; 
00931   EPETRA_TEST_ERR( CrsB->ReplaceDiagonalValues( CrsY ), ierr ) ; 
00932   EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ; 
00933 
00934   EPETRA_TEST_ERR( B->LeftScale( check_y ), ierr ) ; 
00935   EPETRA_TEST_ERR( CrsB->LeftScale( CrsY ), ierr ) ; 
00936   EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ; 
00937 
00938   EPETRA_TEST_ERR( B->RightScale( check_y ), ierr ) ; 
00939   EPETRA_TEST_ERR( CrsB->RightScale( CrsY ), ierr ) ; 
00940   EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ; 
00941 
00942   double B_norm_frob = B->NormFrobenius();
00943   double CrsB_norm_frob = CrsB->NormFrobenius();
00944   //need to use a fairly large tolerance when comparing the frobenius
00945   //norm from a vbr-matrix and a crs-matrix, because the point-entries
00946   //are visited in different orders during the norm calculation, causing
00947   //round-off differences to accumulate. That's why we choose 5.e-5
00948   //instead of a smaller number like 1.e-13 or so.
00949   if (fabs(B_norm_frob-CrsB_norm_frob) > 5.e-5) {
00950     std::cout << "fabs(B_norm_frob-CrsB_norm_frob): "
00951       << fabs(B_norm_frob-CrsB_norm_frob) << std::endl;
00952     std::cout << "VbrMatrix::NormFrobenius test FAILED."<<std::endl;
00953     EPETRA_TEST_ERR(-99, ierr);
00954   }
00955   if (verbose) std::cout << "\n\nVbrMatrix::NormFrobenius OK"<<std::endl<<std::endl;
00956 
00957   if (debug) Comm.Barrier();
00958 
00959   if (verbose) cout << "\n\n*****Testing post construction modifications" << endl<< endl;
00960   if (verbose) cout << "\n\n*****Replace methods should be OK" << endl<< endl;
00961 
00962   //  Check to see if we can restore the matrix to its original value
00963   // Does not work if ExtraBlocks is true
00964 
00965   EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
00966   EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
00967 
00968   if ( ! ExtraBlocks ) 
00969   {
00970     // Add  rows one-at-a-time
00971 
00972     int NumEntries;
00973     int NumMyNonzeros = 0, NumMyEquations = 0;
00974     
00975     for (i=0; i<NumMyElements; i++) {
00976       int CurRow = MyGlobalElements[i];
00977       int RowDim = ElementSizeList[i]-MinSize;
00978       NumMyEquations += BlockEntries[RowDim][RowDim].M();
00979       
00980       if (CurRow==0)
00981   {
00982     Indices[0] = CurRow;
00983     Indices[1] = CurRow+1;
00984     NumEntries = 2;
00985     ColDims[0] = ElementSizeList[i] - MinSize;
00986     ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
00987   }
00988       else if (CurRow == NumGlobalElements-1)
00989   {
00990     Indices[0] = CurRow-1;
00991     Indices[1] = CurRow;
00992     NumEntries = 2;
00993     ColDims[0] = ElementSizeList[i-1] - MinSize;
00994     ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
00995   }
00996       else {
00997   Indices[0] = CurRow-1;
00998   Indices[1] = CurRow;
00999   Indices[2] = CurRow+1;
01000   NumEntries = 3;
01001   if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
01002   else ColDims[0] = ElementSizeList[i-1] - MinSize;
01003   ColDims[1] = ElementSizeList[i] - MinSize;
01004   // ElementSize on MyPID+1
01005   if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
01006   else ColDims[2] = ElementSizeList[i+1] - MinSize;
01007       }
01008       EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
01009       EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
01010       EPETRA_TEST_ERR(!(A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
01011       forierr = 0;
01012       for (j=0; j < NumEntries; j++) {
01013   Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
01014   NumMyNonzeros += AD->M() * AD->N();   
01015   forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
01016       }
01017       EPETRA_TEST_ERR(forierr,ierr);
01018 
01019       A->EndSubmitEntries();
01020     }
01021     EPETRA_TEST_ERR( checkmultiply( false, *A, x, orig_check_y ), ierr );     
01022   }
01023 
01024   //
01025   //  Suminto should cause the matrix to be doubled 
01026   //
01027   if ( ! ExtraBlocks )   {
01028     // Add  rows one-at-a-time
01029 
01030     int NumEntries;
01031     int NumMyNonzeros = 0, NumMyEquations = 0;
01032     
01033     for (i=0; i<NumMyElements; i++) {
01034       int CurRow = MyGlobalElements[i];
01035       int RowDim = ElementSizeList[i]-MinSize;
01036       NumMyEquations += BlockEntries[RowDim][RowDim].M();
01037       
01038       if (CurRow==0)
01039   {
01040     Indices[0] = CurRow;
01041     Indices[1] = CurRow+1;
01042     NumEntries = 2;
01043     ColDims[0] = ElementSizeList[i] - MinSize;
01044     ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
01045   }
01046       else if (CurRow == NumGlobalElements-1)
01047   {
01048     Indices[0] = CurRow-1;
01049     Indices[1] = CurRow;
01050     NumEntries = 2;
01051     ColDims[0] = ElementSizeList[i-1] - MinSize;
01052     ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
01053   }
01054       else {
01055   Indices[0] = CurRow-1;
01056   Indices[1] = CurRow;
01057   Indices[2] = CurRow+1;
01058   NumEntries = 3;
01059   if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
01060   else ColDims[0] = ElementSizeList[i-1] - MinSize;
01061   ColDims[1] = ElementSizeList[i] - MinSize;
01062   // ElementSize on MyPID+1
01063   if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
01064   else ColDims[2] = ElementSizeList[i+1] - MinSize;
01065       }
01066       if ( insertlocal ) {
01067   for ( int ii=0; ii < NumEntries; ii++ ) 
01068     MyIndices[ii] = colmap->LID( Indices[ii] ) ; 
01069   EPETRA_TEST_ERR(!(A->BeginSumIntoMyValues( rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
01070       } else { 
01071   EPETRA_TEST_ERR(!(A->BeginSumIntoGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
01072       }
01073       forierr = 0;
01074       for (j=0; j < NumEntries; j++) {
01075   Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
01076   NumMyNonzeros += AD->M() * AD->N();   
01077   //  This has nothing to do with insertlocal, but that is a convenient bool to key on 
01078   if ( insertlocal ) {
01079     forierr += !(A->SubmitBlockEntry( *AD )==0);
01080   } else { 
01081     forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
01082   }
01083       }
01084       EPETRA_TEST_ERR(forierr,ierr);
01085 
01086       A->EndSubmitEntries();
01087     }
01088     
01089     orig_check_y.Scale(2.0) ;
01090 
01091     //  This will not work with FixedNumEntries unless we fix the fix the above loop to add the sorner elements to the tridi matrix
01092     if ( ! FixedNumEntries ) 
01093       EPETRA_TEST_ERR( checkmultiply( false, *A, x, orig_check_y ), ierr ); 
01094   }
01095 
01096 
01097   {for (int kr=0; kr<SizeRange; kr++) delete [] BlockEntries[kr];}
01098   delete [] BlockEntries;
01099   delete [] ColDims;
01100   delete [] MyIndices ;
01101   delete [] Indices;
01102   delete [] ElementSizeList;
01103 
01104 
01105 
01106 
01107   if (verbose) cout << "\n\n*****Insert methods should not be accepted" << endl<< endl;
01108 
01109   int One = 1;
01110   if (B->MyGRID(0)) EPETRA_TEST_ERR(!(B->BeginInsertGlobalValues(0, 1, &One)==-2),ierr);
01111 
01112   Epetra_Vector checkDiag(B->RowMap());
01113   forierr = 0;
01114   int NumMyEquations1 = B->NumMyRows();
01115   double two1 = 2.0;
01116 
01117   // Test diagonal replacement and extraction methods
01118 
01119   forierr = 0;
01120   for (i=0; i<NumMyEquations1; i++) checkDiag[i]=two1;
01121   EPETRA_TEST_ERR(forierr,ierr);
01122   
01123   EPETRA_TEST_ERR(!(B->ReplaceDiagonalValues(checkDiag)==0),ierr);
01124   
01125   Epetra_Vector checkDiag1(B->RowMap());
01126   EPETRA_TEST_ERR(!(B->ExtractDiagonalCopy(checkDiag1)==0),ierr);
01127   
01128   forierr = 0;
01129   for (i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
01130   EPETRA_TEST_ERR(forierr,ierr);
01131 
01132   if (verbose) cout << "\n\nDiagonal extraction and replacement OK.\n\n" << endl;
01133 
01134   double originfnorm = B->NormInf();
01135   double origonenorm = B->NormOne();
01136   EPETRA_TEST_ERR(!(B->Scale(4.0)==0),ierr);
01137   //  EPETRA_TEST_ERR((B->NormOne()!=origonenorm),ierr);
01138   EPETRA_TEST_ERR(!(B->NormInf()==4.0 * originfnorm),ierr);
01139   EPETRA_TEST_ERR(!(B->NormOne()==4.0 * origonenorm),ierr);
01140 
01141   if (verbose) cout << "\n\nMatrix scale OK.\n\n" << endl;
01142   
01143   // Testing VbrRowMatrix adapter
01144   Epetra_VbrRowMatrix ARowMatrix(A);
01145 
01146   EPETRA_TEST_ERR(checkVbrRowMatrix(*A, ARowMatrix, verbose),ierr);
01147 
01148   if ( PreviousA ) 
01149     delete *PreviousA; 
01150   
01151   //
01152   //  The following code deals with the fact that A has to be delete before graph is, when 
01153   //  A is built with a contructor that takes a graph as input.
01154   //
01155   delete B;
01156   if ( graph ) {
01157     delete A;
01158     delete graph ; 
01159     *PreviousA = 0 ; 
01160   } else { 
01161     *PreviousA = A ; 
01162   }
01163   
01164   delete CrsA;
01165   delete CrsB;
01166   delete OrigCrsA ;
01167 
01168   return ierr; 
01169 
01170 }
01171 
01172 
01173 int main(int argc, char *argv[])
01174 {
01175   int ierr = 0;
01176   bool debug = false;
01177 
01178 #ifdef EPETRA_MPI
01179   MPI_Init(&argc,&argv);
01180   Epetra_MpiComm Comm( MPI_COMM_WORLD );
01181 #else
01182   Epetra_SerialComm Comm;
01183 #endif
01184 
01185   int MyPID = Comm.MyPID();
01186   int NumProc = Comm.NumProc();
01187   bool verbose = false;
01188 
01189   // Check if we should print results to standard out
01190   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
01191 
01192   // char tmp;
01193   //  int rank = Comm.MyPID();
01194   //  if (rank==0) cout << "Press any key to continue..."<< endl;
01195   //  if (rank==0) cin >> tmp;
01196   //  Comm.Barrier();
01197 
01198   //  if ( ! verbose )  
01199   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
01200 
01201   if (verbose && MyPID==0)
01202     cout << Epetra_Version() << endl << endl;
01203 
01204   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
01205               << " is alive."<<endl;
01206 
01207   // Redefine verbose to only print on PE 0
01208   if (verbose && Comm.MyPID()!=0) verbose = false;
01209 
01210 
01211   int NumMyElements = 400;
01212   //  int NumMyElements = 3; 
01213   int MinSize = 2;
01214   int MaxSize = 8;
01215   bool NoExtraBlocks = false; 
01216   bool symmetric = true; 
01217   bool NonSymmetric = false;
01218   bool NoInsertLocal = false ; 
01219   bool InsertLocal = true ; 
01220 
01221   Epetra_VbrMatrix* PreviousA = 0 ; 
01222 
01223   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 1, 1, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );  
01224 
01225   //
01226   //  Check the various constructors
01227   //  
01228   
01229   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
01230 
01231   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
01232 
01233   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );   
01234 
01235   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA ); 
01236 
01237   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_VEPR, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
01238 
01239   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01240 
01241   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_VEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
01242 
01243   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01244 
01245   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, WithGraph, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01246   assert ( PreviousA == 0 ) ; 
01247 
01248 
01249   //
01250   //  Check some various options
01251   //
01252   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01253 
01254   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01255 
01256   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01257   assert ( PreviousA == 0 ) ; 
01258 
01259   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01260 
01261   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01262 
01263   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01264   assert ( PreviousA == 0 ) ; 
01265 
01266 
01267 
01268   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );   
01269 
01270   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, RowMapColMap_FEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );   
01271 
01272 
01273 
01274   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01275 
01276   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 3, 3, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01277 
01278   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01279 
01280   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 5, 5, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01281 
01282   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 6, 6, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01283 
01284   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 7, 7, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01285 
01286   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 8, 8, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01287 
01288   //  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01289 
01290 
01291   delete PreviousA;
01292 
01293   /*
01294   if (verbose) {
01295     // Test ostream << operator (if verbose)
01296     // Construct a Map that puts 2 equations on each PE
01297     
01298     int NumMyElements1 = 2;
01299     int NumMyElements1 = NumMyElements1;
01300     int NumGlobalElements1 = NumMyElements1*NumProc;
01301 
01302     Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
01303     
01304     // Get update list and number of local equations from newly created Map
01305     int * MyGlobalElements1 = new int[Map1.NumMyElements()];
01306     Map1.MyGlobalElements(MyGlobalElements1);
01307     
01308     // Create an integer vector NumNz that is used to build the Petra Matrix.
01309     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
01310     
01311     int * NumNz1 = new int[NumMyElements1];
01312     
01313     // We are building a tridiagonal matrix where each row has (-1 2 -1)
01314     // So we need 2 off-diagonal terms (except for the first and last equation)
01315     
01316     for (i=0; i<NumMyElements1; i++)
01317       if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalElements1-1)
01318   NumNz1[i] = 1;
01319       else
01320   NumNz1[i] = 2;
01321     
01322     // Create a Epetra_Matrix
01323     
01324     Epetra_VbrMatrix A1(Copy, Map1, NumNz1);
01325     
01326     // Add  rows one-at-a-time
01327     // Need some vectors to help
01328     // Off diagonal Values will always be -1
01329     
01330     
01331     int *Indices1 = new int[2];
01332     double two1 = 2.0;
01333     int NumEntries1;
01334 
01335     forierr = 0;
01336     for (i=0; i<NumMyElements1; i++)
01337       {
01338   if (MyGlobalElements1[i]==0)
01339     {
01340       Indices1[0] = 1;
01341       NumEntries1 = 1;
01342     }
01343   else if (MyGlobalElements1[i] == NumGlobalElements1-1)
01344     {
01345       Indices1[0] = NumGlobalElements1-2;
01346       NumEntries1 = 1;
01347     }
01348   else
01349     {
01350       Indices1[0] = MyGlobalElements1[i]-1;
01351       Indices1[1] = MyGlobalElements1[i]+1;
01352       NumEntries1 = 2;
01353     }
01354         forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1)==0);
01355   forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i)>0); // Put in the diagonal entry
01356       }
01357     EPETRA_TEST_ERR(forierr,ierr);
01358     // Finish up
01359     EPETRA_TEST_ERR(!(A1.FillComplete()==0),ierr);
01360     
01361     if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
01362     cout << A1 << endl;
01363     
01364     // Release all objects
01365     delete [] NumNz1;
01366     delete [] Values1;
01367     delete [] Indices1;
01368     delete [] MyGlobalElements1;
01369 
01370   }
01371   */
01372 
01373   EPETRA_TEST_ERR( checkMergeRedundantEntries(Comm, verbose), ierr);
01374 
01375   EPETRA_TEST_ERR( checkExtractMyRowCopy(Comm, verbose), ierr);
01376 
01377   EPETRA_TEST_ERR( checkMatvecSameVectors(Comm, verbose), ierr);
01378 
01379   EPETRA_TEST_ERR( checkEarlyDelete(Comm, verbose), ierr);
01380 
01381   if (verbose) {
01382     if (ierr==0) cout << "All VbrMatrix tests OK" << endl;
01383     else cout << ierr << " errors encountered." << endl;
01384   }
01385 
01386 #ifdef EPETRA_MPI
01387   MPI_Finalize() ;
01388 #endif
01389 
01390 /* end main
01391 */
01392 return ierr ;
01393 }
01394 
01395 int power_method(bool TransA, Epetra_VbrMatrix& A, 
01396      Epetra_MultiVector& q,
01397      Epetra_MultiVector& z, 
01398      Epetra_MultiVector& resid, 
01399      double * lambda, int niters, double tolerance,
01400      bool verbose) {  
01401 
01402   // variable needed for iteration
01403   double normz, residual;
01404 
01405   int ierr = 1;
01406 
01407   for (int iter = 0; iter < niters; iter++)
01408     {
01409       z.Norm2(&normz); // Compute 2-norm of z
01410       q.Scale(1.0/normz, z);
01411       A.Multiply(TransA, q, z); // Compute z = A*q
01412       q.Dot(z, lambda); // Approximate maximum eigenvaluE
01413       if (iter%100==0 || iter+1==niters)
01414   {
01415     resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
01416     resid.Norm2(&residual);
01417     if (verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda 
01418            << "  Residual of A*q - lambda*q = " << residual << endl;
01419   } 
01420       if (residual < tolerance) {
01421   ierr = 0;
01422   break;
01423       }
01424     }
01425   return(ierr);
01426 }
01427 int check(Epetra_VbrMatrix& A, 
01428     int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, 
01429     int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1, 
01430     int * MyGlobalElements, bool verbose) {
01431 
01432   int ierr = 0, forierr = 0;
01433   // Test query functions
01434 
01435   int NumMyRows = A.NumMyRows();
01436   if (verbose) cout << "\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
01437   // TEMP
01438   //if (verbose) cout << "\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
01439   if (verbose) cout << "\n\nNumber of local Rows is = " << NumMyRows << endl<< endl;
01440   if (verbose) cout << "\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
01441 
01442   EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);
01443 
01444   int NumMyNonzeros = A.NumMyNonzeros();
01445   if (verbose) cout << "\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
01446   //if (verbose) cout << "                            Should  = " << NumMyNonzeros1 << endl<< endl;
01447 
01448 
01449   if ( NumMyNonzeros != NumMyNonzeros1 ) {
01450     cout << " MyPID = " << A.Comm().MyPID() 
01451    << " NumMyNonzeros = " << NumMyNonzeros 
01452    << " NumMyNonzeros1 = " << NumMyNonzeros1 << endl ; 
01453   }
01454 
01455   EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);
01456 
01457   int NumGlobalRows = A.NumGlobalRows();
01458   if (verbose) cout << "\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
01459 
01460   EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);
01461 
01462   int NumGlobalNonzeros = A.NumGlobalNonzeros();
01463   if (verbose) cout << "\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
01464 
01465   if ( NumGlobalNonzeros != NumGlobalNonzeros1 ) {
01466     cout << " MyPID = " << A.Comm().MyPID() 
01467    << " NumGlobalNonzeros = " << NumGlobalNonzeros 
01468    << " NumGlobalNonzeros1 = " << NumGlobalNonzeros1 << endl ; 
01469   }
01470   EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);
01471 
01472   int NumMyBlockRows = A.NumMyBlockRows();
01473   if (verbose) cout << "\n\nNumber of local Block Rows = " << NumMyBlockRows << endl<< endl;
01474 
01475   EPETRA_TEST_ERR(!(NumMyBlockRows==NumMyBlockRows1),ierr);
01476 
01477   int NumMyBlockNonzeros = A.NumMyBlockEntries();
01478 
01479   if (verbose) cout << "\n\nNumber of local Nonzero Block entries = " << NumMyBlockNonzeros << endl<< endl;
01480   if (verbose) cout << "\n\nNumber of local Nonzero Block entries 1 = " << NumMyBlockNonzeros1 << endl<< endl;
01481 
01482   EPETRA_TEST_ERR(!(NumMyBlockNonzeros==NumMyBlockNonzeros1),ierr);
01483 
01484   int NumGlobalBlockRows = A.NumGlobalBlockRows();
01485   if (verbose) cout << "\n\nNumber of global Block Rows = " << NumGlobalBlockRows << endl<< endl;
01486 
01487   EPETRA_TEST_ERR(!(NumGlobalBlockRows==NumGlobalBlockRows1),ierr);
01488 
01489   int NumGlobalBlockNonzeros = A.NumGlobalBlockEntries();
01490   if (verbose) cout << "\n\nNumber of global Nonzero Block entries = " << NumGlobalBlockNonzeros << endl<< endl;
01491 
01492   EPETRA_TEST_ERR(!(NumGlobalBlockNonzeros==NumGlobalBlockNonzeros1),ierr);
01493 
01494   
01495   // Test RowMatrix interface implementations
01496   int RowDim, NumBlockEntries, * BlockIndices;
01497   Epetra_SerialDenseMatrix ** Values;
01498   // Get View of last block row
01499   A.ExtractMyBlockRowView(NumMyBlockRows-1, RowDim, NumBlockEntries,
01500         BlockIndices, Values);
01501   int NumMyEntries1 = 0;
01502   {for (int i=0; i < NumBlockEntries; i++) NumMyEntries1 += Values[i]->N();}
01503   int NumMyEntries;
01504   A.NumMyRowEntries(NumMyRows-1, NumMyEntries);
01505   if (verbose) {
01506     cout << "\n\nNumber of nonzero values in last row = "
01507    << NumMyEntries << endl<< endl;
01508   }
01509 
01510   EPETRA_TEST_ERR(!(NumMyEntries==NumMyEntries1),ierr);
01511   
01512   // Other binary tests
01513 
01514   EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
01515   EPETRA_TEST_ERR(!(A.Filled()),ierr);
01516   EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID())),ierr);
01517   EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID())),ierr);
01518   EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID()),ierr);
01519   EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID()),ierr);
01520   EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
01521   EPETRA_TEST_ERR(!(A.MyLRID(NumMyBlockRows-1)),ierr);
01522   EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
01523   EPETRA_TEST_ERR(A.MyLRID(NumMyBlockRows),ierr);
01524 
01525     
01526   int i, j;
01527   int MaxNumBlockEntries = A.MaxNumBlockEntries();
01528 
01529   // Pointer Extraction approach
01530 
01531   //   Local index
01532   int MyPointersRowDim, MyPointersNumBlockEntries;
01533   int * MyPointersBlockIndices = new int[MaxNumBlockEntries];
01534   Epetra_SerialDenseMatrix **MyPointersValuesPointers;
01535   //   Global Index
01536   int GlobalPointersRowDim, GlobalPointersNumBlockEntries;
01537   int * GlobalPointersBlockIndices = new int[MaxNumBlockEntries];
01538   Epetra_SerialDenseMatrix **GlobalPointersValuesPointers;
01539 
01540   // Copy Extraction approach
01541 
01542   //   Local index
01543   int MyCopyRowDim, MyCopyNumBlockEntries;
01544   int * MyCopyBlockIndices = new int[MaxNumBlockEntries];
01545   int * MyCopyColDims = new int[MaxNumBlockEntries];
01546   int * MyCopyLDAs = new int[MaxNumBlockEntries];
01547   int MaxRowDim = A.MaxRowDim();
01548   int MaxColDim = A.MaxColDim();
01549   int MyCopySizeOfValues = MaxRowDim*MaxColDim;
01550   double ** MyCopyValuesPointers = new double*[MaxNumBlockEntries];
01551   for (i=0; i<MaxNumBlockEntries; i++)
01552     MyCopyValuesPointers[i] = new double[MaxRowDim*MaxColDim];
01553 
01554   //   Global Index
01555   int GlobalCopyRowDim, GlobalCopyNumBlockEntries;
01556   int * GlobalCopyBlockIndices = new int[MaxNumBlockEntries];
01557   int * GlobalCopyColDims = new int[MaxNumBlockEntries];
01558   int * GlobalCopyLDAs = new int[MaxNumBlockEntries];
01559   
01560   int GlobalMaxRowDim = A.GlobalMaxRowDim();
01561   int GlobalMaxColDim = A.GlobalMaxColDim();
01562   int GlobalCopySizeOfValues = GlobalMaxRowDim*GlobalMaxColDim;
01563   double ** GlobalCopyValuesPointers = new double*[MaxNumBlockEntries];
01564   for (i=0; i<MaxNumBlockEntries; i++)
01565     GlobalCopyValuesPointers[i] = new double[GlobalMaxRowDim*GlobalMaxColDim];
01566 
01567   // View Extraction approaches
01568 
01569   //   Local index (There is no global view available)
01570   int MyView1RowDim, MyView1NumBlockEntries;
01571   int * MyView1BlockIndices;
01572   Epetra_SerialDenseMatrix **MyView1ValuesPointers = new Epetra_SerialDenseMatrix*[MaxNumBlockEntries];
01573 
01574   //   Local index version 2 (There is no global view available)
01575   int MyView2RowDim, MyView2NumBlockEntries;
01576   int * MyView2BlockIndices;
01577   Epetra_SerialDenseMatrix **MyView2ValuesPointers;
01578 
01579 
01580   // For each row, test six approaches to extracting data from a given local index matrix
01581   forierr = 0;
01582   for (i=0; i<NumMyBlockRows; i++) {
01583     int MyRow = i;
01584     int GlobalRow = A.GRID(i);
01585     // Get a copy of block indices in local index space, pointers to everything else
01586     EPETRA_TEST_ERR( A.ExtractMyBlockRowPointers(MyRow, MaxNumBlockEntries, MyPointersRowDim, 
01587         MyPointersNumBlockEntries, MyPointersBlockIndices,
01588         MyPointersValuesPointers), ierr );
01589     // Get a copy of block indices in local index space, pointers to everything else
01590     EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(GlobalRow, MaxNumBlockEntries, GlobalPointersRowDim, 
01591             GlobalPointersNumBlockEntries, GlobalPointersBlockIndices,
01592             GlobalPointersValuesPointers), ierr ) ;
01593 
01594     // Initiate a copy of block row in local index space.
01595     EPETRA_TEST_ERR( A.BeginExtractMyBlockRowCopy(MyRow, MaxNumBlockEntries, MyCopyRowDim, 
01596          MyCopyNumBlockEntries, MyCopyBlockIndices,
01597          MyCopyColDims), ierr);
01598     // Copy Values
01599     for (j=0; j<MyCopyNumBlockEntries; j++) {
01600       EPETRA_TEST_ERR( A.ExtractEntryCopy(MyCopySizeOfValues, MyCopyValuesPointers[j], MaxRowDim, false), ierr) ;
01601       MyCopyLDAs[j] = MaxRowDim;
01602     }
01603 
01604     // Initiate a copy of block row in global index space.
01605     EPETRA_TEST_ERR( A.BeginExtractGlobalBlockRowCopy(GlobalRow, MaxNumBlockEntries, GlobalCopyRowDim, 
01606             GlobalCopyNumBlockEntries, GlobalCopyBlockIndices,
01607             GlobalCopyColDims), ierr ) ;
01608     // Copy Values
01609     for (j=0; j<GlobalCopyNumBlockEntries; j++) {
01610       EPETRA_TEST_ERR( A.ExtractEntryCopy(GlobalCopySizeOfValues, GlobalCopyValuesPointers[j], GlobalMaxRowDim, false), ierr );
01611       GlobalCopyLDAs[j] = GlobalMaxRowDim;
01612     }
01613 
01614     // Initiate a view of block row in local index space (Version 1)
01615     EPETRA_TEST_ERR( A.BeginExtractMyBlockRowView(MyRow, MyView1RowDim, 
01616          MyView1NumBlockEntries, MyView1BlockIndices), ierr ) ;
01617     // Set pointers to values
01618     for (j=0; j<MyView1NumBlockEntries; j++) 
01619       EPETRA_TEST_ERR ( A.ExtractEntryView(MyView1ValuesPointers[j]), ierr ) ;
01620 
01621 
01622     // Extract a view of block row in local index space (version 2)
01623     EPETRA_TEST_ERR( A.ExtractMyBlockRowView(MyRow, MyView2RowDim, 
01624           MyView2NumBlockEntries, MyView2BlockIndices,
01625           MyView2ValuesPointers), ierr );
01626 
01627     forierr += !(MyPointersNumBlockEntries==GlobalPointersNumBlockEntries);
01628     forierr += !(MyPointersNumBlockEntries==MyCopyNumBlockEntries);
01629     forierr += !(MyPointersNumBlockEntries==GlobalCopyNumBlockEntries);
01630     forierr += !(MyPointersNumBlockEntries==MyView1NumBlockEntries);
01631     forierr += !(MyPointersNumBlockEntries==MyView2NumBlockEntries);
01632     for (j=1; j<MyPointersNumBlockEntries; j++) {
01633       forierr += !(MyCopyBlockIndices[j-1]<MyCopyBlockIndices[j]);
01634       forierr += !(MyView1BlockIndices[j-1]<MyView1BlockIndices[j]);
01635       forierr += !(MyView2BlockIndices[j-1]<MyView2BlockIndices[j]);
01636 
01637       forierr += !(GlobalPointersBlockIndices[j]==A.GCID(MyPointersBlockIndices[j]));
01638       forierr += !(A.LCID(GlobalPointersBlockIndices[j])==MyPointersBlockIndices[j]);
01639       forierr += !(GlobalPointersBlockIndices[j]==GlobalCopyBlockIndices[j]);
01640       
01641       Epetra_SerialDenseMatrix* my = MyPointersValuesPointers[j];
01642       Epetra_SerialDenseMatrix* global = GlobalPointersValuesPointers[j];
01643 
01644       Epetra_SerialDenseMatrix* myview1 = MyView1ValuesPointers[j];
01645       Epetra_SerialDenseMatrix* myview2 = MyView2ValuesPointers[j];
01646 
01647       forierr += !(CompareValues(my->A(), my->LDA(), 
01648          MyPointersRowDim, my->N(), 
01649          global->A(), global->LDA(), 
01650          GlobalPointersRowDim, global->N())==0);
01651       forierr += !(CompareValues(my->A(), my->LDA(), 
01652          MyPointersRowDim, my->N(), 
01653          MyCopyValuesPointers[j], MyCopyLDAs[j], 
01654          MyCopyRowDim, MyCopyColDims[j])==0);
01655       forierr += !(CompareValues(my->A(), my->LDA(), 
01656          MyPointersRowDim, my->N(), 
01657          GlobalCopyValuesPointers[j], GlobalCopyLDAs[j], 
01658          GlobalCopyRowDim, GlobalCopyColDims[j])==0);
01659       forierr += !(CompareValues(my->A(), my->LDA(), 
01660          MyPointersRowDim, my->N(), 
01661          myview1->A(), myview1->LDA(), 
01662          MyView1RowDim, myview1->N())==0);
01663       forierr += !(CompareValues(my->A(), my->LDA(),
01664          MyPointersRowDim, my->N(),
01665          myview2->A(), myview2->LDA(),
01666          MyView2RowDim, myview2->N())==0);
01667     }
01668   }
01669   EPETRA_TEST_ERR(forierr,ierr);
01670 
01671   // GlobalRowView should be illegal (since we have local indices)
01672   EPETRA_TEST_ERR(!(A.BeginExtractGlobalBlockRowView(A.GRID(0), MyView1RowDim, 
01673                  MyView1NumBlockEntries,
01674                  MyView1BlockIndices)==-2),ierr);
01675   
01676   // Extract a view of block row in local index space (version 2)
01677   EPETRA_TEST_ERR(!(A.ExtractGlobalBlockRowView(A.GRID(0), MyView2RowDim, 
01678              MyView2NumBlockEntries, MyView2BlockIndices,
01679              MyView2ValuesPointers)==-2),ierr);
01680   
01681   delete [] MyPointersBlockIndices;
01682   delete [] GlobalPointersBlockIndices;
01683   delete [] MyCopyBlockIndices;
01684   delete [] MyCopyColDims;
01685   delete [] MyCopyLDAs;
01686   for (i=0; i<MaxNumBlockEntries; i++) delete [] MyCopyValuesPointers[i];
01687   delete [] MyCopyValuesPointers;
01688   delete [] GlobalCopyBlockIndices;
01689   delete [] GlobalCopyColDims;
01690   delete [] GlobalCopyLDAs;
01691   for (i=0; i<MaxNumBlockEntries; i++) delete [] GlobalCopyValuesPointers[i];
01692   delete [] GlobalCopyValuesPointers;
01693   delete [] MyView1ValuesPointers;
01694   if (verbose) cout << "\n\nRows sorted check OK" << endl<< endl;
01695   
01696   return ierr;
01697 }
01698 
01699 //=============================================================================
01700 int CompareValues(double * A, int LDA, int NumRowsA, int NumColsA, 
01701       double * B, int LDB, int NumRowsB, int NumColsB) {
01702   
01703   int i, j, ierr = 0, forierr = 0;
01704   double * ptr1 = B;
01705   double * ptr2;
01706   
01707   if (NumRowsA!=NumRowsB) EPETRA_TEST_ERR(-2,ierr);
01708   if (NumColsA!=NumColsB) EPETRA_TEST_ERR(-3,ierr);
01709  
01710 
01711   forierr = 0;
01712   for (j=0; j<NumColsA; j++) {
01713     ptr1 = B + j*LDB;
01714     ptr2 = A + j*LDA;
01715     for (i=0; i<NumRowsA; i++) forierr += (*ptr1++ != *ptr2++);
01716   }
01717   EPETRA_TEST_ERR(forierr,ierr);
01718   return ierr;
01719 }
01720 
01721 int checkMergeRedundantEntries(Epetra_Comm& comm, bool verbose)
01722 {
01723   int numProcs = comm.NumProc();
01724   int localProc = comm.MyPID();
01725 
01726   int myFirstRow = localProc*3;
01727   int myLastRow = myFirstRow+2;
01728   int numMyRows = myLastRow - myFirstRow + 1;
01729   int numGlobalRows = numProcs*numMyRows;
01730   int i,j, ierr;
01731 
01732   //We'll set up a matrix which is globally block-diagonal, i.e., on each
01733   //processor the list of columns == list of rows.
01734   //Set up a list of column-indices which is twice as long as it should be,
01735   //and its contents will be the list of local rows, repeated twice.
01736   int numCols = 2*numMyRows;
01737   int* myCols = new int[numCols];
01738 
01739   int col = myFirstRow;
01740   for(i=0; i<numCols; ++i) {
01741     myCols[i] = col++;
01742     if (col > myLastRow) col = myFirstRow;
01743   }
01744 
01745   int elemSize = 2;
01746   int indexBase = 0;
01747 
01748   Epetra_BlockMap map(numGlobalRows, numMyRows,
01749           elemSize, indexBase, comm);
01750 
01751   Epetra_VbrMatrix A(Copy, map, numCols);
01752 
01753   Epetra_MultiVector x(map, 1), y(map, 1);
01754   x.PutScalar(1.0);
01755 
01756   Epetra_MultiVector x3(map, 3), y3(map, 3);
01757   x.PutScalar(1.0);
01758 
01759   double* coef = new double[elemSize*elemSize];
01760   for(i=0; i<elemSize*elemSize; ++i) {
01761     coef[i] = 0.5;
01762   }
01763 
01764   //we're going to insert each row twice, with coef values of 0.5. So after
01765   //FillComplete, which internally calls MergeRedundantEntries, the
01766   //matrix should contain 1.0 in each entry.
01767 
01768   for(i=myFirstRow; i<=myLastRow; ++i) {
01769     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, numCols, myCols), ierr);
01770 
01771     for(j=0; j<numCols; ++j) {
01772       EPETRA_TEST_ERR( A.SubmitBlockEntry(coef, elemSize,
01773             elemSize, elemSize), ierr);
01774     }
01775 
01776     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
01777   }
01778 
01779   EPETRA_TEST_ERR( A.FillComplete(), ierr);
01780 
01781   delete [] coef;
01782 
01783   if (verbose) cout << "Multiply x"<<endl;
01784   EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr );
01785 
01786 
01787   //Next we're going to extract pointers-to-block-rows and check values to make
01788   //sure that the internal method Epetra_VbrMatrix::mergeRedundantEntries()
01789   //worked correctly. 
01790   //At the same time, we're also going to create another VbrMatrix which will
01791   //be a View of the matrix we've already created. This will serve to provide
01792   //more test coverage of the VbrMatrix code.
01793 
01794   int numBlockEntries = 0;
01795   int RowDim;
01796   int** BlockIndices = new int*[numMyRows];
01797   Epetra_SerialDenseMatrix** Values;
01798   Epetra_VbrMatrix Aview(View, map, numMyRows);
01799 
01800   for(i=myFirstRow; i<=myLastRow; ++i) {
01801     BlockIndices[i-myFirstRow] = new int[numCols];
01802     EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(i, numCols,
01803                  RowDim, numBlockEntries,
01804                  BlockIndices[i-myFirstRow],
01805                  Values), ierr);
01806 
01807     EPETRA_TEST_ERR( Aview.BeginInsertGlobalValues(i, numBlockEntries,
01808                 BlockIndices[i-myFirstRow]), ierr);
01809 
01810     if (numMyRows != numBlockEntries) return(-1);
01811     if (RowDim != elemSize) return(-2);
01812     for(j=0; j<numBlockEntries; ++j) {
01813       if (Values[j]->A()[0] != 1.0) {
01814   cout << "Row " << i << " Values["<<j<<"][0]: "<< Values[j][0]
01815        << " should be 1.0" << endl;
01816   return(-3); //comment-out this return to de-activate this test
01817       }
01818 
01819       EPETRA_TEST_ERR( Aview.SubmitBlockEntry(Values[j]->A(),
01820                 Values[j]->LDA(),
01821                 Values[j]->M(),
01822                 Values[j]->N()), ierr);
01823     }
01824 
01825     EPETRA_TEST_ERR( Aview.EndSubmitEntries(), ierr);
01826   }
01827 
01828   EPETRA_TEST_ERR( Aview.FillComplete(), ierr);
01829 
01830   //So the test appears to have passed for the original matrix A. Now check the
01831   //values of our second "view" of the matrix, 'Aview'.
01832 
01833   for(i=myFirstRow; i<=myLastRow; ++i) {
01834     EPETRA_TEST_ERR( Aview.ExtractGlobalBlockRowPointers(i, numMyRows,
01835                RowDim, numBlockEntries,
01836                BlockIndices[i-myFirstRow],
01837                Values), ierr);
01838 
01839     if (numMyRows != numBlockEntries) return(-1);
01840     if (RowDim != elemSize) return(-2);
01841     for(j=0; j<numBlockEntries; ++j) {
01842       if (Values[j]->A()[0] != 1.0) {
01843   cout << "Aview: Row " << i << " Values["<<j<<"][0]: "<< Values[j][0]
01844        << " should be 1.0" << endl;
01845   return(-3); //comment-out this return to de-activate this test
01846       }
01847     }
01848     
01849     delete [] BlockIndices[i-myFirstRow];
01850   }
01851 
01852 
01853   if (verbose&&localProc==0) {
01854     cout << "checkMergeRedundantEntries, A:" << endl;
01855   }
01856 
01857 
01858   delete [] BlockIndices;
01859   delete [] myCols;
01860 
01861   return(0);
01862 }
01863 
01864 int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose)
01865 {
01866   int numProcs = comm.NumProc();
01867   int localProc = comm.MyPID();
01868 
01869   int myFirstRow = localProc*3;
01870   int myLastRow = myFirstRow+2;
01871   int numMyRows = myLastRow - myFirstRow + 1;
01872   int numGlobalRows = numProcs*numMyRows;
01873   int i,j, ierr;
01874 
01875   int numCols = numMyRows;
01876   int* myCols = new int[numCols];
01877 
01878   int col = myFirstRow;
01879   for(i=0; i<numCols; ++i) {
01880     myCols[i] = col++;
01881     if (col > myLastRow) col = myFirstRow;
01882   }
01883 
01884   int elemSize = 2;
01885   int indexBase = 0;
01886 
01887   Epetra_BlockMap map(numGlobalRows, numMyRows,
01888           elemSize, indexBase, comm);
01889 
01890   Epetra_VbrMatrix A(Copy, map, numCols);
01891 
01892   double* coef = new double[elemSize*elemSize];
01893 
01894   for(i=myFirstRow; i<=myLastRow; ++i) {
01895     int myPointRow = i*elemSize;
01896 
01897     //The coefficients need to be laid out in column-major order. i.e., the
01898     //coefficients in a column occur contiguously.
01899     for(int ii=0; ii<elemSize; ++ii) {
01900       for(int jj=0; jj<elemSize; ++jj) {
01901   double val = (myPointRow+ii)*1.0;
01902   coef[ii+elemSize*jj] = val;
01903       }
01904     }
01905 
01906     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, numCols, myCols), ierr);
01907 
01908     for(j=0; j<numCols; ++j) {
01909       EPETRA_TEST_ERR( A.SubmitBlockEntry(coef, elemSize,
01910             elemSize, elemSize), ierr);
01911     }
01912 
01913     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
01914   }
01915 
01916   EPETRA_TEST_ERR( A.FillComplete(), ierr);
01917 
01918   delete [] coef;
01919   delete [] myCols;
01920 
01921   Epetra_SerialDenseMatrix** blockEntries;
01922   int len = elemSize*numCols, checkLen;
01923   double* values = new double[len];
01924   int* indices = new int[len];
01925   int RowDim, numBlockEntries;
01926 
01927   for(i=myFirstRow; i<=myLastRow; ++i) {
01928     EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(i, numMyRows,
01929                  RowDim, numBlockEntries,
01930                  indices,
01931                  blockEntries), ierr);
01932     if (numMyRows != numBlockEntries) return(-1);
01933     if (RowDim != elemSize) return(-2);
01934 
01935     int myPointRow = i*elemSize - myFirstRow*elemSize;
01936     int ii,jj;
01937     for(ii=0; ii<elemSize; ++ii) {
01938       EPETRA_TEST_ERR( A.ExtractMyRowCopy(myPointRow+ii, len,
01939             checkLen, values, indices), ierr);
01940       if (len != checkLen) return(-3);
01941 
01942       double val = (i*elemSize+ii)*1.0;
01943       double blockvalue = blockEntries[0]->A()[ii];
01944 
01945       for(jj=0; jj<len; ++jj) {
01946   if (values[jj] != val) return(-4);
01947   if (values[jj] != blockvalue) return(-5);
01948       }
01949     }
01950   }
01951 
01952   delete [] values;
01953   delete [] indices;
01954 
01955   return(0);
01956 }
01957 
01958 int checkMatvecSameVectors(Epetra_Comm& comm, bool verbose)
01959 {
01960   int numProcs = comm.NumProc();
01961   int localProc = comm.MyPID();
01962 
01963   int myFirstRow = localProc*3;
01964   int myLastRow = myFirstRow+2;
01965   int numMyRows = myLastRow - myFirstRow + 1;
01966   int numGlobalRows = numProcs*numMyRows;
01967   int i,ierr;
01968 
01969   int elemSize = 2;
01970   int num_off_diagonals = 1;
01971 
01972   epetra_test::matrix_data matdata(numGlobalRows, numGlobalRows,
01973            num_off_diagonals, elemSize);
01974 
01975   Epetra_BlockMap map(numGlobalRows, numMyRows, elemSize, 0, comm);
01976 
01977   Epetra_VbrMatrix A(Copy, map, num_off_diagonals*2+1);
01978 
01979   int* rowlengths = matdata.rowlengths();
01980   int** colindices = matdata.colindices();
01981 
01982   for(i=myFirstRow; i<=myLastRow; ++i) {
01983 
01984     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, rowlengths[i],
01985                                                colindices[i]), ierr);
01986 
01987     for(int j=0; j<rowlengths[i]; ++j) {
01988       EPETRA_TEST_ERR( A.SubmitBlockEntry(matdata.coefs(i,colindices[i][j]), elemSize,
01989                                           elemSize, elemSize), ierr);
01990     }
01991 
01992     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
01993   }
01994 
01995   EPETRA_TEST_ERR( A.FillComplete(), ierr);
01996 
01997   Epetra_Vector x(map), y(map);
01998 
01999   x.PutScalar(1.0);
02000 
02001   A.Multiply(false, x, y);
02002   A.Multiply(false, x, x);
02003 
02004   double* xptr = x.Values();
02005   double* yptr = y.Values();
02006 
02007   for(i=0; i<numMyRows; ++i) {
02008     if (xptr[i] != yptr[i]) {
02009       return(-1);
02010     }
02011   }
02012 
02013   return(0);
02014 }
02015 
02016 int checkEarlyDelete(Epetra_Comm& comm, bool verbose)
02017 {
02018   int localProc = comm.MyPID();
02019   int numProcs = comm.NumProc();
02020   int myFirstRow = localProc*3;
02021   int myLastRow = myFirstRow+2;
02022   int numMyRows = myLastRow - myFirstRow + 1;
02023   int numGlobalRows = numProcs*numMyRows;
02024   int i,ierr;
02025 
02026   int elemSize = 2;
02027   int num_off_diagonals = 1;
02028 
02029   epetra_test::matrix_data matdata(numGlobalRows, numGlobalRows,
02030                                    num_off_diagonals, elemSize);
02031 
02032   Epetra_BlockMap map(numGlobalRows, numMyRows, elemSize, 0, comm);
02033 
02034   Epetra_VbrMatrix* A = new Epetra_VbrMatrix(Copy, map, num_off_diagonals*2+1);
02035 
02036   int* rowlengths = matdata.rowlengths();
02037   int** colindices = matdata.colindices();
02038 
02039   for(i=myFirstRow; i<=myLastRow; ++i) {
02040 
02041     EPETRA_TEST_ERR( A->BeginInsertGlobalValues(i, rowlengths[i],
02042                                                colindices[i]), ierr);
02043 
02044     for(int j=0; j<rowlengths[i]; ++j) {
02045       EPETRA_TEST_ERR( A->SubmitBlockEntry(matdata.coefs(i,colindices[i][j]),
02046                                            elemSize, elemSize, elemSize), ierr);
02047     }
02048 
02049     EPETRA_TEST_ERR( A->EndSubmitEntries(), ierr);
02050   }
02051 
02052   //A call to BeginReplaceMyValues should produce an error at this
02053   //point, since IndicesAreLocal should be false.
02054   int errcode = A->BeginReplaceMyValues(0, 0, 0);
02055   if (errcode == 0) EPETRA_TEST_ERR(-1, ierr);
02056 
02057   EPETRA_TEST_ERR( A->FillComplete(), ierr);
02058 
02059   Epetra_VbrMatrix B(Copy, A->Graph());
02060 
02061   delete A;
02062 
02063   for(i=myFirstRow; i<=myLastRow; ++i) {
02064 
02065     EPETRA_TEST_ERR( B.BeginReplaceGlobalValues(i, rowlengths[i],
02066                                                colindices[i]), ierr);
02067 
02068     for(int j=0; j<rowlengths[i]; ++j) {
02069       EPETRA_TEST_ERR( B.SubmitBlockEntry(matdata.coefs(i,colindices[i][j]),
02070                                            elemSize, elemSize, elemSize), ierr);
02071     }
02072 
02073     EPETRA_TEST_ERR( B.EndSubmitEntries(), ierr);
02074   }
02075 
02076   EPETRA_TEST_ERR( B.FillComplete(), ierr);
02077 
02078   return(0);
02079 }
02080 
02081 int checkVbrRowMatrix(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {  
02082 
02083   if (verbose) cout << "Checking VbrRowMatrix Adapter..." << endl;
02084   int ierr = 0;
02085   EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
02086   EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
02087   EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
02088   EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
02089   // EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
02090   // EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
02091   EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
02092   EPETRA_TEST_ERR(!A.NumGlobalCols()==B.NumGlobalCols(),ierr);
02093   EPETRA_TEST_ERR(!A.NumGlobalDiagonals()==B.NumGlobalDiagonals(),ierr);
02094   EPETRA_TEST_ERR(!A.NumGlobalNonzeros()==B.NumGlobalNonzeros(),ierr);
02095   EPETRA_TEST_ERR(!A.NumGlobalRows()==B.NumGlobalRows(),ierr);
02096   EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
02097   EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
02098   EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
02099   for (int i=0; i<A.NumMyRows(); i++) {
02100     int nA, nB;
02101     A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
02102     EPETRA_TEST_ERR(!nA==nB,ierr);
02103   }
02104   EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
02105   EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
02106   EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
02107   EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
02108   EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
02109   // EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
02110   EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);
02111 
02112   int NumVectors = 5;
02113   { // No transpose case
02114     Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
02115     Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
02116     Epetra_MultiVector YA2(YA1);
02117     Epetra_MultiVector YB1(YA1);
02118     Epetra_MultiVector YB2(YA1);
02119     X.Random();
02120     
02121     bool transA = false;
02122     A.SetUseTranspose(transA);
02123     B.SetUseTranspose(transA);
02124     A.Apply(X,YA1);
02125     A.Multiply(transA, X, YA2);
02126     EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
02127     B.Apply(X,YB1);
02128     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
02129     B.Multiply(transA, X, YB2);
02130     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);
02131     
02132   }
02133   {// transpose case
02134     Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
02135     Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
02136     Epetra_MultiVector YA2(YA1);
02137     Epetra_MultiVector YB1(YA1);
02138     Epetra_MultiVector YB2(YA1);
02139     X.Random();
02140     
02141     bool transA = true;
02142     A.SetUseTranspose(transA);
02143     B.SetUseTranspose(transA);
02144     A.Apply(X,YA1);
02145     A.Multiply(transA, X, YA2);
02146     EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
02147     B.Apply(X,YB1);
02148     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
02149     B.Multiply(transA, X,YB2);
02150     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);
02151     
02152   }
02153 
02154   /*
02155     Epetra_Vector diagA(A.RowMatrixRowMap());
02156     EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
02157     Epetra_Vector diagB(B.RowMatrixRowMap());
02158     EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
02159     EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);
02160   */
02161   Epetra_Vector rowA(A.RowMatrixRowMap());
02162   EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
02163   Epetra_Vector rowB(B.RowMatrixRowMap());
02164   EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
02165   EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);
02166 
02167   Epetra_Vector colA(A.OperatorDomainMap());
02168   EPETRA_TEST_ERR(A.InvColSums(colA),ierr);  // -2 error code
02169   Epetra_Vector colB(B.OperatorDomainMap());
02170   EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
02171   EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr); // 1 error code
02172 
02173   EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
02174   EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);
02175 
02176   EPETRA_TEST_ERR(A.RightScale(colA),ierr);  // -3 error code
02177   EPETRA_TEST_ERR(B.RightScale(colB),ierr);  // -3 error code
02178 
02179 
02180   EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
02181   EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);
02182 
02183 
02184   EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf after scaling", verbose), ierr);
02185   EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne after scaling", verbose),ierr);
02186 
02187   vector<double> valuesA(A.MaxNumEntries());
02188   vector<int> indicesA(A.MaxNumEntries());  
02189   vector<double> valuesB(B.MaxNumEntries());
02190   vector<int> indicesB(B.MaxNumEntries());
02191   //return(0);
02192   for (int i=0; i<A.NumMyRows(); i++) {
02193     int nA, nB;
02194     EPETRA_TEST_ERR(A.ExtractMyRowCopy(i, A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr); 
02195     EPETRA_TEST_ERR(B.ExtractMyRowCopy(i, B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
02196     EPETRA_TEST_ERR(!nA==nB,ierr);
02197     for (int j=0; j<nA; j++) {
02198       double curVal = valuesA[j];
02199       int curIndex = indicesA[j];
02200       bool notfound = true;
02201       int jj = 0;
02202       while (notfound && jj< nB) {
02203   if (!checkValues(curVal, valuesB[jj])) notfound = false;
02204   jj++;
02205       }
02206       EPETRA_TEST_ERR(notfound, ierr);
02207       vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex);  // find curIndex in indicesB
02208       EPETRA_TEST_ERR(p==indicesB.end(), ierr);
02209     }
02210 
02211   }
02212 
02213   if (verbose) {
02214     if (ierr==0)
02215       cout << "RowMatrix Methods check OK" << endl;
02216     else
02217       cout << "ierr = " << ierr << ". RowMatrix Methods check failed." << endl;
02218   }
02219   return (ierr);
02220 }
02221 

Generated on Wed May 12 21:41:04 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7