EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/Transpose/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_GenerateCrsProblem(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 int nsizes = 4;
165 int sizes[] = {4, 6, 5, 3};
166
167 Epetra_VbrMatrix * Avbr;
168 Epetra_BlockMap * bmap;
169
170 Trilinos_Util_GenerateVbrProblem(nx, ny, npoints, xoff, yoff, nsizes, sizes,
171 comm, bmap, Avbr, x, b, xexact);
172
173 if (nx<8)
174 {
175 std::cout << *Avbr << std::endl;
176 std::cout << "X exact = " << std::endl << *xexact << std::endl;
177 std::cout << "B = " << std::endl << *b << std::endl;
178 }
179
180 start = timer.ElapsedTime();
182
183 Epetra_CrsMatrix & transA1 = dynamic_cast<Epetra_CrsMatrix&>(transposer1(*Avbr));
184 if (verbose) std::cout << "\nTime to create transpose matrix = " << timer.ElapsedTime() - start << std::endl;
185
186 // Now test output of transposer by performing matvecs
187;
188 ierr += checkResults(Avbr, &transA1, xexact, verbose);
189
190 // Now change values in original matrix and test update facility of transposer
191 // Scale matrix on the left by rowsums
192
193 Epetra_Vector invRowSums(Avbr->RowMap());
194
195 Avbr->InvRowSums(invRowSums);
196 Avbr->LeftScale(invRowSums);
197
198 start = timer.ElapsedTime();
199 transposer1.fwd();
200 if (verbose) std::cout << "\nTime to update transpose matrix = " << timer.ElapsedTime() - start << std::endl;
201
202 ierr += checkResults(Avbr, &transA1, xexact, verbose);
203
204 delete Avbr;
205 delete b;
206 delete x;
207 delete xexact;
208 delete bmap;
209
210#ifdef EPETRA_MPI
211 MPI_Finalize();
212#endif
213
214 return ierr;
215}
216
218 Epetra_Vector * xexact, bool verbose) {
219
220 int n = A->NumGlobalRows();
221
222 if (n<100) std::cout << "A transpose = " << std::endl << *transA << std::endl;
223
224 Epetra_Vector x1(View,A->OperatorDomainMap(), &((*xexact)[0]));
226
227 A->SetUseTranspose(true);
228
229 Epetra_Time timer(A->Comm());
230 double start = timer.ElapsedTime();
231 A->Apply(x1, b1);
232 if (verbose) std::cout << "\nTime to compute b1: matvec with original matrix using transpose flag = " << timer.ElapsedTime() - start << std::endl;
233
234 if (n<100) std::cout << "b1 = " << std::endl << b1 << std::endl;
235 Epetra_Vector x2(View,transA->OperatorRangeMap(), &((*xexact)[0]));
236 Epetra_Vector b2(transA->OperatorDomainMap());
237 start = timer.ElapsedTime();
238 transA->Multiply(false, x2, b2);
239 if (verbose) std::cout << "\nTime to compute b2: matvec with transpose matrix = " << timer.ElapsedTime() - start << std::endl;
240
241 if (n<100) std::cout << "b1 = " << std::endl << b1 << std::endl;
242
243 double residual;
245
246 resid.Update(1.0, b1, -1.0, b2, 0.0);
247 resid.Norm2(&residual);
248 if (verbose) std::cout << "Norm of b1 - b2 = " << residual << std::endl;
249
250 int ierr = 0;
251
252 if (residual > 1.0e-10) ierr++;
253
254 if (ierr!=0 && verbose) std::cerr << "Status: Test failed" << std::endl;
255 else if (verbose) std::cerr << "Status: Test passed" << std::endl;
256
257 return(ierr);
258}
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 int NumGlobalRows() 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)