Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
PartialFactorization.cpp
Go to the documentation of this file.
1//
2// OUR_CHK_ERR always returns 1 on error.
3//
4#define OUR_CHK_ERR(a) { { int epetra_err = a; \
5 if (epetra_err != 0) { std::cerr << "Amesos ERROR " << epetra_err << ", " \
6 << __FILE__ << ", line " << __LINE__ << std::endl; \
7relresidual=1e15; return(1);} }\
8 }
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"
19//
20// PartialFactorization checks to make sure that we clean up after our mess
21// before everything is done.
22//
23
24const int MaxNumSteps = 7 ;
25
26int PartialFactorizationOneStep( const char* AmesosClass,
27 const Epetra_Comm &Comm,
28 bool transpose,
29 bool verbose,
30 Teuchos::ParameterList ParamList,
31 Epetra_CrsMatrix *& Amat,
32 double Rcond,
33 int Steps )
34{
35
36 assert( Steps >= 0 && Steps < MaxNumSteps ) ;
37
38 int iam = Comm.MyPID() ;
39 int errors = 0 ;
40
41 const Epetra_Map *RangeMap =
42 transpose?&Amat->OperatorDomainMap():&Amat->OperatorRangeMap() ;
43 const Epetra_Map *DomainMap =
44 transpose?&Amat->OperatorRangeMap():&Amat->OperatorDomainMap() ;
45
46 Epetra_Vector xexact(*DomainMap);
47 Epetra_Vector x(*DomainMap);
48
49 Epetra_Vector b(*RangeMap);
50 Epetra_Vector bcheck(*RangeMap);
51
52 Epetra_Vector difference(*DomainMap);
53
54 Epetra_LinearProblem Problem;
55 Amesos_BaseSolver* Abase ;
56 Amesos Afactory;
57
58 Abase = Afactory.Create( AmesosClass, Problem ) ;
59
60 std::string AC = AmesosClass ;
61 if ( AC == "Amesos_Mumps" ) {
62 ParamList.set( "NoDestroy", true );
63 Abase->SetParameters( ParamList ) ;
64 }
65
66 double relresidual = 0 ;
67
68 if ( Steps > 0 ) {
69 //
70 // Phase 1: Compute b = A' A' A xexact
71 //
72 Problem.SetOperator( Amat );
73
74 //
75 // We only set transpose if we have to - this allows valgrind to check
76 // that transpose is set to a default value before it is used.
77 //
78 if ( transpose ) OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
79 // if (verbose) ParamList.set( "DebugLevel", 1 );
80 // if (verbose) ParamList.set( "OutputLevel", 1 );
81 if ( Steps > 1 ) {
82 OUR_CHK_ERR( Abase->SetParameters( ParamList ) );
83 if ( Steps > 2 ) {
84
85 xexact.Random();
86 xexact.PutScalar(1.0);
87
88 //
89 // Compute cAx = A' xexact
90 //
91 Amat->Multiply( transpose, xexact, b ) ; // b = A x2 = A A' A'' xexact
92
93#if 0
94 std::cout << __FILE__ << "::" << __LINE__ << "b = " << std::endl ;
95 b.Print( std::cout ) ;
96 std::cout << __FILE__ << "::" << __LINE__ << "xexact = " << std::endl ;
97 xexact.Print( std::cout ) ;
98 std::cout << __FILE__ << "::" << __LINE__ << "x = " << std::endl ;
99 x.Print( std::cout ) ;
100#endif
101 //
102 // Phase 2: Solve A' A' A x = b
103 //
104 //
105 // Solve A sAAx = b
106 //
107 Problem.SetLHS( &x );
108 Problem.SetRHS( &b );
110 if ( Steps > 2 ) {
112 if ( Steps > 3 ) {
113 OUR_CHK_ERR( Abase->NumericFactorization( ) );
114 if ( Steps > 4 ) {
115 OUR_CHK_ERR( Abase->NumericFactorization( ) );
116 if ( Steps > 5 ) {
117 OUR_CHK_ERR( Abase->Solve( ) );
118 if ( Steps > 6 ) {
119 OUR_CHK_ERR( Abase->Solve( ) );
120
121
122 Amat->Multiply( transpose, x, bcheck ) ; // temp = A" x2
123
124 double norm_diff ;
125 double norm_one ;
126
127 difference.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
128 difference.Norm2( &norm_diff ) ;
129 x.Norm2( &norm_one ) ;
130
131 relresidual = norm_diff / norm_one ;
132
133 if (iam == 0 ) {
134 if ( relresidual * Rcond > 1e-16 ) {
135 if (verbose) std::cout << __FILE__ << "::"<< __LINE__
136 << " norm( x - xexact ) / norm(x) = "
137 << norm_diff /norm_one << std::endl ;
138 errors += 1 ;
139 }
140 }
141 }
142 }
143 }
144 }
145 }
146 }
147 }
148}
149 delete Abase;
150
151 return errors;
152
153}
154
155
156int PartialFactorization( const char* AmesosClass,
157 const Epetra_Comm &Comm,
158 bool transpose,
159 bool verbose,
160 Teuchos::ParameterList ParamList,
161 Epetra_CrsMatrix *& Amat,
162 double Rcond ) {
163
164 for( int i =0 ; i < MaxNumSteps ; i ++ ) {
165 std::string AC = AmesosClass ;
166 int epetra_err = PartialFactorizationOneStep(AmesosClass, Comm, transpose, verbose, ParamList, Amat, Rcond, i);
167 if (epetra_err != 0) {
168 std::cerr << "Amesos ERROR " << epetra_err << ", " << __FILE__ << ", line " << __LINE__ << std::endl;
169 return(1);
170 }
171 }
172 return(0);
173}
static bool verbose
Definition: Amesos.cpp:67
const int MaxNumSteps
int PartialFactorizationOneStep(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond, int Steps)
int PartialFactorization(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond)
#define OUR_CHK_ERR(a)
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
virtual int Solve()=0
Solves A X = B (or AT x = B)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69