Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/VbrMatrix/cxx_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
43 
44 //
45 //
46 // valgrind usage:
47 // valgrind --suppressions=Suppressions --leak-check=yes --show-reachable=yes ./VbrMatrix_test.exe
48 //
49 // mpirun -np 2 valgrind --suppressions=Suppressions --logfile=valg.out --leak-check=yes --show-reachable=yes ./VbrMatrix_test.exe
50 // The file Suppressions can be found in packages/epetra/test/VbrMatrix/Suppressions.in
51 //
52 //
54 #include "Epetra_Map.h"
55 #include "Epetra_Time.h"
56 #include "Epetra_Vector.h"
57 #include "Epetra_Flops.h"
58 #include "Epetra_VbrMatrix.h"
59 #include "Epetra_VbrRowMatrix.h"
60 #include "Epetra_CrsMatrix.h"
61 #include <vector>
62 #ifdef EPETRA_MPI
63 #include "Epetra_MpiComm.h"
64 #include "mpi.h"
65 #else
66 #include "Epetra_SerialComm.h"
67 #endif
68 #include "../epetra_test_err.h"
69 #include "../src/Epetra_matrix_data.h"
70 #include "Epetra_Version.h"
71 
72 // prototypes
73 
74 int checkValues( double x, double y, string message = "", bool verbose = false) {
75  if (fabs((x-y)/x) > 0.01) {
76  return(1);
77  if (verbose) cout << "********** " << message << " check failed.********** " << endl;
78  }
79  else {
80  if (verbose) cout << message << " check OK." << endl;
81  return(0);
82  }
83 }
84 
85 int checkMultiVectors( Epetra_MultiVector & X, Epetra_MultiVector & Y, string message = "", bool verbose = false) {
86  int numVectors = X.NumVectors();
87  int length = Y.MyLength();
88  int badvalue = 0;
89  int globalbadvalue = 0;
90  for (int j=0; j<numVectors; j++)
91  for (int i=0; i< length; i++)
92  if (checkValues(X[j][i], Y[j][i])==1) badvalue = 1;
93  X.Map().Comm().MaxAll(&badvalue, &globalbadvalue, 1);
94 
95  if (verbose) {
96  if (globalbadvalue==0) cout << message << " check OK." << endl;
97  else cout << "********* " << message << " check failed.********** " << endl;
98  }
99  return(globalbadvalue);
100 }
101 
102 int checkVbrMatrixOptimizedGraph(Epetra_Comm& comm, bool verbose);
103 
104 int checkVbrRowMatrix(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose);
105 
106 int CompareValues(double * A, int LDA, int NumRowsA, int NumColsA,
107  double * B, int LDB, int NumRowsB, int NumColsB);
108 
110  int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1,
111  int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1,
112  int * MyGlobalElements, bool verbose);
113 
114 int power_method(bool TransA, Epetra_VbrMatrix& A,
117  Epetra_MultiVector& resid,
118  double * lambda, int niters, double tolerance,
119  bool verbose);
120 
121 int checkMergeRedundantEntries(Epetra_Comm& comm, bool verbose);
122 
123 int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose);
124 
125 int checkMatvecSameVectors(Epetra_Comm& comm, bool verbose);
126 
127 int checkEarlyDelete(Epetra_Comm& comm, bool verbose);
128 
129 //
130 // ConvertVbrToCrs is a crude but effective way to convert a Vbr matrix to a Crs matrix
131 // Caller is responsible for deleting the CrsMatrix CrsOut
132 //
134 
135  const Epetra_Map &E_Vbr_RowMap = VbrIn->RowMatrixRowMap() ;
136  const Epetra_Map &E_Vbr_ColMap = VbrIn->RowMatrixColMap() ;
137  // const Epetra_Map &E_Vbr_RowMap = VbrIn->OperatorRangeMap() ;
138  // const Epetra_Map &E_Vbr_ColMap = VbrIn->OperatorDomainMap() ;
139 
140  CrsOut = new Epetra_CrsMatrix( Copy, E_Vbr_RowMap, E_Vbr_ColMap, 0 );
141  // CrsOut = new Epetra_CrsMatrix( Copy, E_Vbr_RowMap, 0 );
142  int NumMyElements = VbrIn->RowMatrixRowMap().NumMyElements() ;
143  vector<int> MyGlobalElements( NumMyElements );
144  VbrIn->RowMatrixRowMap().MyGlobalElements( &MyGlobalElements[0] ) ;
145 
146  int NumMyColumns = VbrIn->RowMatrixColMap().NumMyElements() ;
147  vector<int> MyGlobalColumns( NumMyColumns );
148  VbrIn->RowMatrixColMap().MyGlobalElements( &MyGlobalColumns[0] ) ;
149 
150  int MaxNumIndices = VbrIn->MaxNumEntries();
151  int NumIndices;
152  vector<int> LocalColumnIndices(MaxNumIndices);
153  vector<int> GlobalColumnIndices(MaxNumIndices);
154  vector<double> MatrixValues(MaxNumIndices);
155 
156  for( int LocalRow=0; LocalRow<NumMyElements; ++LocalRow ) {
157 
158  VbrIn->ExtractMyRowCopy( LocalRow,
159  MaxNumIndices,
160  NumIndices,
161  &MatrixValues[0],
162  &LocalColumnIndices[0] );
163 
164  for (int j = 0 ; j < NumIndices ; j++ )
165  {
166  GlobalColumnIndices[j] = MyGlobalColumns[ LocalColumnIndices[j] ] ;
167  }
168 
169 #if 0
170  if ( CrsOut->InsertGlobalValues( MyGlobalElements[LocalRow],
171  NumIndices,
172  &MatrixValues[0],
173  &GlobalColumnIndices[0] )!=0)abort();
174 #else
175  if ( CrsOut->InsertMyValues( LocalRow,
176  NumIndices,
177  &MatrixValues[0],
178  &LocalColumnIndices[0] )!= 0) abort();
179 #endif
180 
181 
182  }
183  CrsOut->FillComplete();
184 }
185 
186 //
187 // checkmultiply checks to make sure that AX=Y using both multiply
188 // and both Crs Matrices, multiply1
189 //
190 
192 
193  int numerrors = 0 ;
194 
195 #if 1
196  //
197  // If X and Y are Epetra_Vectors, we first test Multiply1
198  //
199  Epetra_Vector *vecY = dynamic_cast<Epetra_Vector *>( &Check_Y );
200  Epetra_Vector *vecX = dynamic_cast<Epetra_Vector *>( &X );
201  assert( (vecX && vecY) || (!vecX && !vecY) ) ;
202 
203  if ( vecX && vecY ) {
204  double normY, NormError;
205  Check_Y.NormInf( &normY ) ;
206  Epetra_Vector Y = *vecX ;
207  Y.PutScalar( -13.0 ) ;
208  Epetra_Vector Error = *vecX ;
209  A.Multiply1( transpose, *vecX, Y ) ;
210 
211  Error.Update( 1.0, Y, -1.0, *vecY, 0.0 ) ;
212 
213  Error.NormInf( &NormError ) ;
214  if ( NormError / normY > 1e-13 ) {
215  numerrors++;
216  //cout << "Y = " << Y << endl;
217  //cout << "vecY " << *vecY << endl;
218  //cout << "Error " << Error << endl;
219  //abort();
220  }
221  //
222  // Check x = Ax
223  //
224  Epetra_Vector Z = *vecX ;
225 
226  A.Multiply1( transpose, Z, Z ) ;
227  Error.Update( 1.0, Z, -1.0, *vecY, 0.0 ) ;
228  // Error.Update( 1.0, Y, -1.0, *vecY, 0.0 ) ;
229 
230  Error.NormInf( &NormError ) ;
231  if ( NormError / normY > 1e-13 ) numerrors++;
232  }
233  //
234  // Here we test Multiply
235  //
236  Epetra_MultiVector Y = X ;
237  Epetra_MultiVector Error = X ;
238  A.Multiply( transpose, X, Y ) ;
239 
240  int NumVecs = X.NumVectors() ;
241  vector<double> NormError(NumVecs);
242  vector<double> Normy(NumVecs);
243 
244  Check_Y.NormInf( &Normy[0] ) ;
245 
246  Error.Update( 1.0, Y, -1.0, Check_Y, 0.0 ) ;
247  Error.NormInf( &NormError[0] ) ;
248 
249  bool LoopError = false ;
250  for ( int ii = 0 ; ii < NumVecs ; ii++ ) {
251  if( NormError[ii] / Normy[ii] > 1e-13 ) {
252  LoopError = true ;
253  }
254  }
255  if ( LoopError ) {
256  numerrors++ ;
257  }
258 
259  //
260  // Check X = AX
261  //
262 
263 #endif
264 
265  return numerrors;
266 }
267 //
268 // TestMatrix contains the bulk of the testing.
269 // MinSize and MaxSize control the range of the block sizes - which are chosen randomly
270 // ConstructWithNumNz controls whether a Column Map is passed to the VbrMatrix contructor
271 // if false, the underlying graph will not be optimized
272 // ExtraBlocks, if true, causes extra blocks to be added to the matrix, further
273 // guaranteeing that the matrix and graph will not have optimal storage.
274 // ExtraBlocks is only supported for fixed block sizes (i.e. when minsize=maxsize)
275 //
276 // If TestMatrix() is called with any of the ConsTypes that use PreviousA, the values used to
277 // build A, i.e. MinSize, MaxSize must be the same as they were on the previous call to
278 // TestMatrix (i.e. the call that built PreviousA). Furthermore, the ConsType must be a
279 // matching ConsType.
280 // The ConsType values that cause TestMatrix to
281 // use PreviousA are:
282 // RowMapColMap_VEPR, RowMapColMap_FEPR, RowMapColMap_NEPR,
283 // WithGraph, CopyConstructor
284 // The matching ConsTypes are:
285 // VariableEntriesPerRow, RowMapColMap_VEPR, NoEntriesPerRow, RowMapColMap_NEPR, WithGraph
286 // FixedEntriesPerRow, RowMapColMap_FEPR
287 //
288 // TestMatrix() when called with ConsType=WithGraph, returns with PreviousA = 0 ;
289 //
290 //
294 
295 int TestMatrix( Epetra_Comm& Comm, bool verbose, bool debug,
296  int NumMyElements, int MinSize, int MaxSize,
297  ConsType ConstructorType, bool ExtraBlocks,
298  bool insertlocal,
299  bool symmetric,
300  Epetra_VbrMatrix** PreviousA ) {
301 
302 
303  if (verbose)
304  cout << "MinSize = " << MinSize << endl
305  << "MaxSize = " << MaxSize << endl
306  << "ConstructorType = " << ConstructorType << endl
307  << "ExtraBlocks = " << ExtraBlocks << endl
308  << "insertlocal = " << insertlocal << endl
309  << "symmetric = " << symmetric << endl
310  << "PreviousA = " << PreviousA << endl;
311 
312  int ierr = 0, forierr = 0;
313  int MyPID = Comm.MyPID();
314  if (MyPID < 3) NumMyElements++;
315  if (NumMyElements<2) NumMyElements = 2; // This value must be greater than one on each processor
316 
317  // Define pseudo-random block sizes using an Epetra Vector of random numbers
318  Epetra_Map randmap(-1, NumMyElements, 0, Comm);
319  Epetra_Vector randvec(randmap);
320  randvec.Random(); // Fill with random numbers
321  int * ElementSizeList = new int[NumMyElements];
322  int SizeRange = MaxSize - MinSize + 1;
323  double DSizeRange = SizeRange;
324 
325  const Epetra_BlockMap* rowmap = 0 ;
326  const Epetra_BlockMap* colmap = 0 ;
327  Epetra_CrsGraph* graph = 0 ;
328  if ( *PreviousA != 0 ) {
329 
330  rowmap = &((*PreviousA)->RowMap());
331  colmap = &((*PreviousA)->ColMap());
332  }
333 
334  ElementSizeList[0] = MaxSize;
335  for (int i=1; i<NumMyElements-1; i++) {
336  int curSize = MinSize + (int) (DSizeRange * fabs(randvec[i]) + .99);
337  ElementSizeList[i] = EPETRA_MAX(MinSize, EPETRA_MIN(MaxSize, curSize));
338  }
339  ElementSizeList[NumMyElements-1] = MaxSize;
340 
341 
342 
343  // Construct a Map
344 
345  int *randMyGlobalElements = randmap.MyGlobalElements();
346 
347  if ( ConstructorType == RowMapColMap_VEPR ||
348  ConstructorType == RowMapColMap_FEPR ||
349  ConstructorType == RowMapColMap_NEPR ||
350  ConstructorType == WithGraph ) {
351  rowmap->ElementSizeList( ElementSizeList ) ;
352  }
353 
354  Epetra_BlockMap Map (-1, NumMyElements, randMyGlobalElements, ElementSizeList, 0, Comm);
355 
356  // Get update list and number of local elements from newly created Map
357  int NumGlobalElements = Map.NumGlobalElements();
358  int * MyGlobalElements = Map.MyGlobalElements();
359 
360  // Create an integer vector NumNz that is used to build the Petra Matrix.
361  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
362 
363  int * NumNz = new int[NumMyElements];
364 
365  // We are building a block tridiagonal matrix
366 
367  for (int i=0; i<NumMyElements; i++)
368  if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
369  NumNz[i] = 2;
370  else
371  NumNz[i] = 3;
372  // Create a Epetra_Matrix
373 
374  bool FixedNumEntries = false ; // If FixedNumEntries == true, we add the upper left and lower right hand corners to
375  // our tri-diagonal matrix so that each row has exactly three entries
376  bool HaveColMap = false ; // Matrices constructed without a column map cannot use insertmy to create the matrix
377  bool FixedBlockSize = ( MinSize == MaxSize ) ;
378  bool HaveGraph = false ;
379 
380 
382 
383  switch( ConstructorType ) {
385  A = new Epetra_VbrMatrix( Copy, Map, NumNz ) ;
386  break;
387  case FixedEntriesPerRow:
388  A = new Epetra_VbrMatrix( Copy, Map, 3 ) ;
389  FixedNumEntries = true;
390  break;
391  case NoEntriesPerRow:
392  A = new Epetra_VbrMatrix( Copy, Map, 0 ) ;
393  break;
394  case RowMapColMap_VEPR:
395  A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, NumNz ) ;
396  HaveColMap = true;
397  break;
398  case RowMapColMap_FEPR:
399  A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, 3 ) ;
400  FixedNumEntries = true;
401  HaveColMap = true;
402  break;
403  case RowMapColMap_NEPR:
404  A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, 0 ) ;
405  HaveColMap = true;
406  break;
407  case WithGraph:
408  graph = new Epetra_CrsGraph( (*PreviousA)->Graph() );
409  A = new Epetra_VbrMatrix( Copy, *graph );
410  HaveGraph = true;
411  HaveColMap = true;
412  break;
413  case CopyConstructor:
414  A = new Epetra_VbrMatrix( (**PreviousA) );
415  HaveColMap = true;
416  break;
417  default:
418  assert(false);
419  }
420 
421  if ( insertlocal ) assert( HaveColMap ); // you can't insert local without a column map
422  if ( ExtraBlocks ) assert( FixedBlockSize ); // ExtraBlocks is only supported for fixed block sizes
423  if ( ExtraBlocks ) assert( ! HaveColMap ); // Matrices constructed with a column map cannot be extended
424  if ( FixedNumEntries) assert( FixedBlockSize ) ; // Can't handle a Fixed Number of Entries and a variable block size
425  if ( insertlocal && HaveGraph ) assert( ! FixedNumEntries ) ; // HaveGraph assumes the standard matrix shape
426  if ( insertlocal && HaveGraph ) assert( ! ExtraBlocks ) ; // HaveGraph assumes the standard matrix shape
427 
428 
429  // WORK Insert/Replace/Suminto MY should fail here when there is no colmap
430 
431 
433  if ( ! HaveGraph ) EPETRA_TEST_ERR(A->IndicesAreLocal(),ierr);
434 
435  // Use an array of Epetra_SerialDenseMatrix objects to build VBR matrix
436 
437  Epetra_SerialDenseMatrix ** BlockEntries = new Epetra_SerialDenseMatrix*[SizeRange];
438 
439  // The array of dense matrices will increase in size from MinSize to MaxSize (defined above)
440  for (int kr=0; kr<SizeRange; kr++) {
441  BlockEntries[kr] = new Epetra_SerialDenseMatrix[SizeRange];
442  int RowDim = MinSize+kr;
443  for (int kc = 0; kc<SizeRange; kc++) {
444  int ColDim = MinSize+kc;
445  Epetra_SerialDenseMatrix * curmat = &(BlockEntries[kr][kc]);
446  curmat->Shape(RowDim,ColDim);
447  for (int j=0; j < ColDim; j++)
448  for (int i=0; i < RowDim; i++) {
449  BlockEntries[kr][kc][j][i] = -1.0;
450  if (i==j && kr==kc) BlockEntries[kr][kc][j][i] = 9.0;
451  else BlockEntries[kr][kc][j][i] = -1.0;
452 
453  if ( ! symmetric ) BlockEntries[kr][kc][j][i] += ((double) j)/10000.0;
454  }
455  }
456  }
457 
458  // Add rows one-at-a-time
459 
460  int *Indices = new int[3];
461  int *MyIndices = new int[3];
462  int *ColDims = new int[3];
463  int NumEntries;
464  int NumMyNonzeros = 0, NumMyEquations = 0;
465 
466 
467 
468  for (int i=0; i<NumMyElements; i++) {
469  int CurRow = MyGlobalElements[i];
470  if ( HaveColMap ) {
471  assert ( i == rowmap->LID( CurRow ) ) ;
472  assert ( CurRow == rowmap->GID( i ) ) ;
473  }
474  int RowDim = ElementSizeList[i]-MinSize;
475  NumMyEquations += BlockEntries[RowDim][RowDim].M();
476 
477  if (CurRow==0)
478  {
479  Indices[0] = CurRow;
480  Indices[1] = CurRow+1;
481  if ( FixedNumEntries ) {
482  Indices[2] = NumGlobalElements-1;
483  ColDims[2] = 0 ;
484  assert( ElementSizeList[i] == MinSize ); // This is actually enforced above as well
485  NumEntries = 3;
486  } else {
487  NumEntries = 2;
488  }
489  ColDims[0] = ElementSizeList[i] - MinSize;
490  ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
491  }
492  else if (CurRow == NumGlobalElements-1)
493  {
494  Indices[0] = CurRow-1;
495  Indices[1] = CurRow;
496  if ( FixedNumEntries ) {
497  Indices[2] = 0;
498  ColDims[2] = 0 ;
499  assert( ElementSizeList[i] == MinSize ); // This is actually enforced above as well
500  NumEntries = 3;
501  } else {
502  NumEntries = 2;
503  }
504  ColDims[0] = ElementSizeList[i-1] - MinSize;
505  ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
506  }
507  else {
508  Indices[0] = CurRow-1;
509  Indices[1] = CurRow;
510  Indices[2] = CurRow+1;
511  NumEntries = 3;
512  if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
513  else ColDims[0] = ElementSizeList[i-1] - MinSize;
514  ColDims[1] = ElementSizeList[i] - MinSize;
515  // ElementSize on MyPID+1
516  if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
517  else ColDims[2] = ElementSizeList[i+1] - MinSize;
518  }
519  if ( insertlocal ) {
520  for ( int ii=0; ii < NumEntries; ii++ )
521  MyIndices[ii] = colmap->LID( Indices[ii] ) ;
522  // Epetra_MpiComm* MComm = dynamic_cast<Epetra_MpiComm*>( &Comm ) ;
523  // MComm->SetTracebackMode(1); // This should enable error traceback reporting
524  if ( HaveGraph ) {
525  EPETRA_TEST_ERR(!(A->BeginReplaceMyValues(rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
526  } else {
527  EPETRA_TEST_ERR(!(A->BeginInsertMyValues(rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
528  }
529  // MComm->SetTracebackMode(0); // This should shut down any error traceback reporting
530  } else {
531  if ( HaveGraph ) {
532  EPETRA_TEST_ERR(!(A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
533  } else {
534  //
535  // I, Ken Stanley, think the following call should return an error since it
536  // makes no sense to insert a value with a local index in the absence of a
537  // map indicating what that index means. Instead, this call returns with an
538  // error code of 0, but problems appear later.
539  //
540  // EPETRA_TEST_ERR((A->BeginInsertMyValues(CurRow, NumEntries, Indices)==0),ierr); // Should fail
541  EPETRA_TEST_ERR(!(A->BeginInsertGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
542  }
543  }
544  forierr = 0;
545  for (int j=0; j < NumEntries; j++) {
546  Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
547  NumMyNonzeros += AD->M() * AD->N();
548  forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
549  }
550  EPETRA_TEST_ERR(forierr,ierr);
551 
552  A->EndSubmitEntries();
553  }
554 
555  int NumMyBlockEntries = 3*NumMyElements;
556  if ( ! FixedNumEntries ) {
557  if (A->LRID(0)>=0) NumMyBlockEntries--; // If I own first global row, then there is one less nonzero
558  if (A->LRID(NumGlobalElements-1)>=0) NumMyBlockEntries--; // If I own last global row, then there is one less nonzero
559  }
560 
561  if ( ExtraBlocks ) {
562 
563 
564  //
565  // Add a block to the matrix on each process.
566  // The i index is chosen from among the block rows that this process owns (i.e. MyGlobalElements[0..NumMyElements-1])
567  // The j index is chosen from among all block columns
568  //
569  // Bugs - does not support non-contiguous matrices
570  //
571  // Adding more than one off diagonal block could have resulted in adding an off diagonal block
572  // twice, resulting in errors in NumMyNonzeros and NumMyBlockEntries
573  //
574  const int NumTries = 100;
575  Epetra_SerialDenseMatrix BlockIindex = Epetra_SerialDenseMatrix( NumTries, 1 );
576  Epetra_SerialDenseMatrix BlockJindex = Epetra_SerialDenseMatrix( NumTries, 1 );
577 
578  BlockIindex.Random();
579  BlockJindex.Random();
580 
581  BlockIindex.Scale( 1.0 * NumMyElements );
582  BlockJindex.Scale( 1.0 * A->NumGlobalBlockCols() );
583  bool OffDiagonalBlockAdded = false ;
584  for ( int ii=0; ii < NumTries && ! OffDiagonalBlockAdded; ii++ ) {
585  int i = (int) BlockIindex[0][ii];
586  int j = (int) BlockJindex[0][ii];
587  if ( i < 0 ) i = - i ;
588  if ( j < 0 ) j = - j ;
589  assert( i >= 0 ) ;
590  assert( j >= 0 ) ;
591  assert( i < NumMyElements ) ;
592  assert( j < A->NumGlobalBlockCols() ) ;
593  int CurRow = MyGlobalElements[i];
594  int IndicesL[1] ;
595  IndicesL[0] = j ;
596  Epetra_SerialDenseMatrix * AD = &(BlockEntries[0][0]);
597  EPETRA_TEST_ERR(!(A->BeginInsertGlobalValues( CurRow, 1, IndicesL)==0),ierr);
598 
599  if ( CurRow < j-1 || CurRow > j+1 ) {
600  OffDiagonalBlockAdded = true ;
601  NumMyNonzeros += AD->M() * AD->N();
602  NumMyBlockEntries++ ;
603  }
604 
605  // EPETRA_TEST_ERR(!(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0), ierr);
606  EPETRA_TEST_ERR(!(A->SubmitBlockEntry(*AD)==0), ierr);
607  A->EndSubmitEntries();
608  }
609  }
610 
611  // Finish up
612  if ( ! HaveGraph && ! insertlocal )
613  EPETRA_TEST_ERR(!(A->IndicesAreGlobal()),ierr);
614  EPETRA_TEST_ERR(!(A->FillComplete()==0),ierr);
615  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
616  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
618 
619  A->OptimizeStorage();
620  if ( FixedBlockSize ) {
621  EPETRA_TEST_ERR(!(A->StorageOptimized()),ierr);
622  } else {
623  // EPETRA_TEST_ERR(A->StorageOptimized(),ierr); // Commented out until I figure out why it occasionally fails on one process
624  }
625  EPETRA_TEST_ERR(A->UpperTriangular(),ierr);
626  EPETRA_TEST_ERR(A->LowerTriangular(),ierr);
627 
628 
629 
630 
631 
632 
633 
634  int NumGlobalBlockEntries ;
635  int NumGlobalNonzeros, NumGlobalEquations;
636  Comm.SumAll(&NumMyBlockEntries, &NumGlobalBlockEntries, 1);
637  Comm.SumAll(&NumMyNonzeros, &NumGlobalNonzeros, 1);
638 
639  Comm.SumAll(&NumMyEquations, &NumGlobalEquations, 1);
640 
641  if (! ExtraBlocks ) {
642  if ( FixedNumEntries ) {
643  EPETRA_TEST_ERR( !( NumGlobalBlockEntries == (3*NumGlobalElements)), ierr );
644  } else {
645  EPETRA_TEST_ERR( !( NumGlobalBlockEntries == (3*NumGlobalElements-2)), ierr );
646  }
647  }
648 
649 
650  EPETRA_TEST_ERR(!(check(*A, NumMyEquations, NumGlobalEquations, NumMyNonzeros, NumGlobalNonzeros,
651  NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries,
652  MyGlobalElements, verbose)==0),ierr);
653 
654  forierr = 0;
655  if ( ! ExtraBlocks ) {
656  if ( FixedNumEntries )
657  for (int i=0; i<NumMyElements; i++) forierr += !(A->NumGlobalBlockEntries(MyGlobalElements[i])==3);
658  else
659  for (int i=0; i<NumMyElements; i++) forierr += !(A->NumGlobalBlockEntries(MyGlobalElements[i])==NumNz[i]);
660  }
661  EPETRA_TEST_ERR(forierr,ierr);
662  forierr = 0;
663  if ( ! ExtraBlocks ) {
664  if ( FixedNumEntries )
665  for (int i=0; i<NumMyElements; i++) forierr += !(A->NumMyBlockEntries(i)==3);
666  else
667  for (int i=0; i<NumMyElements; i++) forierr += !(A->NumMyBlockEntries(i)==NumNz[i]);
668  }
669  EPETRA_TEST_ERR(forierr,ierr);
670 
671  if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
672 
673  delete [] NumNz;
674 
675 
676  // Create vectors for Power method
677 
678  Epetra_Vector q(Map);
679  Epetra_Vector z(Map);
680  Epetra_Vector z_initial(Map);
681  Epetra_Vector resid(Map);
682 
683 
684  // Fill z with random Numbers
685  z_initial.Random();
686 
687  // variable needed for iteration
688  double lambda = 0.0;
689  int niters = 100;
690  // int niters = 200;
691  double tolerance = 1.0e-3;
692 
694 
695  // Iterate
696  Epetra_Flops flopcounter;
697  A->SetFlopCounter(flopcounter);
698  q.SetFlopCounter(*A);
699  z.SetFlopCounter(*A);
700  resid.SetFlopCounter(*A);
701  z = z_initial; // Start with common initial guess
702  Epetra_Time timer(Comm);
703  int ierr1 = power_method(false, *A, q, z, resid, &lambda, niters, tolerance, verbose);
704  double elapsed_time = timer.ElapsedTime();
705  double total_flops = flopcounter.Flops();
706  double MFLOPs = total_flops/elapsed_time/1000000.0;
707 
708  if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
709  if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
710 
712 
713  // Solve transpose problem
714 
715  if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
716  << endl;
717  // Iterate
718  lambda = 0.0;
719  z = z_initial; // Start with common initial guess
720  flopcounter.ResetFlops();
721  timer.ResetStartTime();
722  ierr1 = power_method(true, *A, q, z, resid, &lambda, niters, tolerance, verbose);
723  elapsed_time = timer.ElapsedTime();
724  total_flops = flopcounter.Flops();
725  MFLOPs = total_flops/elapsed_time/1000000.0;
726 
727  if (verbose) cout << "\n\nTotal MFLOPs for transpose solve = " << MFLOPs << endl<< endl;
728  if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
729 
731 
732  Epetra_CrsMatrix* OrigCrsA;
733  ConvertVbrToCrs( A, OrigCrsA ) ;
734 
735  // Increase diagonal dominance
736 
737  if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
738  << endl;
739 
740  double AnormInf = -13 ;
741  double AnormOne = -13 ;
742  AnormInf = A->NormInf( ) ;
743  AnormOne = A->NormOne( ) ;
744 
745  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
746  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
747 
748  if (A->MyGlobalBlockRow(0)) {
749  int numvals = A->NumGlobalBlockEntries(0);
750  Epetra_SerialDenseMatrix ** Rowvals;
751  int* Rowinds = new int[numvals];
752  int RowDim;
753  A->ExtractGlobalBlockRowPointers(0, numvals, RowDim, numvals, Rowinds,
754  Rowvals); // Get A[0,:]
755 
756  for (int i=0; i<numvals; i++) {
757  if (Rowinds[i] == 0) {
758  // Rowvals[i]->A()[0] *= 10.0; // Multiply first diag value by 10.0
759  Rowvals[i]->A()[0] += 1000.0; // Add 1000 to first diag value
760  }
761  }
762  delete [] Rowinds;
763  }
764  //
765  // NormOne() and NormInf() will NOT return cached values
766  // See bug #1151
767  //
768 
769  EPETRA_TEST_ERR( ! (AnormOne != A->NormOne( )), ierr );
770  EPETRA_TEST_ERR( ! (AnormInf != A->NormInf( )), ierr );
771  //
772  // On Process 0, let the class know that NormInf_ and NormOne_ are
773  // out of date.
774  //
775  if ( MyPID == 0 ) {
776  EPETRA_TEST_ERR(!(A->BeginSumIntoGlobalValues( 0, 0, 0 )==0),ierr);
777  EPETRA_TEST_ERR( A->EndSubmitEntries(), ierr );
778  }
779  EPETRA_TEST_ERR( ! (AnormOne != A->NormOne( )), ierr );
780  EPETRA_TEST_ERR( ! (AnormInf != A->NormInf( )), ierr );
781  if ( MyPID == 0 )
782  // Iterate (again)
783  lambda = 0.0;
784  z = z_initial; // Start with common initial guess
785  flopcounter.ResetFlops();
786  timer.ResetStartTime();
787  ierr1 = power_method(false, *A, q, z, resid, &lambda, niters, tolerance, verbose);
788  elapsed_time = timer.ElapsedTime();
789  total_flops = flopcounter.Flops();
790  MFLOPs = total_flops/elapsed_time/1000000.0;
791 
792  if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
793  if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
794 
795 
796 
798 
799  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
800  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
801  // Solve transpose problem
802 
803  if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
804  << endl;
805 
806  // Iterate (again)
807  lambda = 0.0;
808  z = z_initial; // Start with common initial guess
809  flopcounter.ResetFlops();
810  timer.ResetStartTime();
811  ierr1 = power_method(true, *A, q, z, resid, &lambda, niters, tolerance, verbose);
812  elapsed_time = timer.ElapsedTime();
813  total_flops = flopcounter.Flops();
814  MFLOPs = total_flops/elapsed_time/1000000.0;
815 
816  if (verbose) cout << "\n\nTotal MFLOPs for tranpose of second solve = " << MFLOPs << endl<< endl;
817  if (verbose && ierr1==1) cout << "***** Power Method did not converge. *****" << endl << endl;
818 
819 
820  if (debug) Comm.Barrier();
821 
822  if (verbose) cout << "\n\n*****Comparing against CrsMatrix " << endl<< endl;
823 
824  Epetra_CrsMatrix* CrsA;
825  ConvertVbrToCrs( A, CrsA ) ;
826 
827  Epetra_Vector CrsX = Epetra_Vector( A->OperatorDomainMap(), false ) ;
828  Epetra_Vector CrsY = Epetra_Vector( A->OperatorRangeMap(), false ) ;
829  Epetra_Vector OrigCrsY = Epetra_Vector( A->OperatorRangeMap(), false ) ;
830  Epetra_Vector Y_Apply = Epetra_Vector( A->OperatorRangeMap(), false ) ;
831  Epetra_Vector x(Map);
832  Epetra_Vector y(Map);
833  Epetra_Vector orig_check_y(Map);
834  Epetra_Vector Apply_check_y(Map);
835  Epetra_Vector check_y(Map);
836  Epetra_Vector check_ytranspose(Map);
837 
838  x.Random() ;
839  CrsX = x;
840  CrsY = CrsX;
841  CrsA->Multiply( false, CrsX, CrsY ) ;
842  OrigCrsA->Multiply( false, CrsX, OrigCrsY ) ;
843 
844 
845  check_y = CrsY ;
846  orig_check_y = OrigCrsY ;
847  CrsA->Multiply( true, CrsX, CrsY ) ;
848  check_ytranspose = CrsY ;
849 
850  EPETRA_TEST_ERR( checkmultiply( false, *A, x, check_y ), ierr );
851 
852  EPETRA_TEST_ERR( checkmultiply( true, *A, x, check_ytranspose ), ierr );
853 
854  if (! symmetric ) EPETRA_TEST_ERR( !checkmultiply( false, *A, x, check_ytranspose ), ierr ); // Just to confirm that A is not symmetric
855 
856  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
857  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
858 
859  if (verbose) cout << "\n\n*****Try the Apply method " << endl<< endl;
860 
861 
862  A->Apply( CrsX, Y_Apply ) ;
863  Apply_check_y = Y_Apply ;
864  EPETRA_TEST_ERR( checkmultiply( false, *A, x, Apply_check_y ), ierr );
865 
866  if (verbose) cout << "\n\n*****Multiply multivectors " << endl<< endl;
867 
868  const int NumVecs = 4 ;
869 
870  Epetra_Map Amap = A->OperatorDomainMap() ;
871 
872  // Epetra_MultiVector CrsMX = Epetra_MultiVector( A->OperatorDomainMap(), false ) ;
873  Epetra_MultiVector CrsMX( Amap, NumVecs, false ) ;
874  Epetra_MultiVector CrsMY( A->OperatorRangeMap(), NumVecs, false ) ;
875  Epetra_MultiVector mx(Map, NumVecs);
876  Epetra_MultiVector my(Map, NumVecs);
877  Epetra_MultiVector check_my(Map, NumVecs);
878  Epetra_MultiVector check_mytranspose(Map, NumVecs);
879  mx.Random(); // Fill mx with random numbers
880 #if 0
881  CrsMX = mx;
882  CrsA->Multiply( false, CrsMX, CrsMY ) ;
883 #else
884  CrsMY = mx;
885  CrsA->Multiply( false, CrsMY, CrsMY ) ;
886 #endif
887  check_my = CrsMY ;
888  CrsMY = mx;
889  CrsA->Multiply( true, CrsMY, CrsMY ) ;
890  check_mytranspose = CrsMY ;
891 
892 
893  EPETRA_TEST_ERR( checkmultiply( false, *A, mx, check_my ), ierr );
894 
895  EPETRA_TEST_ERR( checkmultiply( true, *A, mx, check_mytranspose ), ierr );
896 
897 
898  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
899  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
900  if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
901 
902  // B was changed to a pointer so that we could delete it before the
903  // underlying graph is deleted. This was necessary before Bug #1116 was
904  // fixed, to avoid seg-faults. Bug #1116 is now fixed so this is no longer
905  // an issue.
907 
908  EPETRA_TEST_ERR(!(check(*B, NumMyEquations, NumGlobalEquations, NumMyNonzeros, NumGlobalNonzeros,
909  NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries,
910  MyGlobalElements, verbose)==0),ierr);
911 
912  EPETRA_TEST_ERR( ! ( A->StorageOptimized() == B->StorageOptimized() ), ierr ) ;
913 
914  EPETRA_TEST_ERR( checkmultiply( false, *B, mx, check_my ), ierr );
915 
916 
917  *B = *B; // Should be harmless - check to make sure that it is
918 
919  EPETRA_TEST_ERR(!(check(*B, NumMyEquations, NumGlobalEquations, NumMyNonzeros, NumGlobalNonzeros,
920  NumMyElements, NumGlobalElements, NumMyBlockEntries, NumGlobalBlockEntries,
921  MyGlobalElements, verbose)==0),ierr);
922 
923  EPETRA_TEST_ERR( ! ( A->StorageOptimized() == B->StorageOptimized() ), ierr ) ;
924 
925  EPETRA_TEST_ERR( checkmultiply( false, *B, mx, check_my ), ierr );
926 
927  AnormInf = A->NormInf( );
928  AnormOne = A->NormOne( );
929  EPETRA_TEST_ERR( ! (AnormOne == B->NormOne( )), ierr );
930  EPETRA_TEST_ERR( ! (AnormInf == B->NormInf( )), ierr );
931 
932 
933  Epetra_CrsMatrix* CrsB;
934  ConvertVbrToCrs( B, CrsB ) ;
935 
936  if (verbose) cout << "\n\n*****Testing PutScalar, LeftScale, RightScale, and ReplaceDiagonalValues" << endl<< endl;
937  //
938  // Check PutScalar,
939  //
940  B->PutScalar( 1.0 ) ;
941  CrsB->PutScalar( 1.0 ) ;
942  EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ;
943 
944 
945  check_y = CrsY ;
946  //
947  EPETRA_TEST_ERR( B->ReplaceDiagonalValues( check_y ), ierr ) ;
948  EPETRA_TEST_ERR( CrsB->ReplaceDiagonalValues( CrsY ), ierr ) ;
949  EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ;
950 
951  EPETRA_TEST_ERR( B->LeftScale( check_y ), ierr ) ;
952  EPETRA_TEST_ERR( CrsB->LeftScale( CrsY ), ierr ) ;
953  EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ;
954 
955  EPETRA_TEST_ERR( B->RightScale( check_y ), ierr ) ;
956  EPETRA_TEST_ERR( CrsB->RightScale( CrsY ), ierr ) ;
957  EPETRA_TEST_ERR(! ( B->NormOne() == CrsB->NormOne() ), ierr ) ;
958 
959  double B_norm_frob = B->NormFrobenius();
960  double CrsB_norm_frob = CrsB->NormFrobenius();
961  //need to use a fairly large tolerance when comparing the frobenius
962  //norm from a vbr-matrix and a crs-matrix, because the point-entries
963  //are visited in different orders during the norm calculation, causing
964  //round-off differences to accumulate. That's why we choose 5.e-5
965  //instead of a smaller number like 1.e-13 or so.
966  if (fabs(B_norm_frob-CrsB_norm_frob) > 5.e-5) {
967  std::cout << "fabs(B_norm_frob-CrsB_norm_frob): "
968  << fabs(B_norm_frob-CrsB_norm_frob) << std::endl;
969  std::cout << "VbrMatrix::NormFrobenius test FAILED."<<std::endl;
970  EPETRA_TEST_ERR(-99, ierr);
971  }
972  if (verbose) std::cout << "\n\nVbrMatrix::NormFrobenius OK"<<std::endl<<std::endl;
973 
974  if (debug) Comm.Barrier();
975 
976  if (verbose) cout << "\n\n*****Testing post construction modifications" << endl<< endl;
977  if (verbose) cout << "\n\n*****Replace methods should be OK" << endl<< endl;
978 
979  // Check to see if we can restore the matrix to its original value
980  // Does not work if ExtraBlocks is true
981 
982  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
983  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
984 
985  if ( ! ExtraBlocks )
986  {
987  // Add rows one-at-a-time
988 
989  NumMyNonzeros = NumMyEquations = 0;
990 
991  for (int i=0; i<NumMyElements; i++) {
992  int CurRow = MyGlobalElements[i];
993  int RowDim = ElementSizeList[i]-MinSize;
994  NumMyEquations += BlockEntries[RowDim][RowDim].M();
995 
996  if (CurRow==0)
997  {
998  Indices[0] = CurRow;
999  Indices[1] = CurRow+1;
1000  NumEntries = 2;
1001  ColDims[0] = ElementSizeList[i] - MinSize;
1002  ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
1003  }
1004  else if (CurRow == NumGlobalElements-1)
1005  {
1006  Indices[0] = CurRow-1;
1007  Indices[1] = CurRow;
1008  NumEntries = 2;
1009  ColDims[0] = ElementSizeList[i-1] - MinSize;
1010  ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
1011  }
1012  else {
1013  Indices[0] = CurRow-1;
1014  Indices[1] = CurRow;
1015  Indices[2] = CurRow+1;
1016  NumEntries = 3;
1017  if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
1018  else ColDims[0] = ElementSizeList[i-1] - MinSize;
1019  ColDims[1] = ElementSizeList[i] - MinSize;
1020  // ElementSize on MyPID+1
1021  if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
1022  else ColDims[2] = ElementSizeList[i+1] - MinSize;
1023  }
1024  EPETRA_TEST_ERR(!(A->IndicesAreLocal()),ierr);
1025  EPETRA_TEST_ERR((A->IndicesAreGlobal()),ierr);
1026  EPETRA_TEST_ERR(!(A->BeginReplaceGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
1027  forierr = 0;
1028  for (int j=0; j < NumEntries; j++) {
1029  Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
1030  NumMyNonzeros += AD->M() * AD->N();
1031  forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
1032  }
1033  EPETRA_TEST_ERR(forierr,ierr);
1034 
1035  A->EndSubmitEntries();
1036  }
1037  EPETRA_TEST_ERR( checkmultiply( false, *A, x, orig_check_y ), ierr );
1038  }
1039 
1040  //
1041  // Suminto should cause the matrix to be doubled
1042  //
1043  if ( ! ExtraBlocks ) {
1044  // Add rows one-at-a-time
1045 
1046  NumMyNonzeros = NumMyEquations = 0;
1047 
1048  for (int i=0; i<NumMyElements; i++) {
1049  int CurRow = MyGlobalElements[i];
1050  int RowDim = ElementSizeList[i]-MinSize;
1051  NumMyEquations += BlockEntries[RowDim][RowDim].M();
1052 
1053  if (CurRow==0)
1054  {
1055  Indices[0] = CurRow;
1056  Indices[1] = CurRow+1;
1057  NumEntries = 2;
1058  ColDims[0] = ElementSizeList[i] - MinSize;
1059  ColDims[1] = ElementSizeList[i+1] - MinSize; // Assumes linear global ordering and > 1 row/proc.
1060  }
1061  else if (CurRow == NumGlobalElements-1)
1062  {
1063  Indices[0] = CurRow-1;
1064  Indices[1] = CurRow;
1065  NumEntries = 2;
1066  ColDims[0] = ElementSizeList[i-1] - MinSize;
1067  ColDims[1] = ElementSizeList[i] - MinSize; // Assumes linear global ordering and > 1 row/proc.
1068  }
1069  else {
1070  Indices[0] = CurRow-1;
1071  Indices[1] = CurRow;
1072  Indices[2] = CurRow+1;
1073  NumEntries = 3;
1074  if (i==0) ColDims[0] = MaxSize - MinSize; // ElementSize on MyPID-1
1075  else ColDims[0] = ElementSizeList[i-1] - MinSize;
1076  ColDims[1] = ElementSizeList[i] - MinSize;
1077  // ElementSize on MyPID+1
1078  if (i==NumMyElements-1) ColDims[2] = MaxSize - MinSize;
1079  else ColDims[2] = ElementSizeList[i+1] - MinSize;
1080  }
1081  if ( insertlocal ) {
1082  for ( int ii=0; ii < NumEntries; ii++ )
1083  MyIndices[ii] = colmap->LID( Indices[ii] ) ;
1084  EPETRA_TEST_ERR(!(A->BeginSumIntoMyValues( rowmap->LID(CurRow), NumEntries, MyIndices)==0),ierr);
1085  } else {
1086  EPETRA_TEST_ERR(!(A->BeginSumIntoGlobalValues(CurRow, NumEntries, Indices)==0),ierr);
1087  }
1088  forierr = 0;
1089  for (int j=0; j < NumEntries; j++) {
1090  Epetra_SerialDenseMatrix * AD = &(BlockEntries[RowDim][ColDims[j]]);
1091  NumMyNonzeros += AD->M() * AD->N();
1092  // This has nothing to do with insertlocal, but that is a convenient bool to key on
1093  if ( insertlocal ) {
1094  forierr += !(A->SubmitBlockEntry( *AD )==0);
1095  } else {
1096  forierr += !(A->SubmitBlockEntry(AD->A(), AD->LDA(), AD->M(), AD->N())==0);
1097  }
1098  }
1099  EPETRA_TEST_ERR(forierr,ierr);
1100 
1101  A->EndSubmitEntries();
1102  }
1103 
1104  orig_check_y.Scale(2.0) ;
1105 
1106  // This will not work with FixedNumEntries unless we fix the fix the above loop to add the sorner elements to the tridi matrix
1107  if ( ! FixedNumEntries )
1108  EPETRA_TEST_ERR( checkmultiply( false, *A, x, orig_check_y ), ierr );
1109  }
1110 
1111 
1112  {for (int kr=0; kr<SizeRange; kr++) delete [] BlockEntries[kr];}
1113  delete [] BlockEntries;
1114  delete [] ColDims;
1115  delete [] MyIndices ;
1116  delete [] Indices;
1117  delete [] ElementSizeList;
1118 
1119 
1120 
1121 
1122  if (verbose) cout << "\n\n*****Insert methods should not be accepted" << endl<< endl;
1123 
1124  int One = 1;
1125  if (B->MyGRID(0)) EPETRA_TEST_ERR(!(B->BeginInsertGlobalValues(0, 1, &One)==-2),ierr);
1126 
1127  Epetra_Vector checkDiag(B->RowMap());
1128  forierr = 0;
1129  int NumMyEquations1 = B->NumMyRows();
1130  double two1 = 2.0;
1131 
1132  // Test diagonal replacement and extraction methods
1133 
1134  forierr = 0;
1135  for (int i=0; i<NumMyEquations1; i++) checkDiag[i]=two1;
1136  EPETRA_TEST_ERR(forierr,ierr);
1137 
1138  EPETRA_TEST_ERR(!(B->ReplaceDiagonalValues(checkDiag)==0),ierr);
1139 
1140  Epetra_Vector checkDiag1(B->RowMap());
1141  EPETRA_TEST_ERR(!(B->ExtractDiagonalCopy(checkDiag1)==0),ierr);
1142 
1143  forierr = 0;
1144  for (int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
1145  EPETRA_TEST_ERR(forierr,ierr);
1146 
1147  if (verbose) cout << "\n\nDiagonal extraction and replacement OK.\n\n" << endl;
1148 
1149  double originfnorm = B->NormInf();
1150  double origonenorm = B->NormOne();
1151  EPETRA_TEST_ERR(!(B->Scale(4.0)==0),ierr);
1152  // EPETRA_TEST_ERR((B->NormOne()!=origonenorm),ierr);
1153  EPETRA_TEST_ERR(!(B->NormInf()==4.0 * originfnorm),ierr);
1154  EPETRA_TEST_ERR(!(B->NormOne()==4.0 * origonenorm),ierr);
1155 
1156  if (verbose) cout << "\n\nMatrix scale OK.\n\n" << endl;
1157 
1158  // Testing VbrRowMatrix adapter
1159  Epetra_VbrRowMatrix ARowMatrix(A);
1160 
1161  EPETRA_TEST_ERR(checkVbrRowMatrix(*A, ARowMatrix, verbose),ierr);
1162 
1163  if ( PreviousA )
1164  delete *PreviousA;
1165 
1166  //
1167  // The following code deals with the fact that A has to be delete before graph is, when
1168  // A is built with a contructor that takes a graph as input.
1169  //
1170  delete B;
1171  if ( graph ) {
1172  delete A;
1173  delete graph ;
1174  *PreviousA = 0 ;
1175  } else {
1176  *PreviousA = A ;
1177  }
1178 
1179  delete CrsA;
1180  delete CrsB;
1181  delete OrigCrsA ;
1182 
1183  return ierr;
1184 
1185 }
1186 
1187 
1188 int main(int argc, char *argv[])
1189 {
1190  int ierr = 0;
1191  bool debug = false;
1192 
1193 #ifdef EPETRA_MPI
1194  MPI_Init(&argc,&argv);
1195  Epetra_MpiComm Comm( MPI_COMM_WORLD );
1196 #else
1197  Epetra_SerialComm Comm;
1198 #endif
1199 
1200  int MyPID = Comm.MyPID();
1201  int NumProc = Comm.NumProc();
1202  bool verbose = false;
1203 
1204  // Check if we should print results to standard out
1205  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
1206 
1207  // char tmp;
1208  // int rank = Comm.MyPID();
1209  // if (rank==0) cout << "Press any key to continue..."<< endl;
1210  // if (rank==0) cin >> tmp;
1211  // Comm.Barrier();
1212 
1213  // if ( ! verbose )
1214  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
1215 
1216  if (verbose && MyPID==0)
1217  cout << Epetra_Version() << endl << endl;
1218 
1219  if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
1220  << " is alive."<<endl;
1221 
1222  // Redefine verbose to only print on PE 0
1223  if (verbose && Comm.MyPID()!=0) verbose = false;
1224 
1225 
1226  int NumMyElements = 400;
1227  // int NumMyElements = 3;
1228  int MinSize = 2;
1229  int MaxSize = 8;
1230  bool NoExtraBlocks = false;
1231  bool symmetric = true;
1232  bool NonSymmetric = false;
1233  bool NoInsertLocal = false ;
1234  bool InsertLocal = true ;
1235 
1236  Epetra_VbrMatrix* PreviousA = 0 ;
1237 
1238  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 1, 1, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1239 
1240  //
1241  // Check the various constructors
1242  //
1243 
1244  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1245 
1246  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, VariableEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1247 
1248  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1249 
1250  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1251 
1252  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_VEPR, NoExtraBlocks, NoInsertLocal, symmetric, &PreviousA );
1253 
1254  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1255 
1256  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_VEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1257 
1258  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1259 
1260  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, WithGraph, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1261  assert ( PreviousA == 0 ) ;
1262 
1263 
1264  //
1265  // Check some various options
1266  //
1267  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1268 
1269  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1270 
1271  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MaxSize, WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1272  assert ( PreviousA == 0 ) ;
1273 
1274  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1275 
1276  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, RowMapColMap_NEPR, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1277 
1278  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, WithGraph, NoExtraBlocks, InsertLocal, NonSymmetric, &PreviousA );
1279  assert ( PreviousA == 0 ) ;
1280 
1281 
1282 
1283  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, FixedEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1284 
1285  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, MinSize, MinSize, RowMapColMap_FEPR, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1286 
1287 
1288 
1289  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1290 
1291  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 3, 3, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1292 
1293  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 4, 4, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1294 
1295  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 5, 5, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1296 
1297  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 6, 6, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1298 
1299  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 7, 7, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1300 
1301  ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 8, 8, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1302 
1303  // ierr +=TestMatrix( Comm, verbose, debug, NumMyElements, 2, 2, NoEntriesPerRow, NoExtraBlocks, NoInsertLocal, NonSymmetric, &PreviousA );
1304 
1305 
1306  delete PreviousA;
1307 
1308  /*
1309  if (verbose) {
1310  // Test ostream << operator (if verbose)
1311  // Construct a Map that puts 2 equations on each PE
1312 
1313  int NumMyElements1 = 2;
1314  int NumMyElements1 = NumMyElements1;
1315  int NumGlobalElements1 = NumMyElements1*NumProc;
1316 
1317  Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
1318 
1319  // Get update list and number of local equations from newly created Map
1320  int * MyGlobalElements1 = new int[Map1.NumMyElements()];
1321  Map1.MyGlobalElements(MyGlobalElements1);
1322 
1323  // Create an integer vector NumNz that is used to build the Petra Matrix.
1324  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
1325 
1326  int * NumNz1 = new int[NumMyElements1];
1327 
1328  // We are building a tridiagonal matrix where each row has (-1 2 -1)
1329  // So we need 2 off-diagonal terms (except for the first and last equation)
1330 
1331  for (int i=0; i<NumMyElements1; i++)
1332  if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalElements1-1)
1333  NumNz1[i] = 1;
1334  else
1335  NumNz1[i] = 2;
1336 
1337  // Create a Epetra_Matrix
1338 
1339  Epetra_VbrMatrix A1(Copy, Map1, NumNz1);
1340 
1341  // Add rows one-at-a-time
1342  // Need some vectors to help
1343  // Off diagonal Values will always be -1
1344 
1345 
1346  int *Indices1 = new int[2];
1347  double two1 = 2.0;
1348  int NumEntries1;
1349 
1350  forierr = 0;
1351  for (int i=0; i<NumMyElements1; i++)
1352  {
1353  if (MyGlobalElements1[i]==0)
1354  {
1355  Indices1[0] = 1;
1356  NumEntries1 = 1;
1357  }
1358  else if (MyGlobalElements1[i] == NumGlobalElements1-1)
1359  {
1360  Indices1[0] = NumGlobalElements1-2;
1361  NumEntries1 = 1;
1362  }
1363  else
1364  {
1365  Indices1[0] = MyGlobalElements1[i]-1;
1366  Indices1[1] = MyGlobalElements1[i]+1;
1367  NumEntries1 = 2;
1368  }
1369  forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1)==0);
1370  forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i)>0); // Put in the diagonal entry
1371  }
1372  EPETRA_TEST_ERR(forierr,ierr);
1373  // Finish up
1374  EPETRA_TEST_ERR(!(A1.FillComplete()==0),ierr);
1375 
1376  if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
1377  cout << A1 << endl;
1378 
1379  // Release all objects
1380  delete [] NumNz1;
1381  delete [] Values1;
1382  delete [] Indices1;
1383  delete [] MyGlobalElements1;
1384 
1385  }
1386  */
1387 
1388  /* Active Issue: 5744: EPETRA_TEST_ERR( checkVbrMatrixOptimizedGraph(Comm, verbose), ierr); */
1389 
1390  EPETRA_TEST_ERR( checkMergeRedundantEntries(Comm, verbose), ierr);
1391 
1392  EPETRA_TEST_ERR( checkExtractMyRowCopy(Comm, verbose), ierr);
1393 
1394  EPETRA_TEST_ERR( checkMatvecSameVectors(Comm, verbose), ierr);
1395 
1396  EPETRA_TEST_ERR( checkEarlyDelete(Comm, verbose), ierr);
1397 
1398  if (verbose) {
1399  if (ierr==0) cout << "All VbrMatrix tests OK" << endl;
1400  else cout << ierr << " errors encountered." << endl;
1401  }
1402 
1403 #ifdef EPETRA_MPI
1404  MPI_Finalize() ;
1405 #endif
1406 
1407 /* end main
1408 */
1409 return ierr ;
1410 }
1411 
1412 int power_method(bool TransA, Epetra_VbrMatrix& A,
1413  Epetra_MultiVector& q,
1414  Epetra_MultiVector& z,
1415  Epetra_MultiVector& resid,
1416  double * lambda, int niters, double tolerance,
1417  bool verbose) {
1418 
1419  // variable needed for iteration
1420  double normz, residual;
1421 
1422  int ierr = 1;
1423 
1424  for (int iter = 0; iter < niters; iter++)
1425  {
1426  z.Norm2(&normz); // Compute 2-norm of z
1427  q.Scale(1.0/normz, z);
1428  A.Multiply(TransA, q, z); // Compute z = A*q
1429  q.Dot(z, lambda); // Approximate maximum eigenvaluE
1430  if (iter%100==0 || iter+1==niters)
1431  {
1432  resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
1433  resid.Norm2(&residual);
1434  if (verbose) cout << "Iter = " << iter << " Lambda = " << *lambda
1435  << " Residual of A*q - lambda*q = " << residual << endl;
1436  }
1437  if (residual < tolerance) {
1438  ierr = 0;
1439  break;
1440  }
1441  }
1442  return(ierr);
1443 }
1445  int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1,
1446  int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1,
1447  int * MyGlobalElements, bool verbose) {
1448 
1449  int ierr = 0, forierr = 0;
1450  // Test query functions
1451 
1452  int NumMyRows = A.NumMyRows();
1453  if (verbose) cout << "\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
1454  // TEMP
1455  //if (verbose) cout << "\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
1456  if (verbose) cout << "\n\nNumber of local Rows is = " << NumMyRows << endl<< endl;
1457  if (verbose) cout << "\n\nNumber of local Rows should = " << NumMyRows1 << endl<< endl;
1458 
1459  EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);
1460 
1461  int NumMyNonzeros = A.NumMyNonzeros();
1462  if (verbose) cout << "\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
1463  //if (verbose) cout << " Should = " << NumMyNonzeros1 << endl<< endl;
1464 
1465 
1466  if ( NumMyNonzeros != NumMyNonzeros1 ) {
1467  cout << " MyPID = " << A.Comm().MyPID()
1468  << " NumMyNonzeros = " << NumMyNonzeros
1469  << " NumMyNonzeros1 = " << NumMyNonzeros1 << endl ;
1470  }
1471 
1472  EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);
1473 
1474  int NumGlobalRows = A.NumGlobalRows();
1475  if (verbose) cout << "\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
1476 
1477  EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);
1478 
1479  int NumGlobalNonzeros = A.NumGlobalNonzeros();
1480  if (verbose) cout << "\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
1481 
1482  if ( NumGlobalNonzeros != NumGlobalNonzeros1 ) {
1483  cout << " MyPID = " << A.Comm().MyPID()
1484  << " NumGlobalNonzeros = " << NumGlobalNonzeros
1485  << " NumGlobalNonzeros1 = " << NumGlobalNonzeros1 << endl ;
1486  }
1487  EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);
1488 
1489  int NumMyBlockRows = A.NumMyBlockRows();
1490  if (verbose) cout << "\n\nNumber of local Block Rows = " << NumMyBlockRows << endl<< endl;
1491 
1492  EPETRA_TEST_ERR(!(NumMyBlockRows==NumMyBlockRows1),ierr);
1493 
1494  int NumMyBlockNonzeros = A.NumMyBlockEntries();
1495 
1496  if (verbose) cout << "\n\nNumber of local Nonzero Block entries = " << NumMyBlockNonzeros << endl<< endl;
1497  if (verbose) cout << "\n\nNumber of local Nonzero Block entries 1 = " << NumMyBlockNonzeros1 << endl<< endl;
1498 
1499  EPETRA_TEST_ERR(!(NumMyBlockNonzeros==NumMyBlockNonzeros1),ierr);
1500 
1501  int NumGlobalBlockRows = A.NumGlobalBlockRows();
1502  if (verbose) cout << "\n\nNumber of global Block Rows = " << NumGlobalBlockRows << endl<< endl;
1503 
1504  EPETRA_TEST_ERR(!(NumGlobalBlockRows==NumGlobalBlockRows1),ierr);
1505 
1506  int NumGlobalBlockNonzeros = A.NumGlobalBlockEntries();
1507  if (verbose) cout << "\n\nNumber of global Nonzero Block entries = " << NumGlobalBlockNonzeros << endl<< endl;
1508 
1509  EPETRA_TEST_ERR(!(NumGlobalBlockNonzeros==NumGlobalBlockNonzeros1),ierr);
1510 
1511 
1512  // Test RowMatrix interface implementations
1513  int RowDim, NumBlockEntries, * BlockIndices;
1514  Epetra_SerialDenseMatrix ** Values;
1515  // Get View of last block row
1516  A.ExtractMyBlockRowView(NumMyBlockRows-1, RowDim, NumBlockEntries,
1517  BlockIndices, Values);
1518  int NumMyEntries1 = 0;
1519  {for (int i=0; i < NumBlockEntries; i++) NumMyEntries1 += Values[i]->N();}
1520  int NumMyEntries;
1521  A.NumMyRowEntries(NumMyRows-1, NumMyEntries);
1522  if (verbose) {
1523  cout << "\n\nNumber of nonzero values in last row = "
1524  << NumMyEntries << endl<< endl;
1525  }
1526 
1527  EPETRA_TEST_ERR(!(NumMyEntries==NumMyEntries1),ierr);
1528 
1529  // Other binary tests
1530 
1531  EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
1532  EPETRA_TEST_ERR(!(A.Filled()),ierr);
1533  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID())),ierr);
1534  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID())),ierr);
1535  EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID()),ierr);
1536  EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID()),ierr);
1537  EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
1538  EPETRA_TEST_ERR(!(A.MyLRID(NumMyBlockRows-1)),ierr);
1539  EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
1540  EPETRA_TEST_ERR(A.MyLRID(NumMyBlockRows),ierr);
1541 
1542 
1543  int MaxNumBlockEntries = A.MaxNumBlockEntries();
1544 
1545  // Pointer Extraction approach
1546 
1547  // Local index
1548  int MyPointersRowDim, MyPointersNumBlockEntries;
1549  int * MyPointersBlockIndices = new int[MaxNumBlockEntries];
1550  Epetra_SerialDenseMatrix **MyPointersValuesPointers;
1551  // Global Index
1552  int GlobalPointersRowDim, GlobalPointersNumBlockEntries;
1553  int * GlobalPointersBlockIndices = new int[MaxNumBlockEntries];
1554  Epetra_SerialDenseMatrix **GlobalPointersValuesPointers;
1555 
1556  // Copy Extraction approach
1557 
1558  // Local index
1559  int MyCopyRowDim, MyCopyNumBlockEntries;
1560  int * MyCopyBlockIndices = new int[MaxNumBlockEntries];
1561  int * MyCopyColDims = new int[MaxNumBlockEntries];
1562  int * MyCopyLDAs = new int[MaxNumBlockEntries];
1563  int MaxRowDim = A.MaxRowDim();
1564  int MaxColDim = A.MaxColDim();
1565  int MyCopySizeOfValues = MaxRowDim*MaxColDim;
1566  double ** MyCopyValuesPointers = new double*[MaxNumBlockEntries];
1567  for (int i=0; i<MaxNumBlockEntries; i++)
1568  MyCopyValuesPointers[i] = new double[MaxRowDim*MaxColDim];
1569 
1570  // Global Index
1571  int GlobalCopyRowDim, GlobalCopyNumBlockEntries;
1572  int * GlobalCopyBlockIndices = new int[MaxNumBlockEntries];
1573  int * GlobalCopyColDims = new int[MaxNumBlockEntries];
1574  int * GlobalCopyLDAs = new int[MaxNumBlockEntries];
1575 
1576  int GlobalMaxRowDim = A.GlobalMaxRowDim();
1577  int GlobalMaxColDim = A.GlobalMaxColDim();
1578  int GlobalCopySizeOfValues = GlobalMaxRowDim*GlobalMaxColDim;
1579  double ** GlobalCopyValuesPointers = new double*[MaxNumBlockEntries];
1580  for (int i=0; i<MaxNumBlockEntries; i++)
1581  GlobalCopyValuesPointers[i] = new double[GlobalMaxRowDim*GlobalMaxColDim];
1582 
1583  // View Extraction approaches
1584 
1585  // Local index (There is no global view available)
1586  int MyView1RowDim, MyView1NumBlockEntries;
1587  int * MyView1BlockIndices;
1588  Epetra_SerialDenseMatrix **MyView1ValuesPointers = new Epetra_SerialDenseMatrix*[MaxNumBlockEntries];
1589 
1590  // Local index version 2 (There is no global view available)
1591  int MyView2RowDim, MyView2NumBlockEntries;
1592  int * MyView2BlockIndices;
1593  Epetra_SerialDenseMatrix **MyView2ValuesPointers;
1594 
1595 
1596  // For each row, test six approaches to extracting data from a given local index matrix
1597  forierr = 0;
1598  for (int i=0; i<NumMyBlockRows; i++) {
1599  int MyRow = i;
1600  int GlobalRow = A.GRID(i);
1601  // Get a copy of block indices in local index space, pointers to everything else
1602  EPETRA_TEST_ERR( A.ExtractMyBlockRowPointers(MyRow, MaxNumBlockEntries, MyPointersRowDim,
1603  MyPointersNumBlockEntries, MyPointersBlockIndices,
1604  MyPointersValuesPointers), ierr );
1605  // Get a copy of block indices in local index space, pointers to everything else
1606  EPETRA_TEST_ERR( A.ExtractGlobalBlockRowPointers(GlobalRow, MaxNumBlockEntries, GlobalPointersRowDim,
1607  GlobalPointersNumBlockEntries, GlobalPointersBlockIndices,
1608  GlobalPointersValuesPointers), ierr ) ;
1609 
1610  // Initiate a copy of block row in local index space.
1611  EPETRA_TEST_ERR( A.BeginExtractMyBlockRowCopy(MyRow, MaxNumBlockEntries, MyCopyRowDim,
1612  MyCopyNumBlockEntries, MyCopyBlockIndices,
1613  MyCopyColDims), ierr);
1614  // Copy Values
1615  for (int j=0; j<MyCopyNumBlockEntries; j++) {
1616  EPETRA_TEST_ERR( A.ExtractEntryCopy(MyCopySizeOfValues, MyCopyValuesPointers[j], MaxRowDim, false), ierr) ;
1617  MyCopyLDAs[j] = MaxRowDim;
1618  }
1619 
1620  // Initiate a copy of block row in global index space.
1621  EPETRA_TEST_ERR( A.BeginExtractGlobalBlockRowCopy(GlobalRow, MaxNumBlockEntries, GlobalCopyRowDim,
1622  GlobalCopyNumBlockEntries, GlobalCopyBlockIndices,
1623  GlobalCopyColDims), ierr ) ;
1624  // Copy Values
1625  for (int j=0; j<GlobalCopyNumBlockEntries; j++) {
1626  EPETRA_TEST_ERR( A.ExtractEntryCopy(GlobalCopySizeOfValues, GlobalCopyValuesPointers[j], GlobalMaxRowDim, false), ierr );
1627  GlobalCopyLDAs[j] = GlobalMaxRowDim;
1628  }
1629 
1630  // Initiate a view of block row in local index space (Version 1)
1631  EPETRA_TEST_ERR( A.BeginExtractMyBlockRowView(MyRow, MyView1RowDim,
1632  MyView1NumBlockEntries, MyView1BlockIndices), ierr ) ;
1633  // Set pointers to values
1634  for (int j=0; j<MyView1NumBlockEntries; j++)
1635  EPETRA_TEST_ERR ( A.ExtractEntryView(MyView1ValuesPointers[j]), ierr ) ;
1636 
1637 
1638  // Extract a view of block row in local index space (version 2)
1639  EPETRA_TEST_ERR( A.ExtractMyBlockRowView(MyRow, MyView2RowDim,
1640  MyView2NumBlockEntries, MyView2BlockIndices,
1641  MyView2ValuesPointers), ierr );
1642 
1643  forierr += !(MyPointersNumBlockEntries==GlobalPointersNumBlockEntries);
1644  forierr += !(MyPointersNumBlockEntries==MyCopyNumBlockEntries);
1645  forierr += !(MyPointersNumBlockEntries==GlobalCopyNumBlockEntries);
1646  forierr += !(MyPointersNumBlockEntries==MyView1NumBlockEntries);
1647  forierr += !(MyPointersNumBlockEntries==MyView2NumBlockEntries);
1648  for (int j=1; j<MyPointersNumBlockEntries; j++) {
1649  forierr += !(MyCopyBlockIndices[j-1]<MyCopyBlockIndices[j]);
1650  forierr += !(MyView1BlockIndices[j-1]<MyView1BlockIndices[j]);
1651  forierr += !(MyView2BlockIndices[j-1]<MyView2BlockIndices[j]);
1652 
1653  forierr += !(GlobalPointersBlockIndices[j]==A.GCID(MyPointersBlockIndices[j]));
1654  forierr += !(A.LCID(GlobalPointersBlockIndices[j])==MyPointersBlockIndices[j]);
1655  forierr += !(GlobalPointersBlockIndices[j]==GlobalCopyBlockIndices[j]);
1656 
1657  Epetra_SerialDenseMatrix* my = MyPointersValuesPointers[j];
1658  Epetra_SerialDenseMatrix* global = GlobalPointersValuesPointers[j];
1659 
1660  Epetra_SerialDenseMatrix* myview1 = MyView1ValuesPointers[j];
1661  Epetra_SerialDenseMatrix* myview2 = MyView2ValuesPointers[j];
1662 
1663  forierr += !(CompareValues(my->A(), my->LDA(),
1664  MyPointersRowDim, my->N(),
1665  global->A(), global->LDA(),
1666  GlobalPointersRowDim, global->N())==0);
1667  forierr += !(CompareValues(my->A(), my->LDA(),
1668  MyPointersRowDim, my->N(),
1669  MyCopyValuesPointers[j], MyCopyLDAs[j],
1670  MyCopyRowDim, MyCopyColDims[j])==0);
1671  forierr += !(CompareValues(my->A(), my->LDA(),
1672  MyPointersRowDim, my->N(),
1673  GlobalCopyValuesPointers[j], GlobalCopyLDAs[j],
1674  GlobalCopyRowDim, GlobalCopyColDims[j])==0);
1675  forierr += !(CompareValues(my->A(), my->LDA(),
1676  MyPointersRowDim, my->N(),
1677  myview1->A(), myview1->LDA(),
1678  MyView1RowDim, myview1->N())==0);
1679  forierr += !(CompareValues(my->A(), my->LDA(),
1680  MyPointersRowDim, my->N(),
1681  myview2->A(), myview2->LDA(),
1682  MyView2RowDim, myview2->N())==0);
1683  }
1684  }
1685  EPETRA_TEST_ERR(forierr,ierr);
1686 
1687  // GlobalRowView should be illegal (since we have local indices)
1688  EPETRA_TEST_ERR(!(A.BeginExtractGlobalBlockRowView(A.GRID(0), MyView1RowDim,
1689  MyView1NumBlockEntries,
1690  MyView1BlockIndices)==-2),ierr);
1691 
1692  // Extract a view of block row in local index space (version 2)
1693  EPETRA_TEST_ERR(!(A.ExtractGlobalBlockRowView(A.GRID(0), MyView2RowDim,
1694  MyView2NumBlockEntries, MyView2BlockIndices,
1695  MyView2ValuesPointers)==-2),ierr);
1696 
1697  delete [] MyPointersBlockIndices;
1698  delete [] GlobalPointersBlockIndices;
1699  delete [] MyCopyBlockIndices;
1700  delete [] MyCopyColDims;
1701  delete [] MyCopyLDAs;
1702  for (int i=0; i<MaxNumBlockEntries; i++) delete [] MyCopyValuesPointers[i];
1703  delete [] MyCopyValuesPointers;
1704  delete [] GlobalCopyBlockIndices;
1705  delete [] GlobalCopyColDims;
1706  delete [] GlobalCopyLDAs;
1707  for (int i=0; i<MaxNumBlockEntries; i++) delete [] GlobalCopyValuesPointers[i];
1708  delete [] GlobalCopyValuesPointers;
1709  delete [] MyView1ValuesPointers;
1710  if (verbose) cout << "\n\nRows sorted check OK" << endl<< endl;
1711 
1712  return ierr;
1713 }
1714 
1715 //=============================================================================
1716 int CompareValues(double * A, int LDA, int NumRowsA, int NumColsA,
1717  double * B, int LDB, int NumRowsB, int NumColsB) {
1718 
1719  int ierr = 0, forierr = 0;
1720  double * ptr1 = B;
1721  double * ptr2;
1722 
1723  if (NumRowsA!=NumRowsB) EPETRA_TEST_ERR(-2,ierr);
1724  if (NumColsA!=NumColsB) EPETRA_TEST_ERR(-3,ierr);
1725 
1726 
1727  forierr = 0;
1728  for (int j=0; j<NumColsA; j++) {
1729  ptr1 = B + j*LDB;
1730  ptr2 = A + j*LDA;
1731  for (int i=0; i<NumRowsA; i++) forierr += (*ptr1++ != *ptr2++);
1732  }
1733  EPETRA_TEST_ERR(forierr,ierr);
1734  return ierr;
1735 }
1736 
1737 int checkMergeRedundantEntries(Epetra_Comm& comm, bool verbose)
1738 {
1739  int numProcs = comm.NumProc();
1740  int localProc = comm.MyPID();
1741 
1742  int myFirstRow = localProc*3;
1743  int myLastRow = myFirstRow+2;
1744  int numMyRows = myLastRow - myFirstRow + 1;
1745  int numGlobalRows = numProcs*numMyRows;
1746  int ierr;
1747 
1748  //We'll set up a matrix which is globally block-diagonal, i.e., on each
1749  //processor the list of columns == list of rows.
1750  //Set up a list of column-indices which is twice as long as it should be,
1751  //and its contents will be the list of local rows, repeated twice.
1752  int numCols = 2*numMyRows;
1753  int* myCols = new int[numCols];
1754 
1755  int col = myFirstRow;
1756  for(int i=0; i<numCols; ++i) {
1757  myCols[i] = col++;
1758  if (col > myLastRow) col = myFirstRow;
1759  }
1760 
1761  int elemSize = 2;
1762  int indexBase = 0;
1763 
1764  Epetra_BlockMap map(numGlobalRows, numMyRows,
1765  elemSize, indexBase, comm);
1766 
1767  Epetra_VbrMatrix A(Copy, map, numCols);
1768 
1769  Epetra_MultiVector x(map, 1), y(map, 1);
1770  x.PutScalar(1.0);
1771 
1772  Epetra_MultiVector x3(map, 3), y3(map, 3);
1773  x.PutScalar(1.0);
1774 
1775  double* coef = new double[elemSize*elemSize];
1776  for(int i=0; i<elemSize*elemSize; ++i) {
1777  coef[i] = 0.5;
1778  }
1779 
1780  //we're going to insert each row twice, with coef values of 0.5. So after
1781  //FillComplete, which internally calls MergeRedundantEntries, the
1782  //matrix should contain 1.0 in each entry.
1783 
1784  for(int i=myFirstRow; i<=myLastRow; ++i) {
1785  EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, numCols, myCols), ierr);
1786 
1787  for(int j=0; j<numCols; ++j) {
1788  EPETRA_TEST_ERR( A.SubmitBlockEntry(coef, elemSize,
1789  elemSize, elemSize), ierr);
1790  }
1791 
1792  EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
1793  }
1794 
1795  EPETRA_TEST_ERR( A.FillComplete(), ierr);
1796 
1797  delete [] coef;
1798 
1799  if (verbose) cout << "Multiply x"<<endl;
1800  EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr );
1801 
1802 
1803  //Next we're going to extract pointers-to-block-rows and check values to make
1804  //sure that the internal method Epetra_VbrMatrix::mergeRedundantEntries()
1805  //worked correctly.
1806  //At the same time, we're also going to create another VbrMatrix which will
1807  //be a View of the matrix we've already created. This will serve to provide
1808  //more test coverage of the VbrMatrix code.
1809 
1810  int numBlockEntries = 0;
1811  int RowDim;
1812  int** BlockIndices = new int*[numMyRows];
1813  Epetra_SerialDenseMatrix** Values;
1814  Epetra_VbrMatrix Aview(View, map, numMyRows);
1815 
1816  for(int i=myFirstRow; i<=myLastRow; ++i) {
1817  BlockIndices[i-myFirstRow] = new int[numCols];
1819  RowDim, numBlockEntries,
1820  BlockIndices[i-myFirstRow],
1821  Values), ierr);
1822 
1823  EPETRA_TEST_ERR( Aview.BeginInsertGlobalValues(i, numBlockEntries,
1824  BlockIndices[i-myFirstRow]), ierr);
1825 
1826  if (numMyRows != numBlockEntries) return(-1);
1827  if (RowDim != elemSize) return(-2);
1828  for(int j=0; j<numBlockEntries; ++j) {
1829  if (Values[j]->A()[0] != 1.0) {
1830  cout << "Row " << i << " Values["<<j<<"][0]: "<< Values[j][0]
1831  << " should be 1.0" << endl;
1832  return(-3); //comment-out this return to de-activate this test
1833  }
1834 
1835  EPETRA_TEST_ERR( Aview.SubmitBlockEntry(Values[j]->A(),
1836  Values[j]->LDA(),
1837  Values[j]->M(),
1838  Values[j]->N()), ierr);
1839  }
1840 
1841  EPETRA_TEST_ERR( Aview.EndSubmitEntries(), ierr);
1842  }
1843 
1844  EPETRA_TEST_ERR( Aview.FillComplete(), ierr);
1845 
1846  //So the test appears to have passed for the original matrix A. Now check the
1847  //values of our second "view" of the matrix, 'Aview'.
1848 
1849  for(int i=myFirstRow; i<=myLastRow; ++i) {
1850  EPETRA_TEST_ERR( Aview.ExtractGlobalBlockRowPointers(i, numMyRows,
1851  RowDim, numBlockEntries,
1852  BlockIndices[i-myFirstRow],
1853  Values), ierr);
1854 
1855  if (numMyRows != numBlockEntries) return(-1);
1856  if (RowDim != elemSize) return(-2);
1857  for(int j=0; j<numBlockEntries; ++j) {
1858  if (Values[j]->A()[0] != 1.0) {
1859  cout << "Aview: Row " << i << " Values["<<j<<"][0]: "<< Values[j][0]
1860  << " should be 1.0" << endl;
1861  return(-3); //comment-out this return to de-activate this test
1862  }
1863  }
1864 
1865  delete [] BlockIndices[i-myFirstRow];
1866  }
1867 
1868 
1869  if (verbose&&localProc==0) {
1870  cout << "checkMergeRedundantEntries, A:" << endl;
1871  }
1872 
1873 
1874  delete [] BlockIndices;
1875  delete [] myCols;
1876 
1877  return(0);
1878 }
1879 
1880 int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose)
1881 {
1882  int numProcs = comm.NumProc();
1883  int localProc = comm.MyPID();
1884 
1885  int myFirstRow = localProc*3;
1886  int myLastRow = myFirstRow+2;
1887  int numMyRows = myLastRow - myFirstRow + 1;
1888  int numGlobalRows = numProcs*numMyRows;
1889  int ierr;
1890 
1891  int numCols = numMyRows;
1892  int* myCols = new int[numCols];
1893 
1894  int col = myFirstRow;
1895  for(int i=0; i<numCols; ++i) {
1896  myCols[i] = col++;
1897  if (col > myLastRow) col = myFirstRow;
1898  }
1899 
1900  int elemSize = 2;
1901  int indexBase = 0;
1902 
1903  Epetra_BlockMap map(numGlobalRows, numMyRows,
1904  elemSize, indexBase, comm);
1905 
1906  Epetra_VbrMatrix A(Copy, map, numCols);
1907 
1908  double* coef = new double[elemSize*elemSize];
1909 
1910  for(int i=myFirstRow; i<=myLastRow; ++i) {
1911  int myPointRow = i*elemSize;
1912 
1913  //The coefficients need to be laid out in column-major order. i.e., the
1914  //coefficients in a column occur contiguously.
1915  for(int ii=0; ii<elemSize; ++ii) {
1916  for(int jj=0; jj<elemSize; ++jj) {
1917  double val = (myPointRow+ii)*1.0;
1918  coef[ii+elemSize*jj] = val;
1919  }
1920  }
1921 
1922  EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, numCols, myCols), ierr);
1923 
1924  for(int j=0; j<numCols; ++j) {
1925  EPETRA_TEST_ERR( A.SubmitBlockEntry(coef, elemSize,
1926  elemSize, elemSize), ierr);
1927  }
1928 
1929  EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
1930  }
1931 
1932  EPETRA_TEST_ERR( A.FillComplete(), ierr);
1933 
1934  delete [] coef;
1935  delete [] myCols;
1936 
1937  Epetra_SerialDenseMatrix** blockEntries;
1938  int len = elemSize*numCols, checkLen;
1939  double* values = new double[len];
1940  int* indices = new int[len];
1941  int RowDim, numBlockEntries;
1942 
1943  for(int i=myFirstRow; i<=myLastRow; ++i) {
1945  RowDim, numBlockEntries,
1946  indices,
1947  blockEntries), ierr);
1948  if (numMyRows != numBlockEntries) return(-1);
1949  if (RowDim != elemSize) return(-2);
1950 
1951  int myPointRow = i*elemSize - myFirstRow*elemSize;
1952  for(int ii=0; ii<elemSize; ++ii) {
1953  EPETRA_TEST_ERR( A.ExtractMyRowCopy(myPointRow+ii, len,
1954  checkLen, values, indices), ierr);
1955  if (len != checkLen) return(-3);
1956 
1957  double val = (i*elemSize+ii)*1.0;
1958  double blockvalue = blockEntries[0]->A()[ii];
1959 
1960  for(int jj=0; jj<len; ++jj) {
1961  if (values[jj] != val) return(-4);
1962  if (values[jj] != blockvalue) return(-5);
1963  }
1964  }
1965  }
1966 
1967  delete [] values;
1968  delete [] indices;
1969 
1970  return(0);
1971 }
1972 
1973 int checkMatvecSameVectors(Epetra_Comm& comm, bool verbose)
1974 {
1975  int numProcs = comm.NumProc();
1976  int localProc = comm.MyPID();
1977 
1978  int myFirstRow = localProc*3;
1979  int myLastRow = myFirstRow+2;
1980  int numMyRows = myLastRow - myFirstRow + 1;
1981  int numGlobalRows = numProcs*numMyRows;
1982  int ierr;
1983 
1984  int elemSize = 2;
1985  int num_off_diagonals = 1;
1986 
1987  epetra_test::matrix_data matdata(numGlobalRows, numGlobalRows,
1988  num_off_diagonals, elemSize);
1989 
1990  Epetra_BlockMap map(numGlobalRows, numMyRows, elemSize, 0, comm);
1991 
1992  Epetra_VbrMatrix A(Copy, map, num_off_diagonals*2+1);
1993 
1994  int* rowlengths = matdata.rowlengths();
1995  int** colindices = matdata.colindices();
1996 
1997  for(int i=myFirstRow; i<=myLastRow; ++i) {
1998 
1999  EPETRA_TEST_ERR( A.BeginInsertGlobalValues(i, rowlengths[i],
2000  colindices[i]), ierr);
2001 
2002  for(int j=0; j<rowlengths[i]; ++j) {
2003  EPETRA_TEST_ERR( A.SubmitBlockEntry(matdata.coefs(i,colindices[i][j]), elemSize,
2004  elemSize, elemSize), ierr);
2005  }
2006 
2007  EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
2008  }
2009 
2010  EPETRA_TEST_ERR( A.FillComplete(), ierr);
2011 
2012  Epetra_Vector x(map), y(map);
2013 
2014  x.PutScalar(1.0);
2015 
2016  A.Multiply(false, x, y);
2017  A.Multiply(false, x, x);
2018 
2019  double* xptr = x.Values();
2020  double* yptr = y.Values();
2021 
2022  for(int i=0; i<numMyRows; ++i) {
2023  if (xptr[i] != yptr[i]) {
2024  return(-1);
2025  }
2026  }
2027 
2028  return(0);
2029 }
2030 
2031 int checkEarlyDelete(Epetra_Comm& comm, bool verbose)
2032 {
2033  int localProc = comm.MyPID();
2034  int numProcs = comm.NumProc();
2035  int myFirstRow = localProc*3;
2036  int myLastRow = myFirstRow+2;
2037  int numMyRows = myLastRow - myFirstRow + 1;
2038  int numGlobalRows = numProcs*numMyRows;
2039  int ierr;
2040 
2041  int elemSize = 2;
2042  int num_off_diagonals = 1;
2043 
2044  epetra_test::matrix_data matdata(numGlobalRows, numGlobalRows,
2045  num_off_diagonals, elemSize);
2046 
2047  Epetra_BlockMap map(numGlobalRows, numMyRows, elemSize, 0, comm);
2048 
2049  Epetra_VbrMatrix* A = new Epetra_VbrMatrix(Copy, map, num_off_diagonals*2+1);
2050 
2051  int* rowlengths = matdata.rowlengths();
2052  int** colindices = matdata.colindices();
2053 
2054  for(int i=myFirstRow; i<=myLastRow; ++i) {
2055 
2056  EPETRA_TEST_ERR( A->BeginInsertGlobalValues(i, rowlengths[i],
2057  colindices[i]), ierr);
2058 
2059  for(int j=0; j<rowlengths[i]; ++j) {
2060  EPETRA_TEST_ERR( A->SubmitBlockEntry(matdata.coefs(i,colindices[i][j]),
2061  elemSize, elemSize, elemSize), ierr);
2062  }
2063 
2064  EPETRA_TEST_ERR( A->EndSubmitEntries(), ierr);
2065  }
2066 
2067  //A call to BeginReplaceMyValues should produce an error at this
2068  //point, since IndicesAreLocal should be false.
2069  int errcode = A->BeginReplaceMyValues(0, 0, 0);
2070  if (errcode == 0) EPETRA_TEST_ERR(-1, ierr);
2071 
2072  EPETRA_TEST_ERR( A->FillComplete(), ierr);
2073 
2074  Epetra_VbrMatrix B(Copy, A->Graph());
2075 
2076  delete A;
2077 
2078  for(int i=myFirstRow; i<=myLastRow; ++i) {
2079 
2080  EPETRA_TEST_ERR( B.BeginReplaceGlobalValues(i, rowlengths[i],
2081  colindices[i]), ierr);
2082 
2083  for(int j=0; j<rowlengths[i]; ++j) {
2084  EPETRA_TEST_ERR( B.SubmitBlockEntry(matdata.coefs(i,colindices[i][j]),
2085  elemSize, elemSize, elemSize), ierr);
2086  }
2087 
2088  EPETRA_TEST_ERR( B.EndSubmitEntries(), ierr);
2089  }
2090 
2091  EPETRA_TEST_ERR( B.FillComplete(), ierr);
2092 
2093  return(0);
2094 }
2095 
2097 {
2098  using std::vector;
2099 
2100  int ierr = 0;
2101 
2102  const int node = comm.MyPID();
2103  const int nodes = comm.NumProc();
2104 
2105  int Ni = 4;
2106  int Nj = 4;
2107  int Gi = 4;
2108  // int Gj = 4;
2109  int i_off = 0;
2110  int j_off = 0;
2111 
2112  int first = 0;
2113  int last = 3;
2114 
2115  Epetra_BlockMap* map;
2116  if (nodes == 1)
2117  {
2118  map = new Epetra_BlockMap(-1, 16, 3, 0, comm);
2119  }
2120  else if (nodes == 2)
2121  {
2122  Ni = 2;
2123  vector<int> l2g(8);
2124  if (node == 1)
2125  {
2126  i_off = 2;
2127  }
2128  for (int j = 0; j < 4; ++j)
2129  {
2130  for (int i = 0; i < 2; ++i)
2131  {
2132  l2g[i + 2 * j] = (i + i_off) + (j + j_off) * Gi;
2133  }
2134  }
2135  map = new Epetra_BlockMap(-1, Ni*Nj, &l2g[0], 3, 0, comm);
2136  }
2137  else if (nodes == 4)
2138  {
2139  Ni = 2;
2140  Nj = 2;
2141  vector<int> l2g(4);
2142  if (node == 1)
2143  {
2144  i_off = 2;
2145  }
2146  else if (node == 2)
2147  {
2148  j_off = 2;
2149  }
2150  else if (node == 3)
2151  {
2152  i_off = 2;
2153  j_off = 2;
2154  }
2155  for (int j = 0; j < 2; ++j)
2156  {
2157  for (int i = 0; i < 2; ++i)
2158  {
2159  l2g[i + 2 * j] = (i + i_off) + (j + j_off) * Gi;
2160  }
2161  }
2162  map = new Epetra_BlockMap(-1, Ni*Nj, &l2g[0], 3, 0, comm);
2163  }
2164  else {
2165  map = NULL;
2166  return 0;
2167  }
2168 
2169  // graph
2170  Epetra_CrsGraph *graph = new Epetra_CrsGraph(Copy, *map, 5);
2171  int indices[5];
2172 
2173  for (int j = 0; j < Nj; ++j)
2174  {
2175  for (int i = 0; i < Ni; ++i)
2176  {
2177  int ctr = 0;
2178  int gi = i + i_off;
2179  int gj = j + j_off;
2180  indices[ctr++] = gi + gj * Gi;
2181  if (gi > first)
2182  indices[ctr++] = (gi - 1) + gj * Gi;
2183  if (gi < last)
2184  indices[ctr++] = (gi + 1) + gj * Gi;
2185  if (gj > first)
2186  indices[ctr++] = gi + (gj - 1) * Gi;
2187  if (gj < last)
2188  indices[ctr++] = gi + (gj + 1) * Gi;
2189  EPETRA_TEST_ERR( ! (ctr <= 5), ierr );
2190  // assign the indices to the graph
2191  graph->InsertGlobalIndices(indices[0], ctr, &indices[0]);
2192  }
2193  }
2194 
2195  // complete the graph
2196  int result = graph->FillComplete();
2197  EPETRA_TEST_ERR( ! result == 0, ierr );
2198  EPETRA_TEST_ERR( ! graph->Filled(), ierr );
2199  EPETRA_TEST_ERR( graph->StorageOptimized(), ierr );
2200  EPETRA_TEST_ERR( ! graph->IndicesAreLocal(), ierr );
2201  result = graph->OptimizeStorage();
2202  EPETRA_TEST_ERR( ! result == 0, ierr );
2203  EPETRA_TEST_ERR( ! graph->StorageOptimized(), ierr );
2204 
2205  EPETRA_TEST_ERR( ! Ni*Nj == graph->NumMyBlockRows(), ierr );
2206  EPETRA_TEST_ERR( ! Ni*Nj*3 == graph->NumMyRows(), ierr );
2207 
2208  Epetra_VbrMatrix *matrix = new Epetra_VbrMatrix(Copy, *graph);
2209  EPETRA_TEST_ERR( ! matrix->IndicesAreLocal(), ierr );
2210 
2212  Epetra_SerialDenseMatrix L(3, 3);
2213  Epetra_SerialDenseMatrix R(3, 3);
2214  EPETRA_TEST_ERR( ! 3 == C.LDA(), ierr );
2215  EPETRA_TEST_ERR( ! 3 == L.LDA(), ierr );
2216  EPETRA_TEST_ERR( ! 3 == R.LDA(), ierr );
2217  std::fill(C.A(), C.A()+9, -4.0);
2218  std::fill(L.A(), L.A()+9, 2.0);
2219  std::fill(R.A(), R.A()+9, 2.0);
2220 
2221  // fill matrix
2222  {
2223  for (int j = 0; j < Nj; ++j)
2224  {
2225  for (int i = 0; i < Ni; ++i)
2226  {
2227  int ctr = 0;
2228  int gi = i + i_off;
2229  int gj = j + j_off;
2230 
2231  int local = i + j * Ni;
2232  int global = gi + gj * Gi;
2233 
2234  int left = (gi - 1) + gj * Gi;
2235  int right = (gi + 1) + gj * Gi;
2236  int bottom = gi + (gj - 1) * Gi;
2237  int top = gi + (gj + 1) * Gi;
2238 
2239  EPETRA_TEST_ERR( ! local == matrix->LCID(global), ierr );
2240 
2241  indices[ctr++] = local;
2242  if (gi > first)
2243  indices[ctr++] = left;
2244  if (gi < last)
2245  indices[ctr++] = right;;
2246  if (gj > first)
2247  indices[ctr++] = bottom;
2248  if (gj < last)
2249  indices[ctr++] = top;
2250 
2251  matrix->BeginReplaceMyValues(local, ctr, &indices[0]);
2252  matrix->SubmitBlockEntry(C);
2253  if (gi > first) matrix->SubmitBlockEntry(L);
2254  if (gi < last) matrix->SubmitBlockEntry(R);
2255  if (gj > first) matrix->SubmitBlockEntry(L);
2256  if (gj < last) matrix->SubmitBlockEntry(R);
2257  matrix->EndSubmitEntries();
2258  }
2259  }
2260  }
2261  matrix->FillComplete();
2262  EPETRA_TEST_ERR( ! matrix->Filled(), ierr );
2263  EPETRA_TEST_ERR( ! matrix->StorageOptimized(), ierr );
2264  EPETRA_TEST_ERR( ! matrix->IndicesAreContiguous(), ierr );
2265  // EPETRA_TEST_ERR( matrix->StorageOptimized(), ierr );
2266  // EPETRA_TEST_ERR( matrix->IndicesAreContiguous(), ierr );
2267 
2268  delete matrix;
2269  delete graph;
2270  delete map;
2271 
2272  return(0);
2273 }
2274 
2275 
2277 
2278  if (verbose) cout << "Checking VbrRowMatrix Adapter..." << endl;
2279  int ierr = 0;
2280  EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
2281  EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
2282  EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
2283  EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
2284  // EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
2285  // EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
2291  EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
2294  for (int i=0; i<A.NumMyRows(); i++) {
2295  int nA, nB;
2296  A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
2297  EPETRA_TEST_ERR(!nA==nB,ierr);
2298  }
2299  EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
2304  // EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
2305  EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);
2306 
2307  int NumVectors = 5;
2308  { // No transpose case
2309  Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
2310  Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
2311  Epetra_MultiVector YA2(YA1);
2312  Epetra_MultiVector YB1(YA1);
2313  Epetra_MultiVector YB2(YA1);
2314  X.Random();
2315 
2316  bool transA = false;
2317  A.SetUseTranspose(transA);
2318  B.SetUseTranspose(transA);
2319  A.Apply(X,YA1);
2320  A.Multiply(transA, X, YA2);
2321  EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
2322  B.Apply(X,YB1);
2323  EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
2324  B.Multiply(transA, X, YB2);
2325  EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);
2326 
2327  }
2328  {// transpose case
2329  Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
2330  Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
2331  Epetra_MultiVector YA2(YA1);
2332  Epetra_MultiVector YB1(YA1);
2333  Epetra_MultiVector YB2(YA1);
2334  X.Random();
2335 
2336  bool transA = true;
2337  A.SetUseTranspose(transA);
2338  B.SetUseTranspose(transA);
2339  A.Apply(X,YA1);
2340  A.Multiply(transA, X, YA2);
2341  EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
2342  B.Apply(X,YB1);
2343  EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
2344  B.Multiply(transA, X,YB2);
2345  EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);
2346 
2347  }
2348 
2349  /*
2350  Epetra_Vector diagA(A.RowMatrixRowMap());
2351  EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
2352  Epetra_Vector diagB(B.RowMatrixRowMap());
2353  EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
2354  EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);
2355  */
2356  Epetra_Vector rowA(A.RowMatrixRowMap());
2357  EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
2358  Epetra_Vector rowB(B.RowMatrixRowMap());
2359  EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
2360  EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);
2361 
2362  Epetra_Vector colA(A.OperatorDomainMap());
2363  EPETRA_TEST_ERR(A.InvColSums(colA),ierr); // -2 error code
2364  Epetra_Vector colB(B.OperatorDomainMap());
2365  EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
2366  EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr); // 1 error code
2367 
2368  EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
2369  EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);
2370 
2371  EPETRA_TEST_ERR(A.RightScale(colA),ierr); // -3 error code
2372  EPETRA_TEST_ERR(B.RightScale(colB),ierr); // -3 error code
2373 
2374 
2375  EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
2376  EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);
2377 
2378 
2379  EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf after scaling", verbose), ierr);
2380  EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne after scaling", verbose),ierr);
2381 
2382  vector<double> valuesA(A.MaxNumEntries());
2383  vector<int> indicesA(A.MaxNumEntries());
2384  vector<double> valuesB(B.MaxNumEntries());
2385  vector<int> indicesB(B.MaxNumEntries());
2386  //return(0);
2387  for (int i=0; i<A.NumMyRows(); i++) {
2388  int nA, nB;
2389  EPETRA_TEST_ERR(A.ExtractMyRowCopy(i, A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr);
2390  EPETRA_TEST_ERR(B.ExtractMyRowCopy(i, B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
2391  EPETRA_TEST_ERR(!nA==nB,ierr);
2392  for (int j=0; j<nA; j++) {
2393  double curVal = valuesA[j];
2394  int curIndex = indicesA[j];
2395  bool notfound = true;
2396  int jj = 0;
2397  while (notfound && jj< nB) {
2398  if (!checkValues(curVal, valuesB[jj])) notfound = false;
2399  jj++;
2400  }
2401  EPETRA_TEST_ERR(notfound, ierr);
2402  vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex); // find curIndex in indicesB
2403  EPETRA_TEST_ERR(p==indicesB.end(), ierr);
2404  }
2405 
2406  }
2407 
2408  if (verbose) {
2409  if (ierr==0)
2410  cout << "RowMatrix Methods check OK" << endl;
2411  else
2412  cout << "ierr = " << ierr << ". RowMatrix Methods check failed." << endl;
2413  }
2414  return (ierr);
2415 }
2416 
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
int NumGlobalElements() const
Number of elements across all processors.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int Random()
Set multi-vector values to random numbers.
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_MultiVector X in Y.
int NumMyBlockRows() const
Returns the number of block matrix rows on this processor.
virtual int RightScale(const Epetra_Vector &x)=0
Scales the Epetra_RowMatrix on the right with a Epetra_Vector x.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
const Epetra_Map & RowMatrixColMap() const
Returns the Epetra_Map object associated with columns of this matrix.
int Random()
Set matrix values to random numbers.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, transpose of this operator will be applied.
virtual double NormOne() const =0
Returns the one norm of the global matrix.
int BeginReplaceGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate replacement of current values with this list of entries for a given global row of the matrix...
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
int BeginInsertMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given local row of the matrix, values are inserted via ...
int Multiply1(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_Vector x in y.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
const Epetra_CrsGraph & Graph() const
Returns a pointer to the Epetra_CrsGraph object associated with this matrix.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
int ExtractEntryCopy(int SizeOfValues, double *Values, int LDA, bool SumInto) const
Extract a copy of an entry from the block row specified by one of the BeginExtract routines...
#define EPETRA_TEST_ERR(a, b)
int PutScalar(double ScalarConstant)
Initialize all values in graph of the matrix with constant value.
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
double NormInf() const
Returns the infinity norm of the global matrix.
int CompareValues(double *A, int LDA, int NumRowsA, int NumColsA, double *B, int LDB, int NumRowsB, int NumColsB)
double ElapsedTime(void) const
Epetra_Time elapsed time function.
bool NoDiagonal() const
If matrix has no diagonal entries based on global row/column index comparisons, this query returns tr...
int checkMatvecSameVectors(Epetra_Comm &comm, bool verbose)
int checkEarlyDelete(Epetra_Comm &comm, bool verbose)
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
bool MyGRID(int GRID_in) const
Returns true if the GRID passed in belongs to the calling processor in this map, otherwise returns fa...
int NumMyRows() const
Returns the number of matrix rows on this processor.
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
int BeginExtractGlobalBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices) const
Initiates a view of the specified global row, only works if matrix indices are in global mode...
int power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
int NumMyBlockRows() const
Returns the number of Block matrix rows owned by the calling processor.
int checkVbrMatrixOptimizedGraph(Epetra_Comm &comm, bool verbose)
double NormOne() const
Returns the one norm of the global matrix.
Epetra_VbrMatrix: A class for the construction and use of real-valued double-precision variable block...
virtual const Epetra_Map & OperatorDomainMap() const =0
Returns the Epetra_Map object associated with the domain of this operator.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int ExtractEntryView(Epetra_SerialDenseMatrix *&entry) const
Returns a pointer to the current block entry.
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A &lt;- ScalarConstant * A)...
#define EPETRA_MIN(x, y)
int BeginSumIntoGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate summing into current values with this list of entries for a given global row of the matrix...
virtual int LeftScale(const Epetra_Vector &x)=0
Scales the Epetra_RowMatrix on the left with a Epetra_Vector x.
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
int ExtractMyBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Initiates a view of the specified local row, only works if matrix indices are in local mode...
int BeginSumIntoMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate summing into current values with this list of entries for a given local row of the matrix...
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
std::string Epetra_Version()
virtual void Barrier() const =0
Epetra_Comm Barrier function.
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
bool MyLRID(int LRID_in) const
Returns true if the LRID passed in belongs to the calling processor in this map, otherwise returns fa...
virtual int InvRowSums(Epetra_Vector &x) const =0
Computes the sum of absolute values of the rows of the Epetra_RowMatrix, results returned in x...
int NumGlobalBlockCols() const
Returns the number of global Block matrix columns.
virtual int NumGlobalNonzeros() const =0
Returns the number of nonzero entries in the global matrix.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
int BeginReplaceMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate replacement of current values with this list of entries for a given local row of the matrix...
int NumVectors() const
Returns the number of vectors in the multi-vector.
int GRID(int LRID_in) const
Returns the global row index for give local row index, returns IndexBase-1 if we don&#39;t have this loca...
virtual int MyPID() const =0
Return my process ID.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
const Epetra_Comm & Comm() const
Fills a matrix with rows from a source matrix based on the specified importer.
virtual int NumMyCols() const =0
Returns the number of matrix columns owned by the calling processor.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
int MaxColDim() const
Returns the maximum column dimension of all block entries on this processor.
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
bool LowerTriangular() const
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
double * A() const
Returns pointer to the this matrix.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
int NumMyElements() const
Number of elements on the calling processor.
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
Returns the number of nonzero entries in MyRow.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
void ConvertVbrToCrs(Epetra_VbrMatrix *VbrIn, Epetra_CrsMatrix *&CrsOut)
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
virtual const Epetra_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int BeginExtractMyBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices) const
Initiates a view of the specified local row, only works if matrix indices are in local mode...
Epetra_VbrRowMatrix: A class for using an existing Epetra_VbrMatrix object as an Epetra_RowMatrix obj...
int checkMergeRedundantEntries(Epetra_Comm &comm, bool verbose)
int GlobalMaxColDim() const
Returns the maximum column dimension of all block entries across all processors.
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
int NumMyNonzeros() const
Returns the number of nonzero entriesowned by the calling processor .
virtual bool UseTranspose() const =0
Returns the current UseTranspose setting.
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the right with a Epetra_Vector x.
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
int NumGlobalBlockEntries() const
Returns the number of nonzero block entries in the global matrix.
virtual double NormInf() const =0
Returns the infinity norm of the global matrix.
int ExtractGlobalBlockRowPointers(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Copy the block indices into user-provided array, set pointers for rest of data for specified global b...
int TestMatrix(Epetra_Comm &Comm, bool verbose, bool debug, int NumMyElements, int MinSize, int MaxSize, ConsType ConstructorType, bool ExtraBlocks, bool insertlocal, bool symmetric, Epetra_VbrMatrix **PreviousA)
virtual int NumMyRows() const =0
Returns the number of matrix rows owned by the calling processor.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
bool MyGlobalBlockRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
virtual bool Filled() const =0
If FillComplete() has been called, this query returns true, otherwise it returns false.
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
int MaxRowDim() const
Returns the maximum row dimension of all block entries on this processor.
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false...
virtual int NumGlobalDiagonals() const =0
Returns the number of global nonzero diagonal entries, based on global row/column index comparisons...
virtual int NumGlobalCols() const =0
Returns the number of global matrix columns.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int GlobalMaxRowDim() const
Returns the maximum row dimension of all block entries across all processors.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the left with a Epetra_Vector x.
int checkmultiply(bool transpose, Epetra_VbrMatrix &A, Epetra_MultiVector &X, Epetra_MultiVector &Check_Y)
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
int MaxMyGID() const
Returns the maximum global ID owned by this processor.
int ExtractMyBlockRowPointers(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Copy the block indices into user-provided array, set pointers for rest of data for specified local bl...
int MinMyGID() const
Returns the minimum global ID owned by this processor.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
const Epetra_Map & RowMatrixRowMap() const
Returns the EpetraMap object associated with the rows of this matrix.
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
int BeginExtractMyBlockRowCopy(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, int *ColDims) const
Initiates a copy of the specified local row in user-provided arrays.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
matrix_data is a very simple data source to be used for filling test matrices.
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
int LDA() const
Returns the leading dimension of the this matrix.
int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
virtual int NumMyDiagonals() const =0
Returns the number of local nonzero diagonal entries, based on global row/column index comparisons...
int NumGlobalRows() const
Returns the number of global matrix rows.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
virtual int NumProc() const =0
Returns total number of processes.
virtual bool HasNormInf() const =0
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator. ...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int BeginExtractGlobalBlockRowCopy(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, int *ColDims) const
Initiates a copy of the specified global row in user-provided arrays.
int FillComplete()
Signal that data entry is complete, perform transformations to local index space. ...
bool UpperTriangular() const
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given global row of the matrix, values are inserted via...
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
Submit a block entry to the indicated block row and column specified in the Begin routine...
int main(int argc, char *argv[])
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
Returns a copy of the specified local row in user-provided arrays.
int GCID(int LCID_in) const
Returns the global column index for give local column index, returns IndexBase-1 if we don&#39;t have thi...
virtual const Epetra_Map & RowMatrixColMap() const =0
Returns the Epetra_Map object associated with the columns of this matrix.
virtual int InvColSums(Epetra_Vector &x) const =0
Computes the sum of absolute values of the columns of the Epetra_RowMatrix, results returned in x...
int N() const
Returns column dimension of system.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
int checkValues(double x, double y, string message="", bool verbose=false)
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the with those in the user-provided vector.
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int InsertMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given local row of the matrix.
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
int NumMyBlockEntries() const
Returns the number of nonzero block entries in the calling processor&#39;s portion of the matrix...
int MaxNumBlockEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
int checkMultiVectors(Epetra_MultiVector &X, Epetra_MultiVector &Y, string message="", bool verbose=false)
int RightScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the right with a Epetra_Vector x.
int M() const
Returns row dimension of system.
int ExtractGlobalBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Initiates a view of the specified global row, only works if matrix indices are in global mode...
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
int NumGlobalBlockRows() const
Returns the number of global Block matrix rows.
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
virtual int NumMyNonzeros() const =0
Returns the number of nonzero entries in the calling processor&#39;s portion of the matrix.
double NormOne() const
Returns the one norm of the global matrix.
#define EPETRA_MAX(x, y)
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the left with a Epetra_Vector x.
int checkVbrRowMatrix(Epetra_RowMatrix &A, Epetra_RowMatrix &B, bool verbose)
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false...
int OptimizeStorage()
Eliminates memory that is used for construction. Make consecutive row index sections contiguous...
int checkExtractMyRowCopy(Epetra_Comm &comm, bool verbose)