Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Test_MultipleSolves/cxx_main.cpp
Go to the documentation of this file.
1#include "Amesos_ConfigDefs.h"
2
3#ifdef HAVE_MPI
4#include "mpi.h"
5#include "Epetra_MpiComm.h"
6#else
7#include "Epetra_SerialComm.h"
8#endif
9#include "Epetra_Map.h"
10#include "Epetra_MultiVector.h"
11#include "Epetra_Time.h"
12#include "Epetra_CrsMatrix.h"
13#include "Amesos.h"
14#include "Teuchos_ParameterList.hpp"
15#include "Teuchos_Array.hpp"
16#include "Galeri_Maps.h"
17#include "Galeri_CrsMatrices.h"
18
19using namespace Galeri;
20
21bool TestAmesos(char ProblemType[],
22 Epetra_RowMatrix& A,
23 int NumVectors)
24{
25
26 const Epetra_BlockMap& Map = A.OperatorDomainMap();
27
28 Epetra_MultiVector x2(Map,NumVectors);
29 Epetra_MultiVector x1(Map,NumVectors);
30 Epetra_MultiVector x(Map,NumVectors);
31 Epetra_MultiVector b(Map,NumVectors);
32 Epetra_MultiVector residual(Map,NumVectors);
33 Epetra_MultiVector temp(Map,NumVectors);
34
35 Teuchos::ParameterList ParamList ;
36 Epetra_LinearProblem Problem;
37 Amesos_BaseSolver* Abase ;
38 Amesos Afactory;
39
40 // Note that Abase is created with an empty Problem, none of A, x or b
41 // have been specified at this point.
42 Abase = Afactory.Create(ProblemType,Problem );
43 assert (Abase != 0);
44
45 // Factor A
46 Problem.SetOperator(&A);
49
50 b.Random();
51
52 // Solve Ax = b
53 //
54 Problem.SetLHS(&x);
55 Problem.SetRHS(&b);
56 EPETRA_CHK_ERR(Abase->Solve());
57
58 // Solve Ax1 = x
59 //
60 Problem.SetLHS(&x1);
61 Problem.SetRHS(&x);
62 EPETRA_CHK_ERR(Abase->Solve());
63
64 // Solve Ax2 = x1
65 //
66 Problem.SetLHS(&x2);
67 Problem.SetRHS(&x1);
68 EPETRA_CHK_ERR(Abase->Solve());
69
70 // Compute the residual: A^3 x2 - b
71 //
72 A.Multiply(false,x2, temp) ; // temp = A x2
73 A.Multiply(false,temp, x2) ; // x2 = A^2 x2
74 A.Multiply(false,x2, temp) ; // temp = A^3 x2
75 residual.Update( 1.0, temp, -1.0, b, 0.0 ) ;
76
77 Teuchos::Array<double> norm_residual(NumVectors);
78 residual.Norm2( norm_residual.getRawPtr() ) ;
79
80 if (A.Comm().MyPID() == 0) {
81 std::cout << "norm2(A^3 x-b) = " << norm_residual << std::endl ;
82 }
83
84 delete Abase;
85
86 for (Teuchos::Array<double>::const_iterator it = norm_residual.begin(); it != norm_residual.end(); ++it) {
87 if (!(*it < (1e-5)))
88 return(false);
89 }
90 return(true);
91}
92
93int main(int argc, char *argv[]) {
94
95#ifdef HAVE_MPI
96 MPI_Init(&argc, &argv);
97 Epetra_MpiComm Comm(MPI_COMM_WORLD);
98#else
99 Epetra_SerialComm Comm;
100#endif
101
102 int NumVectors = 2;
103
104 Teuchos::ParameterList GaleriList;
105 GaleriList.set("nx", 8);
106 GaleriList.set("ny", 8 * Comm.NumProc());
107 GaleriList.set("mx", 1);
108 GaleriList.set("my", Comm.NumProc());
109
110 Epetra_Map* Map = CreateMap("Cartesian2D", Comm, GaleriList);
111 Epetra_CrsMatrix* A = CreateCrsMatrix("Laplace2D", Map, GaleriList);
112
113 Amesos Factory;
114
115 std::vector<std::string> SolverType;
116 // SolverType.push_back("Amesos_Lapack");
117 SolverType.push_back("Amesos_Klu");
118 SolverType.push_back("Amesos_Umfpack");
119 SolverType.push_back("Amesos_Pardiso");
120 SolverType.push_back("Amesos_Taucs");
121 SolverType.push_back("Amesos_Superlu");
122 SolverType.push_back("Amesos_Superludist");
123 SolverType.push_back("Amesos_Mumps");
124 SolverType.push_back("Amesos_Dscpack");
125
126 bool TestPassed = true;
127
128 for (unsigned int i = 0 ; i < SolverType.size() ; ++i) {
129 std::string Solver = SolverType[i];
130 std::cout << Solver << " next " << std::endl;
131 if (Factory.Query((char*)Solver.c_str())) {
132 if (Comm.MyPID() == 0)
133 std::cout << "Testing " << Solver << std::endl;
134 if(TestAmesos((char*)Solver.c_str(), *A, NumVectors) == false) {
135 std::cout << Solver << " Failed " << std::endl;
136 TestPassed = false;
137 } else {
138 std::cout << Solver << " Passed " << std::endl;
139 }
140 } else
141 if (Comm.MyPID() == 0) {
142 std::cerr << std::endl;
143 std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
144 std::cerr << std::endl;
145 }
146 }
147
148 delete A;
149 delete Map;
150
151#ifdef HAVE_MPI
152 MPI_Finalize();
153#endif
154
155 if (TestPassed) {
156 if (Comm.MyPID() == 0)
157 std::cout << "TESTS PASSED!" << std::endl;
158 return( EXIT_SUCCESS );
159 }
160 else {
161 if (Comm.MyPID() == 0)
162 std::cout << "TESTS FAILED!" << std::endl;
163 return( EXIT_FAILURE );
164 }
165
166}
#define EPETRA_CHK_ERR(xxx)
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Definition: TestOptions.cpp:81
int main(int argc, char *argv[])
bool TestAmesos(char ProblemType[], Epetra_RowMatrix &A, int NumVectors)
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.
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
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193