EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
68int 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
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();
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;
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}
View
Transform to form the explicit transpose of a Epetra_RowMatrix.
const Epetra_Map & OperatorRangeMap() const
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
int NumMyRows() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_Map & OperatorDomainMap() const
int NumProc() const
int MyPID() const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Norm2(double *Result) const
static void SetTracebackMode(int TracebackModeValue)
virtual int SetUseTranspose(bool UseTranspose)=0
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Comm & Comm() const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
virtual long long NumGlobalRows64() const=0
double ElapsedTime(void) const
int InvRowSums(Epetra_Vector &x) const
int LeftScale(const Epetra_Vector &x)
const Epetra_BlockMap & RowMap() const
std::string EpetraExt_Version()
int main(int argc, char *argv[])
int checkResults(Epetra_RowMatrix *A, Epetra_CrsMatrix *transA, Epetra_Vector *xexact, bool verbose)