33 #include "Epetra_MpiComm.h" 35 #include "Epetra_SerialComm.h" 39 #include "Epetra_Map.h" 40 #include "Epetra_CrsMatrix.h" 41 #include "Epetra_Vector.h" 42 #include "Epetra_LinearProblem.h" 68 int main(
int argc,
char *argv[])
72 MPI_Init(&argc, &argv);
73 Epetra_MpiComm Comm(MPI_COMM_WORLD);
75 Epetra_SerialComm Comm;
78 int NumGlobalElements = 100;
87 Epetra_Map Map(NumGlobalElements, 0, Comm);
90 Epetra_CrsMatrix
A(
Copy, Map, 0);
93 int NumMyElements = Map.NumMyElements();
106 int* MyGlobalElements = Map.MyGlobalElements();
110 for (
int i = 0; i < NumMyElements; i++)
112 if (MyGlobalElements[i] == 0)
118 else if (MyGlobalElements[i] == NumGlobalElements-1)
120 Indices[0] = NumGlobalElements-1;
121 Indices[1] = NumGlobalElements-2;
126 Indices[0] = MyGlobalElements[i]-1;
127 Indices[1] = MyGlobalElements[i];
128 Indices[2] = MyGlobalElements[i]+1;
133 NumEntries, Values, Indices));
155 Epetra_LinearProblem Problem;
157 Problem.SetOperator(&
A);
169 std::cerr <<
"Selected solver is not available" << std::endl;
180 Solver->SymbolicFactorization();
185 for (
int i = 0; i < NumMyElements; i++)
187 if (MyGlobalElements[i] == 0)
195 else if (MyGlobalElements[i] == NumGlobalElements-1)
197 Indices[0] = NumGlobalElements - 1;
198 Indices[1] = NumGlobalElements - 2;
205 Indices[0] = MyGlobalElements[i] - 1;
206 Indices[1] = MyGlobalElements[i];
207 Indices[2] = MyGlobalElements[i] + 1;
215 NumEntries, Values, Indices));
219 Solver->NumericFactorization();
222 Epetra_Vector b(Map);
224 Epetra_Vector x(Map);
236 Solver->GetTiming( TimingsList );
239 double sfact_time, nfact_time, solve_time;
240 double mtx_conv_time, mtx_redist_time, vec_redist_time;
244 sfact_time = TimingsList.
get(
"Total symbolic factorization time", 0.0 );
248 nfact_time = Teuchos::getParameter<double>( TimingsList,
"Total numeric factorization time" );
252 solve_time = Teuchos::getParameter<double>( TimingsList,
"Total solve time" );
256 mtx_conv_time = Teuchos::getParameter<double>( TimingsList,
"Total solve time" );
259 mtx_redist_time = TimingsList.get(
"Total matrix redistribution time", 0.0 );
262 vec_redist_time = TimingsList.get(
"Total vector redistribution time", 0.0 );
264 std::cout <<
" sfact_time " << sfact_time
265 <<
" nfact_time " << nfact_time
266 <<
" solve_time " << solve_time
267 <<
" mtx_conv_time " << mtx_conv_time
268 <<
" mtx_redist_time " << mtx_redist_time
269 <<
" vec_redist_time " << vec_redist_time
282 Epetra_Vector Ax(b.Map());
283 A.Multiply(
false, x, Ax);
284 Ax.Update(1.0, b, -1.0);
288 std::cout <<
"After AMESOS solution, ||b-Ax||_2 = " << residual << std::endl;
295 return(EXIT_FAILURE);
300 return(EXIT_SUCCESS);
T & get(ParameterList &l, const std::string &name)
int main(int argc, char *argv[])
#define AMESOS_CHK_ERR(a)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...