test/VbrMatrix/cxx_main.cpp

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

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1