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
|
/* ******************************************************************** */
/* See the file COPYRIGHT for a complete copyright notice, contact */
/* person and disclaimer. */
/* ******************************************************************** */
#include "ml_config.h"
#include "ml_common.h"
#if defined(HAVE_ML_MLAPI) && defined(HAVE_ML_GALERI)
#include "MLAPI_Space.h"
#include "MLAPI_Operator.h"
#include "MLAPI_MultiVector.h"
#include "MLAPI_Gallery.h"
#include "MLAPI_Expressions.h"
using namespace MLAPI;
// ============== //
// example driver //
// ============== //
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
#endif
try {
// initialize the MLAPI workspace
Init();
// Global dimension of our problem. It must be a square number
int NumGlobalRows = 16;
// Create the space in which vectors and operators live
Space MySpace(NumGlobalRows);
// Define vectors and operator, based on `MySpace'. The linear system
// matrix is created using the MLAPI Gallery, which is only a
// convenient wrapper of the Galeri package.
MultiVector x(MySpace), y(MySpace);
Operator A = Gallery("Tridiag", MySpace);
// We can now start coding the power method. We want a random vector of
// unitary 2-norm. First, we set random elements in the vector. Then,
// we divide each element by the 2-norm of the vector, computed as
// x * x. For two general vectors, x * y represents the scalar product
// between x and y (defined on the same space).
x.Random();
x = x / sqrt(x * x);
// ==================== //
// loop of power method //
// ==================== //
int MaxIters = 10;
for (int i = 0 ; i < MaxIters ; ++i) {
y = A * x; // matrix-vector product
// note that you need parenthesis in the following expression
// and that x * y is a global operation (therefore all
// processes must execute it)
double RQ = (y * x) / (x * x);
if (GetMyPID() == 0)
std::cout << "iter = " << i << ", RQ = " << RQ << std::endl;
x = y / sqrt(y * y);
}
// finalize the MLAPI workspace
Finalize();
}
catch (const int e) {
std::cout << "Caught integer exception, code = " << e << std::endl;
}
catch (...) {
std::cout << "problems here..." << std::endl;
}
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(0);
}
#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(0);
}
#endif // #if defined(HAVE_ML_MLAPI)
|
