Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Test_Epetra_CrsMatrix/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// I have no idea why we have to include Epetra_Map.h here (and not
9// elsewhere, such as Test_EpteraRowMatrix/cxx_main.cpp)
10// but this was added as a workaround to bug #1869
11#include "Epetra_Map.h"
12#endif
13#include "Epetra_Vector.h"
14#include "Epetra_CrsMatrix.h"
15#include "Amesos.h"
16#include "Teuchos_ParameterList.hpp"
17
18// ======================================================================
19// this function tests two things:
20// - the return code from Amesos;
21// - the true residual.
22// The return value is `true' if both tests are passed, `false' otherwise.
23//
24// \author Marzio Sala, SNL 9214.
25//
26// \date Last updated on 21-Apr-04.
27// ======================================================================
28
29bool TestAmesos(char ProblemType[], Teuchos::ParameterList& AmesosList,
30 bool UseTranspose, Epetra_RowMatrix* A, Epetra_MultiVector* lhs,
31 Epetra_MultiVector* rhs)
32{
33 Epetra_LinearProblem Problem;
34 Amesos A_Factory;
35
36 Amesos_BaseSolver * Solver = A_Factory.Create(ProblemType, Problem);
37 assert (Solver != 0);
38
39 // Both sentences should work
40 Solver->SetUseTranspose(UseTranspose);
41
42 Solver->SetParameters(AmesosList);
43
44 // create a rhs corresponding to lhs or 1's
45 lhs->PutScalar(1.0);
46 A->Multiply(UseTranspose,*lhs,*rhs);
47 lhs->PutScalar(0.0);
48
49 Problem.SetOperator(A);
50
51 if (Solver->SymbolicFactorization()) return(false);
52 if (Solver->NumericFactorization()) return(false);
53
54 // set sol and rhs here
55 Problem.SetLHS(lhs);
56 Problem.SetRHS(rhs);
57
58 if (Solver->Solve()) return(false);
59
60 // compute difference between exact solution and Amesos
61 double d = 0.0, d_tot = 0.0;
62
63 for (int i=0 ; i<lhs->Map().NumMyElements() ; ++i)
64 for (int j = 0 ; j < lhs->NumVectors() ; ++j)
65 d += ((*lhs)[j][i] - 1.0) * ((*lhs)[j][i] - 1.0);
66
67 A->Comm().SumAll(&d,&d_tot,1);
68
69 // compute ||Ax - b||
70 std::vector<double> Norm(rhs->NumVectors());
71
72 Epetra_MultiVector Ax(*rhs);
73 A->Multiply(UseTranspose, *lhs, Ax);
74 Ax.Update(1.0, *rhs, -1.0);
75 Ax.Norm2(&Norm[0]);
76
77 std::string msg = ProblemType;
78
79 if (!A->Comm().MyPID())
80 {
81 std::cout << std::endl;
82 std::cout << msg << " : Using " << A->Comm().NumProc() << " processes, UseTranspose = " << UseTranspose << std::endl;
83 for (int j = 0 ; j < rhs->NumVectors() ; ++j)
84 std::cout << msg << " : eq " << j
85 << ", ||A x - b||_2 = " << Norm[j] << std::endl;
86 std::cout << msg << " : ||x_exact - x||_2 = " << sqrt(d_tot) << std::endl;
87 }
88
89 if (Norm[0] > 1e-9)
90 {
91 std::cerr << std::endl << msg << " WARNING : TEST FAILED!" << std::endl;
92 return(false);
93 }
94
95 delete Solver;
96
97 return(true);
98}
99
100// =========== //
101// main driver //
102// =========== //
103//
104int main(int argc, char *argv[])
105{
106#ifdef HAVE_MPI
107 MPI_Init(&argc, &argv);
108 Epetra_MpiComm Comm(MPI_COMM_WORLD);
109#else
110 Epetra_SerialComm Comm;
111#endif
112
113 int NumGlobalRows = 100 * Comm.NumProc();
114
115 Epetra_Map Map(NumGlobalRows, 0, Comm);
116
117 Epetra_CrsMatrix Matrix(Copy, Map, 0);
118
119 int NumMyRows = Map.NumMyElements();
120 int* MyGlobalElements = Map.MyGlobalElements();
121
122 int Indices[3];
123 double Values[3];
124 int NumEntries;
125
126 for (int LRID = 0 ; LRID < NumMyRows ; ++LRID)
127 {
128 int GRID = MyGlobalElements[LRID];
129
130 Indices[0] = GRID;
131 Values[0] = 2.0;
132 NumEntries = 1;
133
134 if (GRID != 0)
135 {
136 Indices[NumEntries] = GRID - 1;
137 Values[NumEntries] = -1.0;
138 ++NumEntries;
139 }
140 if (GRID != NumGlobalRows - 1)
141 {
142 Indices[NumEntries] = GRID + 1;
143 Values[NumEntries] = -1.0;
144 ++NumEntries;
145 }
146
147 Matrix.InsertGlobalValues(GRID, NumEntries, Values, Indices);
148 }
149
150 Matrix.FillComplete();
151
152 Epetra_MultiVector LHS(Map, 3);
153 Epetra_MultiVector RHS(Map, 3);
154
155 Amesos Factory;
156
157 std::vector<std::string> SolverType;
158 SolverType.push_back("Amesos_Lapack");
159 SolverType.push_back("Amesos_Klu");
160 SolverType.push_back("Amesos_Umfpack");
161 // SolverType.push_back("Amesos_Pardiso"); bug #1994
162 SolverType.push_back("Amesos_Taucs");
163 SolverType.push_back("Amesos_Superlu");
164 SolverType.push_back("Amesos_Superludist");
165 SolverType.push_back("Amesos_Mumps");
166 SolverType.push_back("Amesos_Dscpack");
167 SolverType.push_back("Amesos_Scalapack");
168#ifndef NDEBUG
169 bool res;
170#endif
171 // If a given test fails, than the code stops, due to the assert()
172 // statement.
173 for (unsigned int i = 0 ; i < SolverType.size() ; ++i)
174 {
175 std::string Solver = SolverType[i];
176
177 if (Factory.Query((char*)Solver.c_str()))
178 {
179 if (1) {
180 // solve with matrix
181 Teuchos::ParameterList AmesosList;
182#ifndef NDEBUG
183 res =
184#endif
185 TestAmesos((char*)Solver.c_str(), AmesosList, false,
186 &Matrix, &LHS, &RHS);
187 assert (res == true);
188 }
189 if (1) {
190 // solve transpose with matrix
191 if (Solver != "Amesos_Superludist") {// still not implementes
192 Teuchos::ParameterList AmesosList;
193#ifndef NDEBUG
194 res =
195#endif
196 TestAmesos((char*)Solver.c_str(), AmesosList, true,
197 &Matrix, &LHS, &RHS);
198 assert (res == true);
199 }
200 }
201 }
202 else
203 if (!Comm.MyPID())
204 {
205 std::cerr << std::endl;
206 std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
207 std::cerr << std::endl;
208 }
209 }
210
211#ifdef HAVE_MPI
212 MPI_Finalize();
213#endif
214 return(0);
215}
int main(int argc, char *argv[])
bool TestAmesos(char ProblemType[], Teuchos::ParameterList &AmesosList, bool UseTranspose, Epetra_RowMatrix *A, Epetra_MultiVector *lhs, Epetra_MultiVector *rhs)
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
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193