Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Test_MUMPS/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_Vector.h"
11#include "Epetra_Time.h"
12#include "Epetra_CrsMatrix.h"
13#include "Epetra_Util.h"
14#include "Amesos_Mumps.h"
16#include "Teuchos_ParameterList.hpp"
17// using namespace Trilinos_Util; commented out to resolve bug #1886
18
19//=============================================================================
20bool CheckError(const Epetra_RowMatrix& A,
21 const Epetra_MultiVector& x,
22 const Epetra_MultiVector& b,
23 const Epetra_MultiVector& x_exact)
24{
25 std::vector<double> Norm;
26 int NumVectors = x.NumVectors();
27 Norm.resize(NumVectors);
28 Epetra_MultiVector Ax(x);
29 A.Multiply(false,x,Ax);
30 Ax.Update(1.0,b,-1.0);
31 Ax.Norm2(&Norm[0]);
32 bool TestPassed = false;
33 double TotalNorm = 0.0;
34 for (int i = 0 ; i < NumVectors ; ++i) {
35 TotalNorm += Norm[i];
36 }
37 if (A.Comm().MyPID() == 0)
38 std::cout << "||Ax - b|| = " << TotalNorm << std::endl;
39 if (TotalNorm < 1e-5 )
40 TestPassed = true;
41 else
42 TestPassed = false;
43
44 Ax.Update (1.0,x,-1.0,x_exact,0.0);
45 Ax.Norm2(&Norm[0]);
46 for (int i = 0 ; i < NumVectors ; ++i) {
47 TotalNorm += Norm[i];
48 }
49 if (A.Comm().MyPID() == 0)
50 std::cout << "||Ax - b|| = " << TotalNorm << std::endl;
51#ifdef HAVE_AMESOS_SMUMPS
52 if (TotalNorm < 1e-2 )
53#else
54 if (TotalNorm < 1e-5 )
55#endif
56 TestPassed = true;
57 else
58 TestPassed = false;
59
60 return(TestPassed);
61}
62
63//=============================================================================
64int main(int argc, char *argv[]) {
65
66#ifdef HAVE_MPI
67 MPI_Init(&argc, &argv);
68 Epetra_MpiComm Comm(MPI_COMM_WORLD);
69#else
70 Epetra_SerialComm Comm;
71#endif
72
73 int NumGlobalElements = 100;
74 int NumVectors = 7;
75
76 // =================== //
77 // create a random map //
78 // =================== //
79
80 int* part = new int[NumGlobalElements];
81
82 if (Comm.MyPID() == 0) {
83 Epetra_Util Util;
84
85 for( int i=0 ; i<NumGlobalElements ; ++i ) {
86 unsigned int r = Util.RandomInt();
87 part[i] = r%(Comm.NumProc());
88 }
89 }
90
91 Comm.Broadcast(part,NumGlobalElements,0);
92
93 // count the elements assigned to this proc
94 int NumMyElements = 0;
95 for (int i = 0 ; i < NumGlobalElements ; ++i) {
96 if (part[i] == Comm.MyPID())
97 NumMyElements++;
98 }
99
100 // get the loc2global list
101 int* MyGlobalElements = new int[NumMyElements];
102 int count = 0;
103 for (int i = 0 ; i < NumGlobalElements ; ++i) {
104 if (part[i] == Comm.MyPID() )
105 MyGlobalElements[count++] = i;
106 }
107
108 Epetra_Map Map(NumGlobalElements,NumMyElements,MyGlobalElements,
109 0,Comm);
110
111 delete [] part;
112
113 // ===================== //
114 // Create a dense matrix //
115 // ===================== //
116
117 Epetra_CrsMatrix Matrix(Copy,Map,NumGlobalElements);
118
119 int* Indices = new int[NumGlobalElements];
120 double* Values = new double[NumGlobalElements];
121
122 for (int i = 0 ; i < NumGlobalElements ; ++i)
123 Indices[i] = i;
124
125 for (int i = 0 ; i < NumMyElements ; ++i) {
126 int iGlobal = MyGlobalElements[i];
127 for (int jGlobal = 0 ; jGlobal < NumGlobalElements ; ++jGlobal) {
128 if (iGlobal >= jGlobal)
129 Values[jGlobal] = 1.0 * (jGlobal + 1);
130 else
131 Values[jGlobal] = 1.0 * (iGlobal + 1);
132 }
133 Matrix.InsertGlobalValues(MyGlobalElements[i],
134 NumGlobalElements, Values, Indices);
135
136 }
137
138 Matrix.FillComplete();
139
140 delete [] Indices;
141 delete [] Values;
142
143 // ======================== //
144 // other data for this test //
145 // ======================== //
146
147 Amesos_TestRowMatrix A(&Matrix);
148 Epetra_MultiVector x(Map,NumVectors);
149 Epetra_MultiVector x_exact(Map,NumVectors);
150 Epetra_MultiVector b(Map,NumVectors);
151 x_exact.Random();
152 A.Multiply(false,x_exact,b);
153
154 // =========== //
155 // AMESOS PART //
156 // =========== //
157
158 Epetra_LinearProblem Problem;
159 Amesos_Mumps* Solver = new Amesos_Mumps(Problem);
160
161 Problem.SetOperator(&A);
162 Problem.SetLHS(&x);
163 Problem.SetRHS(&b);
164
165 Teuchos::ParameterList List;
166 List.set("MaxProcs",2);
167 AMESOS_CHK_ERR(Solver->SetParameters(List));
168
171 AMESOS_CHK_ERR(Solver->Solve());
172
173 bool TestPassed = true;
174
175 TestPassed = TestPassed &&
176 CheckError(A,x,b,x_exact);
177
178 delete Solver;
179
180 AMESOS_CHK_ERR( ! TestPassed ) ;
181
182#ifdef HAVE_MPI
183 MPI_Finalize();
184#endif
185
186 return(EXIT_SUCCESS);
187}
#define AMESOS_CHK_ERR(a)
bool CheckError(const Epetra_RowMatrix &A, const Epetra_MultiVector &x, const Epetra_MultiVector &b, const Epetra_MultiVector &x_exact)
int main(int argc, char *argv[])
Amesos_Mumps: An object-oriented wrapper for the double precision version of MUMPS.
Definition: Amesos_Mumps.h:117
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int Solve()
Solves A X = B (or AT x = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.