Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
74int 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
85int 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
102int checkVbrMatrixOptimizedGraph(Epetra_Comm& comm, bool verbose);
103
104int checkVbrRowMatrix(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose);
105
106int 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
114int power_method(bool TransA, Epetra_VbrMatrix& A,
117 Epetra_MultiVector& resid,
118 double * lambda, int niters, double tolerance,
119 bool verbose);
120
121int checkMergeRedundantEntries(Epetra_Comm& comm, bool verbose);
122
123int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose);
124
125int checkMatvecSameVectors(Epetra_Comm& comm, bool verbose);
126
127int 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
295int 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;
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;
395 A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, NumNz ) ;
396 HaveColMap = true;
397 break;
399 A = new Epetra_VbrMatrix( Copy, *rowmap, *colmap, 3 ) ;
400 FixedNumEntries = true;
401 HaveColMap = true;
402 break;
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);
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 }
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);
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);
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);
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);
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);
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
1188int 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*/
1409return ierr ;
1410}
1411
1412int power_method(bool TransA, Epetra_VbrMatrix& A,
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//=============================================================================
1716int 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
1738{
1739 int numProcs = comm.NumProc();
1740 int localProc = comm.MyPID();
1741
1742 int myFirstRow = localProc*3;
1743 unsigned 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 unsigned 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(unsigned 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
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(unsigned 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(unsigned int i=myFirstRow; i<=myLastRow; ++i) {
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
1880int checkExtractMyRowCopy(Epetra_Comm& comm, bool verbose)
1881{
1882 int numProcs = comm.NumProc();
1883 int localProc = comm.MyPID();
1884
1885 int myFirstRow = localProc*3;
1886 unsigned 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 unsigned 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(unsigned 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
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(unsigned 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
1973int 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
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
2031int 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
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
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);
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 */
2357 EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
2359 EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
2360 EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);
2361
2363 EPETRA_TEST_ERR(A.InvColSums(colA),ierr); // -2 error code
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
#define EPETRA_MIN(x, y)
#define EPETRA_MAX(x, y)
@ View
@ Copy
std::string Epetra_Version()
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int MinMyGID() const
Returns the minimum global ID owned by this processor.
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor.
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
int NumGlobalElements() const
Number of elements across all processors.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
int MaxMyGID() const
Returns the maximum global ID owned by this processor.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int NumMyBlockRows() const
Returns the number of block matrix rows on this processor.
int NumMyRows() const
Returns the number of matrix rows on this processor.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the left with a Epetra_Vector x.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the right with a Epetra_Vector x.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
double NormOne() const
Returns the one norm of the global matrix.
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
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.
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.
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 ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
double * Values() const
Get pointer to MultiVector values.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, transpose of this operator will be applied.
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.
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual const Epetra_Map & OperatorDomainMap() const =0
Returns the Epetra_Map object associated with the domain of this operator.
virtual bool HasNormInf() const =0
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
virtual bool UseTranspose() const =0
Returns the current UseTranspose setting.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual int NumMyRows() const =0
Returns the number of matrix rows owned by the calling processor.
virtual int NumMyCols() const =0
Returns the number of matrix columns owned by the calling processor.
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.
virtual int NumGlobalCols() const =0
Returns the number of global matrix columns.
virtual int NumGlobalNonzeros() const =0
Returns the number of nonzero entries in the global matrix.
virtual const Epetra_Map & RowMatrixColMap() const =0
Returns the Epetra_Map object associated with the columns of this matrix.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
virtual int NumMyNonzeros() const =0
Returns the number of nonzero entries in the calling processor's portion of the matrix.
virtual int NumMyDiagonals() const =0
Returns the number of local nonzero diagonal entries, based on global row/column index comparisons.
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
Returns the number of nonzero entries in MyRow.
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
virtual int NumGlobalDiagonals() const =0
Returns the number of global nonzero diagonal entries, based on global row/column index comparisons.
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
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.
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.
virtual int RightScale(const Epetra_Vector &x)=0
Scales the Epetra_RowMatrix on the right with a Epetra_Vector x.
virtual double NormOne() const =0
Returns the one norm of the global matrix.
virtual double NormInf() const =0
Returns the infinity norm of the global matrix.
virtual bool Filled() const =0
If FillComplete() has been called, this query returns true, otherwise it returns false.
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.
virtual int LeftScale(const Epetra_Vector &x)=0
Scales the Epetra_RowMatrix on the left with a Epetra_Vector x.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * A() const
Returns pointer to the this matrix.
int LDA() const
Returns the leading dimension of the this matrix.
int Random()
Set matrix values to random numbers.
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_VbrMatrix: A class for the construction and use of real-valued double-precision variable block...
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false.
int MaxColDim() const
Returns the maximum column dimension of all block entries on this processor.
int NumMyBlockEntries() const
Returns the number of nonzero block entries in the calling processor's portion of the matrix.
int ExtractEntryView(Epetra_SerialDenseMatrix *&entry) const
Returns a pointer to the current block entry.
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 Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
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...
bool UpperTriangular() const
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
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 PutScalar(double ScalarConstant)
Initialize all values in graph of the matrix with constant value.
int NumGlobalBlockEntries() const
Returns the number of nonzero block entries in the global matrix.
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.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
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 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.
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 NumMyNonzeros() const
Returns the number of nonzero entriesowned by the calling processor .
bool NoDiagonal() const
If matrix has no diagonal entries based on global row/column index comparisons, this query returns tr...
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,...
int GlobalMaxRowDim() const
Returns the maximum row dimension of all block entries across all processors.
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...
bool MyGlobalBlockRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
int MaxRowDim() const
Returns the maximum row dimension of all block entries on this processor.
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.
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the with those in the user-provided vector.
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.
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 LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
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 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...
const Epetra_CrsGraph & Graph() const
Returns a pointer to the Epetra_CrsGraph object associated with this matrix.
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
const Epetra_Comm & Comm() const
Fills a matrix with rows from a source matrix based on the specified importer.
int MaxNumBlockEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
int FillComplete()
Signal that data entry is complete, perform transformations to local index space.
int NumMyBlockRows() const
Returns the number of Block matrix rows owned by the calling processor.
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 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 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.
double NormInf() const
Returns the infinity norm of the global matrix.
bool IndicesAreContiguous() const
If matrix indices are packed into single array (done in OptimizeStorage()) return true,...
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.
bool LowerTriangular() const
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
int GCID(int LCID_in) const
Returns the global column index for give local column index, returns IndexBase-1 if we don't have thi...
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 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,...
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the left with a Epetra_Vector x.
int NumGlobalBlockRows() const
Returns the number of global Block matrix rows.
const Epetra_Map & RowMatrixColMap() const
Returns the Epetra_Map object associated with columns of this matrix.
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
int GRID(int LRID_in) const
Returns the global row index for give local row index, returns IndexBase-1 if we don't have this loca...
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the right with a Epetra_Vector x.
int OptimizeStorage()
Eliminates memory that is used for construction. Make consecutive row index sections contiguous.
double NormOne() const
Returns the one norm of the global matrix.
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A <- ScalarConstant * A).
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.
bool MyLRID(int LRID_in) const
Returns true if the LRID passed in belongs to the calling processor in this map, otherwise returns fa...
int NumGlobalBlockCols() const
Returns the number of global Block matrix columns.
int GlobalMaxColDim() const
Returns the maximum column dimension of all block entries across all processors.
const Epetra_Map & RowMatrixRowMap() const
Returns the EpetraMap object associated with the rows of this matrix.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
int NumGlobalRows() const
Returns the number of global matrix rows.
Epetra_VbrRowMatrix: A class for using an existing Epetra_VbrMatrix object as an Epetra_RowMatrix obj...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
matrix_data is a very simple data source to be used for filling test matrices.
#define EPETRA_TEST_ERR(a, b)
int checkMatvecSameVectors(Epetra_Comm &comm, bool verbose)
int checkVbrRowMatrix(Epetra_RowMatrix &A, Epetra_RowMatrix &B, bool verbose)
int main(int argc, char *argv[])
@ VariableEntriesPerRow
void ConvertVbrToCrs(Epetra_VbrMatrix *VbrIn, Epetra_CrsMatrix *&CrsOut)
int checkMergeRedundantEntries(Epetra_Comm &comm, bool verbose)
int checkVbrMatrixOptimizedGraph(Epetra_Comm &comm, bool verbose)
int CompareValues(double *A, int LDA, int NumRowsA, int NumColsA, double *B, int LDB, int NumRowsB, int NumColsB)
int checkExtractMyRowCopy(Epetra_Comm &comm, bool verbose)
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)
int checkEarlyDelete(Epetra_Comm &comm, bool verbose)
int checkmultiply(bool transpose, Epetra_VbrMatrix &A, Epetra_MultiVector &X, Epetra_MultiVector &Check_Y)
int checkValues(double x, double y, string message="", bool verbose=false)
int power_method(bool TransA, Epetra_VbrMatrix &A, Epetra_MultiVector &q, Epetra_MultiVector &z, Epetra_MultiVector &resid, double *lambda, int niters, double tolerance, bool verbose)
int checkMultiVectors(Epetra_MultiVector &X, Epetra_MultiVector &Y, string message="", bool verbose=false)
int check(Epetra_VbrMatrix &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int NumMyBlockRows1, int NumGlobalBlockRows1, int NumMyBlockNonzeros1, int NumGlobalBlockNonzeros1, int *MyGlobalElements, bool verbose)