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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
#ifndef __PYARPACKSERVICES_HPP__
#define __PYARPACKSERVICES_HPP__
#include <vector>
#include <string>
#include <complex>
#include <iostream>
#include <cmath> // sqrt.
#include <Eigen/Sparse>
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace bn = boost::python::numpy;
#define ARPACKSOLVERMEMBER(pyarpackSolver) \
.def_readwrite("symPb", &pyarpackSolver<RC, FD, EM, SLV>::symPb, \
"symmetric problem - default: true") \
.def_readwrite("nbEV", &pyarpackSolver<RC, FD, EM, SLV>::nbEV, \
"number of eigen vectors to find - default: 1") \
.def_readwrite("nbCV", &pyarpackSolver<RC, FD, EM, SLV>::nbCV, \
"number of dimensions of the workspace - default: 3") \
.def_readwrite("tol", &pyarpackSolver<RC, FD, EM, SLV>::tol, \
"tolerance - default: 1.e-6") \
.def_readwrite("sigmaReal", &pyarpackSolver<RC, FD, EM, SLV>::sigmaReal, \
"shift over real axis - default: 0.") \
.def_readwrite("sigmaImag", &pyarpackSolver<RC, FD, EM, SLV>::sigmaImag, \
"shift over imaginary axis - default: 0.") \
.def_readwrite("dumpToFile", &pyarpackSolver<RC, FD, EM, SLV>::dumpToFile, \
"dump eigen vectors to arpackSolver.*.out files - default: false") \
.def_readwrite("restartFromFile", &pyarpackSolver<RC, FD, EM, SLV>::restartFromFile, \
"restart from eigen vectors found in arpackSolver.*.out files - default: false") \
.def_readwrite("mag", &pyarpackSolver<RC, FD, EM, SLV>::mag, \
"magnitude - default: LM") \
.def_readwrite("maxIt", &pyarpackSolver<RC, FD, EM, SLV>::maxIt, \
"maximum number of arpack iterations - default: 100") \
.def_readwrite("schur", &pyarpackSolver<RC, FD, EM, SLV>::schur, \
"compute schur vectors - default: false") \
.def_readwrite("verbose", &pyarpackSolver<RC, FD, EM, SLV>::verbose, \
"verbosity level - default: 0") \
.def_readonly ("stdPb", &pyarpackSolver<RC, FD, EM, SLV>::stdPb, \
"standard or generalised problem - default: true") \
.def_readonly ("val", &pyarpackSolver<RC, FD, EM, SLV>::val, \
"eigen values found") \
.def_readonly ("vec", &pyarpackSolver<RC, FD, EM, SLV>::vec, \
"eigen vectors found") \
.def_readonly ("mode", &pyarpackSolver<RC, FD, EM, SLV>::mode, \
"selected arpack mode (according to input options: std/gen, shift, ...)") \
.def_readonly ("nbIt", &pyarpackSolver<RC, FD, EM, SLV>::nbIt, \
"number of arpack iterations") \
.def_readonly ("imsTime", &pyarpackSolver<RC, FD, EM, SLV>::imsTime, \
"time spent to initialize the mode solver if needed") \
.def_readonly ("rciTime", &pyarpackSolver<RC, FD, EM, SLV>::rciTime, \
"time spent in Reverse Communication Interface") \
.def_readwrite("debug", &pyarpackSolver<RC, FD, EM, SLV>::debug, \
"debug traces (up to 3) - default: 0") \
#define ARPACKSOLVERDEBUGSTAT() \
if (debug > 3) debug = 3; \
debug_c(6, -6, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, \
debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug); \
stat_c(nopx, nbx, nrorth, nitref, nrstrt, tsaupd, tsaup2, \
tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, tnaup2, \
tnaitr, tneigt, tngets, tnapps, tnconv, tcaupd, tcaup2, \
tcaitr, tceigt, tcgets, tcapps, tcconv, tmvopx, tmvbx, \
tgetv0, titref, trvec); \
void pyarpackThrowError(std::string const & msg) {
std::string const info = "Error: " + msg;
std::cerr << info << std::endl;
PyErr_SetString(PyExc_IndexError, info.c_str());
bp::throw_error_already_set();
};
template<typename RC, typename EM>
class pyarpackServices {
// Public methods.
public:
static int buildSparseMatrice(bp::tuple const & T, Eigen::SparseMatrix<RC> & M,
a_int const & debug, std::string const & msg) {
// Get boost data as C++ data.
if (bp::len(T) != 4) {pyarpackThrowError(msg + " must be a 3-tuple"); return 1;}
bp::extract<int> nExt(T[0]);
bp::extract<bn::ndarray> iExt(T[1]);
bp::extract<bn::ndarray> jExt(T[2]);
bp::extract<bn::ndarray> mijExt(T[3]);
if (! nExt.check()) {pyarpackThrowError(msg + "[0] must be an integer" ); return 1;}
if (! iExt.check()) {pyarpackThrowError(msg + "[1] must be numpy.array"); return 1;}
if (! jExt.check()) {pyarpackThrowError(msg + "[2] must be numpy.array"); return 1;}
if (!mijExt.check()) {pyarpackThrowError(msg + "[3] must be numpy.array"); return 1;}
bn::ndarray iArray = iExt();
bn::ndarray jArray = jExt();
bn::ndarray mijArray = mijExt();
if (iArray.get_dtype() != bn::dtype::get_builtin<a_int>()) {pyarpackThrowError(msg + "[1] type is not consistent"); return 1;}
if (jArray.get_dtype() != bn::dtype::get_builtin<a_int>()) {pyarpackThrowError(msg + "[2] type is not consistent"); return 1;}
if (mijArray.get_dtype() != bn::dtype::get_builtin<RC>() ) {pyarpackThrowError(msg + "[3] type is not consistent with arpack type"); return 1;}
a_int iSz = iArray.shape(0);
a_int * iPtr = reinterpret_cast<a_int*>(iArray.get_data());
a_int jSz = jArray.shape(0);
a_int * jPtr = reinterpret_cast<a_int*>(jArray.get_data());
a_int mSz = mijArray.shape(0);
RC * mPtr = reinterpret_cast<RC*>(mijArray.get_data());
if (iSz != jSz) {pyarpackThrowError(msg + "[1] and " + msg + "[2] must have same lenght"); return 1;}
if (iSz != mSz) {pyarpackThrowError(msg + "[1] and " + msg + "[3] must have same lenght"); return 1;}
// Debug on demand: casting value on numpy.append is MANDATORY or C++ won't get the expected type..
for (auto k = 0; debug && k < mSz; k++) {
std::cout << "pyarpackServices::buildSparseMatrice - " << msg << "[" << iPtr[k] << ", " << jPtr[k] << "] = " << mPtr[k] << std::endl;
};
// Build sparse matrice.
a_uint n = nExt();
a_uint iMin = n+1, jMin = n+1;
for (auto k = 0; k < mSz; k++) {
if (iPtr[k] < iMin) iMin = iPtr[k];
if (jPtr[k] < jMin) jMin = jPtr[k];
};
if (iMin != 0 && iMin != 1) {pyarpackThrowError(msg + ": smallest row indice must be 0 or 1"); return 1;}
if (jMin != 0 && jMin != 1) {pyarpackThrowError(msg + ": smallest column indice must be 0 or 1"); return 1;}
a_int iBased = 0, jBased = 0;
if (iMin == 1) iBased = 1;
if (jMin == 1) jBased = 1;
M = Eigen::SparseMatrix<RC>(n, n); // Set matrice dimensions.
std::vector<Eigen::Triplet<RC>> triplets;
a_uint nnz = mSz;
triplets.reserve(nnz);
for (auto k = 0; k < nnz; k++) triplets.emplace_back(iPtr[k] - iBased, jPtr[k] - jBased, mPtr[k]);
M.setFromTriplets(triplets.begin(), triplets.end()); // Set all (i, j, Mij).
return 0;
};
static int buildDenseMatrice(bp::tuple const & T, Eigen::Matrix<RC, Eigen::Dynamic, Eigen::Dynamic> & M,
a_int const & debug, std::string const & msg) {
// Get boost data as C++ data.
if (bp::len(T) != 2) {pyarpackThrowError(msg + " must be a 2-tuple"); return 1;}
bp::extract<bn::ndarray> mijExt(T[0]);
bp::extract<bool> oExt(T[1]);
if (!mijExt.check()) {pyarpackThrowError(msg + " must be numpy.array"); return 1;}
if ( !oExt.check()) {pyarpackThrowError(msg + " must be a boolean"); return 1;}
bn::ndarray mijArray = mijExt();
bool rowOrdered = oExt();
if (mijArray.get_dtype() != bn::dtype::get_builtin<RC>()) {pyarpackThrowError(msg + " type is not consistent with arpack type"); return 1;}
a_int mSz = mijArray.shape(0);
RC * mPtr = reinterpret_cast<RC*>(mijArray.get_data());
a_uint n = std::sqrt(mSz);
if (n*n != mSz) {pyarpackThrowError(msg + " must be a squared matrice"); return 1;}
// Debug on demand: casting value on numpy.append is MANDATORY or C++ won't get the expected type..
for (auto k = 0; debug && k < mSz; k++) {
std::cout << "pyarpackServices::buildDenseMatrice - " << msg << "[" << k << "] = " << mPtr[k] << std::endl;
};
// Build dense matrice.
M = Eigen::Matrix<RC, Eigen::Dynamic, Eigen::Dynamic>(n, n); // Set matrice dimensions.
M.setZero(n, n); // Avoid spurious/random values which may break solves (LU, QR, ...).
if (rowOrdered) {
for (size_t k = 0; k < n; k++) {
for (size_t l = 0; l < n; l++) M(k, l) = mPtr[l+k*n];
}
}
else {
for (size_t l = 0; l < n; l++) {
for (size_t k = 0; k < n; k++) M(k, l) = mPtr[k+l*n];
}
}
return 0;
};
};
#endif
// Local Variables:
// mode: c++
// c-file-style:"stroustrup"
// show-trailing-whitespace: t
// End:
/* vim: set sw=2 ts=2 et smartindent :*/
|
