1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 // SUMMARY : ...
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : Frederic Hecht
21 // E-MAIL : frederic.hecht@sorbonne-unviersite.fr
22
23 /* clang-format off */
24 //ff-c++-LIBRARY-dep: mumps_seq blas libseq fc pthread
25 //ff-c++-cpp-dep:
26 /* clang-format on */
27
28 // F. Hecht december 2011
29 // ----------------------------
30 // file to add MUMPS sequentiel interface for sparce linear solver with dynamic load.
31
32 #include <iostream>
33 using namespace std;
34
35 #include "ff++.hpp"
36
37 #include "mumps_seq/mpi.h"
38 #include "dmumps_c.h"
39 #include "zmumps_c.h"
40
41 const int JOB_INIT = -1;
42 const int JOB_END = -2;
43 const int JOB_ANA = 1;
44 const int JOB_FAC = 2;
45 const int JOB_ANA_FAC = 4;
46 const int JOB_SOLVE = 3;
47 const int USE_COMM_WORLD = -987654;
48
49 template< typename RR >
50 struct MUMPS_STRUC_TRAIT {
51 typedef void MUMPS;
52 typedef void R;
53 };
54
55 template<>
56 struct MUMPS_STRUC_TRAIT< double > {
57 typedef DMUMPS_STRUC_C MUMPS;
58 typedef double R;
59 };
60
61 template<>
62 struct MUMPS_STRUC_TRAIT< Complex > {
63 typedef ZMUMPS_STRUC_C MUMPS;
64 typedef ZMUMPS_COMPLEX R;
65 };
66
mumps_c(DMUMPS_STRUC_C * id)67 void mumps_c(DMUMPS_STRUC_C *id) { dmumps_c(id); }
68
mumps_c(ZMUMPS_STRUC_C * id)69 void mumps_c(ZMUMPS_STRUC_C *id) { zmumps_c(id); }
70
71 // template<typename R>
72 template< class R = double >
73 class SolveMUMPS_seq : public VirtualSolver< int, R > {
74 public:
75 static const int orTypeSol;
76 typedef HashMatrix< int, R > HMat;
77 typedef R K; //
78 HMat &A;
79
80 // typedef double R;
81 long verb;
82 double eps;
83 mutable double epsr;
84 double tgv;
85 int cn, cs;
86 typedef typename MUMPS_STRUC_TRAIT< R >::R MR;
87 mutable typename MUMPS_STRUC_TRAIT< R >::MUMPS id;
88 KN< double > *rinfog;
89 KN< long > *infog;
90
ICNTL(int i) const91 int &ICNTL(int i) const { return id.icntl[i - 1]; }
CNTL(int i) const92 double &CNTL(int i) const { return id.cntl[i - 1]; }
INFO(int i) const93 int &INFO(int i) const { return id.info[i - 1]; }
RINFO(int i) const94 double &RINFO(int i) const { return id.rinfo[i - 1]; }
INFOG(int i) const95 int &INFOG(int i) const { return id.infog[i - 1]; }
RINFOG(int i) const96 double &RINFOG(int i) const { return id.rinfog[i - 1]; }
97
SetVerb() const98 void SetVerb( ) const {
99 ICNTL(1) = 6; // output stream for error messages.
100 ICNTL(2) = 6; // stream for diagnostic printing, statistics, and warning messages.
101 ICNTL(3) = 6; // output stream global information, collected on the host.
102 ICNTL(4) = min(max(verb - 2, 1L), 4L); // the level of printing for error, warning, and diag
103 if (verb == 0) ICNTL(4) = 0;
104 ICNTL(11) = 0; // noerroranalysisisperformed(nostatistics).
105 if (id.job == JOB_SOLVE &&
106 verb > 99) { // computes statistics related to an error analysis of the linear system
107 if (verb > 999)
108 ICNTL(11) = 1; // All Stat (veryexpensive)
109 else
110 ICNTL(11) = 2; // compute main statistics
111 }
112 }
Clean()113 void Clean( ) {
114 delete[] id.irn;
115 delete[] id.jcn;
116 delete[] id.a;
117 id.irn = 0;
118 id.jcn = 0;
119 id.a = 0;
120 }
to_mumps_mat()121 void to_mumps_mat( ) {
122 Clean( );
123
124 id.nrhs = 0; //
125 int n = A.n;
126 int nz = A.nnz;
127 ffassert(A.n == A.m);
128
129 int *irn = new int[nz];
130 int *jcn = new int[nz];
131 R *a = new R[nz];
132 A.CSR( );
133
134 for (int i = 0; i < n; ++i) {
135 for (int k = A.p[i]; k < A.p[i + 1]; ++k) {
136 irn[k] = i + 1;
137 jcn[k] = A.j[k] + 1;
138 a[k] = A.aij[k];
139 }
140 }
141
142 id.n = n;
143 id.nz = nz;
144 id.irn = irn;
145 id.jcn = jcn;
146 id.a = (MR *)(void *)a;
147 id.rhs = 0;
148 ffassert(A.half == id.sym); //
149 ICNTL(5) = 0; // input matrix type
150 ICNTL(7) = 7; // NUMBERING ...
151
152 ICNTL(9) = 1; // 1: A x = b, !1 : tA x = b during slove phase
153 ICNTL(18) = 0;
154 }
Check(const char * msg="mumps_seq")155 void Check(const char *msg = "mumps_seq") {
156 if (INFO(1) != 0) {
157 cout << " Erreur Mumps seq: number " << INFO(1) << endl;
158 cout << " Fatal Erreur " << msg << endl;
159 Clean( );
160 id.job = JOB_END;
161 mumps_c(&id); /* Terminate instance */
162 ErrorExec(msg, INFO(1));
163 }
164 }
CopyInfo()165 void CopyInfo( ) {
166 if (rinfog) {
167 // copy rinfog
168 if (rinfog->N( ) < 40) {
169 rinfog->resize(40);
170 }
171
172 for (int i = 0; i < 40; ++i) {
173 (*rinfog)[i] = RINFOG(i + 1);
174 }
175 }
176
177 if (infog) {
178 // copy ginfo
179 if (infog->N( ) < 40) {
180 infog->resize(40);
181 }
182
183 for (int i = 0; i < 40; ++i) {
184 (*infog)[i] = INFOG(i + 1);
185 }
186 }
187 }
SolveMUMPS_seq(HMat & AA,const Data_Sparse_Solver & ds,Stack stack)188 SolveMUMPS_seq(HMat &AA, const Data_Sparse_Solver &ds, Stack stack)
189 : A(AA), verb(ds.verb), eps(ds.epsilon), epsr(0), tgv(ds.tgv), cn(0), cs(0), rinfog(ds.rinfo),
190 infog(ds.info) {
191
192 int myid = 0;
193
194 id.irn = 0;
195 id.jcn = 0;
196 id.a = 0;
197
198 id.job = JOB_INIT;
199 id.par = 1;
200 id.sym = A.half;
201 id.comm_fortran = USE_COMM_WORLD;
202 SetVerb( );
203 mumps_c(&id);
204
205 Check("MUMPS_seq build/init");
206 if (verbosity > 3) {
207 cout << " -- MUMPS n= " << id.n << ", peak Mem: " << INFOG(22) << " Mb"
208 << " sym: " << id.sym << endl;
209 }
210 }
211
~SolveMUMPS_seq()212 ~SolveMUMPS_seq( ) {
213 Clean( );
214 id.job = JOB_END;
215 SetVerb( );
216 mumps_c(&id); /* Terminate instance */
217 /*int ierr = */
218 }
219
dosolver(K * x,K * b,int N,int trans)220 void dosolver(K *x, K *b, int N, int trans) {
221 if (verbosity > 1) {
222 cout << " -- MUMPS solve, peak Mem : " << INFOG(22) << " Mb, n = " << id.n
223 << " sym =" << id.sym << " trans = " << trans << endl;
224 }
225 ICNTL(9) = trans == 0; // 1: A x = b, !1 : tA x = b during slove phase
226 id.nrhs = N;
227 myscopy(id.n, b, x);
228 id.rhs = (MR *)(void *)(R *)x;
229 id.job = JOB_SOLVE; // performs the analysis. and performs the factorization.
230 SetVerb( );
231 mumps_c(&id);
232 Check("MUMPS_seq dosolver");
233
234 if (verb > 9) {
235
236 for (int j = 0; j < N; ++j) {
237 KN_< R > B(b + j * id.n, id.n);
238 cout << j << " b linfty " << B.linfty( ) << endl;
239 }
240 }
241
242 if (verb > 2) {
243
244 for (int j = 0; j < N; ++j) {
245 KN_< R > B(x + j * id.n, id.n);
246 cout << " x " << j << " linfty " << B.linfty( ) << endl;
247 }
248 }
249 CopyInfo( );
250 }
251
fac_init()252 void fac_init( ) { to_mumps_mat( ); } // n, nzz fixe
fac_symbolic()253 void fac_symbolic( ) {
254 id.job = JOB_ANA;
255 SetVerb( );
256 mumps_c(&id);
257 Check("MUMPS_seq Analyse");
258 CopyInfo( );
259 }
fac_numeric()260 void fac_numeric( ) {
261 id.job = JOB_FAC;
262 SetVerb( );
263 mumps_c(&id);
264 Check("MUMPS_seq Factorize");
265 CopyInfo( );
266 }
UpdateState()267 void UpdateState( ) {
268 if (A.GetReDoNumerics( )) cn++;
269 if (A.GetReDoSymbolic( )) cs++;
270 this->ChangeCodeState(A.n, cs, cn);
271 }
272 };
273 struct InitEnd {
InitEndInitEnd274 InitEnd( ) {
275 cout << "init MUMPS_SEQ: MPI_Init" << endl;
276 int argc = 0;
277 char **argv = 0;
278 // BOF BOF
279 MPI_Init(&argc, &argv);
280 }
~InitEndInitEnd281 ~InitEnd( ) {
282 cout << "close MUMPS_SEQ: MPI_Finalize" << endl;
283 MPI_Finalize( );
284 }
285 };
286 static InitEnd global; // To init SEQ MPI ????
287
288 // 1 unsym , 2 herm, 4 sym, 8 pos , 16 nopos, 32 seq, 64 ompi, 128 mpi
289 template<> const int SolveMUMPS_seq<double>::orTypeSol = 1|2|4|8|16|32;
290 template<> const int SolveMUMPS_seq<std::complex<double>>::orTypeSol = 1|4|8|16|32;
291
Load_Init()292 static void Load_Init( ) {
293 addsolver< SolveMUMPS_seq< double > >("MUMPS", 50, 1);
294 addsolver< SolveMUMPS_seq< Complex > >("MUMPS", 50, 1);
295 addsolver< SolveMUMPS_seq< double > >("MUMPSSEQ", 50, 1);
296 addsolver< SolveMUMPS_seq< Complex > >("MUMPSSEQ", 50, 1);
297 setptrstring(def_solver, "MUMPSSEQ");
298 }
299 LOADFUNC(Load_Init)
300