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