Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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; \
8 relerror = 1.3e15; relresidual=1e15; return(1);} }\
9  }
10 
11 #include "Epetra_Comm.h"
13 #include "Amesos.h"
14 #include "Epetra_CrsMatrix.h"
15 #include "Epetra_Map.h"
16 #include "Epetra_Vector.h"
17 #include "Epetra_LinearProblem.h"
18 #include "PerformOneSolveAndTest.h"
19 #include "PartialFactorization.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 
81 int 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 
95  int ierr = 0;
96 
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  bool OptStorage = MyMat->StorageOptimized();
125  assert( OptStorage) ;
126  break;
127  }
128  bool OptStorage = MyMat->StorageOptimized();
129 
130  Epetra_CrsMatrix* MatPtr = &*MyMat ;
131 
132  const std::string AC = AmesosClass ;
133  if ( ExpectedError == 0 ) {
134  if ( AC != "Amesos_Pardiso" ) {
135  OUR_CHK_ERR ( PartialFactorization( AmesosClass, Comm, transpose, MyVerbose,
136  ParamList, MatPtr, Rcond ) );
137  } else {
138  if (MyVerbose) std::cout << " AC = " << AC << " not tested in Partial Factorization" <<std::endl ; // bug #1915
139  }
140  }
141 
142  if (ParamList.isParameter("AddToDiag") ) {
143  int oldtracebackmode = InMat->GetTracebackMode( ) ;
144 
145  //
146  // If AddToDiag is set, create a matrix which is numerically identical, but structurally
147  // has no missing diagaonal entries. In other words, every diagonal element in MyMayWithDiag
148  // has an entry in the matrix, though that entry will be zero if InMat has no entry for that
149  // particular diagonal element.
150  //
151  // It turns out that this does not actually make sure that all diagonal entries exist because it
152  // does not deal with maps that are missing the diagonal element.
153  //
154  if ( AddToAllDiagonalElements ) {
155  MyMatWithDiag = NewMatNewMap( *InMat, 2, 0, 0, 0, 0 ); // Ensure that all diagonal entries exist ;
156  } else {
157  InMat->SetTracebackMode( EPETRA_MIN(oldtracebackmode,1) ) ; // If the matrix is mimssing diagonal elements, the call to ReplaceDiagonalElements will return 1 indicating missing diagonal elements
158  MyMatWithDiag = NewMatNewMap( *InMat, 0, 0, 0, 0, 0 ); // Leave the matrix unchanged structurally
159  }
160  //
161  // Now add AddToDiag to each diagonal element.
162  //
163  double AddToDiag = ParamList.get("AddToDiag", 0.0 );
164  Epetra_Vector Diag( MyMatWithDiag->RowMap() );
165  Epetra_Vector AddConstVecToDiag( MyMatWithDiag->RowMap() );
166  AddConstVecToDiag.PutScalar( AddToDiag );
167 
168  ierr = MyMatWithDiag->ExtractDiagonalCopy( Diag );
169  assert( ierr == 0 );
170  Diag.Update( 1.0, AddConstVecToDiag, 1.0 ) ;
171  ierr = MyMatWithDiag->ReplaceDiagonalValues( Diag ); // This may return 1 indicating that structurally non-zero elements were left untouched.
172  assert( ierr >= 0 );
173 
174  InMat->SetTracebackMode( oldtracebackmode ) ;
175  } else {
176  MyMatWithDiag = rcp( new Epetra_CrsMatrix( *InMat ) );
177  }
178 
179  if ( MyVerbose ) std::cout << " Partial Factorization complete " << std::endl ;
180 
181  relerror = 0 ;
182  relresidual = 0 ;
183 
184  assert( Levels >= 1 && Levels <= 3 ) ;
185 
186  int iam = Comm.MyPID() ;
187  int errors = 0 ;
188 
189  const Epetra_Map *RangeMap =
190  transpose?&MyMat->OperatorDomainMap():&MyMat->OperatorRangeMap() ;
191  const Epetra_Map *DomainMap =
192  transpose?&MyMat->OperatorRangeMap():&MyMat->OperatorDomainMap() ;
193 
194  Epetra_Vector xexact(*DomainMap);
195  Epetra_Vector x(*DomainMap);
196 
197  Epetra_Vector cAx(*DomainMap);
198  Epetra_Vector sAx(*DomainMap);
199  Epetra_Vector kAx(*DomainMap);
200 
201  Epetra_Vector cAAx(*DomainMap);
202  Epetra_Vector sAAx(*DomainMap);
203  Epetra_Vector kAAx(*DomainMap);
204 
205  Epetra_Vector FixedLHS(*DomainMap);
206  Epetra_Vector FixedRHS(*RangeMap);
207 
208  Epetra_Vector b(*RangeMap);
209  Epetra_Vector bcheck(*RangeMap);
210 
211  Epetra_Vector DomainDiff(*DomainMap);
212  Epetra_Vector RangeDiff(*RangeMap);
213 
214  Epetra_LinearProblem Problem;
216  Amesos Afactory;
217 
218 
219 
220 
221  Abase = Teuchos::rcp(Afactory.Create( AmesosClass, Problem )) ;
222 
223  relerror = 0 ;
224  relresidual = 0 ;
225 
226  if ( Abase == Teuchos::null )
227  assert( false ) ;
228  else {
229 
230  //
231  // Phase 1: Compute b = A' A' A xexact
232  //
233  Problem.SetOperator( MyRowMat );
234  // Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
235 
236  //
237  // We only set transpose if we have to - this allows valgrind to check
238  // that transpose is set to a default value before it is used.
239  //
240  if ( transpose ) OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
241  if (MyVerbose) ParamList.set( "DebugLevel", 1 );
242  if (MyVerbose) ParamList.set( "OutputLevel", 1 );
243  OUR_CHK_ERR( Abase->SetParameters( ParamList ) );
244 
245  if ( TrustMe ) {
246  Problem.SetLHS( &FixedLHS );
247  Problem.SetRHS( &FixedRHS );
248  assert( OptStorage) ;
249  }
250 
251  // bug #2184
252  // Structurally singular matrices are not detected by
253  // Amesos_Klu::SymbolicFactorization() but they are detected by
254  // Amesos_Klu::NumericFactorization()
255  if ( ExpectedError == StructurallySingularMatrixError )
256  ExpectedError = NumericallySingularMatrixError ;
257  if ( ExpectedError == StructurallySingularMatrixError ) {
258  Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
259  int oldtracebackmode = ECM->GetTracebackMode( ) ;
260  ECM->SetTracebackMode(0); // We expect an error, but we don't want it to print out
261 
262  const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
263  ECM->SetTracebackMode(oldtracebackmode);
264  // bug #2245 - Amesos fails to return error consistently across all
265  // processes. When this bug is fixed, remove "iam == 0 &&" from the next line
266  if ( iam == 0 && SymbolicFactorizationReturn != ExpectedError ) {
267  std::cout << " SymbolicFactorization returned " << SymbolicFactorizationReturn
268  << " should be " << ExpectedError << std::endl ;
269  OUR_CHK_ERR( 1 ) ;
270  } else {
271  return 0; // Returned the correct error code for this matrix
272  }
273  } else {
274  const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
275  OUR_CHK_ERR( SymbolicFactorizationReturn ) ;
276  }
277  if ( ExpectedError == NumericallySingularMatrixError ) {
278  Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
279  int oldtracebackmode = ECM->GetTracebackMode( ) ;
280  ECM->SetTracebackMode(0); // We expect an error, but we don't want it to print out
281 
282  const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
283  ECM->SetTracebackMode(oldtracebackmode);
284  // bug #2245 - Amesos fails to return error consistently across all
285  // processes. When this bug is fixed, remove "iam == 0 &&" from the next line
286  if ( iam == 0 && NumericFactorizationReturn != ExpectedError ) {
287  std::cout << " NumericFactorization returned " << NumericFactorizationReturn
288  << " should be " << ExpectedError << std::endl ;
289  OUR_CHK_ERR( 1 ) ;
290  } else {
291  return 0; // Returned the correct error code for this matrix
292  }
293  } else {
294  const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
295  OUR_CHK_ERR( NumericFactorizationReturn ) ;
296  }
297  int ind[1];
298  double val[1];
299  ind[0] = 0;
300  xexact.Random();
301  xexact.PutScalar(1.0);
302 
303  //
304  // Compute cAx = A' xexact
305  //
306  double Value = 1.0 ;
307  if ( Levels == 3 )
308  {
309  val[0] = Value ;
310  if ( MyMatWithDiag->MyGRID( 0 ) ) {
311  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
312  }
313  MyMatWithDiag->Multiply( transpose, xexact, cAx ) ;
314 
315  val[0] = - Value ;
316  if ( MyMatWithDiag->MyGRID( 0 ) )
317  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
318  }
319  else
320  {
321  cAx = xexact ;
322  }
323 
324  //
325  // Compute cAAx = A' cAx
326  //
327  if ( Levels >= 2 )
328  {
329  val[0] = Value ;
330  if ( MyMatWithDiag->MyGRID( 0 ) )
331  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
332  MyMatWithDiag->Multiply( transpose, cAx, cAAx ) ; // x2 = A' x1
333 
334  val[0] = - Value ;
335  if ( MyMatWithDiag->MyGRID( 0 ) )
336  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
337  }
338  else
339  {
340  cAAx = cAx ;
341  }
342 
343  if ( MyVerbose ) std::cout << " Compute b = A x2 = A A' A'' xexact " << std::endl ;
344 
345  MyMatWithDiag->Multiply( transpose, cAAx, b ) ; // b = A x2 = A A' A'' xexact
346 
347 
348  // Phase 2: Solve A' A' A x = b
349  //
350  //
351  // Solve A sAAx = b
352  //
353  if ( TrustMe ) {
354  FixedRHS = b;
355  } else {
356  Problem.SetLHS( &sAAx );
357  Problem.SetRHS( &b );
358  }
359 
360 
361  OUR_CHK_ERR( Abase->SymbolicFactorization( ) );
362  OUR_CHK_ERR( Abase->SymbolicFactorization( ) ); // This should be irrelevant, but should nonetheless be legal
363  OUR_CHK_ERR( Abase->NumericFactorization( ) );
364  OUR_CHK_ERR( Abase->Solve( ) );
365  if ( TrustMe ) sAAx = FixedLHS ;
366 
367  if ( Levels >= 2 )
368  {
369  OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
370  if ( TrustMe ) {
371  FixedRHS = sAAx ;
372  } else {
373  Problem.SetLHS( &sAx );
374  Problem.SetRHS( &sAAx );
375  }
376  val[0] = Value ;
377  if ( MyMat->MyGRID( 0 ) )
378  MyMat->SumIntoMyValues( 0, 1, val, ind ) ;
379  OUR_CHK_ERR( Abase->NumericFactorization( ) );
380  OUR_CHK_ERR( Abase->Solve( ) );
381  if ( TrustMe ) sAx = FixedLHS ;
382 
383  }
384  else
385  {
386  sAx = sAAx ;
387  }
388 
389  if ( Levels >= 3 )
390  {
391  if ( TrustMe ) {
392  FixedRHS = sAx ;
393  } else {
394  Problem.SetLHS( &x );
395  Problem.SetRHS( &sAx );
396  }
397  OUR_CHK_ERR( Abase->Solve( ) );
398  if ( TrustMe ) x = FixedLHS ;
399  }
400  else
401  {
402  x = sAx ;
403  }
404 
405  if ( Levels >= 2 )
406  {
407  val[0] = -Value ;
408  if ( MyMat->MyGRID( 0 ) ) {
409  if ( MyMat->SumIntoMyValues( 0, 1, val, ind ) ) {
410  std::cout << " TestOptions requires a non-zero entry in A(1,1) " << std::endl ;
411  }
412  }
413  }
414 
415  //
416  // Phase 3: Check the residual: bcheck = A' A' A x
417  //
418 
419 
420  if ( Levels >= 3 )
421  {
422  val[0] = Value ;
423  if ( MyMatWithDiag->MyGRID( 0 ) )
424  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
425  MyMatWithDiag->Multiply( transpose, x, kAx ) ;
426  val[0] = -Value ;
427  if ( MyMatWithDiag->MyGRID( 0 ) )
428  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
429  }
430  else
431  {
432  kAx = x ;
433  }
434 
435  if ( Levels >= 2 )
436  {
437  val[0] = Value ;
438  if ( MyMatWithDiag->MyGRID( 0 ) )
439  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
440  MyMatWithDiag->Multiply( transpose, kAx, kAAx ) ;
441  val[0] = -Value ;
442  if ( MyMatWithDiag->MyGRID( 0 ) )
443  MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
444  }
445  else
446  {
447  kAAx = kAx ;
448  }
449 
450  MyMatWithDiag->Multiply( transpose, kAAx, bcheck ) ; // temp = A" x2
451 
452  double norm_diff ;
453  double norm_one ;
454 
455  DomainDiff.Update( 1.0, sAAx, -1.0, cAAx, 0.0 ) ;
456  DomainDiff.Norm2( &norm_diff ) ;
457  sAAx.Norm2( &norm_one ) ;
458  if (MyVerbose) std::cout << __FILE__ << "::" << __LINE__
459  << " norm( sAAx - cAAx ) / norm(sAAx ) = "
460  << norm_diff /norm_one << std::endl ;
461 
462 
463 
464 
465  DomainDiff.Update( 1.0, sAx, -1.0, cAx, 0.0 ) ;
466  DomainDiff.Norm2( &norm_diff ) ;
467  sAx.Norm2( &norm_one ) ;
468  if (MyVerbose) std::cout
469  << __FILE__ << "::" << __LINE__
470  << " norm( sAx - cAx ) / norm(sAx ) = "
471  << norm_diff /norm_one << std::endl ;
472 
473 
474  DomainDiff.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
475  DomainDiff.Norm2( &norm_diff ) ;
476  x.Norm2( &norm_one ) ;
477  if (MyVerbose) std::cout
478  << __FILE__ << "::" << __LINE__
479  << " norm( x - xexact ) / norm(x) = "
480  << norm_diff /norm_one << std::endl ;
481 
482  relerror = norm_diff / norm_one ;
483 
484  DomainDiff.Update( 1.0, sAx, -1.0, kAx, 0.0 ) ;
485  DomainDiff.Norm2( &norm_diff ) ;
486  sAx.Norm2( &norm_one ) ;
487  if (MyVerbose) std::cout
488  << __FILE__ << "::" << __LINE__
489  << " norm( sAx - kAx ) / norm(sAx ) = "
490  << norm_diff /norm_one << std::endl ;
491 
492 
493  DomainDiff.Update( 1.0, sAAx, -1.0, kAAx, 0.0 ) ;
494  DomainDiff.Norm2( &norm_diff ) ;
495  sAAx.Norm2( &norm_one ) ;
496  if (MyVerbose) std::cout
497  << __FILE__ << "::" << __LINE__
498  << " norm( sAAx - kAAx ) / norm(sAAx ) = "
499  << norm_diff /norm_one << std::endl ;
500 
501 
502  RangeDiff.Update( 1.0, bcheck, -1.0, b, 0.0 ) ;
503  RangeDiff.Norm2( &norm_diff ) ;
504  bcheck.Norm2( &norm_one ) ;
505  if (MyVerbose) std::cout
506  << __FILE__ << "::" << __LINE__
507  << " norm( bcheck - b ) / norm(bcheck ) = "
508  << norm_diff /norm_one << std::endl ;
509 
510  relresidual = norm_diff / norm_one ;
511 
512  if (iam == 0 ) {
513  if ( relresidual * Rcond < 1e-16 ) {
514  if (MyVerbose) std::cout << " Test 1 Passed " << std::endl ;
515  } else {
516  std::cout << __FILE__ << "::" << __LINE__ <<
517  " relresidual = " << relresidual <<
518  " TEST FAILED " <<
519  " ParamList = " << ParamList << std::endl ;
520  errors += 1 ;
521  }
522  }
523  }
524 
525  return errors;
526 
527 }
528 
529 
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
bool MyGRID(int GRID_in) const
void SetLHS(Epetra_MultiVector *X)
void SetOperator(Epetra_RowMatrix *A)
bool StorageOptimized() const
virtual int Solve()=0
Solves A X = B (or AT x = B)
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
T & get(ParameterList &l, const std::string &name)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
int PartialFactorization(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond)
static bool verbose
Definition: Amesos.cpp:67
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
#define EPETRA_MIN(x, y)
const Epetra_Map & OperatorDomainMap() const
virtual int MyPID() const =0
const Epetra_Map & RowMap() const
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
void SetRHS(Epetra_MultiVector *B)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
const Epetra_Map & OperatorRangeMap() const
const int NumericallySingularMatrixError
#define OUR_CHK_ERR(a)
const int StructurallySingularMatrixError
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)