1 #include <stdio.h>
2 #include <string.h>
3 #include <math.h>
4 #include "GetDPConfig.h"
5 #include "Sparskit.h"
6 #include "MallocUtils.h"
7 #include "Message.h"
8
9 #if defined(HAVE_UNDERSCORE)
10 #define etime_ etime
11 #define ilut_ ilut
12 #define ilutp_ ilutp
13 #define ilud_ ilud
14 #define iludp_ iludp
15 #define iluk_ iluk
16 #define ilu0_ ilu0
17 #define milu0_ milu0
18 #define cmkreord_ cmkreord
19 #define sortcol_ sortcol
20 #define skit_ skit
21 #define psplot_ psplot
22 #define cg_ cg
23 #define cgnr_ cgnr
24 #define bcg_ bcg
25 #define dbcg_ dbcg
26 #define bcgstab_ bcgstab
27 #define tfqmr_ tfqmr
28 #define fom_ fom
29 #define gmres_ gmres
30 #define fgmres_ fgmres
31 #define dqgmres_ dqgmres
32 #define amux_ amux
33 #define atmux_ atmux
34 #define lusol_ lusol
35 #define lutsol_ lutsol
36 #define csrcoo_ csrcoo
37 #define ma28ad_ ma28ad
38 #define ma28cd_ ma28cd
39 #define dnrm2_ dnrm2
40 #define flu_ flu
41 #define pgmres_ pgmres
42 #define getdia_ getdia
43 #define amudia_ amudia
44 #define diamua_ diamua
45 #define rnrms_ rnrms
46 #define cnrms_ cnrms
47 #endif
48
49 /* Fortran prototypes */
50
51 extern "C" {
52 void ilut_ (int*,double*,int*,int*,int*,double*,
53 sscalar*,int*,int*,int*,double*,int*,int*);
54 void ilutp_ (int*,double*,int*,int*,int*,double*,
55 double*,int*,sscalar*,int*,
56 int*,int*,double*,int*,int*,int*);
57 void ilud_ (int*,double*,int*,int*,double*,
58 double*,sscalar*,int*,
59 int*,int*,double*,int*,int*);
60 void iludp_ (int*,double*,int*,int*,double*,
61 double*,double*,
62 int*,sscalar*,int*,int*,int*,
63 double*,int*,int*,int*);
64 void iluk_ (int*,double*,int*,int*,int*,
65 sscalar*,int*,int*,
66 int*,int*,double*,int*,int*);
67 void ilu0_ (int*,double*,int*,int*,sscalar*,int*,int*,int*,int*);
68 void milu0_ (int*,double*,int*,int*,sscalar*,int*,int*,int*,int*);
69 void cmkreord_ (int*,double*,int*,int*,double*,int*,int*,int*,
70 int*,int*,int*,int*,int*,int*);
71 void sortcol_ (int*,double*,int*,int*,int*,double*);
72 void skit_ (int*,double*,int*,int*,int*,int*,int*);
73 void psplot_ (int*,int*,int*,int*,int*);
74 void cg_ (int*,double*,double*,int*,double*,double*);
75 void cgnr_ (int*,double*,double*,int*,double*,double*);
76 void bcg_ (int*,double*,double*,int*,double*,double*);
77 void dbcg_ (int*,double*,double*,int*,double*,double*);
78 void bcgstab_ (int*,double*,double*,int*,double*,double*);
79 void tfqmr_ (int*,double*,double*,int*,double*,double*);
80 void fom_ (int*,double*,double*,int*,double*,double*);
81 void gmres_ (int*,double*,double*,int*,double*,double*);
82 void fgmres_ (int*,double*,double*,int*,double*,double*);
83 void dqgmres_ (int*,double*,double*,int*,double*,double*);
84 void amux_ (int*,double*,double*,double*,int*,int*);
85 void atmux_ (int*,double*,double*,double*,int*,int*);
86 void lusol_ (int*,double*,double*,sscalar*,int*,int*);
87 void lutsol_ (int*,double*,double*,sscalar*,int*,int*);
88 void csrcoo_ (int*,int*,int*,double*,int*,int*,int*,double*,int*,int*,int*);
89 void ma28ad_ (int*,int*,double*,int*,int*,int*,int*,double*,int*,int*,double*,int*);
90 void ma28cd_ (int*,double*,int*,int*,int*,double*,double*,int*);
91 double dnrm2_ (int*,double*,int*);
92 void flu_ (int*,double*,double*,double*,int*,int*,double*,double*,
93 double*,double*,double*);
94 void pgmres_ (int*,int*,double*,double*,double*,double*,
95 int*,int*,double*,int*,int*,
96 sscalar*,int*,int*,int*);
97 void getdia_ (int*,int*,int*,double*,int*,int*,int*,double*,int*,int*);
98 void diamua_ (int*,int*,double*,int*,int*,double*,double*,int*,int*);
99 void amudia_ (int*,int*,double*,int*,int*,double*,double*,int*,int*);
100 void rnrms_ (int*,int*,double*,int*,int*,double*);
101 void cnrms_ (int*,int*,double*,int*,int*,double*);
102 }
103
104 /* ------------------------------------------------------------------------ */
105 /* s o l v e */
106 /* ------------------------------------------------------------------------ */
107
solve_matrix(Matrix * M,Solver_Params * p,double * b,double * x)108 void solve_matrix (Matrix *M, Solver_Params *p, double *b, double *x){
109 FILE *pf;
110 double fpar[17];
111 double *a, *w, *rhs, *sol, *dx ;
112 double res;
113 int i, j, k, nnz, nnz_ilu, ierr, ipar[17];
114 int *ja, *ia, *jw, *mask, *levels;
115 int its, end, do_permute=0;
116 int zero=0, un=1, deux=2, six=6, douze=12, trente=30, trente_et_un=31;
117 int ROW=0, COLUMN=1;
118 double res1=1.;
119 int TrueNnz=0;
120
121 if (!M->N) {
122 Message::Warning("No equations in linear system");
123 return;
124 }
125
126 for(i=0 ; i<M->N ; i++){
127 if(b[i] != 0.) break;
128 if(i == M->N-1) {
129 Message::Warning("Null right hand side in linear system");
130 /*
131 for(i=0 ; i<M->N ; i++) x[i] = 0. ;
132 return ;
133 */
134 }
135 }
136
137 if(M->T == DENSE){
138
139 if(p->Algorithm == LU){
140
141 Message::Info("Dense LU decomposition");
142 print_matrix_info_DENSE(M->N);
143
144 sol = (double*)Malloc(M->N * sizeof(double));
145 w = (double*)Malloc(M->N * sizeof(double));
146 dx = (double*)Malloc(M->N * sizeof(double));
147
148 ipar[1] = p->Re_Use_LU;
149 ipar[2] = p->Iterative_Improvement;
150 ipar[3] = p->Matrix_Printing;
151 ipar[4] = p->Nb_Iter_Max;
152
153 fpar[1] = p->Stopping_Test;
154
155 flu_(&ipar[1], &fpar[1], M->F.a, M->F.lu, &M->N, &M->N, b, x, dx, sol, w);
156
157 Free(sol);
158 Free(w);
159 Free(dx);
160
161 return ;
162 }
163
164 Message::Info("Dense to sparse matrix conversion") ;
165
166 nnz = M->N * M->N ;
167
168 M->S.a = List_Create(1, 1, sizeof(double));
169 M->S.a->n = nnz ;
170 M->S.a->array = (char*) M->F.a ;
171
172 M->S.jptr = List_Create(M->N+1, M->N, sizeof(int));
173 M->S.ai = List_Create(nnz, M->N, sizeof(int));
174
175 for(i=1 ; i<=nnz ; i+=M->N){
176 List_Add(M->S.jptr, &i) ;
177 for(j=1 ; j<=M->N ; j++){
178 List_Add(M->S.ai, &j) ;
179 }
180 }
181
182 i = nnz + 1 ; List_Add(M->S.jptr, &i) ;
183
184 if(M->changed){
185 do_permute = 1 ;
186 M->changed = 0 ;
187 }
188
189 for (i=0 ; i<nnz ; i++) if (M->F.a[i]) TrueNnz++ ;
190 Message::Info("Number of nonzeros %d/%d (%.4f)",TrueNnz, nnz, (double)TrueNnz/(double)nnz);
191
192 Message::Cpu("");
193
194 } /* if DENSE */
195 else{
196
197 nnz = List_Nbr(M->S.a);
198 if(M->changed){
199 do_permute = 1 ;
200 csr_format (&M->S, M->N);
201 restore_format (&M->S);
202 M->changed = 0 ;
203 }
204
205 } /* if SPARSE */
206
207 a = (double*) M->S.a->array;
208 ia = (int*) M->S.jptr->array;
209 ja = (int*) M->S.ai->array;
210
211 if(p->Scaling != NONE){
212 Message::Info("Scaling system of equations") ;
213 scale_matrix (p->Scaling, M) ;
214 scale_vector (ROW, M, b) ;
215 }
216 else{
217 Message::Info("No scaling of system of equations") ;
218 }
219
220 for (i=1; i<M->N; i++) {
221 if(ia[i]-ia[i-1] <= 0)
222 Message::Error("Zero row in matrix");
223 }
224
225 rhs = (double*) Malloc(M->N * sizeof(double));
226 sol = (double*) Calloc(M->N, sizeof(double));
227
228 /* Renumbering */
229
230 if (!M->ILU_Exists){
231 M->S.permr = (int*) Malloc(M->N * sizeof(int));
232 M->S.rpermr = (int*) Malloc(M->N * sizeof(int));
233 M->S.permp = (int*) Malloc(2 * M->N * sizeof(int));
234 }
235
236 if(do_permute || !M->ILU_Exists){
237 for(i=0 ; i<M->N ; i++) {
238 M->S.permr[i] = M->S.rpermr[i] = M->S.permp[i+M->N] = i+1;
239 }
240 switch (p->Renumbering_Technique){
241 case NONE:
242 Message::Info("No renumbering");
243 break;
244 case RCMK:
245 Message::Info("RCMK algebraic renumbering");
246 if(!M->ILU_Exists){
247 M->S.a_rcmk = (double*) Malloc(nnz * sizeof(double));
248 M->S.ia_rcmk = (int*) Malloc((M->N + 1) * sizeof(int));
249 M->S.ja_rcmk = (int*) Malloc(nnz * sizeof(int));
250 }
251 mask = (int*) Malloc(nnz * sizeof(int));
252 levels = (int*) Malloc((M->N + 1) * sizeof(int));
253 i = j = k = 1;
254 cmkreord_(&M->N, a, ja, ia, M->S.a_rcmk, M->S.ja_rcmk, M->S.ia_rcmk,
255 &i, M->S.permr, mask, &j, &k, M->S.rpermr, levels);
256 w = (double*) Malloc(nnz * sizeof(double));
257 sortcol_(&M->N, M->S.a_rcmk, M->S.ja_rcmk, M->S.ia_rcmk, mask, w);
258 Free(w); Free(mask); Free(levels);
259 break;
260 default :
261 Message::Error("Unknown renumbering technique");
262 break;
263 }
264 print_matrix_info_CSR(M->N, ia, ja);
265 Message::Cpu("");
266 }
267
268 if(p->Renumbering_Technique == RCMK){
269 if (p->Re_Use_ILU && !M->ILU_Exists && !do_permute){
270 /*
271 This is incorrect if M is to be changed during the process, and
272 we still want to keep the same precond.
273
274 Free(M->S.a->array) ;
275 Free(M->S.jptr->array) ;
276 Free(M->S.ai->array) ;
277 M->S.a->array = (char*)M->S.a_rcmk;
278 M->S.jptr->array = (char*)M->S.ia_rcmk;
279 M->S.ai->array = (char*)M->S.ja_rcmk;
280 */
281 }
282 a = M->S.a_rcmk;
283 ia = M->S.ia_rcmk;
284 ja = M->S.ja_rcmk;
285 }
286
287
288 if (p->Matrix_Printing == 1 || p->Matrix_Printing == 3) {
289 Message::Info("Matrix printing");
290 skit_(&M->N, a, ja, ia, &douze, &douze, &ierr);
291 pf = fopen("fort.13","w");
292 for (i=0 ; i<M->N ; i++) fprintf(pf, "%d %22.15E\n", i+1, b[i]);
293 fclose(pf);
294 psplot_(&M->N, ja, ia, &trente, &zero);
295 }
296
297 /* Incomplete factorizations */
298
299 if (!M->ILU_Exists) {
300
301 if (p->Re_Use_ILU) M->ILU_Exists = 1;
302
303 #if defined(HAVE_ILU_FLOAT)
304 #define ILUSTORAGE "Float"
305 #else
306 #define ILUSTORAGE "Double"
307 #endif
308
309 end = 0 ;
310
311 switch (p->Preconditioner){
312
313 case ILUT : Message::Info("ILUT (%s, fill-in = %d)", ILUSTORAGE, p->Nb_Fill);
314 nnz_ilu = 2 * (M->N+1) * (p->Nb_Fill+1); break;
315
316 case ILUTP : Message::Info("ILUTP (%s, fill-in = %d)", ILUSTORAGE, p->Nb_Fill);
317 nnz_ilu = 2 * (M->N+1) * (p->Nb_Fill+1); break;
318
319 case ILUD : Message::Info("ILUD (%s)", ILUSTORAGE);
320 /* first guess */
321 nnz_ilu = List_Nbr(M->S.a); break;
322
323 case ILUDP : Message::Info("ILUDP (%s)", ILUSTORAGE);
324 /* first guess */
325 nnz_ilu = List_Nbr(M->S.a); break ;
326
327 case ILUK : Message::Info("ILU%d (%s)", p->Nb_Fill, ILUSTORAGE);
328 /* exact for nbfill=0, first guess otherwise */
329 nnz_ilu = (p->Nb_Fill+1) * List_Nbr(M->S.a) + (M->N+1);
330 break;
331
332 case ILU0 : Message::Info("ILU0 (%s)", ILUSTORAGE);
333 nnz_ilu = List_Nbr(M->S.a) + (M->N+1); break;
334
335 case MILU0 : Message::Info("MILU0 (%s)", ILUSTORAGE);
336 nnz_ilu = List_Nbr(M->S.a) + (M->N+1); break;
337
338 case DIAGONAL :
339 Message::Info("Diagonal scaling (%s)", ILUSTORAGE);
340 M->S.alu = (sscalar*) Malloc((M->N+1) * sizeof(sscalar));
341 M->S.jlu = (int*) Malloc((M->N+1) * sizeof(int));
342 M->S.ju = (int*) Malloc((M->N+1) * sizeof(int));
343
344 for (i=0 ; i<M->N ; i++) {
345 M->S.alu[i] = 1.0 ;
346 M->S.jlu[i] = M->N+2 ;
347 M->S.ju[i] = M->N+2 ;
348 }
349 M->S.alu[M->N] = 0.0 ;
350 M->S.jlu[M->N] = M->N+2 ;
351 M->S.ju[M->N] = 0 ;
352 end = 1;
353 ierr = 0;
354 break;
355
356 case NONE :
357 Message::Info("No ILU");
358 end = 1;
359 ierr = 0;
360 break ;
361
362 default :
363 Message::Error("Unknown ILU method");
364 break;
365 }
366
367 if(!end){
368 M->S.alu = (sscalar*) Malloc(nnz_ilu * sizeof(sscalar));
369 M->S.jlu = (int*) Malloc(nnz_ilu * sizeof(int));
370 M->S.ju = (int*) Malloc((M->N+1) * sizeof(int));
371 }
372
373 reallocate :
374
375 switch(p->Preconditioner){
376 case ILUT :
377 w = (double*) Malloc((M->N+1) * sizeof(double));
378 jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));
379 ilut_(&M->N, a, ja, ia, &p->Nb_Fill, &p->Dropping_Tolerance,
380 M->S.alu, M->S.jlu, M->S.ju, &nnz_ilu, w, jw, &ierr);
381 Free(w); Free(jw); break;
382
383 case ILUTP :
384 w = (double*) Malloc((M->N+1) * sizeof(double));
385 jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));
386 ilutp_(&M->N, a, ja, ia, &p->Nb_Fill, &p->Dropping_Tolerance,
387 &p->Permutation_Tolerance, &M->N, M->S.alu, M->S.jlu,
388 M->S.ju, &nnz_ilu, w, jw, M->S.permp, &ierr);
389 Free(jw); Free(w); break;
390
391 case ILUD :
392 w = (double*) Malloc((M->N+1) * sizeof(double));
393 jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));
394 ilud_(&M->N, a, ja, ia, &p->Diagonal_Compensation,
395 &p->Dropping_Tolerance, M->S.alu, M->S.jlu,
396 M->S.ju, &nnz_ilu, w, jw, &ierr);
397 Free(w); Free(jw); break;
398
399 case ILUDP :
400 w = (double*) Malloc((M->N+1) * sizeof(double));
401 jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));
402 iludp_(&M->N, a, ja, ia, &p->Diagonal_Compensation,
403 &p->Dropping_Tolerance, &p->Permutation_Tolerance,
404 &M->N, M->S.alu, M->S.jlu, M->S.ju, &nnz_ilu,
405 w, jw, M->S.permp, &ierr);
406 Free(jw); Free(w); break;
407
408 case ILUK :
409 levels = (int*) Malloc(nnz_ilu * sizeof(int));
410 w = (double*) Malloc((M->N+1) * sizeof(double));
411 jw = (int*) Malloc(3 * (M->N+1) * sizeof(int));
412 iluk_(&M->N, a, ja, ia, &p->Nb_Fill,
413 M->S.alu, M->S.jlu, M->S.ju,
414 levels, &nnz_ilu, w, jw, &ierr);
415 Free(levels); Free(w); Free(jw); break;
416
417 case ILU0 :
418 jw = (int*) Malloc((M->N+1) * sizeof(int));
419 ilu0_(&M->N, a, ja, ia, M->S.alu, M->S.jlu, M->S.ju, jw, &ierr);
420 Free(jw); break;
421
422 case MILU0 :
423 jw = (int*) Malloc((M->N+1) * sizeof(int));
424 milu0_(&M->N, a, ja, ia, M->S.alu, M->S.jlu, M->S.ju, jw, &ierr);
425 Free(jw); break;
426
427 }
428
429 switch (ierr){
430 case 0 :
431 break;
432 case -1 :
433 Message::Error("Input matrix may be wrong");
434 break;
435 case -2 : /* Matrix L in ILU overflows work array 'al' */
436 case -3 : /* Matrix U in ILU overflows work array 'alu' */
437 nnz_ilu += nnz_ilu/2 ;
438 Message::Info("Reallocating ILU (NZ: %d)", nnz_ilu);
439 Free(M->S.alu) ;
440 M->S.alu = (sscalar*) Malloc(nnz_ilu * sizeof(sscalar));
441 Free(M->S.jlu) ;
442 M->S.jlu = (int*) Malloc(nnz_ilu * sizeof(int));
443 goto reallocate ;
444 case -4 :
445 Message::Error("Illegal value of nb_fill in ILU");
446 break;
447 case -5 :
448 Message::Error("Zero row encountered in ILU");
449 break;
450 default :
451 Message::Error("Zero pivot on line %d in ILU",ierr);
452 break;
453 }
454
455 if(p->Preconditioner != NONE)
456 print_matrix_info_MSR(M->N, M->S.alu, M->S.jlu);
457
458 if(p->Matrix_Printing == 2 || p->Matrix_Printing == 3){
459 Message::Info("ILU printing");
460 psplot_(&M->N, M->S.jlu, M->S.jlu, &trente_et_un, &deux);
461 }
462
463 Message::Cpu("");
464
465 }
466
467 /* RHS reordering */
468
469 for(i=0;i<M->N;i++){
470 rhs[i] = b[M->S.rpermr[i] - 1];
471 }
472
473 /* Iterations */
474
475 ipar[1] = 0;
476 ipar[2] = (p->Preconditioner == NONE) ? 0 : p->Preconditioner_Position;
477 ipar[3] = 1;
478 ipar[4] = 0;
479 ipar[5] = p->Krylov_Size;
480 ipar[6] = p->Nb_Iter_Max;
481
482 fpar[1] = p->Stopping_Test;
483 fpar[2] = 0.0;
484 fpar[11] = 0.0;
485
486 switch (p->Algorithm){
487 case CG : Message::Info("Conjugate Gradient (CG)");
488 ipar[4] = 5 * M->N; break;
489 case CGNR : Message::Info("CG Normal Residual equation (CGNR)");
490 ipar[4] = 5 * M->N; break;
491 case BCG : Message::Info("Bi-Conjugate Gradient (BCG)");
492 ipar[4] = 7 * M->N; break;
493 case DBCG : Message::Info("BCG with partial pivoting (DBCG)");
494 ipar[4] = 11 * M->N; break;
495 case BCGSTAB : Message::Info("BCG stabilized (BCGSTAB)");
496 ipar[4] = 8 * M->N; break;
497 case TFQMR : Message::Info("Transpose-Free Quasi-Minimum Residual (TFQMR)");
498 ipar[4] = 11 * M->N; break;
499 case FOM : Message::Info("Full Orthogonalization Method (FOM)");
500 ipar[4] = (M->N+3) * (ipar[5]+2) + (ipar[5]+1) * ipar[5]/2; break;
501 case GMRES : Message::Info("Generalized Minimum RESidual (GMRES)");
502 ipar[4] = (M->N+3) * (ipar[5]+2) + (ipar[5]+1) * ipar[5]/2; break;
503 case FGMRES : Message::Info("Flexible version of Generalized Minimum RESidual (FGMRES)");
504 ipar[4] = 2*M->N * (ipar[5]+1) + (ipar[5]+1)*ipar[5]/2 + 3*ipar[5] + 2; break;
505 case DQGMRES : Message::Info("Direct version of Quasi Generalize Minimum RESidual (DQGMRES)");
506 ipar[4] = M->N + (ipar[5]+1) * (2*M->N+4); break;
507 case PGMRES : Message::Info("Alternative Generalized Minimum RESidual (GMRES)");
508 ipar[4] = (M->N+4) * (ipar[5]+2) + (ipar[5]+1) * ipar[5]/2; break;
509 default : Message::Error("Unknown algorithm for sparse matrix solver"); break;
510 }
511
512 w = (double*) Malloc(ipar[4] * sizeof(double));
513
514 its = 0;
515 end = 0;
516 res = 0.0;
517
518 while(1){
519
520 switch(p->Algorithm){
521 case CG : cg_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
522 case CGNR : cgnr_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
523 case BCG : bcg_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
524 case DBCG : dbcg_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
525 case BCGSTAB : bcgstab_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
526 case TFQMR : tfqmr_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
527 case FOM : fom_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
528 case GMRES : gmres_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
529 case FGMRES : fgmres_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
530 case DQGMRES : dqgmres_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
531 case PGMRES : pgmres_ (&M->N, &p->Krylov_Size, rhs, sol, w, &p->Stopping_Test,
532 &p->Nb_Iter_Max, &six, a, ja, ia,
533 M->S.alu, M->S.jlu, M->S.ju, &ierr);
534 end = 1; break;
535 }
536
537 if(!end){
538
539 if(ipar[7] != its){
540 if(its) Message::Info(" %4d %.7e %.7e", its, res, res/res1);
541 its = ipar[7] ;
542 }
543
544 res = fpar[5];
545 if(its==1) res1 = fpar[5] ;
546
547 switch(ipar[1]){
548 case 1 :
549 amux_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], a, ja, ia); break;
550 case 2 :
551 atmux_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], a, ja, ia); break;
552 case 3 : case 5 :
553 lusol_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], M->S.alu, M->S.jlu, M->S.ju); break;
554 case 4 : case 6 :
555 lutsol_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], M->S.alu, M->S.jlu, M->S.ju); break;
556 case 0 :
557 end = 1; break;
558 case -1 :
559 Message::Warning("Iterative solver has iterated too many times"); end = 1; break;
560 case -2 :
561 Message::Warning("Iterative solver was not given enough work space");
562 Message::Warning("The work space should at least have %d elements", ipar[4]);
563 end = 1; break;
564 case -3 :
565 Message::Warning("Iterative solver is facing a break-down"); end = 1; break;
566 default :
567 Message::Warning("Iterative solver terminated (code = %d)", ipar[1]); end = 1; break;
568 }
569
570 }
571 if(end) break;
572
573 }
574
575
576 /* Convergence results monitoring */
577
578 Message::Info(" %4d %.7e %.7e", ipar[7], fpar[6], fpar[6]/res1);
579
580 amux_(&M->N, sol, w, a, ja, ia);
581
582
583 for(i=0 ; i<M->N ; i++){
584 w[M->N+i] = sol[i] - 1.0 ;
585 w[i] -= rhs[i] ;
586 }
587
588
589 Message::Info("%d Iterations / Residual: %g", ipar[7], dnrm2_(&M->N,w,&un));
590 /*
591 Message::Info("Conv. Rate: %g, |Res|: %g, |Err|: %g",
592 fpar[7], dnrm2_(&M->N,w,&un), dnrm2_(&M->N,&w[M->N],&un));
593 */
594 Free(w);
595
596 /* Inverse renumbering */
597
598 for (i=0;i<M->N;i++) {
599 j = M->S.permr[i] - 1;
600 k = M->S.permp[j+M->N] - 1;
601 x[i] = sol[k];
602 }
603
604 /* Free memory */
605
606 Free(rhs);
607 Free(sol);
608
609 if (!M->ILU_Exists){
610 if(p->Preconditioner != NONE) {
611 Free(M->S.alu);
612 Free(M->S.jlu);
613 Free(M->S.ju);
614 }
615 if (p->Renumbering_Technique == RCMK) {
616 Free(M->S.rpermr);
617 Free(M->S.permr);
618 Free(M->S.permp);
619 Free(M->S.a_rcmk);
620 Free(M->S.ia_rcmk);
621 Free(M->S.ja_rcmk);
622 }
623 }
624
625 if(M->T == DENSE){
626 List_Delete(M->S.a);
627 List_Delete(M->S.jptr);
628 List_Delete(M->S.ai);
629 }
630
631 if(p->Scaling)
632 scale_vector (COLUMN, M, x) ;
633
634 }
635
636 /* ------------------------------------------------------------------------ */
637 /* p r i n t */
638 /* ------------------------------------------------------------------------ */
639
print_parametres(Solver_Params * p)640 void print_parametres (Solver_Params *p){
641 printf(" Matrix_Format : %d\n", p->Matrix_Format);
642 printf(" Matrix_Printing : %d\n", p->Matrix_Printing);
643 printf(" Renumbering_Technique : %d\n", p->Renumbering_Technique);
644 printf(" Preconditioner : %d\n", p->Preconditioner);
645 printf(" Preconditioner_Position : %d\n", p->Preconditioner_Position);
646 printf(" Nb_Fill : %d\n", p->Nb_Fill);
647 printf(" Dropping_Tolerance : %g\n", p->Dropping_Tolerance);
648 printf(" Permutation_Tolerance : %g\n", p->Permutation_Tolerance);
649 printf(" Diagonal_Compensation : %g\n", p->Diagonal_Compensation);
650 printf(" Algorithm : %d\n", p->Algorithm);
651 printf(" Krylov_Size : %d\n", p->Krylov_Size);
652 printf(" IC_Acceleration : %g\n", p->IC_Acceleration);
653 printf(" Iterative_Improvement : %d\n", p->Iterative_Improvement);
654 printf(" Nb_Iter_Max : %d\n", p->Nb_Iter_Max);
655 printf(" Stopping_Test : %g\n", p->Stopping_Test);
656 }
657
658 /* ------------------------------------------------------------------------ */
659 /* i n i t */
660 /* ------------------------------------------------------------------------ */
661
init_matrix(int NbLines,Matrix * M,Solver_Params * p)662 void init_matrix (int NbLines, Matrix *M, Solver_Params *p){
663 int i, j=0;
664
665 M->T = p->Matrix_Format;
666 M->N = NbLines;
667 M->changed = 1 ;
668 M->ILU_Exists = 0;
669 M->notranspose = 0 ;
670 M->scaled = 0 ;
671
672 switch (M->T) {
673 case SPARSE :
674 M->S.a = List_Create (NbLines, NbLines, sizeof(double));
675 M->S.ai = List_Create (NbLines, NbLines, sizeof(int));
676 M->S.ptr = List_Create (NbLines, NbLines, sizeof(int));
677 M->S.jptr = List_Create (NbLines+1, NbLines, sizeof(int));
678 /* '+1' indispensable: csr_format ecrit 'nnz+1' dans jptr[NbLine] */
679 for(i=0; i<NbLines; i++) List_Add (M->S.jptr, &j);
680 break;
681 case DENSE :
682 M->F.LU_Exist = 0;
683 /* Tous les algos iteratifs sont programmes pour resoudre
684 A^T x = b... C'est tres con, mais bon. L'algo LU est le seul
685 qui demande la vraie matrice en entree... */
686 if(p->Algorithm == LU){
687 M->F.lu = (double*) Malloc (NbLines * NbLines * sizeof(double));
688 M->notranspose = 1 ;
689 }
690 else
691 M->F.lu = NULL;
692 M->F.a = (double*) Malloc (NbLines * NbLines * sizeof(double));
693 break;
694 default :
695 Message::Error("Unknown type of matrix storage format: %d", M->T);
696 break;
697 }
698 }
699
700
init_vector(int Nb,double ** V)701 void init_vector (int Nb, double **V){
702 *V = (double*) Malloc (Nb * sizeof(double));
703 }
704
705 /* ------------------------------------------------------------------------ */
706 /* f r e e */
707 /* ------------------------------------------------------------------------ */
708
free_matrix(Matrix * M)709 void free_matrix (Matrix *M){
710
711 if(M->scaled){
712 Free(M->rowscal) ;
713 Free(M->colscal) ;
714 }
715
716 switch (M->T) {
717 case SPARSE :
718 List_Delete(M->S.a);
719 List_Delete(M->S.ai);
720 List_Delete(M->S.ptr);
721 List_Delete(M->S.jptr);
722 break;
723 case DENSE :
724 Free(M->F.a);
725 Free(M->F.lu);
726 break;
727 }
728 }
729
730 /* ------------------------------------------------------------------------ */
731 /* z e r o */
732 /* ------------------------------------------------------------------------ */
733
zero_matrix(Matrix * M)734 void zero_matrix (Matrix *M){
735 int i,j=0;
736
737 M->changed = 1 ;
738
739 switch (M->T) {
740 case SPARSE :
741 List_Reset (M->S.a);
742 List_Reset (M->S.ai);
743 List_Reset (M->S.ptr);
744 List_Reset (M->S.jptr);
745 for (i=0; i<M->N; i++) List_Add (M->S.jptr, &j);
746 break;
747 case DENSE :
748 for(i=0; i<(M->N)*(M->N); i++) M->F.a[i] = 0.0;
749 break;
750 }
751 }
752
zero_matrix2(Matrix * M)753 void zero_matrix2 (Matrix *M){
754 int i, iptr;
755 int *jptr, *ptr;
756 double *a;
757
758 M->changed = 1 ;
759 switch (M->T) {
760 case SPARSE :
761 jptr = (int*) M->S.jptr->array;
762 ptr = (int*) M->S.ptr->array;
763 a = (double*) M->S.a->array;
764 for (i=0; i<M->N; i++) {
765 iptr = jptr[i];
766 while (iptr>0) {
767 a[iptr-1]= 0. ;
768 iptr = ptr[iptr-1];
769 }
770 }
771 break;
772 case DENSE :
773 for(i=0; i<(M->N)*(M->N); i++) M->F.a[i] = 0.0;
774 break;
775 }
776 }
777
zero_vector(int Nb,double * V)778 void zero_vector (int Nb, double *V){
779 int i;
780 for(i=0; i<Nb; i++) V[i] = 0.0;
781 }
782
783
784 /* ------------------------------------------------------------------------ */
785 /* c o p y */
786 /* ------------------------------------------------------------------------ */
787
copy_vector(int Nb,double * U,double * V)788 void copy_vector (int Nb, double *U, double *V){
789 int i;
790 for(i=0; i<Nb; i++) V[i] = U[i];
791 }
792
793 /* ------------------------------------------------------------------------ */
794 /* a d d */
795 /* ------------------------------------------------------------------------ */
796
add_vector_vector(int Nb,double * U,double * V)797 void add_vector_vector (int Nb, double *U, double *V){
798 /* u+v -> u */
799 int i;
800 for(i=0; i<Nb; i++) U[i] += V[i];
801 }
802
add_vector_prod_vector_double(int Nb,double * U,double * V,double d)803 void add_vector_prod_vector_double (int Nb, double *U, double *V, double d){
804 /* u+v*d -> u */
805 int i;
806 for(i=0; i<Nb; i++) U[i] += d*V[i];
807 }
808
add_matrix_double(Matrix * M,int ic,int il,double val)809 void add_matrix_double (Matrix *M, int ic, int il, double val){
810 /* attention a la transposition ! */
811 int *ai, *pp, n, iptr, iptr2, jptr, *ptr, zero = 0;
812 double *a;
813
814 M->changed = 1 ;
815
816 switch (M->T) {
817
818 case SPARSE :
819 il--;
820 pp = (int*) M->S.jptr->array;
821 ptr = (int*) M->S.ptr->array;
822 ai = (int*) M->S.ai->array;
823 a = (double*) M->S.a->array;
824
825 iptr = pp[il];
826 iptr2 = iptr-1;
827
828 while(iptr>0){
829 iptr2 = iptr-1;
830 jptr = ai[iptr2];
831 if(jptr == ic){
832 a[iptr2] += val;
833 return;
834 }
835 iptr = ptr[iptr2];
836 }
837
838 List_Add (M->S.a, &val);
839 List_Add (M->S.ai, &ic);
840 List_Add (M->S.ptr, &zero);
841
842 /* Les pointeurs ont pu etre modifies
843 s'il y a eu une reallocation dans List_Add */
844
845 ptr = (int*) M->S.ptr->array;
846 ai = (int*) M->S.ai->array;
847 a = (double*) M->S.a->array;
848
849 n = List_Nbr(M->S.a);
850 if(!pp[il]) pp[il] = n;
851 else ptr[iptr2] = n;
852 break;
853
854 case DENSE :
855 if(M->notranspose)
856 M->F.a[((M->N))*(il-1)+(ic-1)] += val;
857 else
858 M->F.a[((M->N))*(ic-1)+(il-1)] += val;
859 break;
860
861 }
862 }
863
864
add_matrix_matrix(Matrix * M,Matrix * N)865 void add_matrix_matrix (Matrix *M, Matrix *N){
866 /* M+N -> M */
867 int i, *ai, iptr, *jptr, *ptr;
868 double *a;
869
870 switch (M->T) {
871 case SPARSE :
872 jptr = (int*) N->S.jptr->array;
873 ptr = (int*) N->S.ptr->array;
874 a = (double*) N->S.a->array;
875 ai = (int*) N->S.ai->array;
876
877 for (i=0; i<N->N; i++) {
878 iptr = jptr[i];
879 while (iptr>0) {
880 add_matrix_double (M, ai[iptr-1], i+1, a[iptr-1]);
881 /* add_matrix_double transpose, donc pour additionner,
882 il faut transposer une seconde fois */
883 iptr = ptr[iptr-1];
884 }
885 }
886
887 break;
888 case DENSE :
889 for(i=0; i<(M->N)*(M->N); i++) M->F.a[i] += N->F.a[i];
890 break;
891
892 }
893 }
894
add_matrix_prod_matrix_double(Matrix * M,Matrix * N,double d)895 void add_matrix_prod_matrix_double (Matrix *M, Matrix *N, double d){
896 /* M+N*d -> M */
897 int i, *ai, iptr, *jptr, *ptr;
898 double *a;
899
900 switch (M->T) {
901 case SPARSE :
902 jptr = (int*) N->S.jptr->array;
903 ptr = (int*) N->S.ptr->array;
904 a = (double*) N->S.a->array;
905 ai = (int*) N->S.ai->array;
906
907 for (i=0; i<N->N; i++) {
908 iptr = jptr[i];
909 while (iptr>0) {
910 add_matrix_double (M, ai[iptr-1], i+1, d*a[iptr-1]);
911 /* add_matrix_double transpose, donc pour additionner,
912 il faut transposer une seconde fois */
913 iptr = ptr[iptr-1];
914 }
915 }
916
917 break;
918 case DENSE :
919 for(i=0; i<(M->N)*(M->N); i++) M->F.a[i] += d*N->F.a[i];
920 break;
921
922 }
923 }
924
925 /* ------------------------------------------------------------------------ */
926 /* s u b */
927 /* ------------------------------------------------------------------------ */
928
sub_vector_vector_1(int Nb,double * U,double * V)929 void sub_vector_vector_1 (int Nb, double *U, double *V){
930 /* u-v -> u */
931 int i;
932 for(i=0; i<Nb; i++) U[i] -= V[i];
933 }
934
935
sub_vector_vector_2(int Nb,double * U,double * V)936 void sub_vector_vector_2 (int Nb, double *U , double *V ){
937 /* u-v -> v */
938 int i;
939 for(i=0; i<Nb; i++) V[i] = U[i] - V[i];
940 }
941
942
943 /* ------------------------------------------------------------------------ */
944 /* p r o d */
945 /* ------------------------------------------------------------------------ */
946
prod_vector_double(int Nb,double * U,double a)947 void prod_vector_double (int Nb, double *U, double a){
948 /* u*a -> u */
949 int i;
950 for(i=0; i<Nb; i++) U[i] *= a;
951 }
952
prodsc_vector_vector(int Nb,double * U,double * V,double * prosca)953 void prodsc_vector_vector (int Nb, double *U, double *V, double *prosca){
954 /* u*v -> prosca */
955 int i;
956 *prosca = 0.0 ;
957 for (i=0; i<Nb; i++) *prosca += U[i] * V[i];
958 }
959
960
961
scale_matrix(int scaling,Matrix * M)962 void scale_matrix (int scaling, Matrix *M){
963 int i, *ai, *jptr ;
964 double *a, *rowscal = NULL, *colscal = NULL;
965 int job0=0, job1=1, ioff=0, len, *idiag, norm ;
966
967
968 switch (M->T) {
969
970 case SPARSE :
971
972 jptr = (int*) M->S.jptr->array;
973 a = (double*) M->S.a->array;
974 ai = (int*) M->S.ai->array;
975
976 switch (scaling) {
977
978 case DIAG_SCALING :
979
980 rowscal = colscal = (double*)Malloc(M->N * sizeof(double));
981
982 /* extract diagonal */
983 idiag = (int*)Malloc(M->N * sizeof(int));
984 getdia_ (&M->N, &M->N, &job0, a, ai, jptr, &len, rowscal, idiag, &ioff) ;
985 Free (idiag);
986
987 for (i = 0 ; i < M->N ; i++){
988 if (rowscal[i]){
989 rowscal[i] = 1./sqrt(fabs(rowscal[i])) ;
990 /* printf(" %d %e \n", i, rowscal[i] ); */
991 }
992 else {
993 Message::Warning("Diagonal scaling aborted because of zero diagonal element (%d)",i+1) ;
994 Free (rowscal) ;
995 return ;
996 }
997 }
998
999 diamua_ (&M->N, &job1, a, ai, jptr, rowscal, a, ai, jptr) ;
1000 amudia_ (&M->N, &job1, a, ai, jptr, colscal, a, ai, jptr) ;
1001 break ;
1002
1003 case MAX_SCALING : case NORM1_SCALING : case NORM2_SCALING :
1004
1005 switch (scaling) {
1006 case MAX_SCALING : norm = 0 ; break ;
1007 case NORM1_SCALING : norm = 1 ; break ;
1008 case NORM2_SCALING : norm = 2 ; break ;
1009 }
1010
1011 rowscal = (double*)Malloc(M->N * sizeof(double));
1012 rnrms_ (&M->N, &norm, a, ai, jptr, rowscal);
1013 for (i = 0 ; i < M->N ; i++){
1014 /* printf(" %d %e \n", i, rowscal[i] ); */
1015 if (rowscal[i])
1016 rowscal[i] = 1./rowscal[i] ;
1017 else {
1018 Message::Warning("Scaling aborted because of zero row (%d)", i+1) ;
1019 Free (rowscal) ;
1020 return ;
1021 }
1022 }
1023 diamua_ (&M->N, &job1, a, ai, jptr, rowscal, a, ai, jptr) ;
1024
1025 colscal = (double*)Malloc(M->N * sizeof(double));
1026 cnrms_ (&M->N, &norm, a, ai, jptr, colscal);
1027 for (i = 0 ; i < M->N ; i++){
1028 if (colscal[i]){
1029 colscal[i] = 1./colscal[i] ;
1030
1031 /* printf(" %d %e %e \n", i, 1./rowscal[i], 1./colscal[i] ); */
1032 }
1033 else {
1034 Message::Warning("Scaling aborted because of zero column (%d)", i+1) ;
1035 Free (colscal) ;
1036 return ;
1037 }
1038 }
1039 amudia_ (&M->N, &job1, a, ai, jptr, colscal, a, ai, jptr) ;
1040 break;
1041
1042 default :
1043
1044 Message::Error("Unknown type of matrix scaling: %d", scaling);
1045 break;
1046
1047 }
1048
1049 M->scaled = 1 ;
1050 M->rowscal = rowscal ;
1051 M->colscal = colscal ;
1052
1053 break;
1054
1055 case DENSE :
1056 Message::Warning("Scaling is not implemented for dense matrices") ;
1057 break;
1058 }
1059
1060 }
1061
1062
scale_vector(int ROW_or_COLUMN,Matrix * M,double * V)1063 void scale_vector (int ROW_or_COLUMN, Matrix *M, double *V){
1064 double *scal = NULL;
1065 int i;
1066
1067 if (!M->scaled) return ;
1068
1069 switch (ROW_or_COLUMN) {
1070 case 0 : scal = M->rowscal ; break ;
1071 case 1 : scal = M->colscal ; break ;
1072 }
1073
1074 if (scal == NULL) Message::Error("scale_vector : no scaling factors available !") ;
1075
1076 for (i = 0 ; i < M->N ; i++) V[i] *= scal[i] ;
1077 }
1078
1079
prod_matrix_vector(Matrix * M,double * V,double * res)1080 void prod_matrix_vector (Matrix *M, double *V , double *res ){
1081 /* M*V -> res ou M est la transposee!! */
1082 int k, i, j, *ai, *jptr ;
1083 double *a;
1084
1085 switch (M->T) {
1086 case SPARSE :
1087 /* csr_format transpose!
1088 donc la matrice arrivant dans cette routine doit
1089 bel et bien etre la transposee !!! */
1090 if(M->changed){
1091 csr_format (&M->S, M->N);
1092 restore_format (&M->S);
1093 M->changed = 0 ;
1094 }
1095 jptr = (int*) M->S.jptr->array;
1096 a = (double*) M->S.a->array;
1097 ai = (int*) M->S.ai->array;
1098
1099 for(i=0; i<M->N; i++){
1100 res[i] = 0.0 ;
1101 for(k=jptr[i]; k<=jptr[i+1]-1; k++){
1102 res[i] += V[ai[k-1]-1] * a[k-1];
1103 }
1104 }
1105
1106 break;
1107 case DENSE :
1108 if(M->notranspose){
1109 for(i=0; i<M->N; i++){
1110 res[i] = 0.0 ;
1111 for(j=0; j<M->N; j++) res[i] += M->F.a[(M->N)*i+j] * V[j];
1112 }
1113 }
1114 else{
1115 for(i=0; i<M->N; i++){
1116 res[i] = 0.0 ;
1117 for(j=0; j<M->N; j++) res[i] += M->F.a[(M->N)*j+i] * V[j];
1118 }
1119 }
1120 break;
1121 }
1122 }
1123
1124
prod_matrix_double(Matrix * M,double v)1125 void prod_matrix_double (Matrix *M, double v){
1126 /* M*v -> M */
1127 int i;
1128 double *a;
1129
1130 switch (M->T) {
1131 case SPARSE :
1132 a = (double*) M->S.a->array;
1133 for(i=0; i<List_Nbr(M->S.a); i++){
1134 a[i] *= v;
1135 }
1136 break;
1137 case DENSE :
1138 for(i=0; i<(M->N)*(M->N); i++) M->F.a[i] *= v;
1139 break;
1140 }
1141
1142 }
1143
multi_prod_matrix_double(int n,Matrix ** Mat,double * coef,Matrix * MatRes)1144 void multi_prod_matrix_double(int n, Matrix **Mat, double *coef, Matrix *MatRes){
1145 int k;
1146
1147 zero_matrix(MatRes);
1148 for(k=0;k<n;k++){
1149 if(coef[k]){
1150 prod_matrix_double(Mat[k],coef[k]);
1151 add_matrix_matrix(MatRes,Mat[k]);
1152 prod_matrix_double(Mat[k],1.0/coef[k]);
1153 }
1154 }
1155
1156 }
1157
multi_prod_vector_double(int n,int Sizevec,double ** Vec,double * coef,double * VecRes)1158 void multi_prod_vector_double(int n, int Sizevec, double **Vec, double *coef, double *VecRes ){
1159 int k;
1160
1161 zero_vector(Sizevec,VecRes);
1162 for(k=0;k<n;k++){
1163 if (coef[k]){
1164 prod_vector_double(Sizevec,Vec[k],coef[k]);
1165 add_vector_vector(Sizevec,VecRes,Vec[k]);
1166 prod_vector_double(Sizevec,Vec[k],1.0/coef[k]);
1167 }
1168 }
1169 }
1170
multi_prod_matrix_vector(int n,int Sizevec,Matrix ** Mat,double ** Vec,double * VecRes)1171 void multi_prod_matrix_vector(int n, int Sizevec, Matrix **Mat, double **Vec, double *VecRes ){
1172 int k;
1173 double *work;
1174
1175 init_vector(Sizevec,&work);
1176
1177 zero_vector(Sizevec,VecRes);
1178 for(k=0;k<n;k++){
1179 prod_matrix_vector(Mat[k],Vec[k],work);
1180 add_vector_vector(Sizevec,VecRes,work);
1181 }
1182 }
1183
prodsc_vectorconj_vector(int Nb,double * U,double * V,double * proscar,double * proscai)1184 void prodsc_vectorconj_vector (int Nb, double *U , double *V, double *proscar, double *proscai ){
1185 /* uconjugue * v -> proscar + i prodscai */
1186 int i;
1187 *proscar = *proscai = 0.0 ;
1188 for (i=0; i<Nb; i+=2){
1189 *proscar += U[i] * V[i] + U[i+1] * V[i+1] ;
1190 *proscai += U[i] * V[i+1] - U[i+1] * V[i] ;
1191 }
1192 }
1193
1194 /* ------------------------------------------------------------------------ */
1195 /* n o r m */
1196 /* ------------------------------------------------------------------------ */
1197
norm2_vector(int Nb,double * U,double * norm)1198 void norm2_vector(int Nb, double *U, double *norm){
1199 prodsc_vector_vector(Nb,U,U,norm);
1200 *norm = sqrt(*norm);
1201 }
1202
norminf_vector(int Nb,double * U,double * norm)1203 void norminf_vector(int Nb, double *U, double *norm){
1204 int i;
1205 *norm = 0.;
1206 for(i=0; i<Nb; i++)
1207 if(fabs(U[i]) > *norm) *norm = fabs(U[i]);
1208 }
1209
1210 /* ------------------------------------------------------------------------ */
1211 /* i d e n t i t y */
1212 /* ------------------------------------------------------------------------ */
1213
identity_matrix(Matrix * M)1214 void identity_matrix (Matrix *M){
1215 int i;
1216 zero_matrix(M);
1217 for (i=1;i<=M->N;i++) add_matrix_double(M,i,i,1.0);
1218 }
1219
1220
1221 /* ------------------------------------------------------------------------ */
1222 /* w r i t e */
1223 /* ------------------------------------------------------------------------ */
1224
binary_write_matrix(Matrix * M,const char * name,const char * ext)1225 void binary_write_matrix (Matrix *M, const char *name, const char *ext){
1226
1227 int Nb;
1228 FILE *pfile;
1229 char filename[256];
1230
1231 if(!M->N){
1232 Message::Warning("No elements in matrix");
1233 return;
1234 }
1235
1236 strcpy(filename, name);
1237 strcat(filename, ext);
1238 pfile = fopen(filename, "wb") ;
1239
1240 fprintf(pfile,"%d\n",M->T);
1241
1242 switch (M->T) {
1243 case SPARSE :
1244 Nb = List_Nbr(M->S.a) ;
1245
1246 fprintf(pfile,"%d %d\n", M->N, Nb);
1247
1248 fprintf(pfile,"%d %d %d %d %d\n", M->S.ptr->nmax, M->S.ptr->size,
1249 M->S.ptr->incr, M->S.ptr->n, M->S.ptr->isorder);
1250 fprintf(pfile,"%d %d %d %d %d\n", M->S.ai->nmax, M->S.ai->size,
1251 M->S.ai->incr, M->S.ai->n, M->S.ai->isorder);
1252 fprintf(pfile,"%d %d %d %d %d\n", M->S.jptr->nmax, M->S.jptr->size,
1253 M->S.jptr->incr, M->S.jptr->n, M->S.jptr->isorder);
1254 fprintf(pfile,"%d %d %d %d %d\n", M->S.a->nmax, M->S.a->size,
1255 M->S.a->incr, M->S.a->n, M->S.a->isorder);
1256
1257 fwrite(M->S.ptr->array, sizeof(int), Nb, pfile);
1258 fwrite(M->S.ai->array, sizeof(int), Nb, pfile);
1259 fwrite(M->S.jptr->array, sizeof(int), M->N, pfile);
1260 fwrite(M->S.a->array, sizeof(double), Nb, pfile);
1261 break;
1262
1263 case DENSE :
1264 fprintf(pfile,"%d\n", M->N);
1265 fwrite(M->F.a, sizeof(double), M->N*M->N, pfile);
1266 break;
1267 }
1268
1269 fclose(pfile) ;
1270 }
1271
1272
binary_write_vector(int Nb,double * V,const char * name,const char * ext)1273 void binary_write_vector (int Nb, double *V, const char *name, const char *ext){
1274 char filename[256];
1275 FILE *pfile;
1276
1277 strcpy(filename, name);
1278 strcat(filename, ext);
1279 pfile = fopen(filename, "wb") ;
1280
1281 fwrite(V, sizeof(double), Nb, pfile);
1282
1283 fclose(pfile) ;
1284 }
1285
1286
formatted_write_matrix(FILE * pfile,Matrix * M,int style)1287 void formatted_write_matrix (FILE *pfile, Matrix *M, int style){
1288 int *ptr,*ai,i,j,*jptr, *ia, *ja, *ir, nnz, ierr;
1289 int un=1;
1290 double *a;
1291
1292 if(!M->N){
1293 Message::Warning("No element in matrix");
1294 return;
1295 }
1296
1297 switch (M->T) {
1298
1299 case DENSE :
1300 if(M->notranspose)
1301 for(i=0 ; i<M->N ; i++)
1302 for(j=0 ; j<M->N ; j++)
1303 fprintf(pfile,"%d %d %.16g\n", j+1, i+1, M->F.a[i*(M->N)+j]);
1304 else
1305 for(i=0 ; i<M->N ; i++)
1306 for(j=0 ; j<M->N ; j++)
1307 fprintf(pfile,"%d %d %.16g\n", i+1, j+1, M->F.a[i*(M->N)+j]);
1308 break;
1309
1310 case SPARSE :
1311
1312 switch(style){
1313
1314 case ELAP :
1315 fprintf(pfile,"%d\n",M->T);
1316 a = (double*)M->S.a->array;
1317 ai = (int*)M->S.ai->array;
1318 ptr = (int*)M->S.ptr->array;
1319 jptr = (int*)M->S.jptr->array;
1320 fprintf(pfile,"%d\n",M->N);
1321 fprintf(pfile,"%d\n",List_Nbr(M->S.a));
1322 for(i=0;i<M->N;i++)
1323 fprintf(pfile," %d",jptr[i]);
1324 fprintf(pfile,"\n");
1325 for(i=0;i<List_Nbr(M->S.a);i++)
1326 fprintf(pfile,"%d %d %.16g \n",ai[i],ptr[i],a[i]);
1327 break;
1328
1329 case KUL :
1330 csr_format(&M->S, M->N);
1331 a = (double*) M->S.a->array;
1332 ia = (int*) M->S.jptr->array;
1333 ja = (int*) M->S.ptr->array;
1334 nnz = List_Nbr(M->S.a);
1335 ir = (int*) Malloc(nnz * sizeof(int));
1336 csrcoo_(&M->N, &un, &nnz, a, ja, ia, &nnz, a, ir, ja, &ierr);
1337 for(i=0 ; i<nnz ; i++)
1338 fprintf(pfile,"%d %d %.16g\n", ir[i], ja[i], a[i]);
1339 restore_format(&M->S);
1340 break;
1341
1342 default :
1343 Message::Error("Unknown print style for formatted matrix output");
1344 }
1345 break ;
1346
1347 default :
1348 Message::Error("Unknown matrix format for formatted matrix output");
1349
1350 }
1351 }
1352
1353
formatted_write_vector(FILE * pfile,int Nb,double * V,int style)1354 void formatted_write_vector (FILE *pfile, int Nb, double *V, int style){
1355 int i;
1356
1357 /* for(i=0 ; i<Nb ; i++) fprintf(pfile,"%d %.16g\n", i+1, V[i]); */
1358 for(i=0 ; i<Nb ; i++) fprintf(pfile,"%.16g\n", V[i]);
1359 }
1360
1361
1362 /* ------------------------------------------------------------------------ */
1363 /* r e a d */
1364 /* ------------------------------------------------------------------------ */
1365
1366
binary_read_matrix(Matrix * M,const char * name,const char * ext)1367 void binary_read_matrix (Matrix *M, const char *name , const char *ext){
1368
1369 int Nb;
1370 FILE *pfile;
1371 char filename[256];
1372
1373 strcpy(filename, name);
1374 strcat(filename, ext);
1375 pfile = fopen(filename, "rb") ;
1376
1377 if (pfile == NULL) {
1378 Message::Error("Error opening file '%s'", filename);
1379 }
1380
1381 fscanf(pfile,"%d",&M->T);
1382 M->ILU_Exists = 0;
1383
1384 switch (M->T) {
1385 case SPARSE :
1386 fscanf(pfile,"%d %d\n", &M->N, &Nb);
1387
1388 M->S.ptr = List_Create (Nb, 1, sizeof(int));
1389 M->S.ai = List_Create (Nb, 1, sizeof(int));
1390 M->S.jptr = List_Create (M->N, 1, sizeof(int));
1391 M->S.a = List_Create (Nb, 1, sizeof(double));
1392
1393 fscanf(pfile,"%d %d %d %d %d\n", &M->S.ptr->nmax, &M->S.ptr->size,
1394 &M->S.ptr->incr, &M->S.ptr->n, &M->S.ptr->isorder);
1395 fscanf(pfile,"%d %d %d %d %d\n", &M->S.ai->nmax, &M->S.ai->size,
1396 &M->S.ai->incr, &M->S.ai->n, &M->S.ai->isorder);
1397 fscanf(pfile,"%d %d %d %d %d\n", &M->S.jptr->nmax, &M->S.jptr->size,
1398 &M->S.jptr->incr, &M->S.jptr->n, &M->S.jptr->isorder);
1399 fscanf(pfile,"%d %d %d %d %d\n", &M->S.a->nmax, &M->S.a->size,
1400 &M->S.a->incr, &M->S.a->n, &M->S.a->isorder);
1401
1402 fread(M->S.ptr->array, sizeof(int), Nb, pfile);
1403 fread(M->S.ai->array, sizeof(int), Nb, pfile);
1404 fread(M->S.jptr->array, sizeof(int), M->N, pfile);
1405 fread(M->S.a->array, sizeof(double), Nb, pfile);
1406 break;
1407
1408 case DENSE :
1409 fscanf(pfile,"%d\n", &M->N);
1410 M->F.LU_Exist = 0;
1411 M->F.a = (double*) Malloc(M->N * M->N * sizeof(double));
1412 M->F.lu = (double*) Malloc(M->N * M->N * sizeof(double));
1413 fread(M->F.a, sizeof(double), M->N * M->N, pfile);
1414 break ;
1415 }
1416
1417 fclose(pfile) ;
1418
1419 }
1420
1421
binary_read_vector(int Nb,double ** V,const char * name,const char * ext)1422 void binary_read_vector (int Nb, double **V, const char *name, const char *ext){
1423 char filename[256];
1424 FILE *pfile;
1425
1426 strcpy(filename, name);
1427 strcat(filename, ext);
1428 pfile = fopen(filename, "rb") ;
1429
1430 if (pfile == NULL) {
1431 Message::Error("Error opening file %s", filename);
1432 }
1433
1434 init_vector(Nb, V);
1435 fread(*V, sizeof(double), Nb, pfile);
1436
1437 fclose(pfile) ;
1438 }
1439
1440
formatted_read_matrix(Matrix * M,const char * name,const char * ext,int style)1441 void formatted_read_matrix (Matrix *M, const char *name , const char *ext, int style){
1442 int i,nnz,inb,inb2;
1443 double nb;
1444 FILE *pfile;
1445 char filename[256];
1446
1447 strcpy(filename, name);
1448 strcat(filename, ext);
1449 pfile = fopen(filename, "r") ;
1450
1451 if (pfile == NULL) {
1452 Message::Error("Error opening file %s", filename);
1453 }
1454
1455 fscanf(pfile,"%d",&M->T);
1456 switch (M->T) {
1457 case SPARSE :
1458 List_Reset(M->S.jptr);
1459 fscanf(pfile,"%d",&M->N);
1460 fscanf(pfile,"%d",&nnz);
1461 for(i=0;i<M->N;i++){
1462 fscanf(pfile," %d",&inb);
1463 List_Add(M->S.jptr,&inb);
1464 }
1465 for(i=0;i<nnz;i++){
1466 fscanf(pfile,"%d %d %lf \n",&inb,&inb2,&nb);
1467 List_Add(M->S.ai,&inb);
1468 List_Add(M->S.ptr,&inb2);
1469 List_Add(M->S.a,&nb);
1470 }
1471
1472 break;
1473
1474 case DENSE :
1475 fscanf(pfile,"%d",&M->N);
1476 for(i=0;i<(M->N)*(M->N);i++){
1477 fscanf(pfile,"%d %lf ", &inb, &M->F.a[i]);
1478 }
1479 break;
1480 }
1481 fclose(pfile) ;
1482
1483 }
1484
1485
formatted_read_vector(int Nb,double * V,const char * name,const char * ext,int style)1486 void formatted_read_vector (int Nb, double *V, const char *name, const char *ext, int style){
1487 int i;
1488 FILE *pfile;
1489 char filename[256];
1490
1491 strcpy(filename, name);
1492 strcat(filename, ext);
1493 pfile = fopen(filename, "r") ;
1494
1495 if (pfile == NULL) {
1496 Message::Error("Error opening file %s", filename);
1497 }
1498
1499 for(i=0 ; i<Nb ; i++) fscanf(pfile, "%lf", &V[i]);
1500
1501 fclose(pfile) ;
1502 }
1503
1504 /* ------------------------------------------------------------------------ */
1505 /* p r i n t _ m a t r i x _ i n f o _ X X X */
1506 /* ------------------------------------------------------------------------ */
1507
maximum(int a,int b)1508 int maximum (int a, int b) {
1509 if (a>b)
1510 return(a);
1511 else
1512 return(b);
1513 }
1514
1515
print_matrix_info_CSR(int N,int * jptr,int * ai)1516 void print_matrix_info_CSR (int N, int *jptr, int *ai){
1517 int i, j, k, l, m, n;
1518
1519 l = n = 0;
1520 j = jptr[N]-1 ;
1521 for (i=0; i<N; i++) {
1522 k = jptr[i+1] - jptr[i];
1523 m = ai[jptr[i+1]-2] - ai[jptr[i]-1] + 1;
1524 if (l<k) l = k;
1525 if (n<m) n = m;
1526 }
1527
1528 Message::Info("N: %d, NZ: %d, BW max/avg: %d/%d, SW max: %d",
1529 N, j, l, (int)(j/N), n);
1530 }
1531
print_matrix_info_MSR(int N,sscalar * a,int * jptr)1532 void print_matrix_info_MSR (int N, sscalar *a, int *jptr){
1533 int i, j, k, l, m, n;
1534
1535 l = n = 0;
1536 j = jptr[N]-2;
1537 for (i=0; i<N; i++) {
1538 k = jptr[i+1] - jptr[i] + (a[i]?1:0) ;
1539 if((jptr[i+1] - jptr[i]) == 0)
1540 m = (a[i]?1:0);
1541 else
1542 m = maximum ( jptr[jptr[i+1]-2] - jptr[jptr[i]-1] + 1,
1543 maximum (jptr[jptr[i+1]-2] - (i+1) + 1, (i+1) - jptr[jptr[i]-1] + 1) );
1544 if (l<k) l = k;
1545 if (n<m) n = m;
1546 }
1547
1548 Message::Info("N: %d, NZ: %d, BW max/avg: %d/%d, SW max: %d",
1549 N, j, l, (int)(j/N), n);
1550
1551 }
1552
print_matrix_info_DENSE(int N)1553 void print_matrix_info_DENSE (int N){
1554 Message::Info("N: %d", N);
1555 }
1556
1557
1558 /* ------------------------------------------------------------------------ */
1559 /* get _ column _ in _ m a t r i x */
1560 /* ------------------------------------------------------------------------ */
1561
get_column_in_matrix(Matrix * M,int col,double * V)1562 void get_column_in_matrix (Matrix *M, int col, double *V){
1563
1564 int k, i, j, *ai, *jptr ;
1565 double *a;
1566 int found;
1567
1568 switch (M->T) {
1569 case SPARSE :
1570 /* csr_format transpose!
1571 donc la matrice arrivant dans cette routine doit
1572 bel et bien etre la transposee !!! */
1573 if(M->changed){
1574 csr_format (&M->S, M->N);
1575 restore_format (&M->S);
1576 M->changed = 0 ;
1577 }
1578 jptr = (int*) M->S.jptr->array;
1579 a = (double*) M->S.a->array;
1580 ai = (int*) M->S.ai->array;
1581
1582 for(i=0; i<M->N; i++){ /* lignes */
1583 found=0;
1584 for(k=jptr[i]-1;k<jptr[i+1]-1;k++){ /*colonne */
1585 if(ai[k]-1==col) {
1586 V[i]=a[k]; found=1; break;
1587 }
1588 else if (ai[k]-1 > col) {
1589 break;
1590 }
1591 }
1592 if (!found) V[i]=0;
1593 /* printf(" V[%d] = %g \n",i, V[i]); */
1594 }
1595 break;
1596 case DENSE :
1597 if(M->notranspose){
1598 for(j=0; j<M->N; j++) V[j] = M->F.a[(M->N)*col+j];
1599 }
1600 else{
1601 for(i=0; i<M->N; i++){
1602 for(j=0; j<M->N; j++) V[j] = M->F.a[(M->N)*j+col];
1603 }
1604 }
1605 break;
1606 }
1607 }
1608
1609
1610
get_element_in_matrix(Matrix * M,int row,int col,double * V)1611 void get_element_in_matrix (Matrix *M, int row, int col, double *V){
1612
1613 int k, i, *ai, *jptr ;
1614 double *a;
1615 int found;
1616
1617 switch (M->T) {
1618 case SPARSE :
1619 /* csr_format transpose!
1620 donc la matrice arrivant dans cette routine doit
1621 bel et bien etre la transposee !!! */
1622 if(M->changed){
1623 csr_format (&M->S, M->N);
1624 restore_format (&M->S);
1625 M->changed = 0 ;
1626 }
1627 jptr = (int*) M->S.jptr->array;
1628 a = (double*) M->S.a->array;
1629 ai = (int*) M->S.ai->array;
1630
1631 for(i=0; i<M->N; i++){ /* lignes */
1632 found=0;
1633 for(k=jptr[i]-1;k<jptr[i+1]-1;k++){ /*colonne */
1634 if(ai[k]-1==col) {
1635 V[i]=a[k]; found=1; break;
1636 }
1637 else if (ai[k]-1 > col) {
1638 break;
1639 }
1640 }
1641 if (!found) V[i]=0;
1642 /* printf(" V[%d] = %g \n",i, V[i]); */
1643 }
1644 break;
1645 case DENSE :
1646 if(M->notranspose){
1647 *V = M->F.a[(M->N)*col+row];
1648 }
1649 else{
1650 for(i=0; i<M->N; i++){
1651 *V = M->F.a[(M->N)*row+col];
1652 }
1653 }
1654 break;
1655 }
1656 }
1657
1658
1659 /* ------------------------------------------------------------------------ */
1660 /* S o l v e r p a r a m e t e r s */
1661 /* ------------------------------------------------------------------------ */
1662
1663 static char comALGORITHM[] =
1664 "\n%s (Integer): \n\
1665 - 1 CG Conjugate Gradient \n\
1666 - 2 CGNR CG (Normal Residual equation) \n\
1667 - 3 BCG Bi-Conjugate Gradient \n\
1668 - 4 DBCG BCG with partial pivoting \n\
1669 - 5 BCGSTAB BCG stabilized \n\
1670 - 6 TFQMR Transpose-Free Quasi-Minimum Residual \n\
1671 - 7 FOM Full Orthogonalization Method \n\
1672 - 8 GMRES Generalized Minimum RESidual \n\
1673 - 9 FGMRES Flexible version of GMRES \n\
1674 - 10 DQGMRES Direct versions of GMRES \n\
1675 - 11 LU LU Factorization \n\
1676 - 12 PGMRES Alternative version of GMRES \n\
1677 - default : %d\n";
1678
1679 static char comPRECONDITIONER[] =
1680 "\n%s (Integer): \n\
1681 - 0 NONE No Factorization\n\
1682 - 1 ILUT Incomplete LU factorization with dual truncation strategy \n\
1683 - 2 ILUTP ILUT with column pivoting \n\
1684 - 3 ILUD ILU with single dropping + diagonal compensation (~MILUT) \n\
1685 - 4 ILUDP ILUD with column pivoting \n\
1686 - 5 ILUK level-k ILU \n\
1687 - 6 ILU0 simple ILU(0) preconditioning \n\
1688 - 7 MILU0 MILU(0) preconditioning \n\
1689 - 8 DIAGONAL \n\
1690 - default : %d \n";
1691
1692 static char comPRECONDITIONER_POSITION[] =
1693 "\n%s (Integer): \n\
1694 - 0 No Preconditioner \n\
1695 - 1 Left Preconditioner \n\
1696 - 2 Right Preconditioner \n\
1697 - 3 Both Left and Right Preconditioner \n\
1698 - default : %d \n";
1699
1700 static char comRENUMBERING_TECHNIQUE[] =
1701 "\n%s (Integer): \n\
1702 - 0 No renumbering \n\
1703 - 1 Reverse Cuthill-Mc Kee \n\
1704 - default : %d \n";
1705
1706 static char comNB_ITER_MAX[] =
1707 "\n%s (Integer): Maximum number of iterations \n\
1708 - default : %d \n";
1709
1710 static char comMATRIX_FORMAT[] =
1711 "\n%s (Integer): \n\
1712 - 1 Sparse \n\
1713 - 2 Full \n\
1714 - default : %d\n";
1715
1716 static char comMATRIX_PRINTING[] =
1717 "\n%s (Integer): Disk write ('fort.*') \n\
1718 - 1 matrix (csr) \n\
1719 - 2 preconditioner (msr) \n\
1720 - 3 both \n\
1721 - default : %d\n";
1722
1723 static char comMATRIX_STORAGE[] =
1724 "\n%s (Integer): Disk Write or Read in internal format \n\
1725 - 0 none \n\
1726 - 1 write matrix (sparse) \n\
1727 - 2 read matrix (sparse) \n\
1728 - default : %d\n";
1729
1730 static char comNB_FILL[] =
1731 "\n%s (Integer): \n\
1732 - ILUT/ILUTP : maximum number of elements per line \n\
1733 of L and U (except diagonal element) \n\
1734 - ILUK : each element whose fill-in level is greater than NB_FILL \n\
1735 is dropped. \n\
1736 - default : %d\n";
1737
1738 static char comKRYLOV_SIZE[] =
1739 "\n%s (Integer): Krylov subspace size \n\
1740 - default : %d\n";
1741
1742 static char comSTOPPING_TEST[] =
1743 "\n%s (Real): Target relative residual \n\
1744 - default : %g \n";
1745
1746 static char comIC_ACCELERATION[] =
1747 "\n%s (Real): IC accelerator\n\
1748 - default : %g \n";
1749
1750 static char comITERATIVE_IMPROVEMENT[] =
1751 "\n%s (Integer): Iterative improvement of the solution obtained by a LU \n\
1752 - default : %d\n";
1753
1754 static char comDROPPING_TOLERANCE[] =
1755 "\n%s (Real): \n\
1756 - ILUT/ILUTP/ILUK: a(i,j) is dropped if \n\
1757 abs(a(i,j)) < DROPPING_TOLERANCE * abs(diagonal element in U). \n\
1758 - ILUD/ILUDP : a(i,j) is dropped if \n\
1759 abs(a(i,j)) < DROPPING_TOLERANCE * [weighted norm of line i]. \n\
1760 Weighted norm = 1-norm / number of nonzero elements on the line. \n\
1761 - default : %g\n";
1762
1763 static char comPERMUTATION_TOLERANCE[] =
1764 "\n%s (Real): Tolerance for column permutation in ILUTP/ILUDP. \n\
1765 At stage i, columns i and j are permuted if \n\
1766 abs(a(i,j))*PERMUTATION_TOLERANCE > abs(a(i,i)). \n\
1767 - 0 no permutations \n\
1768 - 0.001 -> 0.1 classical \n\
1769 - default : %g\n";
1770
1771 static char comRE_USE_LU[] =
1772 "\n%s (Integer): Reuse LU decomposition\n\
1773 - 0 no \n\
1774 - 1 yes \n\
1775 - default : %d\n";
1776
1777 static char comRE_USE_ILU[] =
1778 "\n%s (Integer): Reuse ILU decomposition (and renumbering if any)\n\
1779 - 0 no \n\
1780 - 1 yes \n\
1781 - default : %d\n";
1782
1783 static char comDIAGONAL_COMPENSATION[] =
1784 "\n%s (Real): ILUD/ILUDP: the term 'DIAGONAL_COMPENSATION * (sum \n\
1785 of all dropped elements of the line)' is added to the diagonal element in U \n\
1786 - 0 ~ ILU with threshold \n\
1787 1 ~ MILU with threshold. \n\
1788 - default : %g\n";
1789
1790 static char comSCALING[] =
1791 "\n%s (Integer): Scale system \n\
1792 - 0 no \n\
1793 - 1 on basis of diagonal elements (no loss of possible symmetry) \n\
1794 - 2 on basis of inf. norm of first rows and then columns (asymmetric) \n\
1795 - 3 on basis of norm 1 of first rows and then columns (asymmetric) \n\
1796 - 4 on basis of norm 2 of first rows and then columns (asymmetric) \n\
1797 - default : %d\n";
1798
1799 /* ------------------------------------------------------------------------ */
1800 /* A c t i o n s */
1801 /* ------------------------------------------------------------------------ */
1802
1803 #define act_ARGS Solver_Params *p, int i, double d
1804
actALGORITHM(act_ARGS)1805 void actALGORITHM (act_ARGS){ p->Algorithm = i; }
actPRECONDITIONER(act_ARGS)1806 void actPRECONDITIONER (act_ARGS){ p->Preconditioner = i; }
actPRECONDITIONER_POSITION(act_ARGS)1807 void actPRECONDITIONER_POSITION (act_ARGS){ p->Preconditioner_Position = i; }
actRENUMBERING_TECHNIQUE(act_ARGS)1808 void actRENUMBERING_TECHNIQUE (act_ARGS){ p->Renumbering_Technique = i; }
actNB_ITER_MAX(act_ARGS)1809 void actNB_ITER_MAX (act_ARGS){ p->Nb_Iter_Max = i; }
actMATRIX_FORMAT(act_ARGS)1810 void actMATRIX_FORMAT (act_ARGS){ p->Matrix_Format = i; }
actMATRIX_PRINTING(act_ARGS)1811 void actMATRIX_PRINTING (act_ARGS){ p->Matrix_Printing = i; }
actMATRIX_STORAGE(act_ARGS)1812 void actMATRIX_STORAGE (act_ARGS){ p->Matrix_Storage = i; }
actNB_FILL(act_ARGS)1813 void actNB_FILL (act_ARGS){ p->Nb_Fill = i; }
actKRYLOV_SIZE(act_ARGS)1814 void actKRYLOV_SIZE (act_ARGS){ p->Krylov_Size = i; }
actSTOPPING_TEST(act_ARGS)1815 void actSTOPPING_TEST (act_ARGS){ p->Stopping_Test = d; }
actIC_ACCELERATION(act_ARGS)1816 void actIC_ACCELERATION (act_ARGS){ p->IC_Acceleration = d; }
actITERATIVE_IMPROVEMENT(act_ARGS)1817 void actITERATIVE_IMPROVEMENT (act_ARGS){ p->Iterative_Improvement = i; }
actRE_USE_LU(act_ARGS)1818 void actRE_USE_LU (act_ARGS){ p->Re_Use_LU = i; }
actDROPPING_TOLERANCE(act_ARGS)1819 void actDROPPING_TOLERANCE (act_ARGS){ p->Dropping_Tolerance = d; }
actPERMUTATION_TOLERANCE(act_ARGS)1820 void actPERMUTATION_TOLERANCE (act_ARGS){ p->Permutation_Tolerance = d; }
actRE_USE_ILU(act_ARGS)1821 void actRE_USE_ILU (act_ARGS){ p->Re_Use_ILU = i; }
actDIAGONAL_COMPENSATION(act_ARGS)1822 void actDIAGONAL_COMPENSATION (act_ARGS){ p->Diagonal_Compensation = d; }
actSCALING(act_ARGS)1823 void actSCALING (act_ARGS){ p->Scaling = i; }
1824
1825
1826
1827 /* ------------------------------------------------------------------------ */
1828 /* P a r a m e t e r s w i t h d e f a u l t v a l u e s */
1829 /* ------------------------------------------------------------------------ */
1830
1831 #define REEL 1
1832 #define ENTIER 2
1833
1834 typedef struct {
1835 const char *str;
1836 int typeinfo;
1837 int defaultint;
1838 double defaultfloat;
1839 const char *com;
1840 void (*action) (Solver_Params *p , int i , double d);
1841 }InfoSolver;
1842
compInfoSolver(const void * a,const void * b)1843 int compInfoSolver(const void *a, const void *b){
1844 return(strcmp(((InfoSolver*)a)->str, ((InfoSolver*)b)->str));
1845 }
1846
1847 static InfoSolver Tab_Params[] =
1848 {
1849 {"Matrix_Format", ENTIER, 1, 0., comMATRIX_FORMAT, actMATRIX_FORMAT},
1850 {"Matrix_Printing", ENTIER, 0, 0., comMATRIX_PRINTING, actMATRIX_PRINTING},
1851 {"Matrix_Storage", ENTIER, 0, 0., comMATRIX_STORAGE, actMATRIX_STORAGE},
1852 {"Scaling", ENTIER, 0, 0., comSCALING, actSCALING},
1853 {"Renumbering_Technique", ENTIER, 1, 0., comRENUMBERING_TECHNIQUE, actRENUMBERING_TECHNIQUE},
1854 {"Preconditioner", ENTIER, 2, 0., comPRECONDITIONER, actPRECONDITIONER},
1855 {"Preconditioner_Position", ENTIER, 2, 0., comPRECONDITIONER_POSITION, actPRECONDITIONER_POSITION},
1856 {"Nb_Fill", ENTIER, 20, 0., comNB_FILL, actNB_FILL},
1857 {"Permutation_Tolerance", REEL, 0, 5.e-2, comPERMUTATION_TOLERANCE, actPERMUTATION_TOLERANCE},
1858 {"Dropping_Tolerance", REEL, 0, 0., comDROPPING_TOLERANCE, actDROPPING_TOLERANCE},
1859 {"Diagonal_Compensation", REEL, 0, 0., comDIAGONAL_COMPENSATION, actDIAGONAL_COMPENSATION},
1860 {"Re_Use_ILU", ENTIER, 0, 0., comRE_USE_ILU, actRE_USE_ILU},
1861 {"Algorithm", ENTIER, 8, 0., comALGORITHM, actALGORITHM},
1862 {"Krylov_Size", ENTIER, 40, 0., comKRYLOV_SIZE, actKRYLOV_SIZE},
1863 {"IC_Acceleration", REEL, 0, 1., comIC_ACCELERATION, actIC_ACCELERATION},
1864 {"Re_Use_LU", ENTIER, 0, 0., comRE_USE_LU, actRE_USE_LU},
1865 {"Iterative_Improvement", ENTIER, 0, 0., comITERATIVE_IMPROVEMENT, actITERATIVE_IMPROVEMENT},
1866 {"Nb_Iter_Max", ENTIER, 1000, 0., comNB_ITER_MAX, actNB_ITER_MAX},
1867 {"Stopping_Test", REEL, 0, 1.e-10,comSTOPPING_TEST, actSTOPPING_TEST}
1868 };
1869
1870
1871 /* ------------------------------------------------------------------------ */
1872 /* i n i t _ s o l v e r */
1873 /* ------------------------------------------------------------------------ */
1874
1875 #define NbInfosSolver (int)(sizeof(Tab_Params)/sizeof(Tab_Params[0]))
1876
Commentaires(FILE * out)1877 void Commentaires (FILE *out){
1878 int i;
1879 InfoSolver *pI;
1880
1881 for(i=0;i<NbInfosSolver;i++){
1882 pI = &Tab_Params[i];
1883 switch(pI->typeinfo){
1884 case REEL :
1885 fprintf(out,pI->com,pI->str,pI->defaultfloat);
1886 break;
1887 case ENTIER :
1888 fprintf(out,pI->com,pI->str,pI->defaultint);
1889 break;
1890 }
1891 }
1892 fprintf(out,"\n");
1893 }
1894
init_solver(Solver_Params * p,const char * name)1895 void init_solver (Solver_Params *p , const char *name){
1896 char buff[128];
1897 FILE *file;
1898 InfoSolver *pI,I;
1899 int i;
1900 double ff;
1901 int ii;
1902
1903 for(i=0;i<NbInfosSolver;i++){
1904 pI = &Tab_Params[i];
1905 (pI->action)(p,pI->defaultint,pI->defaultfloat);
1906 }
1907
1908 if(!(file = fopen(name,"r"))){
1909 file = fopen(name,"w");
1910 if(!file){
1911 Message::Warning("Could not open solver parameter file");
1912 return;
1913 }
1914 fprintf(file,"/*\n");
1915 Commentaires(file);
1916 fprintf(file,"*/\n\n");
1917 Message::Info("Parameter file not found");
1918 Message::Info("Enter parameter values:");
1919 for(i=0;i<NbInfosSolver;i++){
1920 bool error = false;
1921 pI = &Tab_Params[i];
1922 switch(pI->typeinfo){
1923 case REEL :
1924 getfloat :
1925 if(Message::UseSocket() || Message::UseOnelab())
1926 strcpy(buff, "\n");
1927 else{
1928 printf("%25s (Real) [<h>=help, <return>=%g]: ",pI->str,pI->defaultfloat);
1929 if(!fgets(buff, 128, stdin))
1930 error = true;
1931 }
1932 if(!error && !strcmp(buff,"h\n")){
1933 printf(pI->com,pI->str,pI->defaultfloat);
1934 printf("\n");
1935 goto getfloat;
1936 }
1937 if(error || !strcmp(buff,"\n")){
1938 fprintf(file,"%25s %12g\n",pI->str,pI->defaultfloat);
1939 (pI->action)(p,pI->defaultint,pI->defaultfloat);
1940 }
1941 else{
1942 fprintf(file,"%25s %12g\n",pI->str,atof(buff));
1943 (pI->action)(p,pI->defaultint,atof(buff));
1944 }
1945 break;
1946 case ENTIER :
1947 getint :
1948 if(Message::UseSocket() || Message::UseOnelab()){
1949 strcpy(buff, "\n");
1950 }
1951 else{
1952 printf("%25s (Integer) [<h>=help, <return>=%d]: ",pI->str,pI->defaultint);
1953 if(!fgets(buff, 128, stdin))
1954 error = true;
1955 }
1956 if(!error && !strcmp(buff,"h\n")){
1957 printf(pI->com,pI->str,pI->defaultint);
1958 printf("\n");
1959 goto getint;
1960 }
1961 if(error || !strcmp(buff,"\n")){
1962 fprintf(file,"%25s %12d\n",pI->str,pI->defaultint);
1963 (pI->action)(p,pI->defaultint,pI->defaultfloat);
1964 }
1965 else{
1966 fprintf(file,"%25s %12d\n",pI->str,atoi(buff));
1967 (pI->action)(p,atoi(buff),pI->defaultfloat);
1968 }
1969 break;
1970 }
1971 }
1972 }
1973 else {
1974 qsort(Tab_Params, NbInfosSolver, sizeof(InfoSolver), compInfoSolver);
1975 rewind(file);
1976 while (!feof(file)){
1977 fscanf(file,"%s",buff);
1978 I.str = buff;
1979 if(!(pI = (InfoSolver*)bsearch(&I,Tab_Params, NbInfosSolver,
1980 sizeof(InfoSolver),compInfoSolver))){
1981 if(buff[0] == '/' && buff[1] == '*'){
1982 while(1){
1983 if(feof(file)){
1984 Message::Warning("End of comment not detected");
1985 fclose(file);
1986 return;
1987 }
1988 if((getc(file)=='*')&&(getc(file)=='/')){
1989 break;
1990 }
1991 }
1992 }
1993 else{
1994 Message::Warning("Unknown solver parameter '%s'", buff);
1995 fscanf(file,"%s",buff);
1996 }
1997 }
1998 else{
1999 switch(pI->typeinfo){
2000 case REEL :
2001 fscanf(file,"%lf",&ff);
2002 (pI->action)(p,ii,ff);
2003 break;
2004 case ENTIER :
2005 fscanf(file,"%d",&ii);
2006 (pI->action)(p,ii,ff);
2007 break;
2008 }
2009 }
2010 }
2011 }
2012 fclose(file);
2013 }
2014
2015
init_solver_option(Solver_Params * p,const char * name,const char * value)2016 void init_solver_option (Solver_Params *p , const char *name, const char *value){
2017 InfoSolver *pI;
2018 int i, vali;
2019 float valf;
2020
2021 for(i=0;i<NbInfosSolver;i++){
2022 pI = &Tab_Params[i];
2023
2024 if(!strcmp(pI->str, name)){
2025 switch(pI->typeinfo){
2026 case REEL :
2027 valf = atof(value);
2028 (pI->action)(p,pI->defaultint,valf);
2029 Message::Info("Overriding parameter '%s': %g", pI->str, valf);
2030 break;
2031 case ENTIER :
2032 vali = atoi(value);
2033 (pI->action)(p,vali,pI->defaultfloat);
2034 Message::Info("Overriding parameter '%s': %d", pI->str, vali);
2035 break;
2036 }
2037 return;
2038 }
2039 }
2040
2041 Message::Error("Unknown solver parameter '%s'", name);
2042 }
2043
2044
2045 /* ------------------------------------------------------------------------ */
2046 /* dynamic CSR format */
2047 /* ------------------------------------------------------------------------ */
2048
cmpij(int ai,int aj,int bi,int bj)2049 static int cmpij(int ai,int aj,int bi,int bj){
2050 if(ai<bi)return -1;
2051 if(ai>bi)return 1;
2052 if(aj<bj)return -1;
2053 if(aj>bj)return 1;
2054 return 0;
2055 }
2056
alloc_ivec(long nl,long nh)2057 static int *alloc_ivec(long nl, long nh){
2058 int *v;
2059
2060 v=(int *)Malloc((size_t) ((nh-nl+1+1)*sizeof(int)));
2061 return v-nl+1;
2062 }
2063
free_ivec(int * v,long nl,long nh)2064 static void free_ivec(int *v, long nl, long nh){
2065 Free(v+nl-1);
2066 }
2067
2068 #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
2069 #define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
2070 #define M 7
2071 #define NSTACK 50
2072 #define M1 -1
2073
sort2(unsigned long n,double arr[],int ai[],int aj[])2074 static void sort2(unsigned long n, double arr[], int ai[] , int aj []){
2075 unsigned long i,ir=n,j,k,l=1;
2076 int *istack,jstack=0,tempi;
2077 double a,temp;
2078 int b,c;
2079
2080 istack=alloc_ivec(1,NSTACK);
2081 for (;;) {
2082 if (ir-l < M) {
2083 for (j=l+1;j<=ir;j++) {
2084 a=arr[j M1];
2085 b=ai[j M1];
2086 c=aj[j M1];
2087 for (i=j-1;i>=1;i--) {
2088 if (cmpij(ai[i M1],aj[i M1],b,c) <= 0) break;
2089 arr[i+1 M1]=arr[i M1];
2090 ai[i+1 M1]=ai[i M1];
2091 aj[i+1 M1]=aj[i M1];
2092 }
2093 arr[i+1 M1]=a;
2094 ai[i+1 M1]=b;
2095 aj[i+1 M1]=c;
2096 }
2097 if (!jstack) {
2098 free_ivec(istack,1,NSTACK);
2099 return;
2100 }
2101 ir=istack[jstack];
2102 l=istack[jstack-1];
2103 jstack -= 2;
2104 }
2105 else {
2106 k=(l+ir) >> 1;
2107 SWAP(arr[k M1],arr[l+1 M1])
2108 SWAPI(ai[k M1],ai[l+1 M1])
2109 SWAPI(aj[k M1],aj[l+1 M1])
2110 if (cmpij(ai[l+1 M1],aj[l+1 M1],ai[ir M1],aj[ir M1])>0){
2111 SWAP(arr[l+1 M1],arr[ir M1])
2112 SWAPI(ai[l+1 M1],ai[ir M1])
2113 SWAPI(aj[l+1 M1],aj[ir M1])
2114 }
2115 if (cmpij(ai[l M1],aj[l M1],ai[ir M1],aj[ir M1])>0){
2116 SWAP(arr[l M1],arr[ir M1])
2117 SWAPI(ai[l M1],ai[ir M1])
2118 SWAPI(aj[l M1],aj[ir M1])
2119 }
2120 if (cmpij(ai[l+1 M1],aj[l+1 M1],ai[l M1],aj[l M1])>0){
2121 SWAP(arr[l+1 M1],arr[l M1])
2122 SWAPI(ai[l+1 M1],ai[l M1])
2123 SWAPI(aj[l+1 M1],aj[l M1])
2124 }
2125 i=l+1;
2126 j=ir;
2127 a=arr[l M1];
2128 b=ai[l M1];
2129 c=aj[l M1];
2130 for (;;) {
2131 do i++; while (cmpij(ai[i M1],aj[i M1],b,c) < 0);
2132 do j--; while (cmpij(ai[j M1],aj[j M1],b,c) > 0);
2133 if (j < i) break;
2134 SWAP(arr[i M1],arr[j M1])
2135 SWAPI(ai[i M1],ai[j M1])
2136 SWAPI(aj[i M1],aj[j M1])
2137 }
2138 arr[l M1]=arr[j M1];
2139 arr[j M1]=a;
2140 ai[l M1]=ai[j M1];
2141 ai[j M1]=b;
2142 aj[l M1]=aj[j M1];
2143 aj[j M1]=c;
2144 jstack += 2;
2145 if (jstack > NSTACK) {
2146 Message::Error("NSTACK too small in sort2");
2147 }
2148 if (ir-i+1 >= j-l) {
2149 istack[jstack]=ir;
2150 istack[jstack-1]=i;
2151 ir=j-1;
2152 }
2153 else {
2154 istack[jstack]=j-1;
2155 istack[jstack-1]=l;
2156 l=i;
2157 }
2158 }
2159 }
2160 }
2161
2162 #undef M
2163 #undef NSTACK
2164 #undef SWAP
2165 #undef SWAPI
2166 #undef M1
2167
deblign(int nz,int * ptr,int * jptr,int * ai)2168 static void deblign ( int nz , int *ptr , int *jptr , int *ai){
2169 int i,ilign;
2170
2171 ilign = 1;
2172
2173 jptr[0] = 1;
2174 for(i=1; i<nz; i++) {
2175 if (ai[i-1] < ai[i]) {
2176 jptr[ilign++]=i+1;
2177 ai[i-1] = 0;
2178 }
2179 else{
2180 ai[i-1] = i+1;
2181 }
2182 }
2183 ai[nz-1] = 0;
2184 }
2185
csr_format(Sparse_Matrix * MM,int N)2186 void csr_format (Sparse_Matrix *MM, int N){
2187 int i,*ptr,*jptr,*ai,n,iptr,iptr2;
2188 double *a;
2189
2190 if(!N) return;
2191
2192 ptr = (int*)MM->ptr->array;
2193 jptr = (int*)MM->jptr->array;
2194 ai = (int*)MM->ai->array;
2195 a = (double*)MM->a->array;
2196 n = N;
2197 for(i=0;i<n;i++){
2198 iptr = jptr[i];
2199 while(iptr){
2200 iptr2 = iptr - 1;
2201 iptr = ptr[iptr2];
2202 ptr[iptr2] = i+1;
2203 }
2204 }
2205 sort2(List_Nbr(MM->a),a,ai,ptr);
2206 deblign(List_Nbr(MM->a),ptr,jptr,ai);
2207 jptr[N]=List_Nbr(MM->a)+1;
2208 }
2209
2210
restore_format(Sparse_Matrix * MM)2211 void restore_format (Sparse_Matrix *MM){
2212 char *temp;
2213
2214 temp = MM->ptr->array;
2215 MM->ptr->array = MM->ai->array;
2216 MM->ai->array = temp;
2217 }
2218