EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/Transpose_LL/cxx_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - 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 // Epetra_BlockMap Test routine
43 
44 #include "EpetraExt_Version.h"
45 
46 #include "Epetra_CrsMatrix.h"
47 #include "Epetra_VbrMatrix.h"
48 #include "Epetra_Vector.h"
49 #include "Epetra_MultiVector.h"
50 #include "Epetra_LocalMap.h"
51 #include "Epetra_IntVector.h"
52 #include "Epetra_Map.h"
53 
54 #ifdef EPETRA_MPI
55 #include "Epetra_MpiComm.h"
56 #include <mpi.h>
57 #else
58 #include "Epetra_SerialComm.h"
59 #endif
60 
61 #include "Epetra_Time.h"
63 #include "Trilinos_Util.h"
64 
66  Epetra_Vector * xexact, bool verbose);
67 
68 int main(int argc, char *argv[])
69 {
70  int i;
71 
72 #ifdef EPETRA_MPI
73  // Initialize MPI
74  MPI_Init(&argc,&argv);
75  Epetra_MpiComm comm(MPI_COMM_WORLD);
76 #else
77  Epetra_SerialComm comm;
78 #endif
79 
80  // Uncomment to debug in parallel int tmp; if (comm.MyPID()==0) cin >> tmp; comm.Barrier();
81 
82  bool verbose = false;
83 
84  // Check if we should print results to standard out
85  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
86 
87  if (!verbose) comm.SetTracebackMode(0); // This should shut down any error traceback reporting
88 
89  if (verbose) std::cout << comm << std::endl << std::flush;
90 
91  if (verbose) verbose = (comm.MyPID()==0);
92 
93  if (verbose)
94  std::cout << EpetraExt::EpetraExt_Version() << std::endl << std::endl;
95 
96  int nx = 128;
97  int ny = comm.NumProc()*nx; // Scale y grid with number of processors
98 
99  // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
100 
101  // (i-1,j-1) (i-1,j )
102  // (i ,j-1) (i ,j ) (i ,j+1)
103  // (i+1,j-1) (i+1,j )
104 
105  int npoints = 7;
106 
107  int xoff[] = {-1, 0, 1, -1, 0, 1, 0};
108  int yoff[] = {-1, -1, -1, 0, 0, 0, 1};
109 
110  Epetra_Map * map;
112  Epetra_Vector * x, * b, * xexact;
113 
114  Trilinos_Util_GenerateCrsProblem64(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
115 
116  if (nx<8)
117  {
118  std::cout << *A << std::endl;
119  std::cout << "X exact = " << std::endl << *xexact << std::endl;
120  std::cout << "B = " << std::endl << *b << std::endl;
121  }
122 
123  // Construct transposer
124  Epetra_Time timer(comm);
125 
126  double start = timer.ElapsedTime();
127 
128  //bool IgnoreNonLocalCols = false;
130 
131  if (verbose) std::cout << "\nTime to construct transposer = " << timer.ElapsedTime() - start << std::endl;
132 
133  Epetra_CrsMatrix & transA = dynamic_cast<Epetra_CrsMatrix&>(transposer(*A));
134 
135  start = timer.ElapsedTime();
136  if (verbose) std::cout << "\nTime to create transpose matrix = " << timer.ElapsedTime() - start << std::endl;
137 
138  // Now test output of transposer by performing matvecs
139  int ierr = 0;
140  ierr += checkResults(A, &transA, xexact, verbose);
141 
142 
143  // Now change values in original matrix and test update facility of transposer
144  // Add 2 to the diagonal of each row
145  double Value = 2.0;
146  for (i=0; i< A->NumMyRows(); i++)
147  A->SumIntoMyValues(i, 1, &Value, &i);
148 
149  start = timer.ElapsedTime();
150  transposer.fwd();
151 
152  if (verbose) std::cout << "\nTime to update transpose matrix = " << timer.ElapsedTime() - start << std::endl;
153 
154  ierr += checkResults(A, &transA, xexact, verbose);
155 
156  delete A;
157  delete b;
158  delete x;
159  delete xexact;
160  delete map;
161 
162  if (verbose) std::cout << std::endl << "Checking transposer for VbrMatrix objects" << std::endl<< std::endl;
163 
164 // CJ TODO FIXME: Trilinos_Util_GenerateVbrProblem64 not yet defined since vbr matrices not converted to 64 bit.
165 #if 0
166  int nsizes = 4;
167  int sizes[] = {4, 6, 5, 3};
168 
169  Epetra_VbrMatrix * Avbr;
170  Epetra_BlockMap * bmap;
171 
172  Trilinos_Util_GenerateVbrProblem64(nx, ny, npoints, xoff, yoff, nsizes, sizes,
173  comm, bmap, Avbr, x, b, xexact);
174  if (nx<8)
175  {
176  std::cout << *Avbr << std::endl;
177  std::cout << "X exact = " << std::endl << *xexact << std::endl;
178  std::cout << "B = " << std::endl << *b << std::endl;
179  }
180 
181  start = timer.ElapsedTime();
182  EpetraExt::RowMatrix_Transpose transposer1;
183 
184  Epetra_CrsMatrix & transA1 = dynamic_cast<Epetra_CrsMatrix&>(transposer1(*Avbr));
185  if (verbose) std::cout << "\nTime to create transpose matrix = " << timer.ElapsedTime() - start << std::endl;
186 
187  // Now test output of transposer by performing matvecs
188 ;
189  ierr += checkResults(Avbr, &transA1, xexact, verbose);
190 
191  // Now change values in original matrix and test update facility of transposer
192  // Scale matrix on the left by rowsums
193 
194  Epetra_Vector invRowSums(Avbr->RowMap());
195 
196  Avbr->InvRowSums(invRowSums);
197  Avbr->LeftScale(invRowSums);
198 
199  start = timer.ElapsedTime();
200  transposer1.fwd();
201  if (verbose) std::cout << "\nTime to update transpose matrix = " << timer.ElapsedTime() - start << std::endl;
202 
203  ierr += checkResults(Avbr, &transA1, xexact, verbose);
204 
205  delete Avbr;
206  delete bmap;
207  delete b;
208  delete x;
209  delete xexact;
210 #endif
211 
212 #ifdef EPETRA_MPI
213  MPI_Finalize();
214 #endif
215 
216  return ierr;
217 }
218 
220  Epetra_Vector * xexact, bool verbose) {
221 
222  long long n = A->NumGlobalRows64();
223 
224  if (n<100) std::cout << "A transpose = " << std::endl << *transA << std::endl;
225 
226  Epetra_Vector x1(View,A->OperatorDomainMap(), &((*xexact)[0]));
228 
229  A->SetUseTranspose(true);
230 
231  Epetra_Time timer(A->Comm());
232  double start = timer.ElapsedTime();
233  A->Apply(x1, b1);
234  if (verbose) std::cout << "\nTime to compute b1: matvec with original matrix using transpose flag = " << timer.ElapsedTime() - start << std::endl;
235 
236  if (n<100) std::cout << "b1 = " << std::endl << b1 << std::endl;
237  Epetra_Vector x2(View,transA->OperatorRangeMap(), &((*xexact)[0]));
238  Epetra_Vector b2(transA->OperatorDomainMap());
239  start = timer.ElapsedTime();
240  transA->Multiply(false, x2, b2);
241  if (verbose) std::cout << "\nTime to compute b2: matvec with transpose matrix = " << timer.ElapsedTime() - start << std::endl;
242 
243  if (n<100) std::cout << "b1 = " << std::endl << b1 << std::endl;
244 
245  double residual;
246  Epetra_Vector resid(A->OperatorDomainMap());
247 
248  resid.Update(1.0, b1, -1.0, b2, 0.0);
249  resid.Norm2(&residual);
250  if (verbose) std::cout << "Norm of b1 - b2 = " << residual << std::endl;
251 
252  int ierr = 0;
253 
254  if (residual > 1.0e-10) ierr++;
255 
256  if (ierr!=0 && verbose) std::cerr << "Status: Test failed" << std::endl;
257  else if (verbose) std::cerr << "Status: Test passed" << std::endl;
258 
259  return(ierr);
260 }
virtual int SetUseTranspose(bool UseTranspose)=0
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
double ElapsedTime(void) const
std::string EpetraExt_Version()
virtual const Epetra_Map & OperatorDomainMap() const =0
virtual int LeftScale(const Epetra_Vector &x)=0
int MyPID() const
const Epetra_Map & OperatorDomainMap() const
Transform to form the explicit transpose of a Epetra_RowMatrix.
virtual int InvRowSums(Epetra_Vector &x) const =0
int main(int argc, char **argv)
Definition: HDF5_IO.cpp:67
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
virtual const Epetra_Comm & Comm() const =0
int NumMyRows() const
int checkResults(Epetra_RowMatrix *A, Epetra_CrsMatrix *transA, Epetra_Vector *xexact, bool verbose)
int NumProc() const
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_Map & OperatorRangeMap() const
int n
virtual long long NumGlobalRows64() const =0