Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
PerformOneSolveAndTest.cpp
Go to the documentation of this file.
1//
2// OUR_CHK_ERR always returns 1 on error.
3//
4
5#define OUR_CHK_ERR(a) { { int epetra_err = a; \
6 if (epetra_err != 0) { std::cerr << "Amesos ERROR " << epetra_err << ", " \
7 << __FILE__ << ", line " << __LINE__ << std::endl; \
8relerror = 1.3e15; relresidual=1e15; return(1);} }\
9 }
10
11#include "Epetra_Comm.h"
12#include "Teuchos_ParameterList.hpp"
13#include "Amesos.h"
14#include "Epetra_CrsMatrix.h"
15#include "Epetra_Map.h"
16#include "Epetra_Vector.h"
17#include "Epetra_LinearProblem.h"
20#include "NewMatNewMap.h"
21#include "Amesos_TestRowMatrix.h"
22
23//
24// Returns the number of failures.
25// Note: If AmesosClass is not supported, PerformOneSolveAndTest() will
26// always return 0
27//
28// Still have to decide where we are going to check the residual.
29//
30// The following table shows the variable names that we use for
31// each of the three phases:
32// compute - which computes the correct value of b
33// solve - which solves for x in A' A' A x = b
34// check - which computes bcheck = A' A' A x
35//
36// For ill-conditioned matrices we restrict the test to one or two
37// solves, by setting Levels to 1 or 2 on input to this routine.
38// When Levels is less than 3, some of the transformations
39// shown in the table as "->" and "<-" are not performed, instead
40// a direct copy is made.
41//
42// In the absence of roundoff, each item in a given column should
43// be identical.
44//
45// If Levels = 3, we compute and solve A' A' A x = b and hence all
46// of the transformations are performed
47//
48// If Levels = 2, the transformations shown in the first column of
49// transformations (labelled Levels>=3) are replaced by a direct copy.
50//
51// If Levels = 1, only the transformations shown in the third column
52// are performed, the others being replaced by direct copies.
53//
54// Levels>=3 Levels>=2
55// A' A' A
56// compute xexact -> cAx -> cAAx -> b
57// solve x <- sAx <- sAAx <- b
58// check x -> kAx -> kAAx -> bcheck
59//
60// Note that since Levels 2 and 3 use the same A, there is no need to
61// call NumericFactorization() between the second and third call to Solve.
62//
63// On testing AddToDiag
64// ====================
65//
66// When this code was written, our intent was to have AddToDiag add a constant
67// value to the diagonal even for non-existent diagonal elements. For now, we have
68// backed off of that. If we decide to go back to the earlier semantics for
69// AddToDiag (i.e. that it would force all diagonal elements to exist in the
70// matrix, this can be tested by setting AddToAllDiagonalElements to true ;
71//
72// See bug #1928 for further discussion.
73//
74// ExpectedError
75// =============
76//
77// ExpectedError allows us to test for Singular and Structurally Singular matrices
78// ExpectedError = StructurallySingularMatrix or NumericallySingularMatrix
79
80
81int PerformOneSolveAndTest( const char* AmesosClass,
82 int EpetraMatrixType,
83 const Epetra_Comm &Comm,
84 bool transpose,
85 bool verbose,
86 Teuchos::ParameterList ParamList,
87 Epetra_CrsMatrix *& InMat,
88 int Levels,
89 const double Rcond,
90 double& relerror,
91 double& relresidual,
92 int ExpectedError )
93{
94#ifndef NDEBUG
95 int ierr = 0;
96#endif
97 bool AddToAllDiagonalElements = ParamList.get( "AddZeroToDiag", false ) ;
98
99
100 bool TrustMe = ParamList.get( "TrustMe", false );
101
102 bool MyVerbose = false ; // setting this equal to verbose produces too much output and exceeds the test harnesses 1 Megabyte limit
103
104 RCP<Epetra_CrsMatrix> MyMat ;
105 RCP<Epetra_CrsMatrix> MyMatWithDiag ;
106
107 MyMat = rcp( new Epetra_CrsMatrix( *InMat ) );
108
109 Amesos_TestRowMatrix ATRW( &*MyMat ) ;
110
111 Epetra_RowMatrix* MyRowMat = 0;
112
113 assert ( EpetraMatrixType >= 0 && EpetraMatrixType <= 2 );
114 switch ( EpetraMatrixType ) {
115 case 0:
116 MyRowMat = &*MyMat ;
117 break;
118 case 1:
119 MyRowMat = &ATRW;
120 break;
121 case 2:
122 MyMat->OptimizeStorage();
123 MyRowMat = &*MyMat ;
124#ifndef NDEBUG
125 bool OptStorage =
126#endif
127 MyMat->StorageOptimized();
128 assert( OptStorage) ;
129 break;
130 }
131#ifndef NDEBUG
132 bool OptStorage =
133#endif
134 MyMat->StorageOptimized();
135
136 Epetra_CrsMatrix* MatPtr = &*MyMat ;
137
138 const std::string AC = AmesosClass ;
139 if ( ExpectedError == 0 ) {
140 if ( AC != "Amesos_Pardiso" ) {
141 OUR_CHK_ERR ( PartialFactorization( AmesosClass, Comm, transpose, MyVerbose,
142 ParamList, MatPtr, Rcond ) );
143 } else {
144 if (MyVerbose) std::cout << " AC = " << AC << " not tested in Partial Factorization" <<std::endl ; // bug #1915
145 }
146 }
147
148 if (ParamList.isParameter("AddToDiag") ) {
149 int oldtracebackmode = InMat->GetTracebackMode( ) ;
150
151 //
152 // If AddToDiag is set, create a matrix which is numerically identical, but structurally
153 // has no missing diagaonal entries. In other words, every diagonal element in MyMayWithDiag
154 // has an entry in the matrix, though that entry will be zero if InMat has no entry for that
155 // particular diagonal element.
156 //
157 // It turns out that this does not actually make sure that all diagonal entries exist because it
158 // does not deal with maps that are missing the diagonal element.
159 //
160 if ( AddToAllDiagonalElements ) {
161 MyMatWithDiag = NewMatNewMap( *InMat, 2, 0, 0, 0, 0 ); // Ensure that all diagonal entries exist ;
162 } else {
163 InMat->SetTracebackMode( EPETRA_MIN(oldtracebackmode,1) ) ; // If the matrix is mimssing diagonal elements, the call to ReplaceDiagonalElements will return 1 indicating missing diagonal elements
164 MyMatWithDiag = NewMatNewMap( *InMat, 0, 0, 0, 0, 0 ); // Leave the matrix unchanged structurally
165 }
166 //
167 // Now add AddToDiag to each diagonal element.
168 //
169 double AddToDiag = ParamList.get("AddToDiag", 0.0 );
170 Epetra_Vector Diag( MyMatWithDiag->RowMap() );
171 Epetra_Vector AddConstVecToDiag( MyMatWithDiag->RowMap() );
172 AddConstVecToDiag.PutScalar( AddToDiag );
173
174#ifndef NDEBUG
175 ierr =
176#endif
177 MyMatWithDiag->ExtractDiagonalCopy( Diag );
178 assert( ierr == 0 );
179 Diag.Update( 1.0, AddConstVecToDiag, 1.0 ) ;
180#ifndef NDEBUG
181 ierr =
182#endif
183 MyMatWithDiag->ReplaceDiagonalValues( Diag ); // This may return 1 indicating that structurally non-zero elements were left untouched.
184 assert( ierr >= 0 );
185
186 InMat->SetTracebackMode( oldtracebackmode ) ;
187 } else {
188 MyMatWithDiag = rcp( new Epetra_CrsMatrix( *InMat ) );
189 }
190
191 if ( MyVerbose ) std::cout << " Partial Factorization complete " << std::endl ;
192
193 relerror = 0 ;
194 relresidual = 0 ;
195
196 assert( Levels >= 1 && Levels <= 3 ) ;
197
198 int iam = Comm.MyPID() ;
199 int errors = 0 ;
200
201 const Epetra_Map *RangeMap =
202 transpose?&MyMat->OperatorDomainMap():&MyMat->OperatorRangeMap() ;
203 const Epetra_Map *DomainMap =
204 transpose?&MyMat->OperatorRangeMap():&MyMat->OperatorDomainMap() ;
205
206 Epetra_Vector xexact(*DomainMap);
207 Epetra_Vector x(*DomainMap);
208
209 Epetra_Vector cAx(*DomainMap);
210 Epetra_Vector sAx(*DomainMap);
211 Epetra_Vector kAx(*DomainMap);
212
213 Epetra_Vector cAAx(*DomainMap);
214 Epetra_Vector sAAx(*DomainMap);
215 Epetra_Vector kAAx(*DomainMap);
216
217 Epetra_Vector FixedLHS(*DomainMap);
218 Epetra_Vector FixedRHS(*RangeMap);
219
220 Epetra_Vector b(*RangeMap);
221 Epetra_Vector bcheck(*RangeMap);
222
223 Epetra_Vector DomainDiff(*DomainMap);
224 Epetra_Vector RangeDiff(*RangeMap);
225
226 Epetra_LinearProblem Problem;
227 Teuchos::RCP<Amesos_BaseSolver> Abase ;
228 Amesos Afactory;
229
230
231
232
233 Abase = Teuchos::rcp(Afactory.Create( AmesosClass, Problem )) ;
234
235 relerror = 0 ;
236 relresidual = 0 ;
237
238 if ( Abase == Teuchos::null )
239 assert( false ) ;
240 else {
241
242 //
243 // Phase 1: Compute b = A' A' A xexact
244 //
245 Problem.SetOperator( MyRowMat );
246 // Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
247
248 //
249 // We only set transpose if we have to - this allows valgrind to check
250 // that transpose is set to a default value before it is used.
251 //
252 if ( transpose ) OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
253 if (MyVerbose) ParamList.set( "DebugLevel", 1 );
254 if (MyVerbose) ParamList.set( "OutputLevel", 1 );
255 OUR_CHK_ERR( Abase->SetParameters( ParamList ) );
256
257 if ( TrustMe ) {
258 Problem.SetLHS( &FixedLHS );
259 Problem.SetRHS( &FixedRHS );
260 assert( OptStorage) ;
261 }
262
263 // bug #2184
264 // Structurally singular matrices are not detected by
265 // Amesos_Klu::SymbolicFactorization() but they are detected by
266 // Amesos_Klu::NumericFactorization()
267 if ( ExpectedError == StructurallySingularMatrixError )
268 ExpectedError = NumericallySingularMatrixError ;
269 if ( ExpectedError == StructurallySingularMatrixError ) {
270 Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
271 int oldtracebackmode = ECM->GetTracebackMode( ) ;
272 ECM->SetTracebackMode(0); // We expect an error, but we don't want it to print out
273
274 const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
275 ECM->SetTracebackMode(oldtracebackmode);
276 // bug #2245 - Amesos fails to return error consistently across all
277 // processes. When this bug is fixed, remove "iam == 0 &&" from the next line
278 if ( iam == 0 && SymbolicFactorizationReturn != ExpectedError ) {
279 std::cout << " SymbolicFactorization returned " << SymbolicFactorizationReturn
280 << " should be " << ExpectedError << std::endl ;
281 OUR_CHK_ERR( 1 ) ;
282 } else {
283 return 0; // Returned the correct error code for this matrix
284 }
285 } else {
286 const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
287 OUR_CHK_ERR( SymbolicFactorizationReturn ) ;
288 }
289 if ( ExpectedError == NumericallySingularMatrixError ) {
290 Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
291 int oldtracebackmode = ECM->GetTracebackMode( ) ;
292 ECM->SetTracebackMode(0); // We expect an error, but we don't want it to print out
293
294 const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
295 ECM->SetTracebackMode(oldtracebackmode);
296 // bug #2245 - Amesos fails to return error consistently across all
297 // processes. When this bug is fixed, remove "iam == 0 &&" from the next line
298 if ( iam == 0 && NumericFactorizationReturn != ExpectedError ) {
299 std::cout << " NumericFactorization returned " << NumericFactorizationReturn
300 << " should be " << ExpectedError << std::endl ;
301 OUR_CHK_ERR( 1 ) ;
302 } else {
303 return 0; // Returned the correct error code for this matrix
304 }
305 } else {
306 const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
307 OUR_CHK_ERR( NumericFactorizationReturn ) ;
308 }
309 int ind[1];
310 double val[1];
311 ind[0] = 0;
312 xexact.Random();
313 xexact.PutScalar(1.0);
314
315 //
316 // Compute cAx = A' xexact
317 //
318 double Value = 1.0 ;
319 if ( Levels == 3 )
320 {
321 val[0] = Value ;
322 if ( MyMatWithDiag->MyGRID( 0 ) ) {
323 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
324 }
325 MyMatWithDiag->Multiply( transpose, xexact, cAx ) ;
326
327 val[0] = - Value ;
328 if ( MyMatWithDiag->MyGRID( 0 ) )
329 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
330 }
331 else
332 {
333 cAx = xexact ;
334 }
335
336 //
337 // Compute cAAx = A' cAx
338 //
339 if ( Levels >= 2 )
340 {
341 val[0] = Value ;
342 if ( MyMatWithDiag->MyGRID( 0 ) )
343 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
344 MyMatWithDiag->Multiply( transpose, cAx, cAAx ) ; // x2 = A' x1
345
346 val[0] = - Value ;
347 if ( MyMatWithDiag->MyGRID( 0 ) )
348 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
349 }
350 else
351 {
352 cAAx = cAx ;
353 }
354
355 if ( MyVerbose ) std::cout << " Compute b = A x2 = A A' A'' xexact " << std::endl ;
356
357 MyMatWithDiag->Multiply( transpose, cAAx, b ) ; // b = A x2 = A A' A'' xexact
358
359
360 // Phase 2: Solve A' A' A x = b
361 //
362 //
363 // Solve A sAAx = b
364 //
365 if ( TrustMe ) {
366 FixedRHS = b;
367 } else {
368 Problem.SetLHS( &sAAx );
369 Problem.SetRHS( &b );
370 }
371
372
373 OUR_CHK_ERR( Abase->SymbolicFactorization( ) );
374 OUR_CHK_ERR( Abase->SymbolicFactorization( ) ); // This should be irrelevant, but should nonetheless be legal
375 OUR_CHK_ERR( Abase->NumericFactorization( ) );
376 OUR_CHK_ERR( Abase->Solve( ) );
377 if ( TrustMe ) sAAx = FixedLHS ;
378
379 if ( Levels >= 2 )
380 {
381 OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
382 if ( TrustMe ) {
383 FixedRHS = sAAx ;
384 } else {
385 Problem.SetLHS( &sAx );
386 Problem.SetRHS( &sAAx );
387 }
388 val[0] = Value ;
389 if ( MyMat->MyGRID( 0 ) )
390 MyMat->SumIntoMyValues( 0, 1, val, ind ) ;
391 OUR_CHK_ERR( Abase->NumericFactorization( ) );
392 OUR_CHK_ERR( Abase->Solve( ) );
393 if ( TrustMe ) sAx = FixedLHS ;
394
395 }
396 else
397 {
398 sAx = sAAx ;
399 }
400
401 if ( Levels >= 3 )
402 {
403 if ( TrustMe ) {
404 FixedRHS = sAx ;
405 } else {
406 Problem.SetLHS( &x );
407 Problem.SetRHS( &sAx );
408 }
409 OUR_CHK_ERR( Abase->Solve( ) );
410 if ( TrustMe ) x = FixedLHS ;
411 }
412 else
413 {
414 x = sAx ;
415 }
416
417 if ( Levels >= 2 )
418 {
419 val[0] = -Value ;
420 if ( MyMat->MyGRID( 0 ) ) {
421 if ( MyMat->SumIntoMyValues( 0, 1, val, ind ) ) {
422 std::cout << " TestOptions requires a non-zero entry in A(1,1) " << std::endl ;
423 }
424 }
425 }
426
427 //
428 // Phase 3: Check the residual: bcheck = A' A' A x
429 //
430
431
432 if ( Levels >= 3 )
433 {
434 val[0] = Value ;
435 if ( MyMatWithDiag->MyGRID( 0 ) )
436 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
437 MyMatWithDiag->Multiply( transpose, x, kAx ) ;
438 val[0] = -Value ;
439 if ( MyMatWithDiag->MyGRID( 0 ) )
440 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
441 }
442 else
443 {
444 kAx = x ;
445 }
446
447 if ( Levels >= 2 )
448 {
449 val[0] = Value ;
450 if ( MyMatWithDiag->MyGRID( 0 ) )
451 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
452 MyMatWithDiag->Multiply( transpose, kAx, kAAx ) ;
453 val[0] = -Value ;
454 if ( MyMatWithDiag->MyGRID( 0 ) )
455 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
456 }
457 else
458 {
459 kAAx = kAx ;
460 }
461
462 MyMatWithDiag->Multiply( transpose, kAAx, bcheck ) ; // temp = A" x2
463
464 double norm_diff ;
465 double norm_one ;
466
467 DomainDiff.Update( 1.0, sAAx, -1.0, cAAx, 0.0 ) ;
468 DomainDiff.Norm2( &norm_diff ) ;
469 sAAx.Norm2( &norm_one ) ;
470 if (MyVerbose) std::cout << __FILE__ << "::" << __LINE__
471 << " norm( sAAx - cAAx ) / norm(sAAx ) = "
472 << norm_diff /norm_one << std::endl ;
473
474
475
476
477 DomainDiff.Update( 1.0, sAx, -1.0, cAx, 0.0 ) ;
478 DomainDiff.Norm2( &norm_diff ) ;
479 sAx.Norm2( &norm_one ) ;
480 if (MyVerbose) std::cout
481 << __FILE__ << "::" << __LINE__
482 << " norm( sAx - cAx ) / norm(sAx ) = "
483 << norm_diff /norm_one << std::endl ;
484
485
486 DomainDiff.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
487 DomainDiff.Norm2( &norm_diff ) ;
488 x.Norm2( &norm_one ) ;
489 if (MyVerbose) std::cout
490 << __FILE__ << "::" << __LINE__
491 << " norm( x - xexact ) / norm(x) = "
492 << norm_diff /norm_one << std::endl ;
493
494 relerror = norm_diff / norm_one ;
495
496 DomainDiff.Update( 1.0, sAx, -1.0, kAx, 0.0 ) ;
497 DomainDiff.Norm2( &norm_diff ) ;
498 sAx.Norm2( &norm_one ) ;
499 if (MyVerbose) std::cout
500 << __FILE__ << "::" << __LINE__
501 << " norm( sAx - kAx ) / norm(sAx ) = "
502 << norm_diff /norm_one << std::endl ;
503
504
505 DomainDiff.Update( 1.0, sAAx, -1.0, kAAx, 0.0 ) ;
506 DomainDiff.Norm2( &norm_diff ) ;
507 sAAx.Norm2( &norm_one ) ;
508 if (MyVerbose) std::cout
509 << __FILE__ << "::" << __LINE__
510 << " norm( sAAx - kAAx ) / norm(sAAx ) = "
511 << norm_diff /norm_one << std::endl ;
512
513
514 RangeDiff.Update( 1.0, bcheck, -1.0, b, 0.0 ) ;
515 RangeDiff.Norm2( &norm_diff ) ;
516 bcheck.Norm2( &norm_one ) ;
517 if (MyVerbose) std::cout
518 << __FILE__ << "::" << __LINE__
519 << " norm( bcheck - b ) / norm(bcheck ) = "
520 << norm_diff /norm_one << std::endl ;
521
522 relresidual = norm_diff / norm_one ;
523
524 if (iam == 0 ) {
525 if ( relresidual * Rcond < 1e-16 ) {
526 if (MyVerbose) std::cout << " Test 1 Passed " << std::endl ;
527 } else {
528 std::cout << __FILE__ << "::" << __LINE__ <<
529 " relresidual = " << relresidual <<
530 " TEST FAILED " <<
531 " ParamList = " << ParamList << std::endl ;
532 errors += 1 ;
533 }
534 }
535 }
536
537 return errors;
538
539}
540
541
static bool verbose
Definition Amesos.cpp:67
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
int PartialFactorization(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond)
int PerformOneSolveAndTest(const char *AmesosClass, int EpetraMatrixType, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&InMat, int Levels, const double Rcond, double &relerror, double &relresidual, int ExpectedError)
#define OUR_CHK_ERR(a)
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition Amesos.h:44