Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
67 void GenerateCrsProblem(int nx, int ny, int npoints,
68  int * xoff, int * yoff, int nrhs,
69  const Epetra_Comm &comm,
70  Epetra_Map *& map,
72  Epetra_MultiVector *& x,
73  Epetra_MultiVector *& b,
74  Epetra_MultiVector *&xexact);
75 
76 void 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,
81  Epetra_MultiVector *& x,
82  Epetra_MultiVector *& b,
83  Epetra_MultiVector *&xexact);
84 
85 int 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 
99  Epetra_SerialComm Comm;
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;
283  Epetra_MultiVector resid(A->OperatorRangeMap(), 1);
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 
305 void 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,
310  Epetra_MultiVector *& x,
311  Epetra_MultiVector *& b,
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 
370 void 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,
375  Epetra_MultiVector *& x,
376  Epetra_MultiVector *& b,
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 
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:127
int Random()
Set multi-vector values to random numbers.
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...
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, transpose of this operator will be applied.
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.
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer...
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
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.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
void GenerateVbrProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints, int *xoff, int *yoff, int nsizes, int *sizes, const Epetra_Comm &comm, bool verbose, bool summary, Epetra_BlockMap *&map, Epetra_VbrMatrix *&A, Epetra_Vector *&b, Epetra_Vector *&bt, Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly)
Epetra_VbrMatrix: A class for the construction and use of real-valued double-precision variable block...
virtual const Epetra_Map & OperatorDomainMap() const =0
Returns the Epetra_Map object associated with the domain of this operator.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
std::string Epetra_Version()
int * Values() const
Returns a pointer to an array containing the values of this vector.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator. ...
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
int checkResults(bool trans, Epetra_LinearProblemRedistor *redistor, Epetra_LinearProblem *OrigProb, Epetra_LinearProblem *RedistProb, bool verbose)
Definition: bug1_main.cpp:185
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:83
int IndexBase() const
Index base for this map.
int NumMyElements() const
Number of elements on the calling processor.
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_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int UpdateTransposeValues(Epetra_RowMatrix *MatrixWithNewValues)
Update the values of an already-redistributed problem.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:81
void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints, int *xoff, int *yoff, const Epetra_Comm &comm, bool verbose, bool summary, Epetra_Map *&map, Epetra_CrsMatrix *&A, Epetra_Vector *&b, Epetra_Vector *&bt, Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly)
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Epetra_SerialComm: The Epetra Serial Communication Class.
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 Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Epetra_RowMatrixTransposer: A class for transposing an Epetra_RowMatrix object.
int InvRowSums(Epetra_Vector &x) const
Computes the sum of absolute values of the rows of the Epetra_VbrMatrix, results returned in x...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int main(int argc, char *argv[])
int Broadcast(double *MyVals, int Count, int Root) const
Epetra_SerialComm Broadcast function.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
int n
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
#define EPETRA_MAX(x, y)
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the left with a Epetra_Vector x.