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 A;
188  delete b;
189  delete x;
190  delete xexact;
191  delete map;
192 
193 
194  if (verbose) cout << endl << "Checking transposer for VbrMatrix objects" << endl<< endl;
195 
196  int nsizes = 4;
197  int sizes[] = {4, 6, 5, 3};
198 
199  Epetra_VbrMatrix * Avbr;
200  Epetra_BlockMap * bmap;
201 
202  GenerateVbrProblem(nx, ny, npoints, xoff, yoff, nsizes, sizes,
203  1, Comm, bmap, Avbr, x, b, xexact);
204 
205  if (nx<8) {
206  cout << *Avbr << endl;
207  cout << "X exact = " << endl << *xexact << endl;
208  cout << "B = " << endl << *b << endl;
209  }
210 
211  start = timer.ElapsedTime();
212  Epetra_RowMatrixTransposer transposer1(Avbr);
213  if (verbose) cout << "\nTime to construct transposer = "
214  << timer.ElapsedTime() - start << endl;
215 
216  start = timer.ElapsedTime();
217  transposer1.CreateTranspose(MakeDataContiguous, transA);
218  if (verbose) cout << "\nTime to create transpose matrix = "
219  << timer.ElapsedTime() - start << endl;
220 
221 
222  // Now test output of transposer by performing matvecs
223  ierr += checkResults(Avbr, transA, xexact, verbose);
224 
225  // Now change values in original matrix and test update facility of transposer
226  // Scale matrix on the left by rowsums
227 
228  Epetra_Vector invRowSums(Avbr->RowMap());
229 
230  Avbr->InvRowSums(invRowSums);
231  Avbr->LeftScale(invRowSums);
232 
233  start = timer.ElapsedTime();
234  transposer1.UpdateTransposeValues(Avbr);
235  if (verbose) cout << "\nTime to update transpose matrix = "
236  << timer.ElapsedTime() - start << endl;
237 
238  ierr += checkResults(Avbr, transA, xexact, verbose);
239 
240  delete Avbr;
241  delete b;
242  delete x;
243  delete xexact;
244  delete bmap;
245 
246 #ifdef EPETRA_MPI
247  MPI_Finalize() ;
248 #endif
249 
250 /* end main */
251  return ierr ;
252 }
253 
255  Epetra_MultiVector * xexact, bool verbose)
256 {
257  int n = A->NumGlobalRows();
258 
259  if (n<100) cout << "A transpose = " << endl << *transA << endl;
260 
261  Epetra_MultiVector x1(View, A->OperatorRangeMap(), &((*xexact)[0]), 1);
263 
264  A->SetUseTranspose(true);
265 
266  Epetra_Time timer(A->Comm());
267  double start = timer.ElapsedTime();
268  A->Apply(x1, b1);
269  if (verbose) cout << "\nTime to compute b1: matvec with original matrix using transpose flag = " << timer.ElapsedTime() - start << endl;
270 
271  if (n<100) cout << "b1 = " << endl << b1 << endl;
272  Epetra_MultiVector x2(View, transA->OperatorDomainMap(), &((*xexact)[0]), 1);
273  Epetra_MultiVector b2(transA->OperatorRangeMap(), 1);
274  start = timer.ElapsedTime();
275  transA->Multiply(false, x2, b2);
276  if (verbose) cout << "\nTime to compute b2: matvec with transpose matrix = " << timer.ElapsedTime() - start << endl;
277 
278  if (n<100) cout << "b1 = " << endl << b1 << endl;
279 
280  double residual;
281  Epetra_MultiVector resid(A->OperatorRangeMap(), 1);
282 
283  resid.Update(1.0, b1, -1.0, b2, 0.0);
284 #ifndef NDEBUG
285  int ierr0 =
286 #endif
287  resid.Norm2(&residual);
288  assert(ierr0==0);
289  if (verbose) cout << "Norm of b1 - b2 = " << residual << endl;
290 
291  int ierr = 0;
292 
293  if (residual > 1.0e-10) ierr++;
294 
295  if (ierr!=0 && verbose) {
296  cerr << "Status: Test failed" << endl;
297  }
298  else if (verbose) cerr << "Status: Test passed" << endl;
299 
300  return(ierr);
301 }
302 
303 void GenerateCrsProblem(int nx, int ny, int npoints,
304  int * xoff, int * yoff, int nrhs,
305  const Epetra_Comm &comm,
306  Epetra_Map *& map,
307  Epetra_CrsMatrix *& A,
308  Epetra_MultiVector *& x,
309  Epetra_MultiVector *& b,
310  Epetra_MultiVector *&xexact)
311 {
312  // Number of global equations is nx*ny. These will be distributed in a linear fashion
313  int numGlobalEquations = nx*ny;
314  map = new Epetra_Map(numGlobalEquations, 0, comm); // Create map with equal distribution of equations.
315 
316  int numMyEquations = map->NumMyElements();
317 
318  A = new Epetra_CrsMatrix(Copy, *map, 0); // Construct matrix
319 
320  int * indices = new int[npoints];
321  double * values = new double[npoints];
322 
323  double dnpoints = (double) npoints;
324 
325  for (int i=0; i<numMyEquations; i++) {
326 
327  int rowID = map->GID(i);
328  int numIndices = 0;
329 
330  for (int j=0; j<npoints; j++) {
331  int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
332  if (colID>-1 && colID<numGlobalEquations) {
333  indices[numIndices] = colID;
334  double value = - ((double) rand())/ ((double) RAND_MAX);
335  if (colID==rowID)
336  values[numIndices++] = dnpoints - value; // Make diagonal dominant
337  else
338  values[numIndices++] = -value;
339  }
340  }
341 
342  A->InsertGlobalValues(rowID, numIndices, values, indices);
343  }
344 
345  delete [] indices;
346  delete [] values;
347 
348  A->FillComplete();
349 
350  if (nrhs<=1) {
351  x = new Epetra_Vector(*map);
352  b = new Epetra_Vector(*map);
353  xexact = new Epetra_Vector(*map);
354  }
355  else {
356  x = new Epetra_MultiVector(*map, nrhs);
357  b = new Epetra_MultiVector(*map, nrhs);
358  xexact = new Epetra_MultiVector(*map, nrhs);
359  }
360 
361  xexact->Random(); // Fill xexact with random values
362 
363  A->Multiply(false, *xexact, *b);
364 
365  return;
366 }
367 
368 void GenerateVbrProblem(int nx, int ny, int npoints, int * xoff, int * yoff,
369  int nsizes, int * sizes, int nrhs,
370  const Epetra_Comm &comm,
371  Epetra_BlockMap *& map,
372  Epetra_VbrMatrix *& A,
373  Epetra_MultiVector *& x,
374  Epetra_MultiVector *& b,
375  Epetra_MultiVector *&xexact)
376 {
377  int i, j;
378 
379  // Number of global equations is nx*ny. These will be distributed in a linear fashion
380  int numGlobalEquations = nx*ny;
381  Epetra_Map ptMap(numGlobalEquations, 0, comm); // Create map with equal distribution of equations.
382 
383  int numMyElements = ptMap.NumMyElements();
384 
385  Epetra_IntVector elementSizes(ptMap); // This vector will have the list of element sizes
386  for (i=0; i<numMyElements; i++)
387  elementSizes[i] = sizes[ptMap.GID(i)%nsizes]; // cycle through sizes array
388 
389  map = new Epetra_BlockMap(-1, numMyElements, ptMap.MyGlobalElements(), elementSizes.Values(), ptMap.IndexBase(), ptMap.Comm());
390 
391 
392  A = new Epetra_VbrMatrix(Copy, *map, 0); // Construct matrix
393 
394  int * indices = new int[npoints];
395 
396  // This section of code creates a vector of random values that will be used to create
397  // light-weight dense matrices to pass into the VbrMatrix construction process.
398 
399  int maxElementSize = 0;
400  for (i=0; i< nsizes; i++) maxElementSize = EPETRA_MAX(maxElementSize, sizes[i]);
401 
402  Epetra_LocalMap lmap(maxElementSize*maxElementSize, ptMap.IndexBase(), ptMap.Comm());
403  Epetra_Vector randvec(lmap);
404  randvec.Random();
405  randvec.Scale(-1.0); // Make value negative
406 
407 
408  for (i=0; i<numMyElements; i++) {
409  int rowID = map->GID(i);
410  int numIndices = 0;
411  int rowDim = sizes[rowID%nsizes];
412  for (j=0; j<npoints; j++) {
413  int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
414  if (colID>-1 && colID<numGlobalEquations)
415  indices[numIndices++] = colID;
416  }
417 
418  A->BeginInsertGlobalValues(rowID, numIndices, indices);
419 
420  for (j=0; j < numIndices; j++) {
421  int colDim = sizes[indices[j]%nsizes];
422  A->SubmitBlockEntry(&(randvec[0]), rowDim, rowDim, colDim);
423  }
424  A->EndSubmitEntries();
425  }
426 
427  delete [] indices;
428 
429  A->FillComplete();
430 
431  // Compute the InvRowSums of the matrix rows
432  Epetra_Vector invRowSums(A->RowMap());
433  Epetra_Vector rowSums(A->RowMap());
434  A->InvRowSums(invRowSums);
435  rowSums.Reciprocal(invRowSums);
436 
437  // Jam the row sum values into the diagonal of the Vbr matrix (to make it diag dominant)
438  int numBlockDiagonalEntries;
439  int * rowColDims;
440  int * diagoffsets = map->FirstPointInElementList();
441  A->BeginExtractBlockDiagonalView(numBlockDiagonalEntries, rowColDims);
442  for (i=0; i< numBlockDiagonalEntries; i++) {
443  double * diagVals;
444  int diagLDA;
445  A->ExtractBlockDiagonalEntryView(diagVals, diagLDA);
446  int rowDim = map->ElementSize(i);
447  for (j=0; j<rowDim; j++) diagVals[j+j*diagLDA] = rowSums[diagoffsets[i]+j];
448  }
449 
450  if (nrhs<=1) {
451  x = new Epetra_Vector(*map);
452  b = new Epetra_Vector(*map);
453  xexact = new Epetra_Vector(*map);
454  }
455  else {
456  x = new Epetra_MultiVector(*map, nrhs);
457  b = new Epetra_MultiVector(*map, nrhs);
458  xexact = new Epetra_MultiVector(*map, nrhs);
459  }
460 
461  xexact->Random(); // Fill xexact with random values
462 
463  A->Multiply(false, *xexact, *b);
464 
465  return;
466 }
467 
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:119
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:75
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:73
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.