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: block_precon.c
7 *
8 * description: Composition of a block-preconditioner from
9 * preconditioner for the single blocks.
10 *
11 *******************************************************************************
12 *
13 * authors: Claus-Justus Heine
14 * Abteilung fuer Angewandte Mathematik
15 * Albert-Ludwigs-Universitaet Freiburg
16 * Hermann-Herder-Str. 10
17 * D-79104 Freiburg im Breisgau, Germany
18 *
19 * http://www.alberta-fem.de
20 *
21 * (c) by C.-J. Heine (2009)
22 *
23 ******************************************************************************/
24
25 #ifdef HAVE_CONFIG_H
26 # include "config.h"
27 #endif
28
29 #include "alberta_intern.h"
30 #include "alberta.h"
31
32 typedef struct precon_node PRECON_NODE;
33 struct precon_node
34 {
35 const PRECON *precon;
36 OEM_PRECON type;
37 size_t dim;
38 DOF_MATRIX *A; /* unchained flat copy of the original matrix block */
39 DOF_SCHAR_VEC *mask;
40 DOF_REAL_VEC_D *accu; /* for BlkSSORPrecon */
41 DOF_MATRIX *A_row; /* Row without diagonal block for SSORPrecon */
42 DBL_LIST_NODE chain;
43 };
44
45 typedef struct block_precon_data
46 {
47 PRECON precon;
48
49 const DOF_MATRIX *matrix;
50 const DOF_SCHAR_VEC *mask;
51 size_t dim;
52
53 OEM_PRECON block_precon; /* how to combine the individual precons */
54 DBL_LIST_NODE chain; /* the list of the individual precons */
55
56 REAL omega; /* for BlkSSORPrecon */
57 int n_iter; /* for BlkSSORPrecon */
58
59 DOF_REAL_VEC_D *rhs; /* for BlkSSORPrecon */
60 DOF_REAL_VEC_D *r_skel; /* for BLKSSORPrecon */
61
62 SCRATCH_MEM scratch;
63 } BLOCK_PRECON_DATA;
64
65 /* Simplest method: just do block-diagonal precon.
66 *
67 *
68 * P = diag(P_1, P_2, ...) where P_1, ... are listed in precon_chain.
69 *
70 * Maybe we could give block-SSOR a chance:
71 *
72 * A = (a_{ij}) where a_{ij} are matrices, the precons listed in
73 * precon_chain serve as approximate inverse to a_{ii}
74 *
75 */
block_diag_precon(void * pd,int dim,REAL * r)76 static void block_diag_precon(void *pd, int dim, REAL *r)
77 {
78 BLOCK_PRECON_DATA *data = (struct block_precon_data *)pd;
79 struct precon_node *prec;
80
81 CHAIN_FOREACH(prec, data, struct precon_node) {
82 if (prec->type != NoPrecon) {
83 prec->precon->precon(prec->precon->precon_data, prec->dim, r);
84 }
85 r += prec->dim;
86 }
87 }
88
block_diag_exit_precon(void * pd)89 static void block_diag_exit_precon(void *pd)
90 {
91 BLOCK_PRECON_DATA *data = (struct block_precon_data *)pd;
92 struct precon_node *prec;
93
94 CHAIN_FOREACH(prec, data, struct precon_node) {
95 if (prec->type != NoPrecon) {
96 prec->precon->exit_precon(prec->precon->precon_data);
97 }
98 }
99
100 /* Release all our resources. */
101 SCRATCH_MEM_ZAP(data->scratch);
102 }
103
block_diag_init_precon(void * pd)104 static bool block_diag_init_precon(void *pd)
105 {
106 BLOCK_PRECON_DATA *data = (struct block_precon_data *)pd;
107 struct precon_node *prec;
108
109 CHAIN_FOREACH(prec, data, struct precon_node) {
110 if (prec->type != NoPrecon) {
111 update_dof_matrix_sub_chain(prec->A);
112 if (prec->mask) {
113 update_dof_schar_vec_sub_chain(prec->mask);
114 }
115 if (!prec->precon->init_precon(prec->precon->precon_data)) {
116 return false;
117 }
118 prec->dim = dof_real_vec_d_length(prec->A->row_fe_space);
119 }
120 }
121 return true;
122 }
123
124 /* Use a block SSOR method, using preconditioners for the diagonal
125 * elements to approximate the inverse of the diagonal blocks.
126 */
block_SSOR_precon(void * pd,int dim,REAL * r)127 static void block_SSOR_precon(void *pd, int dim, REAL *r)
128 {
129 BLOCK_PRECON_DATA *data = (struct block_precon_data *)pd;
130 struct precon_node *prec;
131 int iter;
132 REAL *rbase = r;
133
134 /* FIXME: do we need to mask out boundary DOFs? */
135 distribute_to_dof_real_vec_d_skel(data->r_skel, r);
136 dof_copy_dow(data->r_skel, data->rhs);
137 dset(dim, 0.0, r, 1);
138
139 for (iter = 0; iter < data->n_iter; iter++) {
140 /* Forward Gauss-Seidel iteration */
141 r = rbase;
142 CHAIN_FOREACH(prec, data, struct precon_node) {
143 data->r_skel = CHAIN_NEXT(data->r_skel, DOF_REAL_VEC_D);
144 dcopy(prec->dim, data->rhs->vec, 1, prec->accu->vec, 1);
145 dof_gemv_dow(NoTranspose,
146 -1.0, prec->A_row, data->mask, data->r_skel,
147 1.0, prec->accu);
148
149 /* Pseudo-invert the diagonal block */
150 if (prec->type != NoPrecon) {
151 prec->precon->precon(
152 prec->precon->precon_data, prec->dim, prec->accu->vec);
153 }
154
155 /* r[idx] = omega A_diag^{-1} accu + (1-omega) r */
156 dscal(prec->dim, data->omega, prec->accu->vec, 1);
157 dxpay(prec->dim, prec->accu->vec, 1, 1.0 - data->omega, r, 1);
158
159 data->rhs = CHAIN_NEXT(data->rhs, DOF_REAL_VEC_D);
160 data->mask = data->mask ? CHAIN_NEXT(data->mask, DOF_SCHAR_VEC) : NULL;
161 r += prec->dim;
162 }
163 /* Backward Gauss-Seidel iteration */
164 r = rbase + dim;
165 CHAIN_FOREACH_REV(prec, data, struct precon_node) {
166 r -= prec->dim;
167 data->rhs = CHAIN_PREV(data->rhs, DOF_REAL_VEC_D);
168 data->mask = data->mask ? CHAIN_PREV(data->mask, DOF_SCHAR_VEC) : NULL;
169 dcopy(prec->dim, data->rhs->vec, 1, prec->accu->vec, 1);
170 dof_gemv_dow(NoTranspose,
171 -1.0, prec->A_row, data->mask, data->r_skel,
172 1.0, prec->accu);
173
174 /* r[idx] = omega A_diag^{-1} accu + (1-omega) r */
175 if (prec->type != NoPrecon) {
176 prec->precon->precon(
177 prec->precon->precon_data, prec->dim, prec->accu->vec);
178 }
179 dscal(prec->dim, data->omega, prec->accu->vec, 1);
180 dxpay(prec->dim, prec->accu->vec, 1, 1.0 - data->omega, r, 1);
181
182 data->r_skel = CHAIN_PREV(data->r_skel, DOF_REAL_VEC_D);
183 }
184 }
185 }
186
block_SSOR_exit_precon(void * pd)187 static void block_SSOR_exit_precon(void *pd)
188 {
189 BLOCK_PRECON_DATA *data = (struct block_precon_data *)pd;
190 struct precon_node *prec;
191
192 CHAIN_FOREACH(prec, data, struct precon_node) {
193 if (prec->type != NoPrecon) {
194 prec->precon->exit_precon(prec->precon->precon_data);
195 }
196 free_dof_real_vec_d(prec->accu);
197 }
198
199 free_dof_real_vec_d(data->rhs);
200
201 /* Release all our resources. */
202 SCRATCH_MEM_ZAP(data->scratch);
203 }
204
block_SSOR_init_precon(void * pd)205 static bool block_SSOR_init_precon(void *pd)
206 {
207 BLOCK_PRECON_DATA *data = (struct block_precon_data *)pd;
208 struct precon_node *prec;
209
210 CHAIN_FOREACH(prec, data, struct precon_node) {
211 update_dof_matrix_sub_chain(prec->A);
212 update_dof_matrix_sub_chain(prec->A_row);
213 if (prec->mask) {
214 update_dof_schar_vec_sub_chain(prec->mask);
215 }
216 if (!prec->precon->init_precon(prec->precon->precon_data)) {
217 return false;
218 }
219 prec->dim = dof_real_vec_d_length(prec->A->row_fe_space);
220 }
221 return true;
222 }
223
_AI_vget_block_diag_precon(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,va_list ap)224 const PRECON *_AI_vget_block_diag_precon(const DOF_MATRIX *A,
225 const DOF_SCHAR_VEC *mask,
226 int info,
227 va_list ap)
228 {
229 FUNCNAME("_AI_vget_block_diag_precon");
230 PRECON_TYPE precon;
231 OEM_PRECON sub_type = NoPrecon;
232 int i, nblocks = COL_CHAIN_LENGTH(A);
233
234 precon.type = BlkDiagPrecon;
235
236 for (i = 0; i < nblocks && sub_type != PreconRepeat; i++) {
237 TEST_EXIT(i < N_BLOCK_PRECON_MAX,
238 "Sorry, only up to %d x %d blocks are supported.\n",
239 N_BLOCK_PRECON_MAX, N_BLOCK_PRECON_MAX);
240 precon.param.BlkDiag.precon[i].type =
241 sub_type = (OEM_PRECON)va_arg(ap, int);
242 if (sub_type == __SSORPrecon) {
243 precon.param.BlkDiag.precon[i].param.__SSOR.omega = va_arg(ap, REAL);
244 precon.param.BlkDiag.precon[i].param.__SSOR.n_iter = va_arg(ap, int);
245 }
246 }
247
248 return _AI_get_block_precon(A, mask, info, &precon);
249 }
250
_AI_vget_block_SSOR_precon(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,va_list ap)251 const PRECON *_AI_vget_block_SSOR_precon(const DOF_MATRIX *A,
252 const DOF_SCHAR_VEC *mask,
253 int info,
254 va_list ap)
255 {
256 FUNCNAME("_AI_vget_block_diag_precon");
257 PRECON_TYPE precon;
258 OEM_PRECON sub_type = NoPrecon;
259 int i, nblocks = COL_CHAIN_LENGTH(A);
260
261 precon.type = BlkSSORPrecon;
262 precon.param.BlkSSOR.omega = va_arg(ap, REAL);
263 precon.param.BlkSSOR.n_iter = va_arg(ap, int);
264
265 for (i = 0; i < nblocks && sub_type != PreconRepeat; i++) {
266 TEST_EXIT(i < N_BLOCK_PRECON_MAX,
267 "Sorry, only up to %d x %d blocks are supported.\n",
268 N_BLOCK_PRECON_MAX, N_BLOCK_PRECON_MAX);
269 precon.param.BlkDiag.precon[i].type =
270 sub_type = (OEM_PRECON)va_arg(ap, int);
271 if (sub_type == __SSORPrecon) {
272 precon.param.BlkDiag.precon[i].param.__SSOR.omega = va_arg(ap, REAL);
273 precon.param.BlkDiag.precon[i].param.__SSOR.n_iter = va_arg(ap, int);
274 }
275 }
276
277 return _AI_get_block_precon(A, mask, info, &precon);
278 }
279
_AI_get_block_diag_precon(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,...)280 const PRECON *_AI_get_block_diag_precon(const DOF_MATRIX *A,
281 const DOF_SCHAR_VEC *mask,
282 int info,
283 ...)
284 {
285 const PRECON *precon;
286 va_list ap;
287
288 va_start(ap, info);
289 precon = _AI_vget_block_diag_precon(A, mask, info, ap);
290 va_end(ap);
291
292 return precon;
293 }
294
_AI_get_block_precon(const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,int info,const PRECON_TYPE * prec_type)295 const PRECON *_AI_get_block_precon(const DOF_MATRIX *A,
296 const DOF_SCHAR_VEC *mask,
297 int info,
298 const PRECON_TYPE *prec_type)
299 {
300 FUNCNAME("_AI_get_block_precon");
301 BLOCK_PRECON_DATA *data;
302 SCRATCH_MEM scrmem;
303 bool repeat = false;
304 OEM_PRECON sub_type = NoPrecon;
305 int i = 0;
306 const FE_SPACE *col_fe_space;
307
308 TEST_EXIT(ROW_CHAIN_LENGTH(A) == COL_CHAIN_LENGTH(A),
309 "Makes sense for quadratic block-matrices only.\n");
310 TEST_EXIT(ROW_CHAIN_LENGTH(A) < N_BLOCK_PRECON_MAX,
311 "Only implemented for up to %d x %d blocks.\n",
312 N_BLOCK_PRECON_MAX, N_BLOCK_PRECON_MAX);
313
314 SCRATCH_MEM_INIT(scrmem);
315 data = SCRATCH_MEM_ALLOC(scrmem, 1, BLOCK_PRECON_DATA);
316 memset(data, 0, sizeof(*data));
317 SCRATCH_MEM_CPY(data->scratch, scrmem);
318
319 CHAIN_INIT(data);
320
321 data->matrix = A;
322 data->mask = mask;
323 data->dim = dof_real_vec_d_length(A->row_fe_space);
324
325 data->precon.precon_data = data;
326
327 col_fe_space = A->col_fe_space ? A->col_fe_space : A->row_fe_space;
328
329 switch (prec_type->type) {
330 case BlkDiagPrecon:
331 data->block_precon = DiagPrecon;
332 data->precon.precon = block_diag_precon;
333 data->precon.init_precon = block_diag_init_precon;
334 data->precon.exit_precon = block_diag_exit_precon;
335 break;
336 case BlkSSORPrecon:
337 data->block_precon = SSORPrecon;
338 data->precon.precon = block_SSOR_precon;
339 data->precon.init_precon = block_SSOR_init_precon;
340 data->precon.exit_precon = block_SSOR_exit_precon;
341 data->omega = prec_type->param.BlkSSOR.omega;
342 data->n_iter = prec_type->param.BlkSSOR.n_iter;
343
344 data->rhs = get_dof_real_vec_d("SSOR rhs", col_fe_space);
345 data->r_skel = init_dof_real_vec_d_skel(
346 SCRATCH_MEM_ALLOC(scrmem, CHAIN_LENGTH(col_fe_space), DOF_REAL_VEC_D),
347 "SSOR r skeleton",
348 col_fe_space);
349
350 break;
351 default:
352 ERROR_EXIT("Precon type %d is not implemented.\n", prec_type->type);
353 data->block_precon = NoPrecon;
354 }
355
356 COL_CHAIN_DO(A, const DOF_MATRIX) {
357 REAL omega = 1.0;
358 int n_iter = 2;
359 int ilu_level = 0;
360 PRECON_NODE *node = SCRATCH_MEM_ALLOC(scrmem, 1, PRECON_NODE);
361 memset(node, 0, sizeof(*node));
362 CHAIN_INIT(node);
363
364 CHAIN_ADD_TAIL(data, node);
365
366 /* Generate flat, unchained copies from the diagonal block and
367 * boundary masks. This way one can use the ordinary
368 * preconditioners with those flat copies.
369 */
370 node->A = dof_matrix_sub_chain(scrmem, A, 0x1, 0x1);
371 node->mask = mask ? dof_schar_vec_sub_chain(scrmem, mask, 0x1) : NULL;
372
373 if (data->block_precon == SSORPrecon) {
374 /* Extract a sub-chain with all elements of the row safe the diagonal */
375 node->A_row = dof_matrix_sub_chain(scrmem, A, 0x1, ~0x1);
376 /* We need some scratch memory for the rhs and so on. */
377 node->accu = get_dof_real_vec_d("SSOR accu", col_fe_space->unchained);
378 }
379
380 /* Compute the size of our fraction of the index set. */
381 node->dim = dof_real_vec_d_length(node->A->row_fe_space);
382
383 if (!repeat && prec_type->param.BlkDiag.precon[i].type == PreconRepeat) {
384 repeat = true;
385 }
386 if (repeat) {
387 node->type = sub_type;
388 } else {
389 node->type = sub_type = prec_type->param.BlkDiag.precon[i].type;
390 }
391 switch (node->type) {
392 case NoPrecon:
393 break;
394 case DiagPrecon:
395 node->precon = get_diag_precon(node->A, node->mask);
396 break;
397 case HBPrecon:
398 node->precon =
399 get_HB_precon(node->A, node->mask, info);
400 break;
401 case BPXPrecon:
402 node->precon =
403 get_BPX_precon(node->A, node->mask, info);
404 break;
405 case __SSORPrecon:
406 if (!repeat) {
407 omega = prec_type->param.BlkDiag.precon[i].param.__SSOR.omega;
408 n_iter = prec_type->param.BlkDiag.precon[i].param.__SSOR.n_iter;
409
410 TEST(0.0 <= omega && omega <= 2.0,
411 "SSORPrecon: omega = %e???\n", omega);
412 TEST(0 <= n_iter && n_iter < 1000,
413 "SSORPrecon: #iter = %d???\n", n_iter);
414 }
415 case SSORPrecon:
416 node->precon = get_SSOR_precon(node->A, node->mask, omega, n_iter);
417 break;
418 case ILUkPrecon:
419 if (!repeat) {
420 ilu_level = prec_type->param.BlkDiag.precon[i].param.ILUk.level;
421 }
422 node->precon = get_ILUk_precon(node->A, node->mask, ilu_level, info);
423 break;
424 default:
425 ERROR("Unknow precon-type %d, ignoring it.\n");
426 node->type = NoPrecon;
427 break;
428 }
429
430 i++;
431
432 A = ROW_CHAIN_NEXT(A, const DOF_MATRIX); /* walk over the diagonal */
433 col_fe_space = CHAIN_NEXT(col_fe_space, const FE_SPACE);
434 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
435
436 return &data->precon;
437 }
438