Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/RowMatrixTransposer/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#include "Epetra_Map.h"
44#include "Epetra_LocalMap.h"
45#include "Epetra_Time.h"
46#include "Epetra_CrsMatrix.h"
47#include "Epetra_VbrMatrix.h"
48#include "Epetra_Vector.h"
49#include "Epetra_Flops.h"
50#ifdef EPETRA_MPI
51#include "Epetra_MpiComm.h"
52#include "mpi.h"
53#else
54#include "Epetra_SerialComm.h"
55#endif
56#include "../epetra_test_err.h"
57#include "Epetra_IntVector.h"
58#include "Epetra_Version.h"
60#include "Epetra_Time.h"
61
63 Epetra_CrsMatrix * transA,
64 Epetra_MultiVector * xexact,
65 bool verbose);
66
67void GenerateCrsProblem(int nx, int ny, int npoints,
68 int * xoff, int * yoff, int nrhs,
69 const Epetra_Comm &comm,
70 Epetra_Map *& map,
74 Epetra_MultiVector *&xexact);
75
76void GenerateVbrProblem(int nx, int ny, int npoints, int * xoff, int * yoff,
77 int nsizes, int * sizes, int nrhs,
78 const Epetra_Comm &comm,
79 Epetra_BlockMap *& map,
83 Epetra_MultiVector *&xexact);
84
85int main(int argc, char *argv[]) {
86
87 int ierr = 0, i;
88
89#ifdef EPETRA_MPI
90
91 // Initialize MPI
92
93 MPI_Init(&argc,&argv);
94
95 Epetra_MpiComm Comm( MPI_COMM_WORLD );
96
97#else
98
100
101#endif
102
103 bool verbose = false;
104
105 // Check if we should print results to standard out
106 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
107
108 int verbose_int = verbose ? 1 : 0;
109 Comm.Broadcast(&verbose_int, 1, 0);
110 verbose = verbose_int==1 ? true : false;
111
112
113 // char tmp;
114 // if (Comm.MyPID()==0) cout << "Press any key to continue..."<< endl;
115 // if (Comm.MyPID()==0) cin >> tmp;
116 // Comm.Barrier();
117
118 Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
119 int MyPID = Comm.MyPID();
120 int NumProc = Comm.NumProc();
121
122 if(verbose && MyPID==0)
123 cout << Epetra_Version() << endl << endl;
124
125 int nx = 128;
126 int ny = NumProc*nx; // Scale y grid with number of processors
127
128 // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
129
130 // (i-1,j-1) (i-1,j )
131 // (i ,j-1) (i ,j ) (i ,j+1)
132 // (i+1,j-1) (i+1,j )
133
134 int npoints = 7;
135
136 int xoff[] = {-1, 0, 1, -1, 0, 1, 0};
137 int yoff[] = {-1, -1, -1, 0, 0, 0, 1};
138
139 Epetra_Map * map;
141 Epetra_MultiVector * x, * b, * xexact;
142
143 GenerateCrsProblem(nx, ny, npoints, xoff, yoff, 1, Comm, map, A, x, b, xexact);
144
145 if (nx<8) {
146 cout << *A << endl;
147 cout << "X exact = " << endl << *xexact << endl;
148 cout << "B = " << endl << *b << endl;
149 }
150 // Construct transposer object
151
152 Epetra_Time timer(Comm);
153
154 double start = timer.ElapsedTime();
155 Epetra_RowMatrixTransposer transposer(A);
156 if (verbose) cout << "\nTime to construct transposer = "
157 << timer.ElapsedTime() - start << endl;
158
159 bool MakeDataContiguous = true;
160 Epetra_CrsMatrix * transA;
161
162 start = timer.ElapsedTime();
163 transposer.CreateTranspose(MakeDataContiguous, transA);
164 if (verbose) cout << "\nTime to create transpose matrix = "
165 << timer.ElapsedTime() - start << endl;
166
167
168 // Now test output of transposer by performing matvecs
169 ierr += checkResults(A, transA, xexact, verbose);
170
171
172 // Now change values in original matrix and test update facility of transposer
173 // Add 2 to the diagonal of each row
174
175 double Value = 2.0;
176 for (i=0; i< A->NumMyRows(); i++)
177 A->SumIntoMyValues(i, 1, &Value, &i);
178
179
180 start = timer.ElapsedTime();
181 transposer.UpdateTransposeValues(A);
182 if (verbose) cout << "\nTime to update transpose matrix = "
183 << timer.ElapsedTime() - start << endl;
184
185 ierr += checkResults(A, transA, xexact, verbose);
186
187 delete transA;
188 delete A;
189 delete b;
190 delete x;
191 delete xexact;
192 delete map;
193
194
195 if (verbose) cout << endl << "Checking transposer for VbrMatrix objects" << endl<< endl;
196
197 int nsizes = 4;
198 int sizes[] = {4, 6, 5, 3};
199
200 Epetra_VbrMatrix * Avbr;
201 Epetra_BlockMap * bmap;
202
203 GenerateVbrProblem(nx, ny, npoints, xoff, yoff, nsizes, sizes,
204 1, Comm, bmap, Avbr, x, b, xexact);
205
206 if (nx<8) {
207 cout << *Avbr << endl;
208 cout << "X exact = " << endl << *xexact << endl;
209 cout << "B = " << endl << *b << endl;
210 }
211
212 start = timer.ElapsedTime();
213 Epetra_RowMatrixTransposer transposer1(Avbr);
214 if (verbose) cout << "\nTime to construct transposer = "
215 << timer.ElapsedTime() - start << endl;
216
217 start = timer.ElapsedTime();
218 transposer1.CreateTranspose(MakeDataContiguous, transA);
219 if (verbose) cout << "\nTime to create transpose matrix = "
220 << timer.ElapsedTime() - start << endl;
221
222
223 // Now test output of transposer by performing matvecs
224 ierr += checkResults(Avbr, transA, xexact, verbose);
225
226 // Now change values in original matrix and test update facility of transposer
227 // Scale matrix on the left by rowsums
228
229 Epetra_Vector invRowSums(Avbr->RowMap());
230
231 Avbr->InvRowSums(invRowSums);
232 Avbr->LeftScale(invRowSums);
233
234 start = timer.ElapsedTime();
235 transposer1.UpdateTransposeValues(Avbr);
236 if (verbose) cout << "\nTime to update transpose matrix = "
237 << timer.ElapsedTime() - start << endl;
238
239 ierr += checkResults(Avbr, transA, xexact, verbose);
240
241 delete transA;
242 delete Avbr;
243 delete b;
244 delete x;
245 delete xexact;
246 delete bmap;
247
248#ifdef EPETRA_MPI
249 MPI_Finalize() ;
250#endif
251
252/* end main */
253 return ierr ;
254}
255
257 Epetra_MultiVector * xexact, bool verbose)
258{
259 int n = A->NumGlobalRows();
260
261 if (n<100) cout << "A transpose = " << endl << *transA << endl;
262
263 Epetra_MultiVector x1(View, A->OperatorRangeMap(), &((*xexact)[0]), 1);
265
266 A->SetUseTranspose(true);
267
268 Epetra_Time timer(A->Comm());
269 double start = timer.ElapsedTime();
270 A->Apply(x1, b1);
271 if (verbose) cout << "\nTime to compute b1: matvec with original matrix using transpose flag = " << timer.ElapsedTime() - start << endl;
272
273 if (n<100) cout << "b1 = " << endl << b1 << endl;
274 Epetra_MultiVector x2(View, transA->OperatorDomainMap(), &((*xexact)[0]), 1);
275 Epetra_MultiVector b2(transA->OperatorRangeMap(), 1);
276 start = timer.ElapsedTime();
277 transA->Multiply(false, x2, b2);
278 if (verbose) cout << "\nTime to compute b2: matvec with transpose matrix = " << timer.ElapsedTime() - start << endl;
279
280 if (n<100) cout << "b1 = " << endl << b1 << endl;
281
282 double residual;
284
285 resid.Update(1.0, b1, -1.0, b2, 0.0);
286#ifndef NDEBUG
287 int ierr0 =
288#endif
289 resid.Norm2(&residual);
290 assert(ierr0==0);
291 if (verbose) cout << "Norm of b1 - b2 = " << residual << endl;
292
293 int ierr = 0;
294
295 if (residual > 1.0e-10) ierr++;
296
297 if (ierr!=0 && verbose) {
298 cerr << "Status: Test failed" << endl;
299 }
300 else if (verbose) cerr << "Status: Test passed" << endl;
301
302 return(ierr);
303}
304
305void GenerateCrsProblem(int nx, int ny, int npoints,
306 int * xoff, int * yoff, int nrhs,
307 const Epetra_Comm &comm,
308 Epetra_Map *& map,
309 Epetra_CrsMatrix *& A,
312 Epetra_MultiVector *&xexact)
313{
314 // Number of global equations is nx*ny. These will be distributed in a linear fashion
315 int numGlobalEquations = nx*ny;
316 map = new Epetra_Map(numGlobalEquations, 0, comm); // Create map with equal distribution of equations.
317
318 int numMyEquations = map->NumMyElements();
319
320 A = new Epetra_CrsMatrix(Copy, *map, 0); // Construct matrix
321
322 int * indices = new int[npoints];
323 double * values = new double[npoints];
324
325 double dnpoints = (double) npoints;
326
327 for (int i=0; i<numMyEquations; i++) {
328
329 int rowID = map->GID(i);
330 int numIndices = 0;
331
332 for (int j=0; j<npoints; j++) {
333 int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
334 if (colID>-1 && colID<numGlobalEquations) {
335 indices[numIndices] = colID;
336 double value = - ((double) rand())/ ((double) RAND_MAX);
337 if (colID==rowID)
338 values[numIndices++] = dnpoints - value; // Make diagonal dominant
339 else
340 values[numIndices++] = -value;
341 }
342 }
343
344 A->InsertGlobalValues(rowID, numIndices, values, indices);
345 }
346
347 delete [] indices;
348 delete [] values;
349
350 A->FillComplete();
351
352 if (nrhs<=1) {
353 x = new Epetra_Vector(*map);
354 b = new Epetra_Vector(*map);
355 xexact = new Epetra_Vector(*map);
356 }
357 else {
358 x = new Epetra_MultiVector(*map, nrhs);
359 b = new Epetra_MultiVector(*map, nrhs);
360 xexact = new Epetra_MultiVector(*map, nrhs);
361 }
362
363 xexact->Random(); // Fill xexact with random values
364
365 A->Multiply(false, *xexact, *b);
366
367 return;
368}
369
370void GenerateVbrProblem(int nx, int ny, int npoints, int * xoff, int * yoff,
371 int nsizes, int * sizes, int nrhs,
372 const Epetra_Comm &comm,
373 Epetra_BlockMap *& map,
374 Epetra_VbrMatrix *& A,
377 Epetra_MultiVector *&xexact)
378{
379 int i, j;
380
381 // Number of global equations is nx*ny. These will be distributed in a linear fashion
382 int numGlobalEquations = nx*ny;
383 Epetra_Map ptMap(numGlobalEquations, 0, comm); // Create map with equal distribution of equations.
384
385 int numMyElements = ptMap.NumMyElements();
386
387 Epetra_IntVector elementSizes(ptMap); // This vector will have the list of element sizes
388 for (i=0; i<numMyElements; i++)
389 elementSizes[i] = sizes[ptMap.GID(i)%nsizes]; // cycle through sizes array
390
391 map = new Epetra_BlockMap(-1, numMyElements, ptMap.MyGlobalElements(), elementSizes.Values(), ptMap.IndexBase(), ptMap.Comm());
392
393
394 A = new Epetra_VbrMatrix(Copy, *map, 0); // Construct matrix
395
396 int * indices = new int[npoints];
397
398 // This section of code creates a vector of random values that will be used to create
399 // light-weight dense matrices to pass into the VbrMatrix construction process.
400
401 int maxElementSize = 0;
402 for (i=0; i< nsizes; i++) maxElementSize = EPETRA_MAX(maxElementSize, sizes[i]);
403
404 Epetra_LocalMap lmap(maxElementSize*maxElementSize, ptMap.IndexBase(), ptMap.Comm());
405 Epetra_Vector randvec(lmap);
406 randvec.Random();
407 randvec.Scale(-1.0); // Make value negative
408
409
410 for (i=0; i<numMyElements; i++) {
411 int rowID = map->GID(i);
412 int numIndices = 0;
413 int rowDim = sizes[rowID%nsizes];
414 for (j=0; j<npoints; j++) {
415 int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
416 if (colID>-1 && colID<numGlobalEquations)
417 indices[numIndices++] = colID;
418 }
419
420 A->BeginInsertGlobalValues(rowID, numIndices, indices);
421
422 for (j=0; j < numIndices; j++) {
423 int colDim = sizes[indices[j]%nsizes];
424 A->SubmitBlockEntry(&(randvec[0]), rowDim, rowDim, colDim);
425 }
426 A->EndSubmitEntries();
427 }
428
429 delete [] indices;
430
431 A->FillComplete();
432
433 // Compute the InvRowSums of the matrix rows
434 Epetra_Vector invRowSums(A->RowMap());
435 Epetra_Vector rowSums(A->RowMap());
436 A->InvRowSums(invRowSums);
437 rowSums.Reciprocal(invRowSums);
438
439 // Jam the row sum values into the diagonal of the Vbr matrix (to make it diag dominant)
440 int numBlockDiagonalEntries;
441 int * rowColDims;
442 int * diagoffsets = map->FirstPointInElementList();
443 A->BeginExtractBlockDiagonalView(numBlockDiagonalEntries, rowColDims);
444 for (i=0; i< numBlockDiagonalEntries; i++) {
445 double * diagVals;
446 int diagLDA;
447 A->ExtractBlockDiagonalEntryView(diagVals, diagLDA);
448 int rowDim = map->ElementSize(i);
449 for (j=0; j<rowDim; j++) diagVals[j+j*diagLDA] = rowSums[diagoffsets[i]+j];
450 }
451
452 if (nrhs<=1) {
453 x = new Epetra_Vector(*map);
454 b = new Epetra_Vector(*map);
455 xexact = new Epetra_Vector(*map);
456 }
457 else {
458 x = new Epetra_MultiVector(*map, nrhs);
459 b = new Epetra_MultiVector(*map, nrhs);
460 xexact = new Epetra_MultiVector(*map, nrhs);
461 }
462
463 xexact->Random(); // Fill xexact with random values
464
465 A->Multiply(false, *xexact, *b);
466
467 return;
468}
469
#define EPETRA_MAX(x, y)
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 GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor.
int IndexBase() const
Index base for this map.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition Epetra_Comm.h:73
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Add this list of entries to existing values for a given local row of the matrix.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
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.
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer.
int * Values() const
Returns a pointer to an array containing the values of this vector.
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
int Broadcast(double *MyVals, int Count, int Root) const
Epetra_MpiComm Broadcast function.
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 Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int Reciprocal(const Epetra_MultiVector &A)
Puts element-wise reciprocal values of input Multi-vector in target.
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.
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 const Epetra_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
Epetra_RowMatrixTransposer: A class for transposing an Epetra_RowMatrix object.
int UpdateTransposeValues(Epetra_RowMatrix *MatrixWithNewValues)
Update the values of an already-redistributed problem.
int CreateTranspose(const bool MakeDataContiguous, Epetra_CrsMatrix *&TransposeMatrix, Epetra_Map *TransposeRowMap=0)
Generate a new Epetra_CrsMatrix as the transpose of an Epetra_RowMatrix passed into the constructor.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Time: The Epetra Timing Class.
Definition Epetra_Time.h:75
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...
int InvRowSums(Epetra_Vector &x) const
Computes the sum of absolute values of the rows of the Epetra_VbrMatrix, results returned in x.
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
Submit a block entry to the indicated block row and column specified in the Begin routine.
int BeginExtractBlockDiagonalView(int &NumBlockDiagonalEntries, int *&RowColDims) const
Initiates a view of the block diagonal entries.
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 FillComplete()
Signal that data entry is complete, perform transformations to local index space.
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_MultiVector X in Y.
int ExtractBlockDiagonalEntryView(double *&Values, int &LDA) const
Extract a view of a block diagonal entry from the matrix.
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the left with a Epetra_Vector x.
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])
void GenerateVbrProblem(int nx, int ny, int npoints, int *xoff, int *yoff, int nsizes, int *sizes, int nrhs, const Epetra_Comm &comm, Epetra_BlockMap *&map, Epetra_VbrMatrix *&A, Epetra_MultiVector *&x, Epetra_MultiVector *&b, Epetra_MultiVector *&xexact)
int checkResults(Epetra_RowMatrix *A, Epetra_CrsMatrix *transA, Epetra_MultiVector *xexact, bool verbose)
void GenerateCrsProblem(int nx, int ny, int npoints, int *xoff, int *yoff, int nrhs, const Epetra_Comm &comm, Epetra_Map *&map, Epetra_CrsMatrix *&A, Epetra_MultiVector *&x, Epetra_MultiVector *&b, Epetra_MultiVector *&xexact)