1 /*******************************************************************************
2 * ALBERTA: an Adaptive multi Level finite element toolbox using
3 * Bisectioning refinement and Error control by Residual
4 * Techniques for scientific Applications
5 *
6 * file: oem_solve.c
7 *
8 * description: Part of the OEM interface which does not depend on the
9 * DIM_OF_WORLD dimensionality of a DOF_REAL[_D]_VEC[_D]
10 *
11 *******************************************************************************
12 *
13 * authors: Alfred Schmidt
14 * Zentrum fuer Technomathematik
15 * Fachbereich 3 Mathematik/Informatik
16 * Universitaet Bremen
17 * Bibliothekstr. 2
18 * D-28359 Bremen, Germany
19 *
20 * Kunibert G. Siebert
21 * Institut fuer Mathematik
22 * Universitaet Augsburg
23 * Universitaetsstr. 14
24 * D-86159 Augsburg, Germany
25 *
26 * Claus-Justus Heine
27 * Abteilung fuer Angewandte Mathematik
28 * Albert-Ludwigs-Universitaet Freiburg
29 * Hermann-Herder-Str. 10
30 * D-79104 Freiburg im Breisgau, Germany
31 *
32 * http://www.alberta-fem.de
33 *
34 * (c) by A. Schmidt and K.G. Siebert (1996-2003),
35 * C.-J. Heine (1998-2009)
36 *
37 ******************************************************************************/
38
39 #ifdef HAVE_CONFIG_H
40 # include "config.h"
41 #endif
42
43 #include "alberta_intern.h"
44 #include "alberta.h"
45
46 struct mv_data
47 {
48 const DOF_MATRIX *matrix;
49 MatrixTranspose transpose;
50 const FE_SPACE *domain_space;
51 const FE_SPACE *range_space;
52 int dim;
53
54 const DOF_SCHAR_VEC *mask;
55
56 DOF_REAL_VEC_D *x_skel;
57 DOF_REAL_VEC_D *y_skel;
58
59 SCRATCH_MEM scratch;
60 };
61
62 /*---8<---------------------------------------------------------------------*/
63 /*--- y <- A x or y <- A^T x ---*/
64 /*--------------------------------------------------------------------->8---*/
65
oem_mat_vec(void * ud,int dim,const REAL * x,REAL * y)66 int oem_mat_vec(void *ud, int dim, const REAL *x, REAL *y)
67 {
68 FUNCNAME("mat_vec_s");
69 struct mv_data *data = (struct mv_data *)ud;
70 DOF_REAL_VEC_D *dof_x = data->x_skel;
71 DOF_REAL_VEC_D *dof_y = data->y_skel;
72
73 TEST_EXIT(dim == data->dim, "argument dim != FE_SPACE dim\n");
74
75 distribute_to_dof_real_vec_d_skel(data->x_skel, x);
76 distribute_to_dof_real_vec_d_skel(data->y_skel, y);
77
78 dof_mv_dow(data->transpose, data->matrix, data->mask, dof_x, dof_y);
79
80 return 0; /* ??? cH asks: what the heck is the meaning of the return value? */
81 }
82
init_oem_mat_vec(void ** datap,MatrixTranspose transpose,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask)83 OEM_MV_FCT init_oem_mat_vec(void **datap,
84 MatrixTranspose transpose, const DOF_MATRIX *A,
85 const DOF_SCHAR_VEC *mask)
86 {
87 /*FUNCNAME("init_oem_mat_vec");*/
88 struct mv_data *data;
89 SCRATCH_MEM scrmem;
90
91 SCRATCH_MEM_INIT(scrmem);
92
93 data = SCRATCH_MEM_ALLOC(scrmem, 1, struct mv_data);
94 memset(data, 0, sizeof(*data));
95 SCRATCH_MEM_CPY(data->scratch, scrmem);
96
97 data->matrix = A;
98 data->transpose = transpose;
99
100 data->mask = mask;
101
102 if (transpose == NoTranspose) {
103 data->range_space = A->row_fe_space;
104 data->domain_space = A->col_fe_space;
105 } else {
106 data->range_space = A->col_fe_space;
107 data->domain_space = A->row_fe_space;
108 }
109
110 data->x_skel = init_dof_real_vec_d_skel(
111 SCRATCH_MEM_ALLOC(scrmem, CHAIN_LENGTH(data->domain_space), DOF_REAL_VEC_D),
112 "x skel",
113 data->domain_space);
114 data->y_skel = init_dof_real_vec_d_skel(
115 SCRATCH_MEM_ALLOC(scrmem, CHAIN_LENGTH(data->range_space), DOF_REAL_VEC_D),
116 "y skel",
117 data->range_space);
118
119 data->dim = dof_real_vec_d_length(data->range_space);
120
121 *datap = data;
122
123 return oem_mat_vec;
124 }
125
exit_oem_mat_vec(void * vdata)126 void exit_oem_mat_vec(void *vdata)
127 {
128 struct mv_data *data = (struct mv_data *)vdata;
129
130 SCRATCH_MEM_ZAP(data->scratch);
131 }
132
vinit_oem_precon(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,OEM_PRECON prec_type,va_list ap)133 const PRECON *vinit_oem_precon(const DOF_MATRIX *A,
134 const DOF_SCHAR_VEC *mask,
135 int info, OEM_PRECON prec_type,
136 va_list ap)
137 {
138 const PRECON *precon;
139 REAL omega = 1.0;
140 int n_iter = 2;
141 int ilu_level = 0;
142
143 if ((!ROW_CHAIN_SINGLE(A) || !COL_CHAIN_SINGLE(A))
144 &&
145 prec_type < BlkDiagPrecon) {
146 precon = _AI_get_block_diag_precon(A, mask, info, prec_type, PreconRepeat);
147 } else {
148 switch(prec_type) {
149 case NoPrecon:
150 precon = NULL;
151 break;
152 case DiagPrecon:
153 precon = get_diag_precon(A, mask);
154 break;
155 case HBPrecon:
156 TEST_EXIT(ROW_CHAIN_SINGLE(A) && COL_CHAIN_SINGLE(A),
157 "HB-preonditioner does not make sense for "
158 "horizontal direct sums\n");
159 precon = get_HB_precon(A, mask, info);
160 break;
161 case BPXPrecon:
162 TEST_EXIT(ROW_CHAIN_SINGLE(A) && COL_CHAIN_SINGLE(A),
163 "BPX-preonditioner does not make sense for "
164 "horizontal direct sums\n");
165 precon = get_BPX_precon(A, mask, info);
166 break;
167 case __SSORPrecon:
168 omega = va_arg(ap, REAL);
169 n_iter = va_arg(ap, int);
170
171 TEST(0.0 <= omega && omega <= 2.0,
172 "SSORPrecon: omega = %e???\n", omega);
173 TEST(0 <= n_iter && n_iter < 10,
174 "SSORPrecon: #iter = %d???\n", n_iter);
175 case SSORPrecon:
176 TEST_EXIT(ROW_CHAIN_SINGLE(A) && COL_CHAIN_SINGLE(A),
177 "SSOR-preconditioner not implemented for horizontal "
178 "direct sums. Very sorry.\n");
179 precon = get_SSOR_precon(A, mask, omega, n_iter);
180 break;
181 case ILUkPrecon:
182 ilu_level = va_arg(ap, int);
183 precon = get_ILUk_precon(A, mask, ilu_level, info);
184 break;
185 case BlkDiagPrecon:
186 precon = _AI_vget_block_diag_precon(A, mask, info, ap);
187 break;
188 case BlkSSORPrecon:
189 precon = _AI_vget_block_SSOR_precon(A, mask, info, ap);
190 break;
191 default:
192 ERROR_EXIT("Unknown precon type: %d\n", prec_type);
193 break;
194 }
195 }
196
197 return precon;
198 }
199
init_oem_precon(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,OEM_PRECON prec_type,...)200 const PRECON *init_oem_precon(const DOF_MATRIX *A,
201 const DOF_SCHAR_VEC *mask,
202 int info, OEM_PRECON prec_type,
203 ... /* ssor_omega, ssor_n_iter etc. */)
204 {
205 const PRECON *precon;
206 va_list ap;
207
208 va_start(ap, prec_type);
209 precon = vinit_oem_precon(A, mask, info, prec_type, ap);
210 va_end(ap);
211
212 return precon;
213 }
214
init_precon_from_type(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,const PRECON_TYPE * prec_type)215 const PRECON *init_precon_from_type(const DOF_MATRIX *A,
216 const DOF_SCHAR_VEC *mask,
217 int info,
218 const PRECON_TYPE *prec_type)
219 {
220 FUNCNAME("init_precon_from_type");
221
222 const PRECON *precon;
223
224 switch (prec_type->type) {
225 case __SSORPrecon:
226 precon = init_oem_precon(A, mask, info,
227 __SSORPrecon,
228 prec_type->param.__SSOR.omega,
229 prec_type->param.__SSOR.n_iter);
230 break;
231 case ILUkPrecon:
232 precon =
233 init_oem_precon(A, mask, info, ILUkPrecon, prec_type->param.ILUk.level);
234 break;
235 case BlkSSORPrecon:
236 case BlkDiagPrecon:
237 precon = _AI_get_block_precon(A, mask, info, prec_type);
238 break;
239 default:
240 precon = init_oem_precon(A, mask, info, prec_type->type, PreconEnd);
241 break;
242 }
243
244 return precon;
245 }
246
init_oem_solve(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,REAL tol,const PRECON * precon,int restart,int max_iter,int info)247 OEM_DATA *init_oem_solve(const DOF_MATRIX *A,
248 const DOF_SCHAR_VEC *mask,
249 REAL tol, const PRECON *precon,
250 int restart, int max_iter, int info)
251 {
252 FUNCNAME("init_oem_solve");
253 OEM_DATA *oem;
254 const MatrixTranspose transpose = NoTranspose;
255
256 TEST_EXIT(FE_SPACE_EQ_P(A->row_fe_space, A->col_fe_space),
257 "Row and column FE_SPACEs don't match!\n");
258
259 oem = MEM_CALLOC(1, OEM_DATA);
260 oem->mat_vec = init_oem_mat_vec(&oem->mat_vec_data, transpose, A, mask);
261
262 if (precon) {
263 if (precon->init_precon && !(*precon->init_precon)(precon->precon_data)) {
264 precon = NULL;
265 MSG("init_precon() failed, disabling preconditioner!\n");
266 } else {
267 oem->left_precon_data = precon->precon_data;
268 oem->left_precon = precon->precon;
269 }
270 }
271
272 oem->ws = NULL; /* Let the solvers handle this point themselves. */
273 oem->tolerance = tol;
274 oem->restart = restart;
275 oem->max_iter = max_iter;
276 oem->info = MAX(0, info);
277
278 return oem;
279 }
280
release_oem_solve(const OEM_DATA * _oem)281 void release_oem_solve(const OEM_DATA *_oem)
282 {
283 OEM_DATA *oem = (OEM_DATA *)_oem;
284 const PRECON *precon;
285
286 exit_oem_mat_vec(oem->mat_vec_data);
287 if ((precon = (PRECON *)oem->left_precon_data) != NULL &&
288 precon && precon->exit_precon) {
289 precon->exit_precon((void *)precon);
290 }
291 MEM_FREE(oem, 1, OEM_DATA);
292 }
293
get_oem_solver(OEM_SOLVER solver)294 OEM_MV_FCT get_oem_solver(OEM_SOLVER solver)
295 {
296
297 switch (solver) {
298 case BiCGStab:
299 return (OEM_MV_FCT)oem_bicgstab;
300 case CG:
301 return (OEM_MV_FCT)oem_cg;
302 case TfQMR:
303 return (OEM_MV_FCT)oem_tfqmr;
304 case GMRes:
305 return (OEM_MV_FCT)oem_gmres;
306 case GMRes_k:
307 return (OEM_MV_FCT)oem_gmres_k;
308 case ODir:
309 return (OEM_MV_FCT)oem_odir;
310 case ORes:
311 return (OEM_MV_FCT)oem_ores;
312 case SymmLQ:
313 return (OEM_MV_FCT)oem_symmlq;
314 default:
315 ERROR_EXIT("unknown OEM solver %d\n", (int) solver);
316 return NULL;
317 }
318 }
319
320 /* scalar problems */
321
call_oem_solve_s(const OEM_DATA * _oem,OEM_SOLVER solver,const DOF_REAL_VEC * f,DOF_REAL_VEC * u)322 int call_oem_solve_s(const OEM_DATA *_oem, OEM_SOLVER solver,
323 const DOF_REAL_VEC *f, DOF_REAL_VEC *u)
324 {
325 int iter, restart, dim;
326 OEM_DATA *oem = (OEM_DATA *)_oem;
327 REAL *uvec, *fvec;
328
329 TEST_EXIT(FE_SPACE_EQ_P(f->fe_space, u->fe_space),
330 "Row and column FE_SPACEs don't match!\n");
331
332 dim = dof_real_vec_length(f->fe_space);
333
334 if (!CHAIN_SINGLE(u)) {
335 uvec = MEM_ALLOC(dim, REAL);
336 fvec = MEM_ALLOC(dim, REAL);
337 copy_from_dof_real_vec(uvec, u);
338 copy_from_dof_real_vec(fvec, f);
339 } else {
340 FOR_ALL_FREE_DOFS(u->fe_space->admin,
341 if (dof < dim) u->vec[dof] = f->vec[dof] = 0.0);
342 fvec = f->vec;
343 uvec = u->vec;
344 }
345
346 switch (solver) {
347 case BiCGStab:
348 iter = oem_bicgstab(oem, dim, fvec, uvec);
349 break;
350 case CG:
351 iter = oem_cg(oem, dim, fvec, uvec);
352 break;
353 case TfQMR:
354 iter = oem_tfqmr(oem, dim, fvec, uvec);
355 break;
356 case GMRes:
357 restart = oem->restart;
358 oem->restart = MAX(0, MIN(oem->restart, dim));
359 iter = oem_gmres(oem, dim, fvec, uvec);
360 oem->restart = restart;
361 break;
362 case GMRes_k:
363 restart = oem->restart;
364 oem->restart = MAX(0, MIN(oem->restart, dim));
365 iter = oem_gmres_k(oem, dim, fvec, uvec);
366 oem->restart = restart;
367 break;
368 case ODir:
369 iter = oem_odir(oem, dim, fvec, uvec);
370 break;
371 case ORes:
372 iter = oem_ores(oem, dim, fvec, uvec);
373 break;
374 case SymmLQ:
375 iter = oem_symmlq(oem, dim, fvec, uvec);
376 break;
377 default:
378 iter = -1; /* make the compiler happy ... -> no-return attribute? */
379 ERROR_EXIT("unknown OEM solver %d\n", (int) solver);
380 }
381
382 if (!CHAIN_SINGLE(u)) {
383 copy_to_dof_real_vec(u, uvec);
384 MEM_FREE(uvec, dim, REAL);
385 MEM_FREE(fvec, dim, REAL);
386 }
387
388 return iter;
389 }
390
oem_solve_s(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * f,DOF_REAL_VEC * u,OEM_SOLVER solver,REAL tol,const PRECON * precon,int restart,int max_iter,int info)391 int oem_solve_s(const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
392 const DOF_REAL_VEC *f, DOF_REAL_VEC *u,
393 OEM_SOLVER solver, REAL tol, const PRECON *precon,
394 int restart, int max_iter, int info)
395 {
396 const OEM_DATA *oem;
397 int iter;
398
399 oem = init_oem_solve(A, mask, tol, precon, restart, max_iter, info);
400 iter = call_oem_solve_s(oem, solver, f, u);
401 release_oem_solve(oem);
402
403 return iter;
404 }
405
406 /* DIM_OF_WORLD interface */
407
call_oem_solve_dow(const OEM_DATA * _oem,OEM_SOLVER solver,const DOF_REAL_VEC_D * f,DOF_REAL_VEC_D * u)408 int call_oem_solve_dow(const OEM_DATA *_oem, OEM_SOLVER solver,
409 const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u)
410 {
411 int iter, restart, dim;
412 OEM_DATA *oem = (OEM_DATA *)_oem;
413 REAL *uvec, *fvec;
414
415 TEST_EXIT(FE_SPACE_EQ_P(f->fe_space, u->fe_space),
416 "Row and column FE_SPACEs don't match!\n");
417
418 dim = dof_real_vec_d_length(f->fe_space);
419
420 if (!CHAIN_SINGLE(u)) {
421 uvec = MEM_ALLOC(dim, REAL);
422 fvec = MEM_ALLOC(dim, REAL);
423 copy_from_dof_real_vec_d(uvec, u);
424 copy_from_dof_real_vec_d(fvec, f);
425 } else {
426 const DOF_ADMIN *admin = u->fe_space->admin;
427
428 fvec = f->vec;
429 uvec = u->vec;
430
431 FOR_ALL_FREE_DOFS(u->fe_space->admin,
432 if (dof >= admin->size_used) {
433 break;
434 }
435 if (u->stride == 1) {
436 u->vec[dof] = f->vec[dof] = 0.0;
437 } else {
438 SET_DOW(0.0, ((REAL_D *)u->vec)[dof]);
439 SET_DOW(0.0, ((REAL_D *)f->vec)[dof]);
440 });
441 }
442
443 switch (solver) {
444 case BiCGStab:
445 iter = oem_bicgstab(oem, dim, fvec, uvec);
446 break;
447 case CG:
448 iter = oem_cg(oem, dim, fvec, uvec);
449 break;
450 case TfQMR:
451 iter = oem_tfqmr(oem, dim, fvec, uvec);
452 break;
453 case GMRes:
454 restart = oem->restart;
455 oem->restart = MAX(0, MIN(oem->restart, dim));
456 iter = oem_gmres(oem, dim, (REAL *) f->vec, (REAL *) u->vec);
457 oem->restart = restart;
458 break;
459 case GMRes_k:
460 restart = oem->restart;
461 oem->restart = MAX(0, MIN(oem->restart, dim));
462 iter = oem_gmres_k(oem, dim, fvec, uvec);
463 oem->restart = restart;
464 break;
465 case ODir:
466 iter = oem_odir(oem, dim, fvec, uvec);
467 break;
468 case ORes:
469 iter = oem_ores(oem, dim, fvec, uvec);
470 break;
471 case SymmLQ:
472 iter = oem_symmlq(oem, dim, fvec, uvec);
473 break;
474 default:
475 iter = -1; /* make the compiler happy ... -> no-return attribute? */
476 ERROR_EXIT("unknown OEM solver %d\n", (int) solver);
477 }
478
479 if (!CHAIN_SINGLE(u)) {
480 copy_to_dof_real_vec_d(u, uvec);
481 MEM_FREE(uvec, dim, REAL);
482 MEM_FREE(fvec, dim, REAL);
483 }
484
485 return iter;
486 }
487
oem_solve_dow(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * f,DOF_REAL_VEC_D * u,OEM_SOLVER solver,REAL tol,const PRECON * precon,int restart,int max_iter,int info)488 int oem_solve_dow(const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
489 const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u,
490 OEM_SOLVER solver, REAL tol, const PRECON *precon,
491 int restart, int max_iter, int info)
492 {
493 const OEM_DATA *oem;
494 int iter;
495
496 oem = init_oem_solve(A, mask, tol, precon, restart, max_iter, info);
497 iter = call_oem_solve_dow(oem, solver, f, u);
498 release_oem_solve(oem);
499
500 return iter;
501 }
502