Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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 
29 bool 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 //
104 int 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 
169  bool res;
170 
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  res = TestAmesos((char*)Solver.c_str(), AmesosList, false,
183  &Matrix, &LHS, &RHS);
184  assert (res == true);
185  }
186  if (1) {
187  // solve transpose with matrix
188  if (Solver != "Amesos_Superludist") {// still not implementes
189  Teuchos::ParameterList AmesosList;
190  res = TestAmesos((char*)Solver.c_str(), AmesosList, true,
191  &Matrix, &LHS, &RHS);
192  assert (res == true);
193  }
194  }
195  }
196  else
197  if (!Comm.MyPID())
198  {
199  std::cerr << std::endl;
200  std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
201  std::cerr << std::endl;
202  }
203  }
204 
205 #ifdef HAVE_MPI
206  MPI_Finalize();
207 #endif
208  return(0);
209 }
void SetLHS(Epetra_MultiVector *X)
void SetOperator(Epetra_RowMatrix *A)
bool TestAmesos(char ProblemType[], Teuchos::ParameterList &AmesosList, bool UseTranspose, Epetra_RowMatrix *A, Epetra_MultiVector *lhs, Epetra_MultiVector *rhs)
virtual int Solve()=0
Solves A X = B (or AT x = B)
int MyGlobalElements(int *MyGlobalElementList) const
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
int MyPID() const
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual const Epetra_Comm & Comm() const =0
virtual const Epetra_BlockMap & Map() const =0
int main(int argc, char *argv[])
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
int NumProc() const
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 Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
void SetRHS(Epetra_MultiVector *B)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
virtual int NumProc() const =0
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193