Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/LinearProblemRedistor/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// Epetra_LinearProblemRedistor Test routine
44
45//
46// Limitations
47// -----------
48//
49// We test only square matrices. Indeed we perform this test only on
50// a single square (but non-symmetric) matrix.
51//
52// We test only one type of redistor: the one needed by
53// SuperLUdist. i.e. ConstructTranspose = true -and-
54// MakeDataContiguous = true.
55//
56// We exercise only one constructor, again the constructor required by
57// SuperLUdist. i.e. OrigProblem, NumProc, Replicate with
58// NumProc = comm.NumProc() and Replicate = true.
59//
60// Tests performed
61// ---------------
62//
63// We create three Epetra_LinearProblems and redistributions:
64// 1) Create an Epetra_LinearProblem and redistribute it
65// 2) Update the right hand side of an Epetra_LinearProblem and
66// update the redistributed version of the right hand side
67// 3) Update the matrix of an Epetra_LinearProblem and
68// update the redistributed version of the matrix
69// For each of these we compare the linear problem and its
70// redistribution as follows:
71//
72// Method
73// ------
74//
75// We compare a linear problem and the redistributed version of same
76// by multiplying the matrix in the linear problem and the redistributed
77// version of the linear problem by the right hand side to yield a
78// left hand side which we can then compare.
79//
80// We perform a matrix vector multiply instead of a solve to avoid
81// dependence on any particular solver. However, this results in
82// computing x = Ab instead of solving Ax =b. This only makes sense
83// if A is square.
84//
85
86
87#ifdef EPETRA_MPI
88#include "Epetra_MpiComm.h"
89#include <mpi.h>
90#endif
91#include "Epetra_CrsMatrix.h"
92#include "Epetra_VbrMatrix.h"
93#include "Epetra_Vector.h"
94#include "Epetra_MultiVector.h"
95#include "Epetra_LocalMap.h"
96#include "Epetra_IntVector.h"
97#include "Epetra_Map.h"
98
99#include "Epetra_SerialComm.h"
100#include "Epetra_Time.h"
101#include "Epetra_LinearProblem.h"
103#include "Trilinos_Util.h"
104#ifdef HAVE_COMM_ASSERT_EQUAL
105#include "Comm_assert_equal.h"
106#endif
107
108int checkResults( bool trans,
110 Epetra_LinearProblem * OrigProb,
111 Epetra_LinearProblem * RedistProb,
112 bool verbose) ;
113
114int main(int argc, char *argv[]) {
115
116 int i;
117
118#ifdef EPETRA_MPI
119 // Initialize MPI
120 MPI_Init(&argc,&argv);
121 Epetra_MpiComm comm(MPI_COMM_WORLD);
122#else
124#endif
125
126 // Uncomment to debug in parallel int tmp; if (comm.MyPID()==0) cin >> tmp; comm.Barrier();
127
128 bool verbose = false;
129 bool veryVerbose = false;
130
131 // Check if we should print results to standard out
132 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
133
134 if (!verbose) comm.SetTracebackMode(0); // This should shut down any error traceback reporting
135
136 if (verbose) cout << comm << endl << flush;
137
138 bool verbose1 = verbose;
139 if (verbose) verbose = (comm.MyPID()==0);
140
141 int nx = 4;
142 int ny = comm.NumProc()*nx; // Scale y grid with number of processors
143
144 // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
145
146 // (i-1,j-1) (i-1,j )
147 // (i ,j-1) (i ,j ) (i ,j+1)
148 // (i+1,j-1) (i+1,j )
149
150 int npoints = 2;
151
152 int xoff[] = {-1, 0, 1, -1, 0, 1, 0};
153 int yoff[] = {-1, -1, -1, 0, 0, 0, 1};
154
155 Epetra_Map * map;
157 Epetra_Vector * x, * b, * xexact;
158
159 Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
160
161 if (verbose)
162 cout << "npoints = " << npoints << " nx = " << nx << " ny = " << ny << endl ;
163
164 if (verbose && nx<6 ) {
165 cout << *A << endl;
166 cout << "B = " << endl << *b << endl;
167 }
168 // Construct linear problem object
169
170 Epetra_LinearProblem origProblem(A, x, b);
171 Epetra_LinearProblem *redistProblem;
172
173 Epetra_Time timer(comm);
174
175 // Construct redistor object, use all processors and replicate full problem on each
176
177 double start = timer.ElapsedTime();
178 Epetra_LinearProblemRedistor redistor(&origProblem, comm.NumProc(), true);
179 if (verbose) cout << "\nTime to construct redistor = "
180 << timer.ElapsedTime() - start << endl;
181
182 bool ConstructTranspose = true;
183 bool MakeDataContiguous = true;
184
185 start = timer.ElapsedTime();
186 redistor.CreateRedistProblem(ConstructTranspose, MakeDataContiguous, redistProblem);
187 if (verbose) cout << "\nTime to create redistributed problem = "
188 << timer.ElapsedTime() - start << endl;
189
190
191 // Now test output of redistor by performing matvecs
192
193 int ierr = 0;
194 ierr += checkResults( ConstructTranspose, &redistor, &origProblem,
195 redistProblem, verbose);
196
197
198 // Now change values in original rhs and test update facility of redistor
199 // Multiply b by 2
200
201 double Value = 2.0;
202
203 b->Scale(Value); // b = 2*b
204
205 redistor.UpdateRedistRHS(b);
206 if (verbose) cout << "\nTime to update redistributed RHS = "
207 << timer.ElapsedTime() - start << endl;
208
209 ierr += checkResults( ConstructTranspose, &redistor,
210 &origProblem, redistProblem, verbose);
211
212 // Now change values in original matrix and test update facility of redistor
213
214#define CREATE_CONST_MATRIX
215#ifdef CREATE_CONST_MATRIX
216 // The easiest way that I could find to change the matrix without EPETRA_CHK_ERRs
217 A->PutScalar(13.0);
218#else
219
220 // This has no effect on matrices, such as when nx = 4, that have no
221 // diagonal entries. However, it does cause many EPETRA_CHK_ERR prints.
222
223 // Add 2 to the diagonal of each row
224 for (i=0; i< A->NumMyRows(); i++) {
225 // for (i=0; i < 1; i++)
226 cout << " i = " << i ;
227 A->SumIntoMyValues(i, 1, &Value, &i);
228 }
229#endif
230
231
232 start = timer.ElapsedTime();
233 redistor.UpdateRedistProblemValues(&origProblem);
234 if (verbose) cout << "\nTime to update redistributed problem = "
235 << timer.ElapsedTime() - start << endl;
236
237 ierr += checkResults(ConstructTranspose, &redistor, &origProblem, redistProblem, verbose);
238
239 delete A;
240 delete b;
241 delete x;
242 delete xexact;
243 delete map;
244
245
246#ifdef EPETRA_MPI
247 MPI_Finalize();
248#endif
249
250 return ierr;
251}
252
253//
254// checkResults computes Ax = b1 and either Rx=b2 or R^T x=b2 (depending on trans) where
255// A = origProblem and R = redistProblem.
256// checkResults returns 0 (OK) if norm(b1-b2)/norm(b1) < size(b1) * 1e-15;
257//
258// I guess that we need the redistor so that we can move x and b2 around.
259//
260const double error_tolerance = 1e-15 ;
261
262int checkResults( bool trans,
266 bool verbose) {
267
268 int m = A->GetRHS()->MyLength();
269 int n = A->GetLHS()->MyLength();
270 assert( m == n ) ;
271
272 Epetra_MultiVector *x = A->GetLHS() ;
273 Epetra_MultiVector x1( *x ) ;
274 // Epetra_MultiVector Difference( x1 ) ;
275 Epetra_MultiVector *b = A->GetRHS();
276 Epetra_RowMatrix *matrixA = A->GetMatrix();
277 assert( matrixA != 0 ) ;
278 int iam = matrixA->Comm().MyPID();
279
280 // Epetra_Time timer(A->Comm());
281 // double start = timer.ElapsedTime();
282
283 matrixA->Multiply(trans, *b, x1) ; // x = Ab
284
285 int M,N,nz;
286 int *ptr, *ind;
287 double *val, *rhs, *lhs;
288 int Nrhs, ldrhs, ldlhs;
289
290 redistor->ExtractHbData( M, N, nz, ptr, ind,
291 val, Nrhs, rhs, ldrhs,
292 lhs, ldlhs);
293
294 assert( M == N ) ;
295 if ( verbose ) {
296 cout << " iam = " << iam
297 << " m = " << m << " n = " << n << " M = " << M << endl ;
298
299 cout << " iam = " << iam << " ptr = " << ptr[0] << " " << ptr[1] << " " << ptr[2] << " " << ptr[3] << " " << ptr[4] << " " << ptr[5] << endl ;
300
301 cout << " iam = " << iam << " ind = " << ind[0] << " " << ind[1] << " " << ind[2] << " " << ind[3] << " " << ind[4] << " " << ind[5] << endl ;
302
303 cout << " iam = " << iam << " val = " << val[0] << " " << val[1] << " " << val[2] << " " << val[3] << " " << val[4] << " " << val[5] << endl ;
304 }
305 // Create a serial map in case we end up needing it
306 // If it is created inside the else block below it would have to
307 // be with a call to new().
308 int NumMyElements_ = 0 ;
309 if (matrixA->Comm().MyPID()==0) NumMyElements_ = n;
310 Epetra_Map SerialMap( n, NumMyElements_, 0, matrixA->Comm() );
311
312 // These are unnecessary and useless
313 // Epetra_Vector serial_A_rhs( SerialMap ) ;
314 // Epetra_Vector serial_A_lhs( SerialMap ) ;
315
316 // Epetra_Export exporter( matrixA->BlockRowMap(), SerialMap) ;
317
318 //
319 // In each process, we will compute Rb putting the answer into LHS
320 //
321
322
323 for ( int k = 0 ; k < Nrhs; k ++ ) {
324 for ( int i = 0 ; i < M ; i ++ ) {
325 lhs[ i + k * ldlhs ] = 0.0;
326 }
327 for ( int i = 0 ; i < M ; i++ ) {
328 for ( int l = ptr[i]; l < ptr[i+1]; l++ ) {
329 int j = ind[l] ;
330 if ( verbose && N < 40 ) {
331 cout << " i = " << i << " j = " << j ;
332 cout << " l = " << l << " val[l] = " << val[l] ;
333 cout << " rhs = " << rhs[ j + k * ldrhs ] << endl ;
334 }
335 lhs[ i + k * ldrhs ] += val[l] * rhs[ j + k * ldrhs ] ;
336 }
337 }
338
339 if ( verbose && N < 40 ) {
340 cout << " lhs = " ;
341 for ( int j = 0 ; j < N ; j++ ) cout << " " << lhs[j] ;
342 cout << endl ;
343 cout << " rhs = " ;
344 for ( int j = 0 ; j < N ; j++ ) cout << " " << rhs[j] ;
345 cout << endl ;
346 }
347
348 const Epetra_Comm &comm = matrixA->Comm() ;
349#ifdef HAVE_COMM_ASSERT_EQUAL
350 //
351 // Here we double check to make sure that lhs and rhs are
352 // replicated.
353 //
354 for ( int j = 0 ; j < N ; j++ ) {
355 assert( Comm_assert_equal( &comm, lhs[ j + k * ldrhs ] ) ) ;
356 assert( Comm_assert_equal( &comm, rhs[ j + k * ldrhs ] ) ) ;
357 }
358#endif
359 }
360
361 //
362 // Now we have to redistribue them back
363 //
364 redistor->UpdateOriginalLHS( A->GetLHS() ) ;
365 //
366 // Now we want to compare x and x1 which have been computed as follows:
367 // x = Rb
368 // x1 = Ab
369 //
370
371 double Norm_x1, Norm_diff ;
372 EPETRA_CHK_ERR( x1.Norm2( &Norm_x1 ) ) ;
373
374 // cout << " x1 = " << x1 << endl ;
375 // cout << " *x = " << *x << endl ;
376
377 x1.Update( -1.0, *x, 1.0 ) ;
378 EPETRA_CHK_ERR( x1.Norm2( &Norm_diff ) ) ;
379
380 // cout << " diff, i.e. updated x1 = " << x1 << endl ;
381
382 int ierr = 0;
383
384 if ( verbose ) {
385 cout << " Norm_diff = " << Norm_diff << endl ;
386 cout << " Norm_x1 = " << Norm_x1 << endl ;
387 }
388
389 if ( Norm_diff / Norm_x1 > n * error_tolerance ) ierr++ ;
390
391 if (ierr!=0 && verbose) cerr << "Status: Test failed" << endl;
392 else if (verbose) cerr << "Status: Test passed" << endl;
393
394 return(ierr);
395}
#define EPETRA_CHK_ERR(a)
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
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 NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
Epetra_LinearProblemRedistor: A class for redistributing an Epetra_LinearProblem object.
int UpdateRedistRHS(Epetra_MultiVector *RHSWithNewValues)
Update the values of an already-redistributed RHS.
int CreateRedistProblem(const bool ConstructTranspose, const bool MakeDataContiguous, Epetra_LinearProblem *&RedistProblem)
Generate a new Epetra_LinearProblem as a redistribution of the one passed into the constructor.
int UpdateRedistProblemValues(Epetra_LinearProblem *ProblemWithNewValues)
Update the values of an already-redistributed problem.
int UpdateOriginalLHS(Epetra_MultiVector *LHS)
Update LHS of original Linear Problem object.
int ExtractHbData(int &M, int &N, int &nz, int *&ptr, int *&ind, double *&val, int &Nrhs, double *&rhs, int &ldrhs, double *&lhs, int &ldlhs) const
Extract the redistributed problem data in a form usable for other codes that require Harwell-Boeing f...
Epetra_LinearProblem: The Epetra Linear Problem Class.
Epetra_MultiVector * GetRHS() const
Get a pointer to the right-hand-side B.
Epetra_RowMatrix * GetMatrix() const
Get a pointer to the matrix A.
Epetra_MultiVector * GetLHS() const
Get a pointer to the left-hand-side X.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
const double error_tolerance
int main(int argc, char *argv[])
int checkResults(bool trans, Epetra_LinearProblemRedistor *redistor, Epetra_LinearProblem *OrigProb, Epetra_LinearProblem *RedistProb, bool verbose)