Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
bug1_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// This file demonstrates a problem in some destructor somewhere - or
47// perhaps the proble could be in this test code.
48//
49// In any case, the symptom is a Segmentation fualt
50//
51
52
53#ifdef EPETRA_MPI
54#include "Epetra_MpiComm.h"
55#include <mpi.h>
56#endif
57#include "Epetra_CrsMatrix.h"
58#include "Epetra_VbrMatrix.h"
59#include "Epetra_Vector.h"
60#include "Epetra_MultiVector.h"
61#include "Epetra_LocalMap.h"
62#include "Epetra_IntVector.h"
63#include "Epetra_Map.h"
64
65#include "Epetra_SerialComm.h"
66#include "Epetra_Time.h"
69#include "Trilinos_Util.h"
70#ifdef HAVE_COMM_ASSERT_EQUAL
71#include "Comm_assert_equal.h"
72#endif
73
74int checkResults( bool trans,
76 Epetra_LinearProblem * OrigProb,
77 Epetra_LinearProblem * RedistProb,
78 bool verbose) ;
79
80int main(int argc, char *argv[]) {
81
82{
83 int i;
84
85#ifdef EPETRA_MPI
86 // Initialize MPI
87 MPI_Init(&argc,&argv);
88 Epetra_MpiComm comm(MPI_COMM_WORLD);
89#else
91#endif
92
93 // Uncomment to debug in parallel int tmp; if (comm.MyPID()==0) cin >> tmp; comm.Barrier();
94
95 bool verbose = false;
96 bool veryVerbose = false;
97
98 // Check if we should print results to standard out
99 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
100
101 if (!verbose) comm.SetTracebackMode(0); // This should shut down any error traceback reporting
102
103 if (verbose) cout << comm << endl << flush;
104
105 bool verbose1 = verbose;
106 if (verbose) verbose = (comm.MyPID()==0);
107
108 int nx = 4;
109 int ny = comm.NumProc()*nx; // Scale y grid with number of processors
110
111 // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
112
113 // (i-1,j-1) (i-1,j )
114 // (i ,j-1) (i ,j ) (i ,j+1)
115 // (i+1,j-1) (i+1,j )
116
117 int npoints = 2;
118
119 int xoff[] = {-1, 0, 1, -1, 0, 1, 0};
120 int yoff[] = {-1, -1, -1, 0, 0, 0, 1};
121
122 Epetra_Map * map;
124 Epetra_Vector * x, * b, * xexact;
125
126 Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
127
128 if (verbose)
129 cout << "npoints = " << npoints << " nx = " << nx << " ny = " << ny << endl ;
130
131 if (verbose && nx<6 ) {
132 cout << *A << endl;
133 cout << "B = " << endl << *b << endl;
134 }
135 // Construct linear problem object
136
137 Epetra_LinearProblem origProblem(A, x, b);
138 Epetra_LinearProblem *redistProblem;
139
140 Epetra_Time timer(comm);
141
142 // Construct redistor object, use all processors and replicate full problem on each
143
144 double start = timer.ElapsedTime();
145 Epetra_LinearProblemRedistor redistor(&origProblem, comm.NumProc(), true);
146 if (verbose) cout << "\nTime to construct redistor = "
147 << timer.ElapsedTime() - start << endl;
148
149 bool ConstructTranspose = true;
150 bool MakeDataContiguous = true;
151
152 //
153 // Commenting out the following define avoids the segmentation fault
154 //
155#define CREATEREDISTPROBLEM
156#ifdef CREATEREDISTPROBLEM
157 start = timer.ElapsedTime();
158 redistor.CreateRedistProblem(ConstructTranspose, MakeDataContiguous, redistProblem);
159 if (verbose) cout << "\nTime to create redistributed problem = "
160 << timer.ElapsedTime() - start << endl;
161
162#endif
163#if 0
164 // Now test output of redistor by performing matvecs
165
166 int ierr = 0;
167 ierr += checkResults( ConstructTranspose, &redistor, &origProblem,
168 redistProblem, verbose);
169#endif
170 cout << " Here we are just before the return() " << endl ;
171}
172 return(0) ;
173
174}
175
176//
177// checkResults computes Ax = b1 and either Rx=b2 or R^T x=b2 (depending on trans) where
178// A = origProblem and R = redistProblem.
179// checkResults returns 0 (OK) if norm(b1-b2)/norm(b1) < size(b1) * 1e-15;
180//
181// I guess that we need the redistor so that we can move x and b2 around.
182//
183const double error_tolerance = 1e-15 ;
184
185int checkResults( bool trans,
189 bool verbose) {
190
191 int m = A->GetRHS()->MyLength();
192 int n = A->GetLHS()->MyLength();
193 assert( m == n ) ;
194
195 Epetra_MultiVector *x = A->GetLHS() ;
196 Epetra_MultiVector x1( *x ) ;
197 // Epetra_MultiVector Difference( x1 ) ;
198 Epetra_MultiVector *b = A->GetRHS();
199 Epetra_RowMatrix *matrixA = A->GetMatrix();
200 assert( matrixA != 0 ) ;
201 int iam = matrixA->Comm().MyPID();
202
203 // Epetra_Time timer(A->Comm());
204 // double start = timer.ElapsedTime();
205
206 matrixA->Multiply(trans, *b, x1) ; // x = Ab
207
208 int M,N,nz;
209 int *ptr, *ind;
210 double *val, *rhs, *lhs;
211 int Nrhs, ldrhs, ldlhs;
212
213 redistor->ExtractHbData( M, N, nz, ptr, ind,
214 val, Nrhs, rhs, ldrhs,
215 lhs, ldlhs);
216
217 assert( M == N ) ;
218 if ( verbose ) {
219 cout << " iam = " << iam
220 << " m = " << m << " n = " << n << " M = " << M << endl ;
221
222 cout << " iam = " << iam << " ptr = " << ptr[0] << " " << ptr[1] << " " << ptr[2] << " " << ptr[3] << " " << ptr[4] << " " << ptr[5] << endl ;
223
224 cout << " iam = " << iam << " ind = " << ind[0] << " " << ind[1] << " " << ind[2] << " " << ind[3] << " " << ind[4] << " " << ind[5] << endl ;
225
226 cout << " iam = " << iam << " val = " << val[0] << " " << val[1] << " " << val[2] << " " << val[3] << " " << val[4] << " " << val[5] << endl ;
227 }
228 // Create a serial map in case we end up needing it
229 // If it is created inside the else block below it would have to
230 // be with a call to new().
231 int NumMyElements_ = 0 ;
232 if (matrixA->Comm().MyPID()==0) NumMyElements_ = n;
233 Epetra_Map SerialMap( n, NumMyElements_, 0, matrixA->Comm() );
234
235 // These are unnecessary and useless
236 // Epetra_Vector serial_A_rhs( SerialMap ) ;
237 // Epetra_Vector serial_A_lhs( SerialMap ) ;
238
239 // Epetra_Export exporter( matrixA->BlockRowMap(), SerialMap) ;
240
241 //
242 // In each process, we will compute Rb putting the answer into LHS
243 //
244
245
246 for ( int k = 0 ; k < Nrhs; k ++ ) {
247 for ( int i = 0 ; i < M ; i ++ ) {
248 lhs[ i + k * ldlhs ] = 0.0;
249 }
250 for ( int i = 0 ; i < M ; i++ ) {
251 for ( int l = ptr[i]; l < ptr[i+1]; l++ ) {
252 int j = ind[l] ;
253 if ( verbose && N < 40 ) {
254 cout << " i = " << i << " j = " << j ;
255 cout << " l = " << l << " val[l] = " << val[l] ;
256 cout << " rhs = " << rhs[ j + k * ldrhs ] << endl ;
257 }
258 lhs[ i + k * ldrhs ] += val[l] * rhs[ j + k * ldrhs ] ;
259 }
260 }
261
262 if ( verbose && N < 40 ) {
263 cout << " lhs = " ;
264 for ( int j = 0 ; j < N ; j++ ) cout << " " << lhs[j] ;
265 cout << endl ;
266 cout << " rhs = " ;
267 for ( int j = 0 ; j < N ; j++ ) cout << " " << rhs[j] ;
268 cout << endl ;
269 }
270
271 const Epetra_Comm &comm = matrixA->Comm() ;
272#ifdef HAVE_COMM_ASSERT_EQUAL
273 //
274 // Here we double check to make sure that lhs and rhs are
275 // replicated.
276 //
277 for ( int j = 0 ; j < N ; j++ ) {
278 assert( Comm_assert_equal( &comm, lhs[ j + k * ldrhs ] ) ) ;
279 assert( Comm_assert_equal( &comm, rhs[ j + k * ldrhs ] ) ) ;
280 }
281#endif
282 }
283
284 //
285 // Now we have to redistribue them back
286 //
287 redistor->UpdateOriginalLHS( A->GetLHS() ) ;
288 //
289 // Now we want to compare x and x1 which have been computed as follows:
290 // x = Rb
291 // x1 = Ab
292 //
293
294 double Norm_x1, Norm_diff ;
295 EPETRA_CHK_ERR( x1.Norm2( &Norm_x1 ) ) ;
296
297 // cout << " x1 = " << x1 << endl ;
298 // cout << " *x = " << *x << endl ;
299
300 x1.Update( -1.0, *x, 1.0 ) ;
301 EPETRA_CHK_ERR( x1.Norm2( &Norm_diff ) ) ;
302
303 // cout << " diff, i.e. updated x1 = " << x1 << endl ;
304
305 int ierr = 0;
306
307 if ( verbose ) {
308 cout << " Norm_diff = " << Norm_diff << endl ;
309 cout << " Norm_x1 = " << Norm_x1 << endl ;
310 }
311
312 if ( Norm_diff / Norm_x1 > n * error_tolerance ) ierr++ ;
313
314 if (ierr!=0 && verbose) cerr << "Status: Test failed" << endl;
315 else if (verbose) cerr << "Status: Test passed" << endl;
316
317 return(ierr);
318}
#define EPETRA_CHK_ERR(a)
const double error_tolerance
int main(int argc, char *argv[])
Definition bug1_main.cpp:80
int checkResults(bool trans, Epetra_LinearProblemRedistor *redistor, Epetra_LinearProblem *OrigProb, Epetra_LinearProblem *RedistProb, bool verbose)
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...
Epetra_LinearProblemRedistor: A class for redistributing an Epetra_LinearProblem object.
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 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 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.