1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
/* ******************************************************************** */
/* See the file COPYRIGHT for a complete copyright notice, contact */
/* person and disclaimer. */
/* ******************************************************************** */
#include "ml_common.h"
#if defined(HAVE_ML_MLAPI) && defined(HAVE_ML_GALERI) && defined(HAVE_ML_AZTECOO)
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
// required to build the example matrix
#include "Galeri_Maps.h"
#include "Galeri_CrsMatrices.h"
// required by the linear system solver
#include "AztecOO.h"
#include "MLAPI_Space.h"
#include "MLAPI_Operator.h"
#include "MLAPI_MultiLevelSA.h"
#include "MLAPI_EpetraBaseOperator.h"
#include "MLAPI_Workspace.h"
using namespace Teuchos;
using namespace MLAPI;
using namespace Galeri;
// ============== //
// example driver //
// ============== //
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
// get the epetra matrix from the Gallery
int nx = 8;
int ny = 8 * Comm.NumProc();
ParameterList GaleriList;
GaleriList.set("nx", nx);
GaleriList.set("ny", ny);
GaleriList.set("mx", 1);
GaleriList.set("my", Comm.NumProc());
Epetra_Map* Map = CreateMap("Cartesian2D", Comm, GaleriList);
Epetra_CrsMatrix* A = CreateCrsMatrix("Laplace2D", Map, GaleriList);
Epetra_Vector LHS(*Map); LHS.Random();
Epetra_Vector RHS(*Map); RHS.PutScalar(0.0);
Epetra_LinearProblem Problem(A, &LHS, &RHS);
AztecOO solver(Problem);
Init();
try {
Space S(-1, A->NumMyRows(), A->RowMatrixRowMap().MyGlobalElements());
Operator A_MLAPI(S, S, A, false);
Teuchos::ParameterList MLList;
MLList.set("max levels",3);
MLList.set("increasing or decreasing","increasing");
MLList.set("aggregation: type", "Uncoupled");
MLList.set("aggregation: damping factor", 0.0);
MLList.set("smoother: type","symmetric Gauss-Seidel");
MLList.set("smoother: sweeps",1);
MLList.set("smoother: damping factor",1.0);
MLList.set("coarse: max size",3);
MLList.set("smoother: pre or post", "both");
MLList.set("coarse: type","Amesos-KLU");
MultiLevelSA Prec_MLAPI(A_MLAPI, MLList);
Epetra_Operator* Prec = new EpetraBaseOperator(A->RowMatrixRowMap(), Prec_MLAPI);
solver.SetPrecOperator(Prec);
solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
solver.SetAztecOption(AZ_output, 32);
solver.Iterate(500, 1e-12);
// destroy the preconditioner
delete Prec;
}
catch (const int e) {
cerr << "Caught exception, code = " << e << endl;
}
catch (...) {
cerr << "Caught exception..." << endl;
}
Finalize();
// compute the real residual
double residual;
LHS.Norm2(&residual);
if( Comm.MyPID()==0 ) {
cout << "||b-Ax||_2 = " << residual << endl;
}
delete A;
delete Map;
// for testing purposes
if (residual > 1e-5)
exit(EXIT_FAILURE);
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(EXIT_SUCCESS);
}
#else
#include "ml_include.h"
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
#endif
puts("This MLAPI example requires the following configuration options:");
puts("\t--enable-epetra");
puts("\t--enable-teuchos");
puts("\t--enable-ifpack");
puts("\t--enable-amesos");
puts("\t--enable-galeri");
puts("Please check your configure line.");
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(EXIT_SUCCESS);
}
#endif
|
