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