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 2001 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, i, j, 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 (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 (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 (j=0; j < ColDim; j++)
00446   for (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 (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 (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 Indices[1] ; 
00593       Indices[0] = j ; 
00594       Epetra_SerialDenseMatrix * AD = &(BlockEntries[0][0]);
00595       EPETRA_TEST_ERR(!(A->BeginInsertGlobalValues( CurRow, 1, Indices)==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 (i=0; i<NumMyElements; i++) forierr += !(A->NumGlobalBlockEntries(MyGlobalElements[i])==3);
00656     else
00657       for (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 (i=0; i<NumMyElements; i++) forierr += !(A->NumMyBlockEntries(i)==3);
00664     else
00665       for (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 (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     int NumEntries;
00987     int NumMyNonzeros = 0, NumMyEquations = 0;
00988     
00989     for (i=0; i<NumMyElements; i++) {
00990       int CurRow = MyGlobalElements[i];
00991       int RowDim = ElementSizeList[i]-MinSize;
00992       NumMyEquations += BlockEntries[RowDim][RowDim].M();
00993       
00994       if (CurRow==0)
00995   {
00996     Indices[0] = CurRow;
00997     Indices[1] = CurRow+1;
00998     NumEntries = 2;
00999     ColDims[0] = ElementSizeList[i] - MinSize;
01000     ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
01001   }
01002       else if (CurRow == NumGlobalElements-1)
01003   {
01004     Indices[0] = CurRow-1;
01005     Indices[1] = CurRow;
01006     NumEntries = 2;
01007     ColDims[0] = ElementSizeList[i-1] - MinSize;
01008     ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
01009   }
01010       else {
01011   Indices[0] = CurRow-1;
01012   Indices[1] = CurRow;
01013   Indices[2] = CurRow+1;
01014   NumEntries = 3;
01015   if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
01016   else ColDims[0] = ElementSizeList[i-1] - MinSize;
01017   ColDims[1] = ElementSizeList[i] - MinSize;
01018   // ElementSize on MyPID+1
01019   if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
01020   else ColDims[2] = ElementSizeList[i+1] - MinSize;
01021       }
01022       EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
01023       EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
01024       EPETRA_TEST_ERR(!(A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
01025       forierr = 0;
01026       for (j=0; j < NumEntries; j++) {
01027   Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
01028   NumMyNonzeros += AD->M() * AD->N();   
01029   forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
01030       }
01031       EPETRA_TEST_ERR(forierr,ierr);
01032 
01033       A->EndSubmitEntries();
01034     }
01035     EPETRA_TEST_ERR( checkmultiply( false, *A, x, orig_check_y ), ierr );     
01036   }
01037 
01038   //
01039   //  Suminto should cause the matrix to be doubled 
01040   //
01041   if ( ! ExtraBlocks )   {
01042     // Add  rows one-at-a-time
01043 
01044     int NumEntries;
01045     int NumMyNonzeros = 0, NumMyEquations = 0;
01046     
01047     for (i=0; i<NumMyElements; i++) {
01048       int CurRow = MyGlobalElements[i];
01049       int RowDim = ElementSizeList[i]-MinSize;
01050       NumMyEquations += BlockEntries[RowDim][RowDim].M();
01051       
01052       if (CurRow==0)
01053   {
01054     Indices[0] = CurRow;
01055     Indices[1] = CurRow+1;
01056     NumEntries = 2;
01057     ColDims[0] = ElementSizeList[i] - MinSize;
01058     ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
01059   }
01060       else if (CurRow == NumGlobalElements-1)
01061   {
01062     Indices[0] = CurRow-1;
01063     Indices[1] = CurRow;
01064     NumEntries = 2;
01065     ColDims[0] = ElementSizeList[i-1] - MinSize;
01066     ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
01067   }
01068       else {
01069   Indices[0] = CurRow-1;
01070   Indices[1] = CurRow;
01071   Indices[2] = CurRow+1;
01072   NumEntries = 3;
01073   if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
01074   else ColDims[0] = ElementSizeList[i-1] - MinSize;
01075   ColDims[1] = ElementSizeList[i] - MinSize;
01076   // ElementSize on MyPID+1
01077   if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
01078   else ColDims[2] = ElementSizeList[i+1] - MinSize;
01079       }
01080       if ( insertlocal ) {
01081   for ( int ii=0; ii < NumEntries; ii++ ) 
01082     MyIndices[ii] = colmap->LID( Indices[ii] ) ; 
01083   EPETRA_TEST_ERR(!(A->BeginSumIntoMyValues( rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
01084       } else { 
01085   EPETRA_TEST_ERR(!(A->BeginSumIntoGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
01086       }
01087       forierr = 0;
01088       for (j=0; j < NumEntries; j++) {
01089   Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
01090   NumMyNonzeros += AD->M() * AD->N();   
01091   //  This has nothing to do with insertlocal, but that is a convenient bool to key on 
01092   if ( insertlocal ) {
01093     forierr += !(A->SubmitBlockEntry( *AD )==0);
01094   } else { 
01095     forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
01096   }
01097       }
01098       EPETRA_TEST_ERR(forierr,ierr);
01099 
01100       A->EndSubmitEntries();
01101     }
01102     
01103     orig_check_y.Scale(2.0) ;
01104 
01105     //  This will not work with FixedNumEntries unless we fix the fix the above loop to add the sorner elements to the tridi matrix
01106     if ( ! FixedNumEntries ) 
01107       EPETRA_TEST_ERR( checkmultiply( false, *A, x, orig_check_y ), ierr ); 
01108   }
01109 
01110 
01111   {for (int kr=0; kr<SizeRange; kr++) delete [] BlockEntries[kr];}
01112   delete [] BlockEntries;
01113   delete [] ColDims;
01114   delete [] MyIndices ;
01115   delete [] Indices;
01116   delete [] ElementSizeList;
01117 
01118 
01119 
01120 
01121   if (verbose) cout << "\n\n*****Insert methods should not be accepted" << endl<< endl;
01122 
01123   int One = 1;
01124   if (B->MyGRID(0)) EPETRA_TEST_ERR(!(B->BeginInsertGlobalValues(0, 1, &One)==-2),ierr);
01125 
01126   Epetra_Vector checkDiag(B->RowMap());
01127   forierr = 0;
01128   int NumMyEquations1 = B->NumMyRows();
01129   double two1 = 2.0;
01130 
01131   // Test diagonal replacement and extraction methods
01132 
01133   forierr = 0;
01134   for (i=0; i<NumMyEquations1; i++) checkDiag[i]=two1;
01135   EPETRA_TEST_ERR(forierr,ierr);
01136   
01137   EPETRA_TEST_ERR(!(B->ReplaceDiagonalValues(checkDiag)==0),ierr);
01138   
01139   Epetra_Vector checkDiag1(B->RowMap());
01140   EPETRA_TEST_ERR(!(B->ExtractDiagonalCopy(checkDiag1)==0),ierr);
01141   
01142   forierr = 0;
01143   for (i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
01144   EPETRA_TEST_ERR(forierr,ierr);
01145 
01146   if (verbose) cout << "\n\nDiagonal extraction and replacement OK.\n\n" << endl;
01147 
01148   double originfnorm = B->NormInf();
01149   double origonenorm = B->NormOne();
01150   EPETRA_TEST_ERR(!(B->Scale(4.0)==0),ierr);
01151   //  EPETRA_TEST_ERR((B->NormOne()!=origonenorm),ierr);
01152   EPETRA_TEST_ERR(!(B->NormInf()==4.0 * originfnorm),ierr);
01153   EPETRA_TEST_ERR(!(B->NormOne()==4.0 * origonenorm),ierr);
01154 
01155   if (verbose) cout << "\n\nMatrix scale OK.\n\n" << endl;
01156   
01157   // Testing VbrRowMatrix adapter
01158   Epetra_VbrRowMatrix ARowMatrix(A);
01159 
01160   EPETRA_TEST_ERR(checkVbrRowMatrix(*A, ARowMatrix, verbose),ierr);
01161 
01162   if ( PreviousA ) 
01163     delete *PreviousA; 
01164   
01165   //
01166   //  The following code deals with the fact that A has to be delete before graph is, when 
01167   //  A is built with a contructor that takes a graph as input.
01168   //
01169   delete B;
01170   if ( graph ) {
01171     delete A;
01172     delete graph ; 
01173     *PreviousA = 0 ; 
01174   } else { 
01175     *PreviousA = A ; 
01176   }
01177   
01178   delete CrsA;
01179   delete CrsB;
01180   delete OrigCrsA ;
01181 
01182   return ierr; 
01183 
01184 }
01185 
01186 
01187 int main(int argc, char *argv[])
01188 {
01189   int ierr = 0;
01190   bool debug = false;
01191 
01192 #ifdef EPETRA_MPI
01193   MPI_Init(&argc,&argv);
01194   Epetra_MpiComm Comm( MPI_COMM_WORLD );
01195 #else
01196   Epetra_SerialComm Comm;
01197 #endif
01198 
01199   int MyPID = Comm.MyPID();
01200   int NumProc = Comm.NumProc();
01201   bool verbose = false;
01202 
01203   // Check if we should print results to standard out
01204   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
01205 
01206   // char tmp;
01207   //  int rank = Comm.MyPID();
01208   //  if (rank==0) cout << "Press any key to continue..."<< endl;
01209   //  if (rank==0) cin >> tmp;
01210   //  Comm.Barrier();
01211 
01212   //  if ( ! verbose )  
01213   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
01214 
01215   if (verbose && MyPID==0)
01216     cout << Epetra_Version() << endl << endl;
01217 
01218   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
01219               << " is alive."<<endl;
01220 
01221   // Redefine verbose to only print on PE 0
01222   if (verbose && Comm.MyPID()!=0) verbose = false;
01223 
01224 
01225   int NumMyElements = 400;
01226   //  int NumMyElements = 3; 
01227   int MinSize = 2;
01228   int MaxSize = 8;
01229   bool NoExtraBlocks = false; 
01230   bool symmetric = true; 
01231   bool NonSymmetric = false;
01232   bool NoInsertLocal = false ; 
01233   bool InsertLocal = true ; 
01234 
01235   Epetra_VbrMatrix* PreviousA = 0 ; 
01236 
01237   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 1, 1, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );  
01238 
01239   //
01240   //  Check the various constructors
01241   //  
01242   
01243   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
01244 
01245   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
01246 
01247   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );   
01248 
01249   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA ); 
01250 
01251   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_VEPR, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
01252 
01253   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01254 
01255   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_VEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
01256 
01257   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01258 
01259   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, WithGraph, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01260   assert ( PreviousA == 0 ) ; 
01261 
01262 
01263   //
01264   //  Check some various options
01265   //
01266   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01267 
01268   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01269 
01270   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01271   assert ( PreviousA == 0 ) ; 
01272 
01273   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01274 
01275   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01276 
01277   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA ); 
01278   assert ( PreviousA == 0 ) ; 
01279 
01280 
01281 
01282   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );   
01283 
01284   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, RowMapColMap_FEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );   
01285 
01286 
01287 
01288   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01289 
01290   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 3, 3, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01291 
01292   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01293 
01294   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 5, 5, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01295 
01296   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 6, 6, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01297 
01298   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 7, 7, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01299 
01300   ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 8, 8, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01301 
01302   //  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA ); 
01303 
01304 
01305   delete PreviousA;
01306 
01307   /*
01308   if (verbose) {
01309     // Test ostream << operator (if verbose)
01310     // Construct a Map that puts 2 equations on each PE
01311     
01312     int NumMyElements1 = 2;
01313     int NumMyElements1 = NumMyElements1;
01314     int NumGlobalElements1 = NumMyElements1*NumProc;
01315 
01316     Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
01317     
01318     // Get update list and number of local equations from newly created Map
01319     int * MyGlobalElements1 = new int[Map1.NumMyElements()];
01320     Map1.MyGlobalElements(MyGlobalElements1);
01321     
01322     // Create an integer vector NumNz that is used to build the Petra Matrix.
01323     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
01324     
01325     int * NumNz1 = new int[NumMyElements1];
01326     
01327     // We are building a tridiagonal matrix where each row has (-1 2 -1)
01328     // So we need 2 off-diagonal terms (except for the first and last equation)
01329     
01330     for (i=0; i<NumMyElements1; i++)
01331       if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalElements1-1)
01332   NumNz1[i] = 1;
01333       else
01334   NumNz1[i] = 2;
01335     
01336     // Create a Epetra_Matrix
01337     
01338     Epetra_VbrMatrix A1(Copy, Map1, NumNz1);
01339     
01340     // Add  rows one-at-a-time
01341     // Need some vectors to help
01342     // Off diagonal Values will always be -1
01343     
01344     
01345     int *Indices1 = new int[2];
01346     double two1 = 2.0;
01347     int NumEntries1;
01348 
01349     forierr = 0;
01350     for (i=0; i<NumMyElements1; i++)
01351       {
01352   if (MyGlobalElements1[i]==0)
01353     {
01354       Indices1[0] = 1;
01355       NumEntries1 = 1;
01356     }
01357   else if (MyGlobalElements1[i] == NumGlobalElements1-1)
01358     {
01359       Indices1[0] = NumGlobalElements1-2;
01360       NumEntries1 = 1;
01361     }
01362   else
01363     {
01364       Indices1[0] = MyGlobalElements1[i]-1;
01365       Indices1[1] = MyGlobalElements1[i]+1;
01366       NumEntries1 = 2;
01367     }
01368         forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1)==0);
01369   forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i)>0); // Put in the diagonal entry
01370       }
01371     EPETRA_TEST_ERR(forierr,ierr);
01372     // Finish up
01373     EPETRA_TEST_ERR(!(A1.FillComplete()==0),ierr);
01374     
01375     if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
01376     cout << A1 << endl;
01377     
01378     // Release all objects
01379     delete [] NumNz1;
01380     delete [] Values1;
01381     delete [] Indices1;
01382     delete [] MyGlobalElements1;
01383 
01384   }
01385   */
01386 
01387   EPETRA_TEST_ERR( checkMergeRedundantEntries(Comm, verbose), ierr);
01388 
01389   EPETRA_TEST_ERR( checkExtractMyRowCopy(Comm, verbose), ierr);
01390 
01391   EPETRA_TEST_ERR( checkMatvecSameVectors(Comm, verbose), ierr);
01392 
01393   EPETRA_TEST_ERR( checkEarlyDelete(Comm, verbose), ierr);
01394 
01395   if (verbose) {
01396     if (ierr==0) cout << "All VbrMatrix tests OK" << endl;
01397     else cout << ierr << " errors encountered." << endl;
01398   }
01399 
01400 #ifdef EPETRA_MPI
01401   MPI_Finalize() ;
01402 #endif
01403 
01404 /* end main
01405 */
01406 return ierr ;
01407 }
01408 
01409 int power_method(bool TransA, Epetra_VbrMatrix& A, 
01410      Epetra_MultiVector& q,
01411      Epetra_MultiVector& z, 
01412      Epetra_MultiVector& resid, 
01413      double * lambda, int niters, double tolerance,
01414      bool verbose) {  
01415 
01416   // variable needed for iteration
01417   double normz, residual;
01418 
01419   int ierr = 1;
01420 
01421   for (int iter = 0; iter < niters; iter++)
01422     {
01423       z.Norm2(&normz); // Compute 2-norm of z
01424       q.Scale(1.0/normz, z);
01425       A.Multiply(TransA, q, z); // Compute z = A*q
01426       q.Dot(z, lambda); // Approximate maximum eigenvaluE
01427       if (iter%100==0 || iter+1==niters)
01428   {
01429     resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
01430     resid.Norm2(&residual);
01431     if (verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda 
01432            << "  Residual of A*q - lambda*q = " << residual << endl;
01433   } 
01434       if (residual < tolerance) {
01435   ierr = 0;
01436   break;
01437       }
01438     }
01439   return(ierr);
01440 }
01441 int check(Epetra_VbrMatrix& A, 
01442     int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, 
01443     int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1, 
01444     int * MyGlobalElements, bool verbose) {
01445 
01446   int ierr = 0, forierr = 0;
01447   // Test query functions
01448 
01449   int NumMyRows = A.NumMyRows();
01450   if (verbose) cout << "\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
01451   // TEMP
01452   //if (verbose) cout << "\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
01453   if (verbose) cout << "\n\nNumber of local Rows is = " << NumMyRows << endl<< endl;
01454   if (verbose) cout << "\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
01455 
01456   EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);
01457 
01458   int NumMyNonzeros = A.NumMyNonzeros();
01459   if (verbose) cout << "\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
01460   //if (verbose) cout << "                            Should  = " << NumMyNonzeros1 << endl<< endl;
01461 
01462 
01463   if ( NumMyNonzeros != NumMyNonzeros1 ) {
01464     cout << " MyPID = " << A.Comm().MyPID() 
01465    << " NumMyNonzeros = " << NumMyNonzeros 
01466    << " NumMyNonzeros1 = " << NumMyNonzeros1 << endl ; 
01467   }
01468 
01469   EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);
01470 
01471   int NumGlobalRows = A.NumGlobalRows();
01472   if (verbose) cout << "\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
01473 
01474   EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);
01475 
01476   int NumGlobalNonzeros = A.NumGlobalNonzeros();
01477   if (verbose) cout << "\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
01478 
01479   if ( NumGlobalNonzeros != NumGlobalNonzeros1 ) {
01480     cout << " MyPID = " << A.Comm().MyPID() 
01481    << " NumGlobalNonzeros = " << NumGlobalNonzeros 
01482    << " NumGlobalNonzeros1 = " << NumGlobalNonzeros1 << endl ; 
01483   }
01484   EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);
01485 
01486   int NumMyBlockRows = A.NumMyBlockRows();
01487   if (verbose) cout << "\n\nNumber of local Block Rows = " << NumMyBlockRows << endl<< endl;
01488 
01489   EPETRA_TEST_ERR(!(NumMyBlockRows==NumMyBlockRows1),ierr);
01490 
01491   int NumMyBlockNonzeros = A.NumMyBlockEntries();
01492 
01493   if (verbose) cout << "\n\nNumber of local Nonzero Block entries = " << NumMyBlockNonzeros << endl<< endl;
01494   if (verbose) cout << "\n\nNumber of local Nonzero Block entries 1 = " << NumMyBlockNonzeros1 << endl<< endl;
01495 
01496   EPETRA_TEST_ERR(!(NumMyBlockNonzeros==NumMyBlockNonzeros1),ierr);
01497 
01498   int NumGlobalBlockRows = A.NumGlobalBlockRows();
01499   if (verbose) cout << "\n\nNumber of global Block Rows = " << NumGlobalBlockRows << endl<< endl;
01500 
01501   EPETRA_TEST_ERR(!(NumGlobalBlockRows==NumGlobalBlockRows1),ierr);
01502 
01503   int NumGlobalBlockNonzeros = A.NumGlobalBlockEntries();
01504   if (verbose) cout << "\n\nNumber of global Nonzero Block entries = " << NumGlobalBlockNonzeros << endl<< endl;
01505 
01506   EPETRA_TEST_ERR(!(NumGlobalBlockNonzeros==NumGlobalBlockNonzeros1),ierr);
01507 
01508   
01509   // Test RowMatrix interface implementations
01510   int RowDim, NumBlockEntries, * BlockIndices;
01511   Epetra_SerialDenseMatrix ** Values;
01512   // Get View of last block row
01513   A.ExtractMyBlockRowView(NumMyBlockRows-1, RowDim, NumBlockEntries,
01514         BlockIndices, Values);
01515   int NumMyEntries1 = 0;
01516   {for (int i=0; i < NumBlockEntries; i++) NumMyEntries1 += Values[i]->N();}
01517   int NumMyEntries;
01518   A.NumMyRowEntries(NumMyRows-1, NumMyEntries);
01519   if (verbose) {
01520     cout << "\n\nNumber of nonzero values in last row = "
01521    << NumMyEntries << endl<< endl;
01522   }
01523 
01524   EPETRA_TEST_ERR(!(NumMyEntries==NumMyEntries1),ierr);
01525   
01526   // Other binary tests
01527 
01528   EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
01529   EPETRA_TEST_ERR(!(A.Filled()),ierr);
01530   EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID())),ierr);
01531   EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID())),ierr);
01532   EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID()),ierr);
01533   EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID()),ierr);
01534   EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
01535   EPETRA_TEST_ERR(!(A.MyLRID(NumMyBlockRows-1)),ierr);
01536   EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
01537   EPETRA_TEST_ERR(A.MyLRID(NumMyBlockRows),ierr);
01538 
01539     
01540   int i, j;
01541   int MaxNumBlockEntries = A.MaxNumBlockEntries();
01542 
01543   // Pointer Extraction approach
01544 
01545   //   Local index
01546   int MyPointersRowDim, MyPointersNumBlockEntries;
01547   int * MyPointersBlockIndices = new int[MaxNumBlockEntries];
01548   Epetra_SerialDenseMatrix **MyPointersValuesPointers;
01549   //   Global Index
01550   int GlobalPointersRowDim, GlobalPointersNumBlockEntries;
01551   int * GlobalPointersBlockIndices = new int[MaxNumBlockEntries];
01552   Epetra_SerialDenseMatrix **GlobalPointersValuesPointers;
01553 
01554   // Copy Extraction approach
01555 
01556   //   Local index
01557   int MyCopyRowDim, MyCopyNumBlockEntries;
01558   int * MyCopyBlockIndices = new int[MaxNumBlockEntries];
01559   int * MyCopyColDims = new int[MaxNumBlockEntries];
01560   int * MyCopyLDAs = new int[MaxNumBlockEntries];
01561   int MaxRowDim = A.MaxRowDim();
01562   int MaxColDim = A.MaxColDim();
01563   int MyCopySizeOfValues = MaxRowDim*MaxColDim;
01564   double ** MyCopyValuesPointers = new double*[MaxNumBlockEntries];
01565   for (i=0; i<MaxNumBlockEntries; i++)
01566     MyCopyValuesPointers[i] = new double[MaxRowDim*MaxColDim];
01567 
01568   //   Global Index
01569   int GlobalCopyRowDim, GlobalCopyNumBlockEntries;
01570   int * GlobalCopyBlockIndices = new int[MaxNumBlockEntries];
01571   int * GlobalCopyColDims = new int[MaxNumBlockEntries];
01572   int * GlobalCopyLDAs = new int[MaxNumBlockEntries];
01573   
01574   int GlobalMaxRowDim = A.GlobalMaxRowDim();
01575   int GlobalMaxColDim = A.GlobalMaxColDim();
01576   int GlobalCopySizeOfValues = GlobalMaxRowDim*GlobalMaxColDim;
01577   double ** GlobalCopyValuesPointers = new double*[MaxNumBlockEntries];
01578   for (i=0; i<MaxNumBlockEntries; i++)
01579     GlobalCopyValuesPointers[i] = new double[GlobalMaxRowDim*GlobalMaxColDim];
01580 
01581   // View Extraction approaches
01582 
01583   //   Local index (There is no global view available)
01584   int MyView1RowDim, MyView1NumBlockEntries;
01585   int * MyView1BlockIndices;
01586   Epetra_SerialDenseMatrix **MyView1ValuesPointers = new Epetra_SerialDenseMatrix*[MaxNumBlockEntries];
01587 
01588   //   Local index version 2 (There is no global view available)
01589   int MyView2RowDim, MyView2NumBlockEntries;
01590   int * MyView2BlockIndices;
01591   Epetra_SerialDenseMatrix **MyView2ValuesPointers;
01592 
01593 
01594   // For each row, test six approaches to extracting data from a given local index matrix
01595   forierr = 0;
01596   for (i=0; i<NumMyBlockRows; i++) {
01597     int MyRow = i;
01598     int GlobalRow = A.GRID(i);
01599     // Get a copy of block indices in local index space, pointers to everything else
01600     EPETRA_TEST_ERR( A.ExtractMyBlockRowPointers(MyRow, MaxNumBlockEntries, MyPointersRowDim, 
01601         MyPointersNumBlockEntries, MyPointersBlockIndices,
01602         MyPointersValuesPointers), ierr );
01603     // Get a copy of block indices in local index space, pointers to everything else
01604     EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(GlobalRow, MaxNumBlockEntries, GlobalPointersRowDim, 
01605             GlobalPointersNumBlockEntries, GlobalPointersBlockIndices,
01606             GlobalPointersValuesPointers), ierr ) ;
01607 
01608     // Initiate a copy of block row in local index space.
01609     EPETRA_TEST_ERR( A.BeginExtractMyBlockRowCopy(MyRow, MaxNumBlockEntries, MyCopyRowDim, 
01610          MyCopyNumBlockEntries, MyCopyBlockIndices,
01611          MyCopyColDims), ierr);
01612     // Copy Values
01613     for (j=0; j<MyCopyNumBlockEntries; j++) {
01614       EPETRA_TEST_ERR( A.ExtractEntryCopy(MyCopySizeOfValues, MyCopyValuesPointers[j], MaxRowDim, false), ierr) ;
01615       MyCopyLDAs[j] = MaxRowDim;
01616     }
01617 
01618     // Initiate a copy of block row in global index space.
01619     EPETRA_TEST_ERR( A.BeginExtractGlobalBlockRowCopy(GlobalRow, MaxNumBlockEntries, GlobalCopyRowDim, 
01620             GlobalCopyNumBlockEntries, GlobalCopyBlockIndices,
01621             GlobalCopyColDims), ierr ) ;
01622     // Copy Values
01623     for (j=0; j<GlobalCopyNumBlockEntries; j++) {
01624       EPETRA_TEST_ERR( A.ExtractEntryCopy(GlobalCopySizeOfValues, GlobalCopyValuesPointers[j], GlobalMaxRowDim, false), ierr );
01625       GlobalCopyLDAs[j] = GlobalMaxRowDim;
01626     }
01627 
01628     // Initiate a view of block row in local index space (Version 1)
01629     EPETRA_TEST_ERR( A.BeginExtractMyBlockRowView(MyRow, MyView1RowDim, 
01630          MyView1NumBlockEntries, MyView1BlockIndices), ierr ) ;
01631     // Set pointers to values
01632     for (j=0; j<MyView1NumBlockEntries; j++) 
01633       EPETRA_TEST_ERR ( A.ExtractEntryView(MyView1ValuesPointers[j]), ierr ) ;
01634 
01635 
01636     // Extract a view of block row in local index space (version 2)
01637     EPETRA_TEST_ERR( A.ExtractMyBlockRowView(MyRow, MyView2RowDim, 
01638           MyView2NumBlockEntries, MyView2BlockIndices,
01639           MyView2ValuesPointers), ierr );
01640 
01641     forierr += !(MyPointersNumBlockEntries==GlobalPointersNumBlockEntries);
01642     forierr += !(MyPointersNumBlockEntries==MyCopyNumBlockEntries);
01643     forierr += !(MyPointersNumBlockEntries==GlobalCopyNumBlockEntries);
01644     forierr += !(MyPointersNumBlockEntries==MyView1NumBlockEntries);
01645     forierr += !(MyPointersNumBlockEntries==MyView2NumBlockEntries);
01646     for (j=1; j<MyPointersNumBlockEntries; j++) {
01647       forierr += !(MyCopyBlockIndices[j-1]<MyCopyBlockIndices[j]);
01648       forierr += !(MyView1BlockIndices[j-1]<MyView1BlockIndices[j]);
01649       forierr += !(MyView2BlockIndices[j-1]<MyView2BlockIndices[j]);
01650 
01651       forierr += !(GlobalPointersBlockIndices[j]==A.GCID(MyPointersBlockIndices[j]));
01652       forierr += !(A.LCID(GlobalPointersBlockIndices[j])==MyPointersBlockIndices[j]);
01653       forierr += !(GlobalPointersBlockIndices[j]==GlobalCopyBlockIndices[j]);
01654       
01655       Epetra_SerialDenseMatrix* my = MyPointersValuesPointers[j];
01656       Epetra_SerialDenseMatrix* global = GlobalPointersValuesPointers[j];
01657 
01658       Epetra_SerialDenseMatrix* myview1 = MyView1ValuesPointers[j];
01659       Epetra_SerialDenseMatrix* myview2 = MyView2ValuesPointers[j];
01660 
01661       forierr += !(CompareValues(my->A(), my->LDA(), 
01662          MyPointersRowDim, my->N(), 
01663          global->A(), global->LDA(), 
01664          GlobalPointersRowDim, global->N())==0);
01665       forierr += !(CompareValues(my->A(), my->LDA(), 
01666          MyPointersRowDim, my->N(), 
01667          MyCopyValuesPointers[j], MyCopyLDAs[j], 
01668          MyCopyRowDim, MyCopyColDims[j])==0);
01669       forierr += !(CompareValues(my->A(), my->LDA(), 
01670          MyPointersRowDim, my->N(), 
01671          GlobalCopyValuesPointers[j], GlobalCopyLDAs[j], 
01672          GlobalCopyRowDim, GlobalCopyColDims[j])==0);
01673       forierr += !(CompareValues(my->A(), my->LDA(), 
01674          MyPointersRowDim, my->N(), 
01675          myview1->A(), myview1->LDA(), 
01676          MyView1RowDim, myview1->N())==0);
01677       forierr += !(CompareValues(my->A(), my->LDA(),
01678          MyPointersRowDim, my->N(),
01679          myview2->A(), myview2->LDA(),
01680          MyView2RowDim, myview2->N())==0);
01681     }
01682   }
01683   EPETRA_TEST_ERR(forierr,ierr);
01684 
01685   // GlobalRowView should be illegal (since we have local indices)
01686   EPETRA_TEST_ERR(!(A.BeginExtractGlobalBlockRowView(A.GRID(0), MyView1RowDim, 
01687                  MyView1NumBlockEntries,
01688                  MyView1BlockIndices)==-2),ierr);
01689   
01690   // Extract a view of block row in local index space (version 2)
01691   EPETRA_TEST_ERR(!(A.ExtractGlobalBlockRowView(A.GRID(0), MyView2RowDim, 
01692              MyView2NumBlockEntries, MyView2BlockIndices,
01693              MyView2ValuesPointers)==-2),ierr);
01694   
01695   delete [] MyPointersBlockIndices;
01696   delete [] GlobalPointersBlockIndices;
01697   delete [] MyCopyBlockIndices;
01698   delete [] MyCopyColDims;
01699   delete [] MyCopyLDAs;
01700   for (i=0; i<MaxNumBlockEntries; i++) delete [] MyCopyValuesPointers[i];
01701   delete [] MyCopyValuesPointers;
01702   delete [] GlobalCopyBlockIndices;
01703   delete [] GlobalCopyColDims;
01704   delete [] GlobalCopyLDAs;
01705   for (i=0; i<MaxNumBlockEntries; i++) delete [] GlobalCopyValuesPointers[i];
01706   delete [] GlobalCopyValuesPointers;
01707   delete [] MyView1ValuesPointers;
01708   if (verbose) cout << "\n\nRows sorted check OK" << endl<< endl;
01709   
01710   return ierr;
01711 }
01712 
01713 //=============================================================================
01714 int CompareValues(double * A, int LDA, int NumRowsA, int NumColsA, 
01715       double * B, int LDB, int NumRowsB, int NumColsB) {
01716   
01717   int i, j, ierr = 0, forierr = 0;
01718   double * ptr1 = B;
01719   double * ptr2;
01720   
01721   if (NumRowsA!=NumRowsB) EPETRA_TEST_ERR(-2,ierr);
01722   if (NumColsA!=NumColsB) EPETRA_TEST_ERR(-3,ierr);
01723  
01724 
01725   forierr = 0;
01726   for (j=0; j<NumColsA; j++) {
01727     ptr1 = B + j*LDB;
01728     ptr2 = A + j*LDA;
01729     for (i=0; i<NumRowsA; i++) forierr += (*ptr1++ != *ptr2++);
01730   }
01731   EPETRA_TEST_ERR(forierr,ierr);
01732   return ierr;
01733 }
01734 
01735 int checkMergeRedundantEntries(Epetra_Comm& comm, bool verbose)
01736 {
01737   int numProcs = comm.NumProc();
01738   int localProc = comm.MyPID();
01739 
01740   int myFirstRow = localProc*3;
01741   int myLastRow = myFirstRow+2;
01742   int numMyRows = myLastRow - myFirstRow + 1;
01743   int numGlobalRows = numProcs*numMyRows;
01744   int i,j, ierr;
01745 
01746   //We'll set up a matrix which is globally block-diagonal, i.e., on each
01747   //processor the list of columns == list of rows.
01748   //Set up a list of column-indices which is twice as long as it should be,
01749   //and its contents will be the list of local rows, repeated twice.
01750   int numCols = 2*numMyRows;
01751   int* myCols = new int[numCols];
01752 
01753   int col = myFirstRow;
01754   for(i=0; i<numCols; ++i) {
01755     myCols[i] = col++;
01756     if (col > myLastRow) col = myFirstRow;
01757   }
01758 
01759   int elemSize = 2;
01760   int indexBase = 0;
01761 
01762   Epetra_BlockMap map(numGlobalRows, numMyRows,
01763           elemSize, indexBase, comm);
01764 
01765   Epetra_VbrMatrix A(Copy, map, numCols);
01766 
01767   Epetra_MultiVector x(map, 1), y(map, 1);
01768   x.PutScalar(1.0);
01769 
01770   Epetra_MultiVector x3(map, 3), y3(map, 3);
01771   x.PutScalar(1.0);
01772 
01773   double* coef = new double[elemSize*elemSize];
01774   for(i=0; i<elemSize*elemSize; ++i) {
01775     coef[i] = 0.5;
01776   }
01777 
01778   //we're going to insert each row twice, with coef values of 0.5. So after
01779   //FillComplete, which internally calls MergeRedundantEntries, the
01780   //matrix should contain 1.0 in each entry.
01781 
01782   for(i=myFirstRow; i<=myLastRow; ++i) {
01783     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, numCols, myCols), ierr);
01784 
01785     for(j=0; j<numCols; ++j) {
01786       EPETRA_TEST_ERR( A.SubmitBlockEntry(coef, elemSize,
01787             elemSize, elemSize), ierr);
01788     }
01789 
01790     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
01791   }
01792 
01793   EPETRA_TEST_ERR( A.FillComplete(), ierr);
01794 
01795   delete [] coef;
01796 
01797   if (verbose) cout << "Multiply x"<<endl;
01798   EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr );
01799 
01800 
01801   //Next we're going to extract pointers-to-block-rows and check values to make
01802   //sure that the internal method Epetra_VbrMatrix::mergeRedundantEntries()
01803   //worked correctly. 
01804   //At the same time, we're also going to create another VbrMatrix which will
01805   //be a View of the matrix we've already created. This will serve to provide
01806   //more test coverage of the VbrMatrix code.
01807 
01808   int numBlockEntries = 0;
01809   int RowDim;
01810   int** BlockIndices = new int*[numMyRows];
01811   Epetra_SerialDenseMatrix** Values;
01812   Epetra_VbrMatrix Aview(View, map, numMyRows);
01813 
01814   for(i=myFirstRow; i<=myLastRow; ++i) {
01815     BlockIndices[i-myFirstRow] = new int[numCols];
01816     EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(i, numCols,
01817                  RowDim, numBlockEntries,
01818                  BlockIndices[i-myFirstRow],
01819                  Values), ierr);
01820 
01821     EPETRA_TEST_ERR( Aview.BeginInsertGlobalValues(i, numBlockEntries,
01822                 BlockIndices[i-myFirstRow]), ierr);
01823 
01824     if (numMyRows != numBlockEntries) return(-1);
01825     if (RowDim != elemSize) return(-2);
01826     for(j=0; j<numBlockEntries; ++j) {
01827       if (Values[j]->A()[0] != 1.0) {
01828   cout << "Row " << i << " Values["<<j<<"][0]: "<< Values[j][0]
01829        << " should be 1.0" << endl;
01830   return(-3); //comment-out this return to de-activate this test
01831       }
01832 
01833       EPETRA_TEST_ERR( Aview.SubmitBlockEntry(Values[j]->A(),
01834                 Values[j]->LDA(),
01835                 Values[j]->M(),
01836                 Values[j]->N()), ierr);
01837     }
01838 
01839     EPETRA_TEST_ERR( Aview.EndSubmitEntries(), ierr);
01840   }
01841 
01842   EPETRA_TEST_ERR( Aview.FillComplete(), ierr);
01843 
01844   //So the test appears to have passed for the original matrix A. Now check the
01845   //values of our second "view" of the matrix, 'Aview'.
01846 
01847   for(i=myFirstRow; i<=myLastRow; ++i) {
01848     EPETRA_TEST_ERR( Aview.ExtractGlobalBlockRowPointers(i, numMyRows,
01849                RowDim, numBlockEntries,
01850                BlockIndices[i-myFirstRow],
01851                Values), ierr);
01852 
01853     if (numMyRows != numBlockEntries) return(-1);
01854     if (RowDim != elemSize) return(-2);
01855     for(j=0; j<numBlockEntries; ++j) {
01856       if (Values[j]->A()[0] != 1.0) {
01857   cout << "Aview: Row " << i << " Values["<<j<<"][0]: "<< Values[j][0]
01858        << " should be 1.0" << endl;
01859   return(-3); //comment-out this return to de-activate this test
01860       }
01861     }
01862     
01863     delete [] BlockIndices[i-myFirstRow];
01864   }
01865 
01866 
01867   if (verbose&&localProc==0) {
01868     cout << "checkMergeRedundantEntries, A:" << endl;
01869   }
01870 
01871 
01872   delete [] BlockIndices;
01873   delete [] myCols;
01874 
01875   return(0);
01876 }
01877 
01878 int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose)
01879 {
01880   int numProcs = comm.NumProc();
01881   int localProc = comm.MyPID();
01882 
01883   int myFirstRow = localProc*3;
01884   int myLastRow = myFirstRow+2;
01885   int numMyRows = myLastRow - myFirstRow + 1;
01886   int numGlobalRows = numProcs*numMyRows;
01887   int i,j, ierr;
01888 
01889   int numCols = numMyRows;
01890   int* myCols = new int[numCols];
01891 
01892   int col = myFirstRow;
01893   for(i=0; i<numCols; ++i) {
01894     myCols[i] = col++;
01895     if (col > myLastRow) col = myFirstRow;
01896   }
01897 
01898   int elemSize = 2;
01899   int indexBase = 0;
01900 
01901   Epetra_BlockMap map(numGlobalRows, numMyRows,
01902           elemSize, indexBase, comm);
01903 
01904   Epetra_VbrMatrix A(Copy, map, numCols);
01905 
01906   double* coef = new double[elemSize*elemSize];
01907 
01908   for(i=myFirstRow; i<=myLastRow; ++i) {
01909     int myPointRow = i*elemSize;
01910 
01911     //The coefficients need to be laid out in column-major order. i.e., the
01912     //coefficients in a column occur contiguously.
01913     for(int ii=0; ii<elemSize; ++ii) {
01914       for(int jj=0; jj<elemSize; ++jj) {
01915   double val = (myPointRow+ii)*1.0;
01916   coef[ii+elemSize*jj] = val;
01917       }
01918     }
01919 
01920     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, numCols, myCols), ierr);
01921 
01922     for(j=0; j<numCols; ++j) {
01923       EPETRA_TEST_ERR( A.SubmitBlockEntry(coef, elemSize,
01924             elemSize, elemSize), ierr);
01925     }
01926 
01927     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
01928   }
01929 
01930   EPETRA_TEST_ERR( A.FillComplete(), ierr);
01931 
01932   delete [] coef;
01933   delete [] myCols;
01934 
01935   Epetra_SerialDenseMatrix** blockEntries;
01936   int len = elemSize*numCols, checkLen;
01937   double* values = new double[len];
01938   int* indices = new int[len];
01939   int RowDim, numBlockEntries;
01940 
01941   for(i=myFirstRow; i<=myLastRow; ++i) {
01942     EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(i, numMyRows,
01943                  RowDim, numBlockEntries,
01944                  indices,
01945                  blockEntries), ierr);
01946     if (numMyRows != numBlockEntries) return(-1);
01947     if (RowDim != elemSize) return(-2);
01948 
01949     int myPointRow = i*elemSize - myFirstRow*elemSize;
01950     int ii,jj;
01951     for(ii=0; ii<elemSize; ++ii) {
01952       EPETRA_TEST_ERR( A.ExtractMyRowCopy(myPointRow+ii, len,
01953             checkLen, values, indices), ierr);
01954       if (len != checkLen) return(-3);
01955 
01956       double val = (i*elemSize+ii)*1.0;
01957       double blockvalue = blockEntries[0]->A()[ii];
01958 
01959       for(jj=0; jj<len; ++jj) {
01960   if (values[jj] != val) return(-4);
01961   if (values[jj] != blockvalue) return(-5);
01962       }
01963     }
01964   }
01965 
01966   delete [] values;
01967   delete [] indices;
01968 
01969   return(0);
01970 }
01971 
01972 int checkMatvecSameVectors(Epetra_Comm& comm, bool verbose)
01973 {
01974   int numProcs = comm.NumProc();
01975   int localProc = comm.MyPID();
01976 
01977   int myFirstRow = localProc*3;
01978   int myLastRow = myFirstRow+2;
01979   int numMyRows = myLastRow - myFirstRow + 1;
01980   int numGlobalRows = numProcs*numMyRows;
01981   int i,ierr;
01982 
01983   int elemSize = 2;
01984   int num_off_diagonals = 1;
01985 
01986   epetra_test::matrix_data matdata(numGlobalRows, numGlobalRows,
01987            num_off_diagonals, elemSize);
01988 
01989   Epetra_BlockMap map(numGlobalRows, numMyRows, elemSize, 0, comm);
01990 
01991   Epetra_VbrMatrix A(Copy, map, num_off_diagonals*2+1);
01992 
01993   int* rowlengths = matdata.rowlengths();
01994   int** colindices = matdata.colindices();
01995 
01996   for(i=myFirstRow; i<=myLastRow; ++i) {
01997 
01998     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, rowlengths[i],
01999                                                colindices[i]), ierr);
02000 
02001     for(int j=0; j<rowlengths[i]; ++j) {
02002       EPETRA_TEST_ERR( A.SubmitBlockEntry(matdata.coefs(i,colindices[i][j]), elemSize,
02003                                           elemSize, elemSize), ierr);
02004     }
02005 
02006     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
02007   }
02008 
02009   EPETRA_TEST_ERR( A.FillComplete(), ierr);
02010 
02011   Epetra_Vector x(map), y(map);
02012 
02013   x.PutScalar(1.0);
02014 
02015   A.Multiply(false, x, y);
02016   A.Multiply(false, x, x);
02017 
02018   double* xptr = x.Values();
02019   double* yptr = y.Values();
02020 
02021   for(i=0; i<numMyRows; ++i) {
02022     if (xptr[i] != yptr[i]) {
02023       return(-1);
02024     }
02025   }
02026 
02027   return(0);
02028 }
02029 
02030 int checkEarlyDelete(Epetra_Comm& comm, bool verbose)
02031 {
02032   int localProc = comm.MyPID();
02033   int numProcs = comm.NumProc();
02034   int myFirstRow = localProc*3;
02035   int myLastRow = myFirstRow+2;
02036   int numMyRows = myLastRow - myFirstRow + 1;
02037   int numGlobalRows = numProcs*numMyRows;
02038   int i,ierr;
02039 
02040   int elemSize = 2;
02041   int num_off_diagonals = 1;
02042 
02043   epetra_test::matrix_data matdata(numGlobalRows, numGlobalRows,
02044                                    num_off_diagonals, elemSize);
02045 
02046   Epetra_BlockMap map(numGlobalRows, numMyRows, elemSize, 0, comm);
02047 
02048   Epetra_VbrMatrix* A = new Epetra_VbrMatrix(Copy, map, num_off_diagonals*2+1);
02049 
02050   int* rowlengths = matdata.rowlengths();
02051   int** colindices = matdata.colindices();
02052 
02053   for(i=myFirstRow; i<=myLastRow; ++i) {
02054 
02055     EPETRA_TEST_ERR( A->BeginInsertGlobalValues(i, rowlengths[i],
02056                                                colindices[i]), ierr);
02057 
02058     for(int j=0; j<rowlengths[i]; ++j) {
02059       EPETRA_TEST_ERR( A->SubmitBlockEntry(matdata.coefs(i,colindices[i][j]),
02060                                            elemSize, elemSize, elemSize), ierr);
02061     }
02062 
02063     EPETRA_TEST_ERR( A->EndSubmitEntries(), ierr);
02064   }
02065 
02066   //A call to BeginReplaceMyValues should produce an error at this
02067   //point, since IndicesAreLocal should be false.
02068   int errcode = A->BeginReplaceMyValues(0, 0, 0);
02069   if (errcode == 0) EPETRA_TEST_ERR(-1, ierr);
02070 
02071   EPETRA_TEST_ERR( A->FillComplete(), ierr);
02072 
02073   Epetra_VbrMatrix B(Copy, A->Graph());
02074 
02075   delete A;
02076 
02077   for(i=myFirstRow; i<=myLastRow; ++i) {
02078 
02079     EPETRA_TEST_ERR( B.BeginReplaceGlobalValues(i, rowlengths[i],
02080                                                colindices[i]), ierr);
02081 
02082     for(int j=0; j<rowlengths[i]; ++j) {
02083       EPETRA_TEST_ERR( B.SubmitBlockEntry(matdata.coefs(i,colindices[i][j]),
02084                                            elemSize, elemSize, elemSize), ierr);
02085     }
02086 
02087     EPETRA_TEST_ERR( B.EndSubmitEntries(), ierr);
02088   }
02089 
02090   EPETRA_TEST_ERR( B.FillComplete(), ierr);
02091 
02092   return(0);
02093 }
02094 
02095 int checkVbrRowMatrix(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {  
02096 
02097   if (verbose) cout << "Checking VbrRowMatrix Adapter..." << endl;
02098   int ierr = 0;
02099   EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
02100   EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
02101   EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
02102   EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
02103   // EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
02104   // EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
02105   EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
02106   EPETRA_TEST_ERR(!A.NumGlobalCols()==B.NumGlobalCols(),ierr);
02107   EPETRA_TEST_ERR(!A.NumGlobalDiagonals()==B.NumGlobalDiagonals(),ierr);
02108   EPETRA_TEST_ERR(!A.NumGlobalNonzeros()==B.NumGlobalNonzeros(),ierr);
02109   EPETRA_TEST_ERR(!A.NumGlobalRows()==B.NumGlobalRows(),ierr);
02110   EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
02111   EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
02112   EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
02113   for (int i=0; i<A.NumMyRows(); i++) {
02114     int nA, nB;
02115     A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
02116     EPETRA_TEST_ERR(!nA==nB,ierr);
02117   }
02118   EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
02119   EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
02120   EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
02121   EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
02122   EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
02123   // EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
02124   EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);
02125 
02126   int NumVectors = 5;
02127   { // No transpose case
02128     Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
02129     Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
02130     Epetra_MultiVector YA2(YA1);
02131     Epetra_MultiVector YB1(YA1);
02132     Epetra_MultiVector YB2(YA1);
02133     X.Random();
02134     
02135     bool transA = false;
02136     A.SetUseTranspose(transA);
02137     B.SetUseTranspose(transA);
02138     A.Apply(X,YA1);
02139     A.Multiply(transA, X, YA2);
02140     EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
02141     B.Apply(X,YB1);
02142     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
02143     B.Multiply(transA, X, YB2);
02144     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);
02145     
02146   }
02147   {// transpose case
02148     Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
02149     Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
02150     Epetra_MultiVector YA2(YA1);
02151     Epetra_MultiVector YB1(YA1);
02152     Epetra_MultiVector YB2(YA1);
02153     X.Random();
02154     
02155     bool transA = true;
02156     A.SetUseTranspose(transA);
02157     B.SetUseTranspose(transA);
02158     A.Apply(X,YA1);
02159     A.Multiply(transA, X, YA2);
02160     EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
02161     B.Apply(X,YB1);
02162     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
02163     B.Multiply(transA, X,YB2);
02164     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);
02165     
02166   }
02167 
02168   /*
02169     Epetra_Vector diagA(A.RowMatrixRowMap());
02170     EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
02171     Epetra_Vector diagB(B.RowMatrixRowMap());
02172     EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
02173     EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);
02174   */
02175   Epetra_Vector rowA(A.RowMatrixRowMap());
02176   EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
02177   Epetra_Vector rowB(B.RowMatrixRowMap());
02178   EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
02179   EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);
02180 
02181   Epetra_Vector colA(A.OperatorDomainMap());
02182   EPETRA_TEST_ERR(A.InvColSums(colA),ierr);  // -2 error code
02183   Epetra_Vector colB(B.OperatorDomainMap());
02184   EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
02185   EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr); // 1 error code
02186 
02187   EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
02188   EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);
02189 
02190   EPETRA_TEST_ERR(A.RightScale(colA),ierr);  // -3 error code
02191   EPETRA_TEST_ERR(B.RightScale(colB),ierr);  // -3 error code
02192 
02193 
02194   EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
02195   EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);
02196 
02197 
02198   EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf after scaling", verbose), ierr);
02199   EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne after scaling", verbose),ierr);
02200 
02201   vector<double> valuesA(A.MaxNumEntries());
02202   vector<int> indicesA(A.MaxNumEntries());  
02203   vector<double> valuesB(B.MaxNumEntries());
02204   vector<int> indicesB(B.MaxNumEntries());
02205   //return(0);
02206   for (int i=0; i<A.NumMyRows(); i++) {
02207     int nA, nB;
02208     EPETRA_TEST_ERR(A.ExtractMyRowCopy(i, A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr); 
02209     EPETRA_TEST_ERR(B.ExtractMyRowCopy(i, B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
02210     EPETRA_TEST_ERR(!nA==nB,ierr);
02211     for (int j=0; j<nA; j++) {
02212       double curVal = valuesA[j];
02213       int curIndex = indicesA[j];
02214       bool notfound = true;
02215       int jj = 0;
02216       while (notfound && jj< nB) {
02217   if (!checkValues(curVal, valuesB[jj])) notfound = false;
02218   jj++;
02219       }
02220       EPETRA_TEST_ERR(notfound, ierr);
02221       vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex);  // find curIndex in indicesB
02222       EPETRA_TEST_ERR(p==indicesB.end(), ierr);
02223     }
02224 
02225   }
02226 
02227   if (verbose) {
02228     if (ierr==0)
02229       cout << "RowMatrix Methods check OK" << endl;
02230     else
02231       cout << "ierr = " << ierr << ". RowMatrix Methods check failed." << endl;
02232   }
02233   return (ierr);
02234 }
02235 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines