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: dof_admin.c */
7 /* */
8 /* description: implementation of the administration of DOFs */
9 /* */
10 /*--------------------------------------------------------------------------*/
11 /* */
12 /* authors: Alfred Schmidt */
13 /* Zentrum fuer Technomathematik */
14 /* Fachbereich 3 Mathematik/Informatik */
15 /* Universitaet Bremen */
16 /* Bibliothekstr. 2 */
17 /* D-28359 Bremen, Germany */
18 /* */
19 /* Kunibert G. Siebert */
20 /* Institut fuer Mathematik */
21 /* Universitaet Augsburg */
22 /* Universitaetsstr. 14 */
23 /* D-86159 Augsburg, Germany */
24 /* */
25 /* Daniel Koester */
26 /* Institut fuer Mathematik */
27 /* Universitaet Augsburg */
28 /* Universitaetsstr. 14 */
29 /* D-86159 Augsburg, Germany */
30 /* */
31 /* Claus-Justus Heine */
32 /* Abteilung fuer Angewandte Mathematik */
33 /* Universitaet Freiburg */
34 /* Hermann-Herder-Strasse 10 */
35 /* 79104 Freiburg, Germany */
36 /* */
37 /* http://www.alberta-fem.de */
38 /* */
39 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
40 /* D. Koester (2002-2007???), */
41 /* C.-J. Heine (2002-2007). */
42 /* */
43 /* */
44 /*--------------------------------------------------------------------------*/
45
46 #include "alberta_intern.h"
47 #include "alberta.h"
48
49 /*--------------------------------------------------------------------------*/
50 /* use a bit vector to indicate used/unused dofs */
51 /* storage needed: one bit per dof */
52 /*--------------------------------------------------------------------------*/
53
54 #include "dof_free_bit.h"
55
56 /* We use this in get_dof_index(). Even without assembler stuff this
57 * version is O(log(DOF_FREE_SIZE)) instead of O(DOF_FREE_SIZE) which
58 * is the effort of the naive linear-search effort.
59 */
ffs_dfu(DOF_FREE_UNIT unit)60 static inline int ffs_dfu(DOF_FREE_UNIT unit)
61 {
62 #if USE_ASSEMBLER && \
63 !HAVE_FFS_DFU && SIZEOF_DOF_FREE_UNIT == SIZEOF_LONG && __GNUC__ && __amd64
64 # define HAVE_FFS_DFU 1
65 DOF_FREE_UNIT result;
66
67 __asm__("bsfq\t%1,%0\n\t"
68 "jnz\t1f\n\t"
69 "mov\t$-1,%0\n"
70 "1:"
71 :"=r" (result)
72 :"rm" (unit));
73 return result;
74 #endif
75 #if USE_ASSEMBLER && \
76 !HAVE_FFS_DFU && SIZEOF_DOF_FREE_UNIT == SIZEOF_LONG && __GNUC__ && __i386
77 # define HAVE_FFS_DFU 1
78 DOF_FREE_UNIT result;
79
80 __asm__("bsfl\t%1,%0\n\t"
81 "jnz\t1f\n\t"
82 "mov\t$-1,%0\n"
83 "1:"
84 :"=r" (result)
85 :"rm" (unit));
86 return result;
87 #endif
88 /* You are invited to add more stuff here! */
89 #if !HAVE_FFS_DFU && SIZEOF_DOF_FREE_UNIT == SIZEOF_LONG && HAVE_FFSL
90 # define HAVE_FFS_DFU 1
91 # if !HAVE_DECL_FFSL
92 extern int ffsl(long int i);
93 # endif
94 return ffsl(unit)-1;
95 #endif
96 #if !HAVE_FFS_DFU
97 # define HAVE_FFS_DFU 1
98 int r = 0;
99
100 if (!unit) {
101 return -1;
102 }
103
104 # if SIZEOF_DOF_FREE_UNIT >= 16
105 if (!(x & 0xFFFFFFFFFFFFFFFF)) {
106 unit >>= 64;
107 r += 64;
108 }
109 # endif
110 # if SIZEOF_DOF_FREE_UNIT >= 8
111 if (!(unit & 0xFFFFFFFF)) {
112 unit >>= 32;
113 r += 32;
114 }
115 # endif
116 # if SIZEOF_DOF_FREE_UNIT >= 4
117 if (!(unit & 0xFFFF)) {
118 unit >>= 16;
119 r += 16;
120 }
121 # endif
122 # if SIZEOF_DOF_FREE_UNIT >= 2
123 if (!(unit & 0xFF)) {
124 unit >>= 8;
125 r += 8;
126 }
127 # endif
128 # if SIZEOF_DOF_FREE_UNIT >= 1
129 if (!(unit & 0xF)) {
130 unit >>= 4;
131 r += 4;
132 }
133 if (!(unit & 0x3)) {
134 unit >>= 2;
135 r += 2;
136 }
137 if (!(unit & 0x1)) {
138 unit >>= 1;
139 r += 1;
140 }
141 # endif
142 return r;
143 #endif
144 }
145
146 /*--------------------------------------------------------------------------*/
147 /* SIZE_INCREMENT: */
148 /* default increment for DOF_VECs etc. in enlarge_dof_lists */
149 /*--------------------------------------------------------------------------*/
150
151 #define SIZE_INCREMENT (DOF_FREE_SIZE * 32)
152
153 /*--------------------------------------------------------------------------*/
154
free_dof_index(DOF_ADMIN * admin,int dof)155 void free_dof_index(DOF_ADMIN *admin, int dof)
156 {
157 FUNCNAME("free_dof_index");
158 #if 0
159 int i, j;
160 int col, col2;
161 #endif
162 unsigned int iunit, ibit;
163 DOF_MATRIX *matrix;
164 MATRIX_ROW *row, *row2;
165
166 DEBUG_TEST_EXIT(admin, "no admin\n");
167 DEBUG_TEST_EXIT(admin->used_count > 0, "no dofss sin use\n");
168 DEBUG_TEST_EXIT((dof >= 0) && (dof < admin->size),
169 "invalid DOF index %d!\n",dof);
170
171 iunit = dof / DOF_FREE_SIZE;
172 ibit = dof % DOF_FREE_SIZE;
173 if ((admin->dof_free[iunit] & dof_free_bit[ibit]) != 0) {
174 /* unused DOF, probably a double free of a periodic DOF index */
175 ERROR_EXIT("Double free of DOF index.\n");
176 }
177
178 for (matrix = admin->dof_matrix; matrix; matrix = matrix->next) {
179 if (matrix->matrix_row == NULL) {
180 continue;
181 }
182 if ((row = matrix->matrix_row[dof])) {
183 do
184 {
185 #if 0
186 /* This operation is expensive - and leads to wrong results, if fe_space */
187 /* and col_fe_space don't match! (DK) */
188 /* Besides, this sort of cleanup is not done for vectors either. */
189 for (i=0; i<ROW_LENGTH; i++)
190 {
191 col = row->col[i];
192 if (ENTRY_USED(col))
193 {
194 /* remove symmetric entry if exists */
195 if (matrix->col_fe_space->admin == admin && col != dof)
196 {
197 row2 = matrix->matrix_row[col];
198 while (row2)
199 {
200 for (j=0; j<ROW_LENGTH; j++)
201 {
202 col2 = row2->col[j];
203 if (col2 == dof)
204 {
205 row2->col[j] = UNUSED_ENTRY;
206 }
207 else if (col2 == NO_MORE_ENTRIES)
208 {
209 break;
210 }
211 }
212 row2 = row2->next;
213 };
214 }
215 }
216 else if (col == NO_MORE_ENTRIES)
217 {
218 break;
219 }
220 }
221 #endif
222 row2 = row;
223 row = row->next;
224 free_matrix_row(matrix->row_fe_space, row2);
225
226 } while(row);
227
228 matrix->matrix_row[dof] = NULL;
229 }
230 }
231
232 iunit = dof / DOF_FREE_SIZE;
233 ibit = dof % DOF_FREE_SIZE;
234 admin->dof_free[iunit] |= dof_free_bit[ibit];
235 if (admin->first_hole > iunit) {
236 admin->first_hole = iunit;
237 }
238
239 admin->used_count--;
240 admin->hole_count++;
241 }
242
243 /*--------------------------------------------------------------------------*/
244
get_dof_index(DOF_ADMIN * admin)245 DOF get_dof_index(DOF_ADMIN *admin)
246 {
247 FUNCNAME("get_dof_index");
248 DOF dof = 0, i, ibit;
249
250 DEBUG_TEST_EXIT(admin, "no admin\n");
251
252 #if 1
253 if (admin->first_hole < admin->dof_free_size) {
254 ibit = ffs_dfu(admin->dof_free[admin->first_hole]);
255 DEBUG_TEST_EXIT(ibit < DOF_FREE_SIZE, "no free bit in first_hole ?\n");
256 admin->dof_free[admin->first_hole] &= ~(1UL << ibit);
257 dof = DOF_FREE_SIZE * admin->first_hole + ibit;
258 if (admin->dof_free[admin->first_hole] == 0) {
259 for (i = admin->first_hole+1; i < (int)admin->dof_free_size; i++) {
260 if (admin->dof_free[i]) {
261 break;
262 }
263 }
264 admin->first_hole = i;
265 }
266 }
267 #else
268 if (admin->first_hole < admin->dof_free_size) {
269 for (ibit = 0; ibit < DOF_FREE_SIZE; ibit++) {
270 if (admin->dof_free[admin->first_hole] & dof_free_bit[ibit]) {
271 admin->dof_free[admin->first_hole] ^= dof_free_bit[ibit];
272 dof = DOF_FREE_SIZE * admin->first_hole + ibit;
273 if (admin->dof_free[admin->first_hole] == 0) {
274 for (i = admin->first_hole+1; i < admin->dof_free_size; i++) {
275 if (admin->dof_free[i]) {
276 break;
277 }
278 }
279 admin->first_hole = i;
280 }
281 break;
282 }
283 }
284 DEBUG_TEST_EXIT(ibit < DOF_FREE_SIZE, "no free bit in first_hole ?\n");
285 }
286 #endif
287 else {
288 enlarge_dof_lists(admin, 0);
289 DEBUG_TEST_EXIT(admin->first_hole < admin->dof_free_size,
290 "no free entry after enlarge_dof_lists\n");
291 DEBUG_TEST_EXIT(admin->dof_free[admin->first_hole] & dof_free_bit[0],
292 "no free bit 0\n");
293 admin->dof_free[admin->first_hole] ^= dof_free_bit[0];
294 dof = DOF_FREE_SIZE * admin->first_hole;
295 }
296
297 admin->used_count++;
298 if (admin->hole_count > 0) admin->hole_count--;
299 admin->size_used = MAX(admin->size_used, dof+1);
300
301 return(dof);
302 }
303
304 /* Unconditionally set the number of used DOFs to used_count, and make
305 * sure that admin->dof_free[] has these DOFs allocated; the resulting
306 * DOF_ADMIN has no holes.
307 */
_AI_allocate_n_dofs(DOF_ADMIN * admin,int used_count)308 void _AI_allocate_n_dofs(DOF_ADMIN *admin, int used_count)
309 {
310 int i;
311
312 enlarge_dof_lists(admin, used_count);
313 for (i = 0; i < used_count / DOF_FREE_SIZE; i++) {
314 admin->dof_free[i] = ~DOF_UNIT_ALL_FREE;
315 }
316 admin->dof_free[i++] = DOF_UNIT_ALL_FREE << (used_count % DOF_FREE_SIZE);
317 admin->used_count = used_count;
318 admin->size_used = used_count;
319 admin->hole_count = 0;
320 admin->first_hole = used_count / DOF_FREE_SIZE;
321 }
322
323 /*--------------------------------------------------------------------------*/
324
enlarge_dof_lists(DOF_ADMIN * admin,int minsize)325 void enlarge_dof_lists(DOF_ADMIN *admin, int minsize)
326 {
327 FUNCNAME("enlarge_dof_lists");
328 int old_size, new_size, i, new_free_size;
329 DOF_INT_VEC *iv;
330 DOF_DOF_VEC *dv;
331 DOF_UCHAR_VEC *uv;
332 DOF_SCHAR_VEC *sv;
333 DOF_REAL_VEC *rv;
334 DOF_REAL_D_VEC *rdv;
335 DOF_REAL_DD_VEC *rddv;
336 DOF_PTR_VEC *pv;
337 DOF_MATRIX *mat;
338
339 DEBUG_TEST_EXIT(admin, "no admin\n");
340
341 old_size = admin->size;
342 if (minsize > 0) {
343 if (old_size > minsize) return;
344 }
345
346 new_size = MAX(minsize, admin->size + SIZE_INCREMENT);
347
348 new_size += (DOF_FREE_SIZE - (new_size % DOF_FREE_SIZE)) % DOF_FREE_SIZE;
349 admin->size = new_size;
350
351 new_free_size = new_size / DOF_FREE_SIZE;
352 admin->dof_free = MEM_REALLOC(admin->dof_free, admin->dof_free_size,
353 new_free_size, DOF_FREE_UNIT);
354 for (i = admin->dof_free_size; i < new_free_size; i++)
355 admin->dof_free[i] = DOF_UNIT_ALL_FREE;
356 admin->first_hole = admin->dof_free_size;
357 admin->dof_free_size = new_free_size;
358
359 /* enlarge all vectors and matrices */
360 /* but int_dof_vecs don't have to be changed */
361
362 for (iv = admin->dof_int_vec; iv; iv = iv->next) {
363 if (iv->size < new_size) {
364 iv->vec = MEM_REALLOC(iv->vec, iv->size, new_size, int);
365 for (i = iv->size; i < new_size; i++) iv->vec[i] = 0;
366 iv->size = new_size;
367 }
368 }
369
370 for (dv = admin->dof_dof_vec; dv; dv = dv->next) {
371 if (dv->size < new_size) {
372 dv->vec = MEM_REALLOC(dv->vec, dv->size, new_size, DOF);
373 for (i = dv->size; i < new_size; i++) dv->vec[i] = -1;
374 dv->size = new_size;
375 }
376 }
377
378 for (uv = admin->dof_uchar_vec; uv; uv = uv->next) {
379 if (uv->size < new_size) {
380 uv->vec = MEM_REALLOC(uv->vec, uv->size, new_size, U_CHAR);
381 for (i = uv->size; i < new_size; i++) uv->vec[i] = 0;
382 uv->size = new_size;
383 }
384 }
385
386 for (sv = admin->dof_schar_vec; sv; sv = sv->next) {
387 if (sv->size < new_size) {
388 sv->vec = MEM_REALLOC(sv->vec, old_size, new_size, S_CHAR);
389 for (i = sv->size; i < new_size; i++) sv->vec[i] = 0;
390 sv->size = new_size;
391 }
392 }
393
394 for (rv = admin->dof_real_vec; rv; rv = rv->next) {
395 if (rv->size < new_size) {
396 rv->vec = MEM_REALLOC(rv->vec, rv->size, new_size, REAL);
397 for (i = rv->size; i < new_size; i++) rv->vec[i] = 0.0;
398 rv->size = new_size;
399 }
400 }
401
402 for (rdv = admin->dof_real_d_vec; rdv; rdv = rdv->next) {
403 if (rdv->size < new_size) {
404 rdv->vec = MEM_REALLOC(rdv->vec, rdv->size, new_size, REAL_D);
405 for (i = rdv->size; i < new_size; i++)
406 SET_DOW(0.0, rdv->vec[i]);
407 rdv->size = new_size;
408 }
409 }
410
411 for (rddv = admin->dof_real_dd_vec; rddv; rddv = rddv->next) {
412 if (rddv->size < new_size) {
413 rddv->vec = MEM_REALLOC(rddv->vec, rddv->size, new_size, REAL_DD);
414 for (i = rddv->size; i < new_size; i++)
415 MSET_DOW(0.0, rddv->vec[i]);
416 rddv->size = new_size;
417 }
418 }
419
420 for (pv = admin->dof_ptr_vec; pv; pv = pv->next) {
421 if (pv->size < new_size) {
422 pv->vec = MEM_REALLOC(pv->vec, pv->size, new_size, void *);
423 for (i = pv->size; i < new_size; i++) pv->vec[i] = NULL;
424 pv->size = new_size;
425 }
426 }
427
428 for (mat = admin->dof_matrix; mat; mat = mat->next) {
429 if (!mat->is_diagonal) {
430 if (mat->size < new_size) {
431 mat->matrix_row = MEM_REALLOC(mat->matrix_row, mat->size, new_size,
432 MATRIX_ROW *);
433 for (i = mat->size; i < new_size; i++) {
434 mat->matrix_row[i] = NULL;
435 }
436 mat->size = new_size;
437 }
438 } else {
439 mat->size = new_size;
440 }
441 }
442
443 }
444
445
446 /*--------------------------------------------------------------------------*/
447 /* dof_compress: remove holes in dof vectors */
448 /*--------------------------------------------------------------------------*/
449
450 typedef struct dof_admin_traverse_data1 {
451 int *new_dof;
452 int *g_n_dof;
453 int *g_n0_dof;
454 int *g_node;
455 } DOF_ADMIN_TRAVERSE_DATA1;
456
457 /*--------------------------------------------------------------------------*/
458 /* ATTENTION: */
459 /* new_dof_fct() destroys new_dof !!!!!!!!!! */
460 /* should be used only at the end of dof_compress()!!!!! */
461 /*--------------------------------------------------------------------------*/
462
463 /* CHANGE_DOFS_1 changes old dofs to NEGATIVE new dofs. The index is shifted */
464 /* by - 2 so that unused element DOFs (-1) are preserved. */
465
466 #define CHANGE_DOFS_1(el) \
467 dof0 = el->dof[n0+i]; \
468 dof = dof0 + nd0; \
469 if(dof0) \
470 for (j = 0; j < nd; j++) \
471 if ((k = dof[j]) >= 0) \
472 /* do it only once! (dofs are visited more than once) */ \
473 dof[j] = - ud->new_dof[k] - 2;
474
475 /* CHANGE_DOFS_2 changes NEGATIVE new dofs to POSITIVE. The index shift above*/
476 /* is undone. */
477
478 #define CHANGE_DOFS_2(el) \
479 dof0 = el->dof[n0+i]; \
480 dof = dof0 + nd0; \
481 if(dof0) \
482 for (j = 0; j < nd; j++) \
483 if ((k = dof[j]) < -1) \
484 /* do it only once! (dofs are visited more than once) */ \
485 dof[j] = - k - 2;
486
new_dof_fct1(const EL_INFO * el_info,void * data)487 static inline void new_dof_fct1(const EL_INFO *el_info, void *data)
488 {
489 DOF_ADMIN_TRAVERSE_DATA1 *ud = (DOF_ADMIN_TRAVERSE_DATA1 *)data;
490 EL *el = el_info->el;
491 int i, j, k, n0, nd, nd0;
492 int dim = el_info->mesh->dim;
493 DOF *dof0, *dof;
494
495 if ((nd = ud->g_n_dof[VERTEX])) {
496 nd0 = ud->g_n0_dof[VERTEX];
497 n0 = ud->g_node[VERTEX];
498 for (i = 0; i < N_VERTICES(dim); i++) {
499 CHANGE_DOFS_1(el);
500 }
501 }
502
503 if(dim > 1 && (nd = ud->g_n_dof[EDGE])) {
504 nd0 = ud->g_n0_dof[EDGE];
505 n0 = ud->g_node[EDGE];
506 for (i = 0; i < N_EDGES(dim); i++) {
507 CHANGE_DOFS_1(el);
508 }
509 }
510
511 if(dim == 3 && (nd = ud->g_n_dof[FACE])) {
512 nd0 = ud->g_n0_dof[FACE];
513 n0 = ud->g_node[FACE];
514 for (i = 0; i < N_FACES_3D; i++) {
515 CHANGE_DOFS_1(el);
516 }
517 }
518
519 if ((nd = ud->g_n_dof[CENTER])) {
520 nd0 = ud->g_n0_dof[CENTER];
521 n0 = ud->g_node[CENTER];
522 i = 0; /* only one center */
523 CHANGE_DOFS_1(el);
524 }
525
526 return;
527 }
528
529
new_dof_fct2(const EL_INFO * el_info,void * data)530 static inline void new_dof_fct2(const EL_INFO *el_info, void *data)
531 {
532 DOF_ADMIN_TRAVERSE_DATA1 *ud = (DOF_ADMIN_TRAVERSE_DATA1 *)data;
533 EL *el = el_info->el;
534 int i, j, k, n0, nd, nd0;
535 int dim = el_info->mesh->dim;
536 DOF *dof0, *dof;
537
538 if ((nd = ud->g_n_dof[VERTEX])) {
539 nd0 = ud->g_n0_dof[VERTEX];
540 n0 = ud->g_node[VERTEX];
541 for (i = 0; i < N_VERTICES(dim); i++) {
542 CHANGE_DOFS_2(el);
543 }
544 }
545
546 if(dim > 1 && (nd = ud->g_n_dof[EDGE])) {
547 nd0 = ud->g_n0_dof[EDGE];
548 n0 = ud->g_node[EDGE];
549 for (i = 0; i < N_EDGES(dim); i++)
550 {
551 CHANGE_DOFS_2(el);
552 }
553 }
554
555 if(dim == 3 && (nd = ud->g_n_dof[FACE])) {
556 nd0 = ud->g_n0_dof[FACE];
557 n0 = ud->g_node[FACE];
558 for (i = 0; i < N_FACES_3D; i++) {
559 CHANGE_DOFS_2(el);
560 }
561 }
562
563 if ((nd = ud->g_n_dof[CENTER])) {
564 nd0 = ud->g_n0_dof[CENTER];
565 n0 = ud->g_node[CENTER];
566 i = 0; /* only one center */
567 CHANGE_DOFS_2(el);
568 }
569
570 return;
571 }
572
573 #undef CHANGE_DOFS_1
574 #undef CHANGE_DOFS_2
575
576 /*--------------------------------------------------------------------------*/
577
dof_compress(MESH * mesh)578 void dof_compress(MESH *mesh)
579 {
580 FUNCNAME("dof_compress");
581 DOF_ADMIN_TRAVERSE_DATA1 td[1] = {{0}};
582 DOF_ADMIN *compress_admin, *search_admin;
583 int i, j, k, n_rows, size, first, last=0, col, iadmin, jadmin;
584 FLAGS fill_flag;
585 DOF_INT_VEC *div;
586 DOF_DOF_VEC *ddv;
587 DOF_UCHAR_VEC *duv;
588 DOF_SCHAR_VEC *dsv;
589 DOF_REAL_D_VEC *drdv;
590 DOF_REAL_DD_VEC *drddv;
591 DOF_REAL_VEC *drv;
592 DOF_PTR_VEC *dpv;
593 DOF_MATRIX *dm;
594 MATRIX_ROW *row, *row_next;
595 DBL_LIST_NODE *pos;
596 MESH_MEM_INFO *minfo;
597 bool did_compress = false;
598
599 TEST_EXIT(mesh, "No mesh given!\n");
600
601 td->g_node = mesh->node;
602
603 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++)
604 {
605 compress_admin = mesh->dof_admin[iadmin];
606
607 DEBUG_TEST_EXIT(compress_admin, "no dof_admin[%d] in mesh\n", iadmin);
608
609 if ((size = compress_admin->size) < 1) continue;
610 if (compress_admin->used_count < 1) continue;
611 if (compress_admin->hole_count < 1) continue;
612
613 did_compress = true;
614
615 td->g_n_dof = compress_admin->n_dof;
616 td->g_n0_dof = compress_admin->n0_dof;
617
618 td->new_dof = MEM_ALLOC(size, int);
619 for (i = 0; i < size; i++) { /* mark everything unused */
620 td->new_dof[i] = -1;
621 }
622
623 FOR_ALL_DOFS(compress_admin, td->new_dof[dof] = 1);
624
625 n_rows = 0;
626 for (i = 0; i < size; i++) { /* create a MONOTONE compress */
627 if (td->new_dof[i] == 1) {
628 td->new_dof[i] = n_rows++;
629 last = i;
630 }
631 }
632
633 DEBUG_TEST_EXIT(n_rows == compress_admin->used_count,
634 "count %d != used_count %d\n",
635 n_rows, compress_admin->used_count);
636
637 {
638 int n1 = n_rows / DOF_FREE_SIZE;
639
640 for (i = 0; i < n1; i++)
641 compress_admin->dof_free[i] = 0;
642 if (n1 < (int)compress_admin->dof_free_size)
643 compress_admin->dof_free[n1] = (~0UL) << (n_rows % DOF_FREE_SIZE);
644 for (i = n1+1; i < (int)compress_admin->dof_free_size; i++)
645 compress_admin->dof_free[i] = DOF_UNIT_ALL_FREE;
646 compress_admin->first_hole = n1;
647 }
648
649 compress_admin->hole_count = 0;
650 compress_admin->size_used = n_rows;
651
652 first = 0;
653 for (i=0; i<size; i++) {
654 if ((td->new_dof[i] < i) && (td->new_dof[i] >= 0)) {
655 first = i;
656 break;
657 }
658 }
659 if (last >= first) {
660 last++; /* for (i = first; i < last; i++) ... */
661
662 /* compress all vectors associated to admin */
663
664 for (div = compress_admin->dof_int_vec; div; div = div->next) {
665 for (i=first; i<last; i++) {
666 if ((k = td->new_dof[i]) >= 0)
667 div->vec[k] = div->vec[i];
668 }
669 }
670
671 for (ddv = compress_admin->dof_dof_vec; ddv; ddv = ddv->next) {
672 for (i=first; i<last; i++) {
673 if ((k = td->new_dof[i]) >= 0)
674 ddv->vec[k] = ddv->vec[i];
675 }
676 for (i=0; i < n_rows; i++) {
677 if ((k = td->new_dof[ddv->vec[i]]) >= 0)
678 ddv->vec[i] = k;
679 }
680 }
681
682 for (ddv = compress_admin->int_dof_vec; ddv; ddv = ddv->next) {
683 for (i=0; i < ddv->size; i++) {
684 if (ddv->vec[i] >= 0)
685 if ((k = td->new_dof[ddv->vec[i]]) >= 0)
686 ddv->vec[i] = k;
687 }
688 }
689
690 for (duv = compress_admin->dof_uchar_vec; duv; duv = duv->next) {
691 for (i=first; i<last; i++) {
692 if ((k = td->new_dof[i]) >= 0)
693 duv->vec[k] = duv->vec[i];
694 }
695 }
696
697 for (dsv = compress_admin->dof_schar_vec; dsv; dsv = dsv->next) {
698 for (i=first; i<last; i++) {
699 if ((k = td->new_dof[i]) >= 0)
700 dsv->vec[k] = dsv->vec[i];
701 }
702 }
703
704 for (drv = compress_admin->dof_real_vec; drv; drv = drv->next) {
705 for (i=first; i<last; i++) {
706 if ((k = td->new_dof[i]) >= 0)
707 drv->vec[k] = drv->vec[i];
708 }
709 }
710
711 for (drdv = compress_admin->dof_real_d_vec; drdv; drdv = drdv->next) {
712 for (i=first; i<last; i++) {
713 if ((k = td->new_dof[i]) >= 0)
714 COPY_DOW(drdv->vec[i], drdv->vec[k]);
715 }
716 }
717
718 for (drddv = compress_admin->dof_real_dd_vec;
719 drddv; drddv = drddv->next) {
720 for (i = first; i < last; i++) {
721 if ((k = td->new_dof[i]) >= 0)
722 MCOPY_DOW((const REAL_D *)drddv->vec[i], drddv->vec[k]);
723 }
724 }
725
726 for (dpv = compress_admin->dof_ptr_vec; dpv; dpv = dpv->next) {
727 for (i=first; i<last; i++) {
728 if ((k = td->new_dof[i]) >= 0)
729 dpv->vec[k] = dpv->vec[i];
730 }
731 }
732
733 /* Row corrections */
734 for (dm = compress_admin->dof_matrix; dm; dm = dm->next) {
735 if (dm->is_diagonal) {
736 /* probably a diagonal matrix, the row/column corrections
737 * already have been performed on the vector holding the
738 * diagonal entries.
739 */
740 continue;
741 }
742 for (i = first; i < last; i++) {
743 if ((k = td->new_dof[i]) >= 0) {
744 /* free dm->matrix_row[k]; */
745 for (row = dm->matrix_row[k]; row; row = row_next) {
746 /* Should this really happen ???? */
747 row_next = row->next;
748 free_matrix_row(dm->row_fe_space, row);
749 }
750 dm->matrix_row[k] = dm->matrix_row[i];
751 dm->matrix_row[i] = NULL;
752 }
753 }
754 }
755
756 /* Column corrections, we search all matrices for ones with */
757 /* col_fe_space->admin == compress_admin. */
758
759 for (jadmin = 0; jadmin < mesh->n_dof_admin; jadmin++) {
760 search_admin = mesh->dof_admin[jadmin];
761
762 DEBUG_TEST_EXIT(search_admin,
763 "no dof_admin[%d] in mesh\n", jadmin);
764
765 for (dm = search_admin->dof_matrix; dm; dm = dm->next) {
766 if ((dm->col_fe_space == NULL &&
767 dm->row_fe_space->admin == compress_admin) ||
768 (dm->col_fe_space != NULL &&
769 dm->col_fe_space->admin == compress_admin)) {
770 int n_rows = dm->row_fe_space->admin->size_used;
771 if (dm->is_diagonal) {
772 /* probably a diagonal matrix, correct the column indices */
773 for (i=0; i < n_rows; i++) { /* change columns */
774 DOF col_dof = dm->diag_cols->vec[i];
775 if (ENTRY_USED(col_dof)) {
776 dm->diag_cols->vec[i] = td->new_dof[col_dof];
777 }
778 }
779 } else {
780 for (i=0; i < n_rows; i++) { /* change columns */
781 for (row = dm->matrix_row[i]; row; row = row->next) {
782 for (j=0; j<ROW_LENGTH; j++) {
783 col = row->col[j];
784 if (ENTRY_USED(col)) row->col[j] = td->new_dof[col];
785 }
786 }
787 }
788 }
789 }
790 }
791 }
792
793
794 /* Call custom compress hooks */
795 for (pos = compress_admin->compress_hooks.next;
796 pos != &compress_admin->compress_hooks;
797 pos = pos->next) {
798 DOF_COMP_HOOK *hook = dbl_list_entry(pos, DOF_COMP_HOOK, node);
799 hook->handler(first, last, td->new_dof, hook->application_data);
800 }
801
802 /* now, change dofs in MESH's nodes */
803 /* new_dof_fct destroys new_dof[] !!! */
804
805 fill_flag = CALL_EVERY_EL_PREORDER | FILL_NOTHING;
806
807 TRAVERSE_FIRST(mesh, -1, fill_flag) {
808 new_dof_fct1(el_info, td);
809 } TRAVERSE_NEXT();
810 TRAVERSE_FIRST(mesh, -1, fill_flag) {
811 new_dof_fct2(el_info, td);
812 } TRAVERSE_NEXT();
813 }
814
815 MEM_FREE(td->new_dof, size, int);
816 td->new_dof = NULL;
817 }
818
819 mesh->cookie += did_compress;
820
821 minfo = (MESH_MEM_INFO *)mesh->mem_info;
822 for (i = 0; i < minfo->n_slaves; i++) {
823 dof_compress(minfo->slaves[i]);
824 }
825 }
826
add_dof_compress_hook(const DOF_ADMIN * admin,DOF_COMP_HOOK * hook)827 void add_dof_compress_hook(const DOF_ADMIN *admin, DOF_COMP_HOOK *hook)
828 {
829 dbl_list_add_head((DBL_LIST_NODE *)&admin->compress_hooks, &hook->node);
830 }
831
del_dof_compress_hook(DOF_COMP_HOOK * hook)832 void del_dof_compress_hook(DOF_COMP_HOOK *hook)
833 {
834 dbl_list_del(&hook->node);
835 }
836
837 /*--------------------------------------------------------------------------*/
838 /* matrix administration */
839 /*--------------------------------------------------------------------------*/
840 /* AI_add_dof_matrix_to_admin(): enter new matrix to DOF_ADMIN list */
841 /* (re)alloc matrix_row, el_* if necessary*/
842 /* AI_add_dof_int_vec_to_admin(): enter new vec to DOF_ADMIN list */
843 /* AI_add_int_dof_vec_to_admin(): enter new vec to DOF_ADMIN list */
844 /* AI_add_dof_dof_vec_to_admin(): enter new vec to DOF_ADMIN list */
845 /* AI_add_dof_uchar_vec_to_admin(): enter new vec to DOF_ADMIN list */
846 /* AI_add_dof_schar_vec_to_admin(): enter new vec to DOF_ADMIN list */
847 /* AI_add_dof_real_vec_to_admin(): enter new vec to DOF_ADMIN list */
848 /* AI_add_dof_real_d_vec_to_admin(): enter new vec to DOF_ADMIN list */
849 /* AI_add_dof_real_dd_vec_to_admin(): enter new vec to DOF_ADMIN list */
850 /* AI_add_dof_ptr_vec_to_admin(): enter new vec to DOF_ADMIN list */
851 /* (re)alloc vector if necessary */
852 /* */
853 /* update_dof_matrix(): add/substract el_mat to/from matrix */
854 /*--------------------------------------------------------------------------*/
855
856 #define CHECK_IF_PRESENT(admin, TYPE, list, vec) \
857 { \
858 TYPE *v; \
859 for (v=admin->list; v; v = v->next) \
860 if (v==vec) { \
861 ERROR_EXIT("dof_vec %s already associated to admin %s\n", \
862 NAME(vec), NAME(admin)); \
863 return; \
864 } \
865 }
866
867
868 #define DEFUN_ADD_DOF_OBJ_TO_ADMIN(TYPE, list, enlarge) \
869 void add_##list##_to_admin(TYPE *obj, DOF_ADMIN *admin) \
870 { \
871 FUNCNAME("add_"#list"_to_admin"); \
872 \
873 if (!obj) { \
874 MSG("no obj\n"); \
875 return; \
876 } \
877 \
878 CHECK_IF_PRESENT(admin, TYPE, list, obj); \
879 \
880 if (obj->size < admin->size) { \
881 enlarge; \
882 obj->size = admin->size; \
883 } \
884 \
885 obj->next = admin->list; \
886 admin->list = obj; \
887 /* obj->admin = admin; */ \
888 } \
889 struct _AI_semicolon_dummy
890
891 #define DEFUN_ADD_DOF_VEC_TO_ADMIN(TYPE, type, list) \
892 DEFUN_ADD_DOF_OBJ_TO_ADMIN(TYPE, list, { \
893 obj->vec = MEM_REALLOC(obj->vec, obj->size, admin->size, type); \
894 })
895
896 static void refine_diag_cols(DOF_INT_VEC *vec, RC_LIST_EL *rclist, int n);
897
898 #define DEFUN_ADD_DOF_MAT_TO_ADMIN(TYPE, type, list) \
899 DEFUN_ADD_DOF_OBJ_TO_ADMIN(TYPE, list, { \
900 if (!obj->is_diagonal) { \
901 int i; \
902 \
903 obj->matrix_row = MEM_REALLOC(obj->matrix_row, obj->size, \
904 admin->size, type *); \
905 for (i = obj->size; i < admin->size; i++) { \
906 obj->matrix_row[i] = NULL; \
907 } \
908 } else { \
909 obj->diag_cols = get_dof_int_vec( \
910 "diag cols", obj->row_fe_space->unchained); \
911 obj->diag_cols->refine_interpol = refine_diag_cols; \
912 FOR_ALL_DOFS(admin, obj->diag_cols->vec[dof] = UNUSED_ENTRY); \
913 } \
914 obj->size = admin->size; \
915 })
916
917 #define DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(TYPE, list) \
918 void remove_##list##_from_admin(TYPE *obj) \
919 { \
920 FUNCNAME("remove_"#list"_from_admin"); \
921 DOF_ADMIN *admin; \
922 TYPE *obj_last; \
923 \
924 if (obj->fe_space && (admin = (DOF_ADMIN *)obj->fe_space->admin)) { \
925 if (admin->list == obj) { \
926 admin->list = obj->next; \
927 } else { \
928 obj_last = admin->list; \
929 while (obj_last && obj_last->next != obj) \
930 obj_last = obj_last->next; \
931 \
932 if (!obj_last) { \
933 ERROR_EXIT(""#list" %s not in list of dof admin %s found\n", \
934 NAME(obj), NAME(admin)); \
935 } else { \
936 obj_last->next = obj->next; \
937 } \
938 } \
939 } \
940 } \
941 struct _AI_semicolon_dummy
942
943 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_INT_VEC, int, dof_int_vec);
944 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_INT_VEC, dof_int_vec);
945
946 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_DOF_VEC, DOF, dof_dof_vec);
947 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_DOF_VEC, dof_dof_vec);
948
949 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_DOF_VEC, DOF, int_dof_vec);
950 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_DOF_VEC, int_dof_vec);
951
952 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_UCHAR_VEC, U_CHAR, dof_uchar_vec);
953 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_UCHAR_VEC, dof_uchar_vec);
954
955 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_SCHAR_VEC, S_CHAR, dof_schar_vec);
956 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_SCHAR_VEC, dof_schar_vec);
957
958 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_REAL_VEC, REAL, dof_real_vec);
959 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_REAL_VEC, dof_real_vec);
960
961 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_REAL_D_VEC, REAL_D, dof_real_d_vec);
962 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_REAL_D_VEC, dof_real_d_vec);
963
964 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_REAL_DD_VEC, REAL_DD, dof_real_dd_vec);
965 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_REAL_DD_VEC, dof_real_dd_vec);
966
967 DEFUN_ADD_DOF_VEC_TO_ADMIN(DOF_PTR_VEC, void *, dof_ptr_vec);
968 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_PTR_VEC, dof_ptr_vec);
969
970 #define fe_space row_fe_space
971 DEFUN_ADD_DOF_MAT_TO_ADMIN(DOF_MATRIX, MATRIX_ROW, dof_matrix);
972 DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN(DOF_MATRIX, dof_matrix);
973 #undef fe_space
974
975 #undef CHECK_IF_PRESENT
976 #undef DEFUN_ADD_DOF_OBJ_TO_ADMIN
977 #undef DEFUN_ADD_DOF_VEC_TO_ADMIN
978 #undef DEFUN_ADD_DOF_MAT_TO_ADMIN
979 #undef DEFUN_REMOVE_DOF_OBJ_FROM_ADMIN
980
981 /*--------------------------------------------------------------------------*/
982 /* some BLAS level 1 routines: */
983 /* --------------------------- */
984 /* dof_asum: asum = ||X||_l1 */
985 /* dof_nrm2: nrm2 = ||X||_l2 */
986 /* dof_set: X = alpha */
987 /* dof_scal: X = alpha * X */
988 /* dof_copy: Y = X */
989 /* dof_dot: dot = X . Y */
990 /* dof_axpy: Y = Y + alpha * X */
991 /* */
992 /* some BLAS level 2 routines: */
993 /* --------------------------- */
994 /* dof_gemv: y = alpha*A*x + beta*y or y = alpha*A'*x + beta*y */
995 /* */
996 /* some non-BLAS level 2 routines: */
997 /* --------------------------- */
998 /* dof_mv: y = A*x or y = A'*x */
999 /* dof_min: min of x */
1000 /* dof_max: max of x */
1001 /*--------------------------------------------------------------------------*/
1002
1003 static inline
__dof_nrm2(const DOF_REAL_VEC * x)1004 REAL __dof_nrm2(const DOF_REAL_VEC *x)
1005 {
1006 FUNCNAME("dof_nrm2");
1007 REAL nrm;
1008 const DOF_ADMIN *admin = NULL;
1009
1010 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1011 "pointer is NULL: %p, %p\n", x, admin);
1012 TEST_EXIT(x->size >= admin->size_used,
1013 "x->size = %d too small: admin->size_used = %d\n", x->size,
1014 admin->size_used);
1015
1016 nrm = 0.0;
1017 FOR_ALL_DOFS(admin, nrm += x->vec[dof] * x->vec[dof] );
1018
1019 return nrm;
1020 }
1021
1022 static inline
__dof_asum(const DOF_REAL_VEC * x)1023 REAL __dof_asum(const DOF_REAL_VEC *x)
1024 {
1025 FUNCNAME("dof_asum");
1026 REAL nrm;
1027 const DOF_ADMIN *admin = NULL;
1028
1029 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1030 "pointer is NULL: %p, %p\n", x, admin);
1031 TEST_EXIT(x->size >= admin->size_used,
1032 "x->size = %d too small: admin->size_used = %d\n", x->size,
1033 admin->size_used);
1034
1035 nrm = 0.0;
1036 FOR_ALL_DOFS(admin, nrm += ABS(x->vec[dof]) );
1037
1038 return(nrm);
1039 }
1040
1041 static inline
__dof_set(REAL alpha,DOF_REAL_VEC * x)1042 void __dof_set(REAL alpha, DOF_REAL_VEC *x)
1043 {
1044 FUNCNAME("dof_set");
1045 const DOF_ADMIN *admin = NULL;
1046
1047 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1048 "pointer is NULL: %p, %p\n", x, admin);
1049 TEST_EXIT(x->size >= admin->size_used,
1050 "x->size = %d too small: admin->size_used = %d\n", x->size,
1051 admin->size_used);
1052
1053 FOR_ALL_DOFS(admin, x->vec[dof] = alpha );
1054 }
1055
1056 static inline
__dof_scal(REAL alpha,DOF_REAL_VEC * x)1057 void __dof_scal(REAL alpha, DOF_REAL_VEC *x)
1058 {
1059 FUNCNAME("dof_scal");
1060 const DOF_ADMIN *admin = NULL;
1061
1062 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1063 "pointer is NULL: %p, %p\n", x, admin);
1064 TEST_EXIT(x->size >= admin->size_used,
1065 "x->size = %d too small: admin->size_used = %d\n", x->size,
1066 admin->size_used);
1067
1068 FOR_ALL_DOFS(admin, x->vec[dof] *= alpha);
1069 }
1070
1071 static inline
__dof_dot(const DOF_REAL_VEC * x,const DOF_REAL_VEC * y)1072 REAL __dof_dot(const DOF_REAL_VEC *x, const DOF_REAL_VEC *y)
1073 {
1074 FUNCNAME("dof_dot");
1075 REAL dot;
1076 const DOF_ADMIN *admin = NULL;
1077
1078 TEST_EXIT(x && y, "pointer is NULL: %p, %p\n", x,y);
1079 TEST_EXIT(x->fe_space && y->fe_space,
1080 "fe_space is NULL: %p, %p\n", x->fe_space,y->fe_space);
1081 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1082 "no admin or different admins: %p, %p\n",
1083 x->fe_space->admin, y->fe_space->admin);
1084 TEST_EXIT(x->size >= admin->size_used,
1085 "x->size = %d too small: admin->size_used = %d\n", x->size,
1086 admin->size_used);
1087 TEST_EXIT(y->size >= admin->size_used,
1088 "y->size = %d too small: admin->size_used = %d\n", y->size,
1089 admin->size_used);
1090
1091 dot = 0.0;
1092 FOR_ALL_DOFS(admin, dot += x->vec[dof] * y->vec[dof]);
1093
1094 return dot;
1095 }
1096
1097 static inline
__dof_copy(const DOF_REAL_VEC * x,DOF_REAL_VEC * y)1098 void __dof_copy(const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
1099 {
1100 FUNCNAME("dof_copy");
1101 REAL *xvec, *yvec;
1102 const DOF_ADMIN *admin = NULL;
1103
1104 TEST_EXIT(x && y, "pointer is NULL: %p, %p\n", x,y);
1105 TEST_EXIT(x->fe_space && y->fe_space,
1106 "fe_space is NULL: %p, %p\n", x->fe_space,y->fe_space);
1107 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1108 "no admin or different admins: %p, %p\n",
1109 x->fe_space->admin, y->fe_space->admin);
1110 TEST_EXIT(x->size >= admin->size_used,
1111 "x->size = %d too small: admin->size_used = %d\n", x->size,
1112 admin->size_used);
1113 TEST_EXIT(y->size >= admin->size_used,
1114 "y->size = %d too small: admin->size_used = %d\n", y->size,
1115 admin->size_used);
1116 xvec = x->vec;
1117 yvec = y->vec;
1118
1119 FOR_ALL_DOFS(admin, yvec[dof] = xvec[dof]);
1120 }
1121
1122 static inline
__dof_axpy(REAL alpha,const DOF_REAL_VEC * x,DOF_REAL_VEC * y)1123 void __dof_axpy(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
1124 {
1125 FUNCNAME("dof_axpy");
1126 REAL *xvec, *yvec;
1127 const DOF_ADMIN *admin;
1128
1129 TEST_EXIT(x && y, "pointer is NULL: %p, %p\n", x,y);
1130 TEST_EXIT(x->fe_space && y->fe_space,
1131 "fe_space is NULL: %p, %p\n", x->fe_space,y->fe_space);
1132 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1133 "no admin or different admins: %p, %p\n",
1134 x->fe_space->admin, y->fe_space->admin);
1135 TEST_EXIT(x->size >= admin->size_used,
1136 "x->size = %d too small: admin->size = %d\n",
1137 x->size, admin->size_used);
1138 TEST_EXIT(y->size >= admin->size_used,
1139 "y->size = %d too small: admin->size = %d\n",
1140 y->size, admin->size_used);
1141
1142 xvec = x->vec;
1143 yvec = y->vec;
1144
1145 FOR_ALL_DOFS(admin, yvec[dof] += alpha * xvec[dof]);
1146 }
1147
1148 static inline
__dof_xpay(REAL alpha,const DOF_REAL_VEC * x,DOF_REAL_VEC * y)1149 void __dof_xpay(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
1150 {
1151 FUNCNAME("dof_axpy");
1152 REAL *xvec, *yvec;
1153 const DOF_ADMIN *admin;
1154
1155 TEST_EXIT(x && y, "pointer is NULL: %p, %p\n", x,y);
1156 TEST_EXIT(x->fe_space && y->fe_space,
1157 "fe_space is NULL: %p, %p\n", x->fe_space,y->fe_space);
1158 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1159 "no admin or different admins: %p, %p\n",
1160 x->fe_space->admin, y->fe_space->admin);
1161 TEST_EXIT(x->size >= admin->size_used,
1162 "x->size = %d too small: admin->size_used = %d\n",
1163 x->size, admin->size_used);
1164 TEST_EXIT(y->size >= admin->size_used,
1165 "y->size = %d too small: admin->size_used = %d\n",
1166 y->size, admin->size_used);
1167
1168 xvec = x->vec;
1169 yvec = y->vec;
1170
1171 FOR_ALL_DOFS(admin, yvec[dof] = xvec[dof] + alpha*yvec[dof]);
1172 }
1173
1174 static inline
__dof_min(const DOF_REAL_VEC * x)1175 REAL __dof_min(const DOF_REAL_VEC *x)
1176 {
1177 FUNCNAME("dof_min");
1178 REAL m;
1179 const DOF_ADMIN *admin = NULL;
1180
1181 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1182 "pointer is NULL: %p, %p\n", x, admin);
1183 TEST_EXIT(x->size >= admin->size_used,
1184 "x->size = %d too small: admin->size_used = %d\n", x->size,
1185 admin->size_used);
1186
1187 m = REAL_MAX;
1188 FOR_ALL_DOFS(admin, m = MIN(m, x->vec[dof]));
1189
1190 return m;
1191 }
1192
1193 static inline
__dof_max(const DOF_REAL_VEC * x)1194 REAL __dof_max(const DOF_REAL_VEC *x)
1195 {
1196 FUNCNAME("dof_max");
1197 REAL m;
1198 const DOF_ADMIN *admin = NULL;
1199
1200 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1201 "pointer is NULL: %p, %p\n", x, admin);
1202 TEST_EXIT(x->size >= admin->size_used,
1203 "x->size = %d too small: admin->size_used = %d\n", x->size,
1204 admin->size_used);
1205
1206 m = REAL_MIN;
1207 FOR_ALL_DOFS(admin, m = MAX(m, x->vec[dof]));
1208
1209 return m;
1210 }
1211
dof_nrm2(const DOF_REAL_VEC * x)1212 REAL dof_nrm2(const DOF_REAL_VEC *x)
1213 {
1214 REAL nrm2 = 0.0;
1215
1216 CHAIN_DO(x, const DOF_REAL_VEC) {
1217 nrm2 += __dof_nrm2(x);
1218 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1219
1220 return sqrt(nrm2);
1221 }
1222
dof_asum(const DOF_REAL_VEC * x)1223 REAL dof_asum(const DOF_REAL_VEC *x)
1224 {
1225 REAL asum = 0.0;
1226
1227 CHAIN_DO(x, const DOF_REAL_VEC) {
1228 asum += __dof_asum(x);
1229 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1230
1231 return asum;
1232 }
1233
dof_set(REAL alpha,DOF_REAL_VEC * x)1234 void dof_set(REAL alpha, DOF_REAL_VEC *x)
1235 {
1236 CHAIN_DO(x, DOF_REAL_VEC) {
1237 __dof_set(alpha, x);
1238 } CHAIN_WHILE(x, DOF_REAL_VEC);
1239 }
1240
dof_scal(REAL alpha,DOF_REAL_VEC * x)1241 void dof_scal(REAL alpha, DOF_REAL_VEC *x)
1242 {
1243 CHAIN_DO(x, DOF_REAL_VEC) {
1244 __dof_scal(alpha, x);
1245 } CHAIN_WHILE(x, DOF_REAL_VEC);
1246 }
1247
dof_dot(const DOF_REAL_VEC * x,const DOF_REAL_VEC * y)1248 REAL dof_dot(const DOF_REAL_VEC *x, const DOF_REAL_VEC *y)
1249 {
1250 REAL dot = 0.0;
1251
1252 CHAIN_DO(x, const DOF_REAL_VEC) {
1253 dot += __dof_dot(x, y);
1254 y = CHAIN_NEXT(y, const DOF_REAL_VEC);
1255 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1256
1257 return dot;
1258 }
1259
dof_copy(const DOF_REAL_VEC * x,DOF_REAL_VEC * y)1260 void dof_copy(const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
1261 {
1262 CHAIN_DO(x, const DOF_REAL_VEC) {
1263 __dof_copy(x, y);
1264 y = CHAIN_NEXT(y, DOF_REAL_VEC);
1265 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1266 }
1267
dof_axpy(REAL alpha,const DOF_REAL_VEC * x,DOF_REAL_VEC * y)1268 void dof_axpy(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
1269 {
1270 CHAIN_DO(x, const DOF_REAL_VEC) {
1271 __dof_axpy(alpha, x, y);
1272 y = CHAIN_NEXT(y, DOF_REAL_VEC);
1273 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1274 }
1275
dof_xpay(REAL alpha,const DOF_REAL_VEC * x,DOF_REAL_VEC * y)1276 void dof_xpay(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
1277 {
1278 CHAIN_DO(x, const DOF_REAL_VEC) {
1279 __dof_xpay(alpha, x, y);
1280 y = CHAIN_NEXT(y, DOF_REAL_VEC);
1281 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1282 }
1283
dof_min(const DOF_REAL_VEC * x)1284 REAL dof_min(const DOF_REAL_VEC *x)
1285 {
1286 REAL m = REAL_MAX, m2;
1287 CHAIN_DO(x, const DOF_REAL_VEC) {
1288 m2 = __dof_min(x);
1289 m = MIN(m, m2);
1290 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1291 return m;
1292 }
1293
dof_max(const DOF_REAL_VEC * x)1294 REAL dof_max(const DOF_REAL_VEC *x)
1295 {
1296 REAL m = REAL_MIN, m2;
1297 CHAIN_DO(x, const DOF_REAL_VEC) {
1298 m2 = __dof_max(x);
1299 m = MAX(m, m2);
1300 } CHAIN_WHILE(x, const DOF_REAL_VEC);
1301
1302 return m;
1303 }
1304
1305 /*--------------------------------------------------------------------------*/
1306 /* now the same routines for REAL_D vectors */
1307 /*--------------------------------------------------------------------------*/
1308
1309 static inline
__dof_set_d(REAL alpha,DOF_REAL_D_VEC * x)1310 void __dof_set_d(REAL alpha, DOF_REAL_D_VEC *x)
1311 {
1312 FUNCNAME("dof_set_d");
1313 const DOF_ADMIN *admin = NULL;
1314
1315 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1316 "pointer is NULL: x: %p, x->fe_space: %p, x->fe_space->admin :%p\n",
1317 x, x->fe_space, admin);
1318 TEST_EXIT(x->size >= admin->size_used,
1319 "x->size = %d too small: admin->size_used = %d\n",
1320 x->size, admin->size_used);
1321
1322 FOR_ALL_DOFS(admin, SET_DOW(alpha, x->vec[dof]));
1323 }
1324
1325 static inline
__dof_scal_d(REAL alpha,DOF_REAL_D_VEC * x)1326 void __dof_scal_d(REAL alpha, DOF_REAL_D_VEC *x)
1327 {
1328 FUNCNAME("dof_scal_d");
1329 const DOF_ADMIN *admin = NULL;
1330
1331 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1332 "pointer is NULL: x: %p, x->fe_space: %p, x->fe_space->admin :%p\n",
1333 x, admin);
1334 TEST_EXIT(x->size >= admin->size_used,
1335 "x->size = %d too small: admin->size_used = %d\n",
1336 x->size, admin->size_used);
1337
1338 FOR_ALL_DOFS(admin, SCAL_DOW(alpha, x->vec[dof]));
1339 }
1340
1341 static inline
__dof_dot_d(const DOF_REAL_D_VEC * x,const DOF_REAL_D_VEC * y)1342 REAL __dof_dot_d(const DOF_REAL_D_VEC *x, const DOF_REAL_D_VEC *y)
1343 {
1344 FUNCNAME("dof_dot_d");
1345 REAL dot;
1346 const DOF_ADMIN *admin = NULL;
1347
1348 TEST_EXIT(x && y,
1349 "pointer to DOF_REAL_D_VEC is NULL: x: %p, y: %p\n",
1350 x, y);
1351 TEST_EXIT(x->fe_space && y->fe_space,
1352 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1353 x->fe_space, y->fe_space);
1354
1355 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1356 "no admin or admins: x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1357 x->fe_space->admin, y->fe_space->admin);
1358
1359 TEST_EXIT(x->size >= admin->size_used,
1360 "x->size = %d too small: admin->size_used = %d\n",
1361 x->size, admin->size_used);
1362 TEST_EXIT(y->size >= admin->size_used,
1363 "y->size = %d too small: admin->size_used = %d\n",
1364 y->size, admin->size_used);
1365
1366 dot = 0.0;
1367 FOR_ALL_DOFS(admin, dot += SCP_DOW(x->vec[dof], y->vec[dof]));
1368
1369 return(dot);
1370 }
1371
1372 static inline
__dof_copy_d(const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)1373 void __dof_copy_d(const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
1374 {
1375 FUNCNAME("dof_copy_d");
1376 const DOF_ADMIN *admin = NULL;
1377
1378 TEST_EXIT(x && y,"pointer to DOF_REAL_D_VEC is NULL: x: %p, y: %p\n",
1379 x, y);
1380 TEST_EXIT(x->fe_space && y->fe_space,
1381 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1382 x->fe_space, y->fe_space);
1383
1384 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1385 "admin == NULL or admins differ: "
1386 "x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1387 x->fe_space->admin, y->fe_space->admin);
1388 TEST_EXIT(x->size >= admin->size_used,
1389 "x->size = %d too small: admin->size_used = %d\n",
1390 x->size, admin->size_used);
1391 TEST_EXIT(y->size >= admin->size_used,
1392 "y->size = %d too small: admin->size_used = %d\n",
1393 y->size, admin->size_used);
1394
1395 FOR_ALL_DOFS(admin, COPY_DOW(x->vec[dof], y->vec[dof]));
1396
1397 }
1398
1399 static inline
__dof_axpy_d(REAL alpha,const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)1400 void __dof_axpy_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
1401 {
1402 FUNCNAME("dof_axpy_d");
1403 const DOF_ADMIN *admin;
1404
1405 TEST_EXIT(x && y,"pointer to DOF_REAL_D_VEC is NULL: x: %p, y: %p\n",
1406 x, y);
1407 TEST_EXIT(x->fe_space && y->fe_space,
1408 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1409 x->fe_space, y->fe_space);
1410
1411 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1412 "no admin or admins: x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1413 x->fe_space->admin, y->fe_space->admin);
1414 TEST_EXIT(x->size >= admin->size_used,
1415 "x->size = %d too small: admin->size_used = %d\n",
1416 x->size, admin->size_used);
1417 TEST_EXIT(y->size >= admin->size_used,
1418 "y->size = %d too small: admin->size_used = %d\n",
1419 y->size, admin->size_used);
1420
1421 FOR_ALL_DOFS(admin, AXPY_DOW(alpha, x->vec[dof], y->vec[dof]));
1422 }
1423
1424 static inline
__dof_nrm2_d(const DOF_REAL_D_VEC * x)1425 REAL __dof_nrm2_d(const DOF_REAL_D_VEC *x)
1426 {
1427 FUNCNAME("dof_nrm2_d");
1428 REAL nrm;
1429 const DOF_ADMIN *admin = NULL;
1430
1431 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1432 "pointer is NULL: %p, %p\n", x, admin);
1433 TEST_EXIT(x->size >= admin->size_used,
1434 "x->size = %d too small: admin->size_used = %d\n", x->size,
1435 admin->size_used);
1436
1437 nrm = 0.0;
1438 FOR_ALL_DOFS(admin, nrm += NRM2_DOW(x->vec[dof]));
1439
1440 return(sqrt(nrm));
1441 }
1442
1443 static inline
__dof_asum_d(const DOF_REAL_D_VEC * x)1444 REAL __dof_asum_d(const DOF_REAL_D_VEC *x)
1445 {
1446 FUNCNAME("dof_nrm2_d");
1447 REAL nrm;
1448 const DOF_ADMIN *admin = NULL;
1449
1450 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1451 "pointer is NULL: %p, %p\n", x, admin);
1452 TEST_EXIT(x->size >= admin->size_used,
1453 "x->size = %d too small: admin->size_used = %d\n", x->size,
1454 admin->size_used);
1455
1456 nrm = 0.0;
1457 FOR_ALL_DOFS(admin, nrm += NORM1_DOW(x->vec[dof]));
1458
1459 return nrm;
1460 }
1461
1462
1463 static inline
__dof_xpay_d(REAL alpha,const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)1464 void __dof_xpay_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
1465 {
1466 FUNCNAME("dof_xpay_d");
1467 const DOF_ADMIN *admin;
1468 int n;
1469
1470 TEST_EXIT(x && y,"pointer to DOF_REAL_D_VEC is NULL: x: %p, y: %p\n",
1471 x, y);
1472 TEST_EXIT(x->fe_space && y->fe_space,
1473 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1474 x->fe_space, y->fe_space);
1475
1476 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1477 "no admin or admins: x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1478 x->fe_space->admin, y->fe_space->admin);
1479 TEST_EXIT(x->size >= admin->size_used,
1480 "x->size = %d too small: admin->size_used = %d\n",
1481 x->size, admin->size_used);
1482 TEST_EXIT(y->size >= admin->size_used,
1483 "y->size = %d too small: admin->size_used = %d\n",
1484 y->size, admin->size_used);
1485
1486 FOR_ALL_DOFS(admin,
1487 for (n = 0; n < DIM_OF_WORLD; n++)
1488 y->vec[dof][n] = x->vec[dof][n] + alpha*y->vec[dof][n]);
1489 }
1490
1491 static inline
__dof_min_d(const DOF_REAL_D_VEC * x)1492 REAL __dof_min_d(const DOF_REAL_D_VEC *x)
1493 {
1494 FUNCNAME("dof_min_d");
1495 REAL m, norm;
1496 const DOF_ADMIN *admin = NULL;
1497
1498 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1499 "pointer is NULL: %p, %p\n", x, admin);
1500 TEST_EXIT(x->size >= admin->size_used,
1501 "x->size = %d too small: admin->size_used = %d\n", x->size,
1502 admin->size_used);
1503
1504 m = 1.0E30;
1505 FOR_ALL_DOFS(admin,
1506 norm = NORM_DOW(x->vec[dof]);
1507 m = MIN(m, norm));
1508
1509 return(m);
1510 }
1511
1512 static inline
__dof_max_d(const DOF_REAL_D_VEC * x)1513 REAL __dof_max_d(const DOF_REAL_D_VEC *x)
1514 {
1515 FUNCNAME("dof_max_d");
1516 REAL m, norm;
1517 const DOF_ADMIN *admin = NULL;
1518
1519 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1520 "pointer is NULL: %p, %p\n", x, admin);
1521 TEST_EXIT(x->size >= admin->size_used,
1522 "x->size = %d too small: admin->size_used = %d\n", x->size,
1523 admin->size_used);
1524
1525 m = 0.0;
1526 FOR_ALL_DOFS(admin,
1527 norm = NORM_DOW(x->vec[dof]);
1528 m = MAX(m, norm));
1529
1530 return(m);
1531 }
1532
dof_nrm2_d(const DOF_REAL_D_VEC * x)1533 REAL dof_nrm2_d(const DOF_REAL_D_VEC *x)
1534 {
1535 REAL nrm2 = 0.0;
1536
1537 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1538 nrm2 += __dof_nrm2_d(x);
1539 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1540
1541 return sqrt(nrm2);
1542 }
1543
dof_asum_d(const DOF_REAL_D_VEC * x)1544 REAL dof_asum_d(const DOF_REAL_D_VEC *x)
1545 {
1546 REAL asum = 0.0;
1547 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1548 asum += __dof_asum_d(x);
1549 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1550
1551 return asum;
1552 }
1553
dof_set_d(REAL alpha,DOF_REAL_D_VEC * x)1554 void dof_set_d(REAL alpha, DOF_REAL_D_VEC *x)
1555 {
1556 CHAIN_DO(x, DOF_REAL_D_VEC) {
1557 __dof_set_d(alpha, x);
1558 } CHAIN_WHILE(x, DOF_REAL_D_VEC);
1559 }
1560
dof_scal_d(REAL alpha,DOF_REAL_D_VEC * x)1561 void dof_scal_d(REAL alpha, DOF_REAL_D_VEC *x)
1562 {
1563 CHAIN_DO(x, DOF_REAL_D_VEC) {
1564 __dof_scal_d(alpha, x);
1565 } CHAIN_WHILE(x, DOF_REAL_D_VEC);
1566 }
1567
dof_dot_d(const DOF_REAL_D_VEC * x,const DOF_REAL_D_VEC * y)1568 REAL dof_dot_d(const DOF_REAL_D_VEC *x, const DOF_REAL_D_VEC *y)
1569 {
1570 REAL dot = 0.0;
1571
1572 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1573 dot += __dof_dot_d(x, y);
1574 y = CHAIN_NEXT(y, const DOF_REAL_D_VEC);
1575 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1576
1577 return dot;
1578 }
1579
dof_copy_d(const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)1580 void dof_copy_d(const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
1581 {
1582 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1583 __dof_copy_d(x, y);
1584 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
1585 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1586 }
1587
dof_axpy_d(REAL alpha,const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)1588 void dof_axpy_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
1589 {
1590 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1591 __dof_axpy_d(alpha, x, y);
1592 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
1593 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1594 }
1595
dof_xpay_d(REAL alpha,const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)1596 void dof_xpay_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
1597 {
1598 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1599 __dof_xpay_d(alpha, x, y);
1600 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
1601 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1602 }
1603
dof_min_d(const DOF_REAL_D_VEC * x)1604 REAL dof_min_d(const DOF_REAL_D_VEC *x)
1605 {
1606 REAL m = REAL_MAX, m2;
1607 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1608 m2 = __dof_min_d(x);
1609 m = MIN(m, m2);
1610 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1611 return m;
1612 }
1613
dof_max_d(const DOF_REAL_D_VEC * x)1614 REAL dof_max_d(const DOF_REAL_D_VEC *x)
1615 {
1616 REAL m = REAL_MIN, m2;
1617 CHAIN_DO(x, const DOF_REAL_D_VEC) {
1618 m2 = __dof_max_d(x);
1619 m = MAX(m, m2);
1620 } CHAIN_WHILE(x, const DOF_REAL_D_VEC);
1621
1622 return m;
1623 }
1624
1625 /*******************************************************************************
1626 *
1627 * now the same routines for REAL_DDD vectors
1628 *
1629 ******************************************************************************/
1630
1631 static inline
__dof_set_dd(REAL alpha,DOF_REAL_DD_VEC * x)1632 void __dof_set_dd(REAL alpha, DOF_REAL_DD_VEC *x)
1633 {
1634 FUNCNAME("dof_set_d");
1635 const DOF_ADMIN *admin = NULL;
1636
1637 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1638 "pointer is NULL: x: %p, x->fe_space: %p, x->fe_space->admin :%p\n",
1639 x, x->fe_space, admin);
1640 TEST_EXIT(x->size >= admin->size_used,
1641 "x->size = %d too small: admin->size_used = %d\n",
1642 x->size, admin->size_used);
1643
1644 FOR_ALL_DOFS(admin, MSET_DOW(alpha, x->vec[dof]));
1645 }
1646
1647 static inline
__dof_scal_dd(REAL alpha,DOF_REAL_DD_VEC * x)1648 void __dof_scal_dd(REAL alpha, DOF_REAL_DD_VEC *x)
1649 {
1650 FUNCNAME("dof_scal_d");
1651 const DOF_ADMIN *admin = NULL;
1652
1653 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1654 "pointer is NULL: x: %p, x->fe_space: %p, x->fe_space->admin :%p\n",
1655 x, admin);
1656 TEST_EXIT(x->size >= admin->size_used,
1657 "x->size = %d too small: admin->size_used = %d\n",
1658 x->size, admin->size_used);
1659
1660 FOR_ALL_DOFS(admin, MSCAL_DOW(alpha, x->vec[dof]));
1661 }
1662
1663 static inline
__dof_dot_dd(const DOF_REAL_DD_VEC * x,const DOF_REAL_DD_VEC * y)1664 REAL __dof_dot_dd(const DOF_REAL_DD_VEC *x, const DOF_REAL_DD_VEC *y)
1665 {
1666 FUNCNAME("dof_dot_d");
1667 REAL dot;
1668 const DOF_ADMIN *admin = NULL;
1669
1670 TEST_EXIT(x && y,
1671 "pointer to DOF_REAL_DD_VEC is NULL: x: %p, y: %p\n",
1672 x, y);
1673 TEST_EXIT(x->fe_space && y->fe_space,
1674 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1675 x->fe_space, y->fe_space);
1676
1677 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1678 "no admin or admins: x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1679 x->fe_space->admin, y->fe_space->admin);
1680
1681 TEST_EXIT(x->size >= admin->size_used,
1682 "x->size = %d too small: admin->size_used = %d\n",
1683 x->size, admin->size_used);
1684 TEST_EXIT(y->size >= admin->size_used,
1685 "y->size = %d too small: admin->size_used = %d\n",
1686 y->size, admin->size_used);
1687
1688 dot = 0.0;
1689 FOR_ALL_DOFS(admin, dot += MSCP_DOW((const REAL_D *)x->vec[dof],
1690 (const REAL_D *)y->vec[dof]));
1691
1692 return dot;
1693 }
1694
1695 static inline
__dof_copy_dd(const DOF_REAL_DD_VEC * x,DOF_REAL_DD_VEC * y)1696 void __dof_copy_dd(const DOF_REAL_DD_VEC *x, DOF_REAL_DD_VEC *y)
1697 {
1698 FUNCNAME("dof_copy_d");
1699 const DOF_ADMIN *admin = NULL;
1700
1701 TEST_EXIT(x && y,"pointer to DOF_REAL_DD_VEC is NULL: x: %p, y: %p\n",
1702 x, y);
1703 TEST_EXIT(x->fe_space && y->fe_space,
1704 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1705 x->fe_space, y->fe_space);
1706
1707 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1708 "admin == NULL or admins differ: "
1709 "x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1710 x->fe_space->admin, y->fe_space->admin);
1711 TEST_EXIT(x->size >= admin->size_used,
1712 "x->size = %d too small: admin->size_used = %d\n",
1713 x->size, admin->size_used);
1714 TEST_EXIT(y->size >= admin->size_used,
1715 "y->size = %d too small: admin->size_used = %d\n",
1716 y->size, admin->size_used);
1717
1718 FOR_ALL_DOFS(admin, MCOPY_DOW((const REAL_D *)x->vec[dof], y->vec[dof]));
1719
1720 }
1721
1722 static inline
__dof_axpy_dd(REAL alpha,const DOF_REAL_DD_VEC * x,DOF_REAL_DD_VEC * y)1723 void __dof_axpy_dd(REAL alpha, const DOF_REAL_DD_VEC *x, DOF_REAL_DD_VEC *y)
1724 {
1725 FUNCNAME("dof_axpy_d");
1726 const DOF_ADMIN *admin;
1727
1728 TEST_EXIT(x && y,"pointer to DOF_REAL_DD_VEC is NULL: x: %p, y: %p\n",
1729 x, y);
1730 TEST_EXIT(x->fe_space && y->fe_space,
1731 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1732 x->fe_space, y->fe_space);
1733
1734 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1735 "no admin or admins: x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1736 x->fe_space->admin, y->fe_space->admin);
1737 TEST_EXIT(x->size >= admin->size_used,
1738 "x->size = %d too small: admin->size_used = %d\n",
1739 x->size, admin->size_used);
1740 TEST_EXIT(y->size >= admin->size_used,
1741 "y->size = %d too small: admin->size_used = %d\n",
1742 y->size, admin->size_used);
1743
1744 FOR_ALL_DOFS(admin,
1745 MAXPY_DOW(alpha, (const REAL_D *)x->vec[dof], y->vec[dof]));
1746 }
1747
1748 static inline
__dof_nrm2_dd(const DOF_REAL_DD_VEC * x)1749 REAL __dof_nrm2_dd(const DOF_REAL_DD_VEC *x)
1750 {
1751 FUNCNAME("dof_nrm2_d");
1752 REAL nrm;
1753 const DOF_ADMIN *admin = NULL;
1754
1755 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1756 "pointer is NULL: %p, %p\n", x, admin);
1757 TEST_EXIT(x->size >= admin->size_used,
1758 "x->size = %d too small: admin->size_used = %d\n", x->size,
1759 admin->size_used);
1760
1761 nrm = 0.0;
1762 FOR_ALL_DOFS(admin, nrm += MNRM2_DOW((const REAL_D *)x->vec[dof]));
1763
1764 return sqrt(nrm);
1765 }
1766
1767 static inline
__dof_asum_dd(const DOF_REAL_DD_VEC * x)1768 REAL __dof_asum_dd(const DOF_REAL_DD_VEC *x)
1769 {
1770 FUNCNAME("dof_nrm2_d");
1771 REAL nrm;
1772 const DOF_ADMIN *admin = NULL;
1773
1774 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1775 "pointer is NULL: %p, %p\n", x, admin);
1776 TEST_EXIT(x->size >= admin->size_used,
1777 "x->size = %d too small: admin->size_used = %d\n", x->size,
1778 admin->size_used);
1779
1780 nrm = 0.0;
1781 FOR_ALL_DOFS(admin, nrm += MNORM1_DOW((const REAL_D *)x->vec[dof]));
1782
1783 return nrm;
1784 }
1785
1786
1787 static inline
__dof_xpay_dd(REAL alpha,const DOF_REAL_DD_VEC * x,DOF_REAL_DD_VEC * y)1788 void __dof_xpay_dd(REAL alpha, const DOF_REAL_DD_VEC *x, DOF_REAL_DD_VEC *y)
1789 {
1790 FUNCNAME("dof_xpay_d");
1791 const DOF_ADMIN *admin;
1792
1793 TEST_EXIT(x && y,"pointer to DOF_REAL_DD_VEC is NULL: x: %p, y: %p\n",
1794 x, y);
1795 TEST_EXIT(x->fe_space && y->fe_space,
1796 "pointer to FE_SPACE is NULL: x->fe_space: %p, y->fe_space: %p\n",
1797 x->fe_space, y->fe_space);
1798
1799 TEST_EXIT((admin = x->fe_space->admin) && (admin == y->fe_space->admin),
1800 "no admin or admins: x->fe_space->admin: %p, y->fe_space->admin: %p\n",
1801 x->fe_space->admin, y->fe_space->admin);
1802 TEST_EXIT(x->size >= admin->size_used,
1803 "x->size = %d too small: admin->size_used = %d\n",
1804 x->size, admin->size_used);
1805 TEST_EXIT(y->size >= admin->size_used,
1806 "y->size = %d too small: admin->size_used = %d\n",
1807 y->size, admin->size_used);
1808
1809 FOR_ALL_DOFS(admin,
1810 MAXPBY_DOW(1.0, (const REAL_D *)x->vec[dof],
1811 alpha, (const REAL_D *)y->vec[dof], y->vec[dof]));
1812 }
1813
1814 static inline
__dof_min_dd(const DOF_REAL_DD_VEC * x)1815 REAL __dof_min_dd(const DOF_REAL_DD_VEC *x)
1816 {
1817 FUNCNAME("dof_min_d");
1818 REAL m, norm;
1819 const DOF_ADMIN *admin = NULL;
1820
1821 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1822 "pointer is NULL: %p, %p\n", x, admin);
1823 TEST_EXIT(x->size >= admin->size_used,
1824 "x->size = %d too small: admin->size_used = %d\n", x->size,
1825 admin->size_used);
1826
1827 m = 1.0E30;
1828 FOR_ALL_DOFS(admin,
1829 norm = MNORM8_DOW((const REAL_D *)x->vec[dof]);
1830 m = MIN(m, norm));
1831
1832 return m;
1833 }
1834
1835 static inline
__dof_max_dd(const DOF_REAL_DD_VEC * x)1836 REAL __dof_max_dd(const DOF_REAL_DD_VEC *x)
1837 {
1838 FUNCNAME("dof_max_d");
1839 REAL m, norm;
1840 const DOF_ADMIN *admin = NULL;
1841
1842 TEST_EXIT(x && x->fe_space && (admin = x->fe_space->admin),
1843 "pointer is NULL: %p, %p\n", x, admin);
1844 TEST_EXIT(x->size >= admin->size_used,
1845 "x->size = %d too small: admin->size_used = %d\n", x->size,
1846 admin->size_used);
1847
1848 m = 0.0;
1849 FOR_ALL_DOFS(admin,
1850 norm = MNORM8_DOW((const REAL_D *)x->vec[dof]);
1851 m = MAX(m, norm));
1852
1853 return(m);
1854 }
1855
dof_nrm2_dd(const DOF_REAL_DD_VEC * x)1856 REAL dof_nrm2_dd(const DOF_REAL_DD_VEC *x)
1857 {
1858 REAL nrm2 = 0.0;
1859
1860 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1861 nrm2 += __dof_nrm2_dd(x);
1862 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1863
1864 return sqrt(nrm2);
1865 }
1866
dof_asum_dd(const DOF_REAL_DD_VEC * x)1867 REAL dof_asum_dd(const DOF_REAL_DD_VEC *x)
1868 {
1869 REAL asum = 0.0;
1870 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1871 asum += __dof_asum_dd(x);
1872 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1873
1874 return asum;
1875 }
1876
dof_set_dd(REAL alpha,DOF_REAL_DD_VEC * x)1877 void dof_set_dd(REAL alpha, DOF_REAL_DD_VEC *x)
1878 {
1879 CHAIN_DO(x, DOF_REAL_DD_VEC) {
1880 __dof_set_dd(alpha, x);
1881 } CHAIN_WHILE(x, DOF_REAL_DD_VEC);
1882 }
1883
dof_scal_dd(REAL alpha,DOF_REAL_DD_VEC * x)1884 void dof_scal_dd(REAL alpha, DOF_REAL_DD_VEC *x)
1885 {
1886 CHAIN_DO(x, DOF_REAL_DD_VEC) {
1887 __dof_scal_dd(alpha, x);
1888 } CHAIN_WHILE(x, DOF_REAL_DD_VEC);
1889 }
1890
dof_dot_dd(const DOF_REAL_DD_VEC * x,const DOF_REAL_DD_VEC * y)1891 REAL dof_dot_dd(const DOF_REAL_DD_VEC *x, const DOF_REAL_DD_VEC *y)
1892 {
1893 REAL dot = 0.0;
1894
1895 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1896 dot += __dof_dot_dd(x, y);
1897 y = CHAIN_NEXT(y, const DOF_REAL_DD_VEC);
1898 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1899
1900 return dot;
1901 }
1902
dof_copy_dd(const DOF_REAL_DD_VEC * x,DOF_REAL_DD_VEC * y)1903 void dof_copy_dd(const DOF_REAL_DD_VEC *x, DOF_REAL_DD_VEC *y)
1904 {
1905 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1906 __dof_copy_dd(x, y);
1907 y = CHAIN_NEXT(y, DOF_REAL_DD_VEC);
1908 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1909 }
1910
dof_axpy_dd(REAL alpha,const DOF_REAL_DD_VEC * x,DOF_REAL_DD_VEC * y)1911 void dof_axpy_dd(REAL alpha, const DOF_REAL_DD_VEC *x, DOF_REAL_DD_VEC *y)
1912 {
1913 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1914 __dof_axpy_dd(alpha, x, y);
1915 y = CHAIN_NEXT(y, DOF_REAL_DD_VEC);
1916 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1917 }
1918
dof_xpay_dd(REAL alpha,const DOF_REAL_DD_VEC * x,DOF_REAL_DD_VEC * y)1919 void dof_xpay_dd(REAL alpha, const DOF_REAL_DD_VEC *x, DOF_REAL_DD_VEC *y)
1920 {
1921 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1922 __dof_xpay_dd(alpha, x, y);
1923 y = CHAIN_NEXT(y, DOF_REAL_DD_VEC);
1924 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1925 }
1926
dof_min_dd(const DOF_REAL_DD_VEC * x)1927 REAL dof_min_dd(const DOF_REAL_DD_VEC *x)
1928 {
1929 REAL m = REAL_MAX, m2;
1930 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1931 m2 = __dof_min_dd(x);
1932 m = MIN(m, m2);
1933 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1934 return m;
1935 }
1936
dof_max_dd(const DOF_REAL_DD_VEC * x)1937 REAL dof_max_dd(const DOF_REAL_DD_VEC *x)
1938 {
1939 REAL m = REAL_MIN, m2;
1940 CHAIN_DO(x, const DOF_REAL_DD_VEC) {
1941 m2 = __dof_max_dd(x);
1942 m = MAX(m, m2);
1943 } CHAIN_WHILE(x, const DOF_REAL_DD_VEC);
1944
1945 return m;
1946 }
1947
1948 /* DOF_REAL_VEC_D versions (vector valued basis functions) */
1949
1950 #define __DRV DOF_REAL_VEC
1951 #define __DRDV DOF_REAL_D_VEC
1952
dof_nrm2_dow(const DOF_REAL_VEC_D * x)1953 REAL dof_nrm2_dow(const DOF_REAL_VEC_D *x)
1954 {
1955 REAL nrm2 = 0.0;
1956
1957 CHAIN_DO(x, const DOF_REAL_VEC_D) {
1958 nrm2 += x->stride != 1 ? __dof_nrm2_d((__DRDV *)x) : __dof_nrm2((__DRV *)x);
1959 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
1960
1961 return sqrt(nrm2);
1962 }
1963
dof_asum_dow(const DOF_REAL_VEC_D * x)1964 REAL dof_asum_dow(const DOF_REAL_VEC_D *x)
1965 {
1966 REAL asum = 0.0;
1967
1968 CHAIN_DO(x, const DOF_REAL_VEC_D) {
1969 asum += x->stride != 1 ? __dof_asum_d((__DRDV *)x) : __dof_asum((__DRV *)x);
1970 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
1971
1972 return asum;
1973 }
1974
dof_set_dow(REAL alpha,DOF_REAL_VEC_D * x)1975 void dof_set_dow(REAL alpha, DOF_REAL_VEC_D *x)
1976 {
1977 CHAIN_DO(x, DOF_REAL_VEC_D) {
1978 if (x->stride != 1) {
1979 __dof_set_d(alpha, (__DRDV *)x);
1980 } else {
1981 __dof_set(alpha, (__DRV *)x);
1982 }
1983 } CHAIN_WHILE(x, DOF_REAL_VEC_D);
1984 }
1985
dof_scal_dow(REAL alpha,DOF_REAL_VEC_D * x)1986 void dof_scal_dow(REAL alpha, DOF_REAL_VEC_D *x)
1987 {
1988 CHAIN_DO(x, DOF_REAL_VEC_D) {
1989 if (x->stride != 1) {
1990 __dof_scal_d(alpha, (__DRDV *)x);
1991 } else {
1992 __dof_scal(alpha, (__DRV *)x);
1993 }
1994 } CHAIN_WHILE(x, DOF_REAL_VEC_D);
1995 }
1996
dof_dot_dow(const DOF_REAL_VEC_D * x,const DOF_REAL_VEC_D * y)1997 REAL dof_dot_dow(const DOF_REAL_VEC_D *x, const DOF_REAL_VEC_D *y)
1998 {
1999 REAL dot = 0.0;
2000
2001 CHAIN_DO(x, const DOF_REAL_VEC_D) {
2002 if (x->stride != 1) {
2003 dot += __dof_dot_d((__DRDV *)x, (__DRDV *)y);
2004 } else {
2005 dot += __dof_dot((__DRV *)x, (__DRV *)y);
2006 }
2007 y = CHAIN_NEXT(y, const DOF_REAL_VEC_D);
2008 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
2009
2010 return dot;
2011 }
2012
dof_copy_dow(const DOF_REAL_VEC_D * x,DOF_REAL_VEC_D * y)2013 void dof_copy_dow(const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y)
2014 {
2015 CHAIN_DO(x, const DOF_REAL_VEC_D) {
2016 if (x->stride != 1) {
2017 __dof_copy_d((__DRDV *)x, (__DRDV *)y);
2018 } else {
2019 __dof_copy((__DRV *)x, (__DRV *)y);
2020 }
2021 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2022 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
2023 }
2024
dof_axpy_dow(REAL alpha,const DOF_REAL_VEC_D * x,DOF_REAL_VEC_D * y)2025 void dof_axpy_dow(REAL alpha, const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y)
2026 {
2027 CHAIN_DO(x, const DOF_REAL_VEC_D) {
2028 if (x->stride != 1) {
2029 __dof_axpy_d(alpha, (__DRDV *)x, (__DRDV *)y);
2030 } else {
2031 __dof_axpy(alpha, (__DRV *)x, (__DRV *)y);
2032 }
2033 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2034 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
2035 }
2036
dof_xpay_dow(REAL alpha,const DOF_REAL_VEC_D * x,DOF_REAL_VEC_D * y)2037 void dof_xpay_dow(REAL alpha, const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y)
2038 {
2039 CHAIN_DO(x, const DOF_REAL_VEC_D) {
2040 if (x->stride != 1) {
2041 __dof_xpay_d(alpha, (__DRDV *)x, (__DRDV *)y);
2042 } else {
2043 __dof_xpay(alpha, (__DRV *)x, (__DRV *)y);
2044 }
2045 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2046 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
2047 }
2048
dof_min_dow(const DOF_REAL_VEC_D * x)2049 REAL dof_min_dow(const DOF_REAL_VEC_D *x)
2050 {
2051 REAL m = REAL_MAX, m2;
2052 CHAIN_DO(x, const DOF_REAL_VEC_D) {
2053 if (x->stride != 1) {
2054 m2 = __dof_min_d((__DRDV *)x);
2055 } else {
2056 m2 = __dof_min((__DRV *)x);
2057 }
2058 m = MIN(m, m2);
2059 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
2060 return m;
2061 }
2062
dof_max_dow(const DOF_REAL_VEC_D * x)2063 REAL dof_max_dow(const DOF_REAL_VEC_D *x)
2064 {
2065 REAL m = REAL_MIN, m2;
2066 CHAIN_DO(x, const DOF_REAL_VEC_D) {
2067 if (x->stride != 1) {
2068 m2 = __dof_max_d((__DRDV *)x);
2069 } else {
2070 m2 = __dof_max((__DRV *)x);
2071 }
2072 m = MAX(m, m2);
2073 } CHAIN_WHILE(x, const DOF_REAL_VEC_D);
2074
2075 return m;
2076 }
2077
2078 #undef __DRV
2079 #undef __DRDV
2080
2081 /* Copy one DOF_MATRIX to another, e.g. to speed up the assembly for
2082 * operator splitting schemes.
2083 */
2084 static inline
_AI_matrix_row_copy_single(MATRIX_ROW * dst,const MATRIX_ROW * src)2085 void _AI_matrix_row_copy_single(MATRIX_ROW *dst, const MATRIX_ROW *src)
2086 {
2087 MATRIX_ROW *dst_next = dst->next;
2088
2089 DEBUG_TEST_EXIT(dst->type == src->type, "matrix types do not match");
2090
2091 switch (dst->type) {
2092 case MATENT_REAL:
2093 memcpy(dst, src, sizeof(MATRIX_ROW_REAL));
2094 break;
2095 case MATENT_REAL_D:
2096 memcpy(dst, src, sizeof(MATRIX_ROW_REAL_D));
2097 break;
2098 case MATENT_REAL_DD:
2099 memcpy(dst, src, sizeof(MATRIX_ROW_REAL_DD));
2100 break;
2101 case MATENT_NONE:
2102 ERROR_EXIT("Uninitialized DOF_MATRIX.\n");
2103 }
2104 dst->next = dst_next;
2105 }
2106
2107 static inline void _AI_clear_dof_matrix_single(DOF_MATRIX *matrix);
2108
2109 FLATTEN_ATTR
2110 static inline
_AI_dof_matrix_copy_single(DOF_MATRIX * dst,const DOF_MATRIX * src)2111 void _AI_dof_matrix_copy_single(DOF_MATRIX *dst, const DOF_MATRIX *src)
2112 {
2113 const FE_SPACE *row_fe_space = dst->row_fe_space;
2114 const DOF_ADMIN *m_admin = row_fe_space->admin;
2115 int dof;
2116
2117 #if ALBERTA_DEBUG
2118 const FE_SPACE *col_fe_space =
2119 dst->col_fe_space ? dst->col_fe_space : row_fe_space;
2120 DEBUG_TEST_EXIT(dst->row_fe_space->admin == src->row_fe_space->admin &&
2121 (src->col_fe_space == NULL ||
2122 col_fe_space->admin == src->col_fe_space->admin),
2123 "Attempt to copy onto incompatible DOF_MATRIX.\n");
2124
2125 #endif
2126
2127 if (dst->type != src->type) {
2128 _AI_clear_dof_matrix_single(dst);
2129 dst->type = src->type;
2130 }
2131 BNDRY_FLAGS_CPY(dst->dirichlet_bndry, src->dirichlet_bndry);
2132 /* remaining stuff has to be handled by calling function */
2133
2134 if (src->is_diagonal) {
2135 dof_matrix_set_diagonal(dst, true);
2136 FOR_ALL_DOFS(src->row_fe_space->admin,
2137 dst->diag_cols->vec[dof] = src->diag_cols->vec[dof]);
2138 switch (src->type) {
2139 case MATENT_REAL:
2140 if (dst->diagonal.real == NULL) {
2141 dst->diagonal.real =
2142 get_dof_real_vec("matrix diagonal", dst->row_fe_space->unchained);
2143 }
2144 dof_copy(src->diagonal.real, dst->diagonal.real);
2145 break;
2146 case MATENT_REAL_D:
2147 if (dst->diagonal.real_d == NULL) {
2148 dst->diagonal.real_d =
2149 get_dof_real_d_vec("matrix diagonal", dst->row_fe_space->unchained);
2150 }
2151 dof_copy_d(src->diagonal.real_d, dst->diagonal.real_d);
2152 break;
2153 case MATENT_REAL_DD:
2154 if (dst->diagonal.real_dd == NULL) {
2155 dst->diagonal.real_dd =
2156 get_dof_real_dd_vec("matrix diagonal", dst->row_fe_space->unchained);
2157 }
2158 dof_copy_dd(src->diagonal.real_dd, dst->diagonal.real_dd);
2159 break;
2160 default:
2161 break;
2162 }
2163 } else {
2164 dof_matrix_set_diagonal(dst, false);
2165 for (dof = 0; dof < m_admin->size_used; ++dof) {
2166 MATRIX_ROW *src_row, **dst_row_ptr, *dst_row;
2167
2168 for (src_row = src->matrix_row[dof],
2169 dst_row_ptr = &dst->matrix_row[dof];
2170 src_row != NULL;
2171 src_row = src_row->next, dst_row_ptr = &(*dst_row_ptr)->next) {
2172 if (*dst_row_ptr == NULL) {
2173 *dst_row_ptr = get_matrix_row(row_fe_space, dst->type);
2174 }
2175 _AI_matrix_row_copy_single(*dst_row_ptr, src_row);
2176 }
2177 /* zap remaining part of row */
2178 dst_row = *dst_row_ptr;
2179 *dst_row_ptr = NULL;
2180 while (dst_row) {
2181 MATRIX_ROW *free_row = dst_row;
2182 dst_row = free_row->next;
2183 free_matrix_row(row_fe_space, free_row);
2184 }
2185 }
2186 }
2187 }
2188
2189 FLATTEN_ATTR
dof_matrix_copy(DOF_MATRIX * dst,const DOF_MATRIX * src)2190 void dof_matrix_copy(DOF_MATRIX *dst, const DOF_MATRIX *src)
2191 {
2192 DEBUG_TEST_EXIT(fe_space_is_eq(dst->row_fe_space, src->row_fe_space) &&
2193 (dst->col_fe_space == NULL ||
2194 fe_space_is_eq(dst->col_fe_space,
2195 src->col_fe_space == NULL
2196 ? src->row_fe_space : src->col_fe_space)),
2197 "Attempt to copy onto incompatible DOF_MATRIX\n");
2198 ROW_CHAIN_DO(dst, DOF_MATRIX) {
2199 COL_CHAIN_DO(dst, DOF_MATRIX) {
2200 _AI_dof_matrix_copy_single(dst, src);
2201 src = COL_CHAIN_NEXT(src, const DOF_MATRIX);
2202 } COL_CHAIN_WHILE(dst, DOF_MATRIX);
2203 src = ROW_CHAIN_NEXT(src, const DOF_MATRIX);
2204 } ROW_CHAIN_WHILE(dst, DOF_MATRIX);
2205 }
2206
2207 /* dof_mv() and friends:
2208 *
2209 * matrix type domain type range type interpretation function
2210 * ============== =========== ========== ============== =========
2211 * MATENT_REAL REAL_VEC REAL_VEC ordinary matrix dof_mv()
2212 * MATENT_REAL REAL_D_VEC REAL_VEC 1xDOW (scalar blocks) dof_mv_rrd()
2213 * MATENT_REAL REAL_VEC REAL_D_VEC DOWx1 (scalar blocks) dof_mv_rdr()
2214 * MATENT_REAL REAL_D_VEC REAL_D_VEC scalar matrix blocks dof_mv_d()
2215 * MATENT_REAL_D REAL_D_VEC REAL_VEC 1xDOW rectangular dof_mv_rrd()
2216 * MATENT_REAL_D REAL_VEC REAL_D_VEC DOWx1 rectangular dof_mv_rdr()
2217 * MATENT_REAL_D REAL_D_VEC REAL_D_VEC diagonal matrix (?) dof_mv_d()
2218 * MATENT_REAL_DD REAL_VEC ---------- does not make sense --------
2219 * MATENT_REAL_DD REAL_D_VEC REAL_D_VEC DOWxDOW square blocks dof_mv_d()
2220 *
2221 * The suffix codes the type of the domain and range vectors. "rdr"
2222 * "REAL_D x REAL" and "rrd" means "REAL x REAL_D". We do not need
2223 * _dowb() versions any longer, at least.
2224 */
2225
2226 /* Used in all mv routines: */
2227 #define MV_TEST_AND_GET_ADMINS(m, x, y, m_admin, x_admin, y_admin) \
2228 { \
2229 DEBUG_TEST_EXIT(m && x && y,"pointer is NULL: %p, %p, %p\n", m, x, y); \
2230 \
2231 DEBUG_TEST_EXIT(m->row_fe_space && m->col_fe_space && \
2232 x->fe_space && y->fe_space, \
2233 "fe_space is NULL: %p, %p, %p, %p\n", \
2234 m->row_fe_space, m->col_fe_space, \
2235 x->fe_space, y->fe_space); \
2236 \
2237 m_admin = m->row_fe_space->admin; \
2238 DEBUG_TEST_EXIT(m_admin != NULL, \
2239 "no matrix row-admin: %p.\n", m->row_fe_space->admin); \
2240 x_admin = x->fe_space->admin; \
2241 DEBUG_TEST_EXIT(x_admin != NULL, \
2242 "no admin for x: %p.\n", x->fe_space->admin); \
2243 y_admin = y->fe_space->admin; \
2244 DEBUG_TEST_EXIT(y_admin != NULL, \
2245 "no admin for y: %p.\n", y->fe_space->admin); \
2246 \
2247 DEBUG_TEST_EXIT(x->size >= x_admin->size_used, \
2248 "x->size = %d too small: admin->size_used = %d\n", \
2249 x->size, x_admin->size_used); \
2250 DEBUG_TEST_EXIT(y->size >= y_admin->size_used, \
2251 "y->size = %d too small: admin->size_used = %d\n", \
2252 y->size, y_admin->size_used); \
2253 DEBUG_TEST_EXIT(a->size >= m_admin->size_used, \
2254 "a->size = %d too small: admin->size_used = %d\n", \
2255 a->size, m_admin->size_used); \
2256 }
2257
2258 /* <<< REAL x REAL */
2259
2260 FORCE_INLINE_ATTR
2261 static inline
__dof_gemv(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,REAL beta,DOF_REAL_VEC * y)2262 void __dof_gemv(MatrixTranspose transpose, REAL alpha,
2263 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2264 const DOF_REAL_VEC *x, REAL beta, DOF_REAL_VEC *y)
2265 {
2266 FUNCNAME("dof_gemv");
2267 int ysize, dof;
2268 REAL sum, ax;
2269 REAL *xvec, *yvec;
2270 const DOF_ADMIN *m_admin, *y_admin, *x_admin;
2271
2272 MV_TEST_AND_GET_ADMINS(a, x, y, m_admin, x_admin, y_admin);
2273
2274 DEBUG_TEST_EXIT(a->type == MATENT_REAL, "incompatible block-matrix type");
2275
2276 xvec = x->vec;
2277 yvec = y->vec;
2278
2279 ysize = y->size;
2280
2281 FOR_ALL_FREE_DOFS(y_admin, if (dof < ysize) yvec[dof] = 0.0);
2282
2283 if (a->is_diagonal) {
2284 if (x_admin == y_admin) {
2285 FOR_ALL_DOFS(m_admin, {
2286 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2287 yvec[dof] =
2288 beta * yvec[dof] + alpha * a->diagonal.real->vec[dof] * xvec[dof];
2289 } else {
2290 yvec[dof] *= beta;
2291 }
2292 });
2293 } else if (transpose == NoTranspose) {
2294 DOF *col_dofs = a->diag_cols->vec;
2295 FOR_ALL_DOFS(m_admin, {
2296 DOF col_dof = col_dofs[dof];
2297 if (ENTRY_USED(col_dof) &&
2298 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
2299 yvec[dof] =
2300 beta * yvec[dof]
2301 +
2302 alpha * a->diagonal.real->vec[dof] * xvec[col_dof];
2303 } else {
2304 yvec[dof] *= beta;
2305 }
2306 });
2307 } else {
2308 DOF *col_dofs = a->diag_cols->vec;
2309 FOR_ALL_DOFS(m_admin, {
2310 DOF col_dof = col_dofs[dof];
2311 if (ENTRY_USED(col_dof) &&
2312 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
2313 yvec[col_dof] =
2314 beta * yvec[col_dof]
2315 +
2316 alpha * a->diagonal.real->vec[dof] * xvec[dof];
2317 } else {
2318 yvec[dof] *= beta;
2319 }
2320 });
2321 }
2322 } else if (transpose == NoTranspose) { // not diagonal
2323
2324 DEBUG_TEST_EXIT(m_admin == y_admin,
2325 "matrix- and y-admins do not match: %p %p.\n",
2326 m_admin, y_admin);
2327
2328 for (dof = 0; dof < m_admin->size_used; dof++) {
2329 sum = 0.0;
2330 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2331 FOR_ALL_MAT_COLS(REAL, a->matrix_row[dof],
2332 sum += row->entry[col_idx] * xvec[col_dof]);
2333 }
2334 yvec[dof] = beta * yvec[dof] + alpha * sum;
2335 }
2336 } else if (transpose == Transpose) {
2337
2338 TEST_EXIT(m_admin == x_admin,
2339 "matrix- and x-admins do not match: %p %p.\n",
2340 m_admin, x_admin);
2341
2342 FOR_ALL_DOFS(y_admin, yvec[dof] *= beta);
2343 for (dof = 0; dof < m_admin->size_used; dof++) {
2344 ax = alpha * xvec[dof];
2345 FOR_ALL_MAT_COLS(REAL, a->matrix_row[dof],
2346 if (mask == NULL || mask->vec[col_dof] < DIRICHLET) {
2347 yvec[col_dof] += ax * row->entry[col_idx];
2348 });
2349 }
2350 } else {
2351 ERROR_EXIT("transpose=%d\n", transpose);
2352 }
2353 }
2354
2355 FLATTEN_ATTR
dof_gemv(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,REAL beta,DOF_REAL_VEC * y)2356 void dof_gemv(MatrixTranspose transpose, REAL alpha,
2357 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
2358 const DOF_REAL_VEC *x,
2359 REAL beta, DOF_REAL_VEC *y)
2360 {
2361 const DOF_MATRIX *A_chain;
2362
2363 if (transpose == NoTranspose) {
2364 COL_CHAIN_DO(A, const DOF_MATRIX) {
2365 if (mask) {
2366 __dof_gemv(transpose, alpha, A, mask, x, beta, y);
2367 } else {
2368 __dof_gemv(transpose, alpha, A, NULL, x, beta, y);
2369 }
2370 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2371 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2372 if (mask) {
2373 __dof_gemv(transpose, alpha, A_chain, mask, x, 1.0, y);
2374 } else {
2375 __dof_gemv(transpose, alpha, A_chain, NULL, x, 1.0, y);
2376 }
2377 }
2378 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2379 y = CHAIN_NEXT(y, DOF_REAL_VEC);
2380 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2381 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
2382 } else {
2383 ROW_CHAIN_DO(A, const DOF_MATRIX) {
2384 if (mask) {
2385 __dof_gemv(transpose, alpha, A, mask, x, beta, y);
2386 } else {
2387 __dof_gemv(transpose, alpha, A, NULL, x, beta, y);
2388 }
2389 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2390 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2391 if (mask) {
2392 __dof_gemv(transpose, alpha, A_chain, mask, x, 1.0, y);
2393 } else {
2394 __dof_gemv(transpose, alpha, A_chain, NULL, x, 1.0, y);
2395 }
2396 }
2397 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2398 y = CHAIN_NEXT(y, DOF_REAL_VEC);
2399 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2400 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
2401 }
2402 }
2403
2404 /*--------------------------------------------------------------------------*/
2405
2406 FORCE_INLINE_ATTR
2407 static inline
__dof_mv(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,DOF_REAL_VEC * y)2408 void __dof_mv(MatrixTranspose transpose,
2409 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2410 const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
2411 {
2412 FUNCNAME("dof_mv");
2413 REAL sum, ax, *xvec, *yvec;
2414 const DOF_ADMIN *m_admin, *y_admin, *x_admin;
2415 int dof, ysize;
2416
2417 MV_TEST_AND_GET_ADMINS(a, x, y, m_admin, x_admin, y_admin);
2418
2419 DEBUG_TEST_EXIT(a->type == MATENT_REAL, "incompatible block-matrix type");
2420
2421 xvec = x->vec;
2422 yvec = y->vec;
2423
2424 ysize = y->size;
2425 FOR_ALL_FREE_DOFS(y_admin, if (dof < ysize) yvec[dof] = 0.0);
2426
2427 if (a->is_diagonal) {
2428 if (x_admin == y_admin) {
2429 FOR_ALL_DOFS(m_admin, {
2430 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2431 yvec[dof] = a->diagonal.real->vec[dof] * xvec[dof];
2432 }
2433 });
2434 } else if (transpose == NoTranspose) {
2435 DOF *col_dofs = a->diag_cols->vec;
2436 FOR_ALL_DOFS(m_admin, {
2437 if (ENTRY_USED(col_dofs[dof]) &&
2438 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
2439 yvec[dof] = a->diagonal.real->vec[dof] * xvec[col_dofs[dof]];
2440 }
2441 });
2442 } else {
2443 DOF *col_dofs = a->diag_cols->vec;
2444 FOR_ALL_DOFS(m_admin, {
2445 if (ENTRY_USED(col_dofs[dof]) &&
2446 (mask == NULL || mask->vec[col_dofs[dof]] < DIRICHLET)) {
2447 yvec[col_dofs[dof]] = a->diagonal.real->vec[dof] * xvec[dof];
2448 }
2449 });
2450 }
2451 } else if (transpose == NoTranspose) {
2452
2453 DEBUG_TEST_EXIT(m_admin == y_admin,
2454 "matrix- and y-admins do not match: %p %p.\n",
2455 m_admin, y_admin);
2456
2457 for (dof = 0; dof < m_admin->size_used; dof++) {
2458 sum = 0.0;
2459 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2460 FOR_ALL_MAT_COLS(REAL, a->matrix_row[dof],
2461 sum += row->entry[col_idx] * xvec[col_dof]);
2462 }
2463 yvec[dof] = sum;
2464 }
2465 } else if (transpose == Transpose) {
2466
2467 DEBUG_TEST_EXIT(m_admin == x_admin,
2468 "matrix- and x-admins do not match: %p %p.\n",
2469 m_admin, x_admin);
2470
2471 FOR_ALL_DOFS(y_admin, yvec[dof] = 0.0);
2472
2473 for (dof = 0; dof < m_admin->size_used; dof++) {
2474 ax = xvec[dof];
2475 FOR_ALL_MAT_COLS(REAL, a->matrix_row[dof],
2476 if (mask == NULL || mask->vec[col_dof] < DIRICHLET) {
2477 yvec[col_dof] += ax * row->entry[col_idx];
2478 });
2479 }
2480 } else {
2481 ERROR_EXIT("transpose=%d\n", transpose);
2482 }
2483 }
2484
2485 FLATTEN_ATTR
dof_mv(MatrixTranspose transpose,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,DOF_REAL_VEC * y)2486 void dof_mv(MatrixTranspose transpose,
2487 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
2488 const DOF_REAL_VEC *x, DOF_REAL_VEC *y)
2489 {
2490 const DOF_MATRIX *A_chain;
2491
2492 if (transpose == NoTranspose) {
2493 COL_CHAIN_DO(A, const DOF_MATRIX) {
2494 if (mask) {
2495 __dof_mv(transpose, A, mask, x, y);
2496 } else {
2497 __dof_mv(transpose, A, NULL, x, y);
2498 }
2499 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2500 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2501 if (mask) {
2502 __dof_gemv(transpose, 1.0, A_chain, mask, x, 1.0, y);
2503 } else {
2504 __dof_gemv(transpose, 1.0, A_chain, NULL, x, 1.0, y);
2505 }
2506 }
2507 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2508 y = CHAIN_NEXT(y, DOF_REAL_VEC);
2509 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2510 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
2511 } else {
2512 ROW_CHAIN_DO(A, const DOF_MATRIX) {
2513 if (mask) {
2514 __dof_mv(transpose, A, mask, x, y);
2515 } else {
2516 __dof_mv(transpose, A, NULL, x, y);
2517 }
2518 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2519 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2520 if (mask) {
2521 __dof_gemv(transpose, 1.0, A_chain, mask, x, 1.0, y);
2522 } else {
2523 __dof_gemv(transpose, 1.0, A_chain, NULL, x, 1.0, y);
2524 }
2525 }
2526 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2527 y = CHAIN_NEXT(y, DOF_REAL_VEC);
2528 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2529 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
2530 }
2531 }
2532
2533 /* >>> */
2534
2535 /* <<< REAL_D x REAL */
2536
2537 /* <<< fixed stride */
2538
2539 /* Matrix-Vector for DOF_RDR_MATRIX matrices. */
2540 FORCE_INLINE_ATTR
2541 static inline
__dof_gemv_rdr(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,REAL beta,DOF_REAL_D_VEC * y)2542 void __dof_gemv_rdr(MatrixTranspose transpose,
2543 REAL alpha,
2544 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2545 const DOF_REAL_VEC *x, REAL beta, DOF_REAL_D_VEC *y)
2546 {
2547 FUNCNAME("dof_gemv_rdr");
2548 int ysize, dof;
2549 REAL_D sum;
2550 REAL *xvec;
2551 REAL_D *yvec;
2552 const DOF_ADMIN *m_admin, *y_admin, *x_admin;
2553
2554 MV_TEST_AND_GET_ADMINS(a, x, y, m_admin, x_admin, y_admin);
2555
2556 DEBUG_TEST_EXIT(a->type == MATENT_REAL || a->type == MATENT_REAL_D,
2557 "incompatible block-matrix type");
2558
2559 xvec = x->vec;
2560 yvec = y->vec;
2561
2562 ysize = y->size;
2563
2564 FOR_ALL_FREE_DOFS(y_admin,
2565 if (dof < ysize) {
2566 SET_DOW(0.0, yvec[dof]);
2567 });
2568
2569 if (a->is_diagonal) {
2570 if (a->type == MATENT_REAL_D) {
2571 if (x_admin == y_admin) {
2572 FOR_ALL_DOFS(m_admin, {
2573 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2574 AXPBY_DOW(beta, yvec[dof],
2575 alpha * xvec[dof], a->diagonal.real_d->vec[dof], yvec[dof]);
2576 } else {
2577 SCAL_DOW(beta, yvec[dof]);
2578 }
2579 });
2580 } else if (transpose == NoTranspose) {
2581 DOF *col_dofs = a->diag_cols->vec;
2582 FOR_ALL_DOFS(m_admin, {
2583 DOF col_dof = col_dofs[dof];
2584 if (ENTRY_USED(col_dof) &&
2585 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
2586 AXPBY_DOW(beta, yvec[dof],
2587 alpha * xvec[col_dof], a->diagonal.real_d->vec[dof], yvec[dof]);
2588 } else {
2589 SCAL_DOW(beta, yvec[dof]);
2590 }
2591 });
2592 } else {
2593 DOF *col_dofs = a->diag_cols->vec;
2594 FOR_ALL_DOFS(m_admin, {
2595 DOF col_dof = col_dofs[dof];
2596 if (ENTRY_USED(col_dof) &&
2597 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
2598 AXPBY_DOW(beta, yvec[col_dof],
2599 alpha * xvec[dof], a->diagonal.real_d->vec[dof], yvec[col_dof]);
2600 } else {
2601 SCAL_DOW(beta, yvec[col_dof]);
2602 }
2603 });
2604 }
2605 } else {
2606 // Matrix has REAL entries, a block-wise scalar matrix
2607 if (x_admin == y_admin) {
2608 FOR_ALL_DOFS(m_admin, {
2609 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2610 DMDMSCMAXPBY_DOW(beta, yvec[dof],
2611 alpha * xvec[dof], a->diagonal.real->vec[dof],
2612 yvec[dof]);
2613 } else {
2614 SCAL_DOW(beta, yvec[dof]);
2615 }
2616 });
2617 } else if (transpose == NoTranspose) {
2618 DOF *col_dofs = a->diag_cols->vec;
2619 FOR_ALL_DOFS(m_admin, {
2620 DOF col_dof = col_dofs[dof];
2621 if (ENTRY_USED(col_dof) &&
2622 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
2623 DMDMSCMAXPBY_DOW(beta, yvec[dof],
2624 alpha * xvec[col_dof], a->diagonal.real->vec[dof],
2625 yvec[dof]);
2626 } else {
2627 SCAL_DOW(beta, yvec[dof]);
2628 }
2629 });
2630 } else {
2631 DOF *col_dofs = a->diag_cols->vec;
2632 FOR_ALL_DOFS(m_admin, {
2633 DOF col_dof = col_dofs[dof];
2634 if (ENTRY_USED(col_dof) &&
2635 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
2636 DMDMSCMAXPBY_DOW(beta, yvec[col_dof],
2637 alpha * xvec[dof], a->diagonal.real->vec[dof],
2638 yvec[col_dof]);
2639 } else {
2640 SCAL_DOW(beta, yvec[col_dof]);
2641 }
2642 });
2643 }
2644 }
2645 } else if (transpose == NoTranspose) {
2646 TEST_EXIT(m_admin == y_admin,
2647 "matrix- and y-admins do not match: %p %p.\n",
2648 m_admin, y_admin);
2649
2650 /* Handle some very special cases efficiently (s.t. the compiler
2651 * can optimize them).
2652 */
2653 #define RDR_GEMV_LOOP(mat_type, mat_pfx, add_insn) \
2654 for (dof = 0; dof < m_admin->size_used; dof++) { \
2655 SET_DOW(0.0, sum); \
2656 if (mask == NULL || mask->vec[dof] < DIRICHLET) { \
2657 FOR_ALL_MAT_COLS(mat_type, a->matrix_row[dof], \
2658 DM##mat_pfx##AXPY_DOW(xvec[col_dof], \
2659 row->entry[col_idx], sum)); \
2660 } \
2661 add_insn; \
2662 }
2663 #define RDR_GEMV_CASES(type, pfx) \
2664 if (beta == 0.0) { \
2665 if (alpha == 1.0) { \
2666 RDR_GEMV_LOOP(type, pfx, COPY_DOW(sum, yvec[dof])); \
2667 } else { \
2668 RDR_GEMV_LOOP(type, pfx, AXEY_DOW(alpha, sum, yvec[dof])); \
2669 } \
2670 } else if (beta == 1.0) { \
2671 if (alpha == 1.0) { \
2672 RDR_GEMV_LOOP(type, pfx, AXPY_DOW(1.0, sum, yvec[dof])); \
2673 } else { \
2674 RDR_GEMV_LOOP(type, pfx, AXPY_DOW(alpha, sum, yvec[dof])); \
2675 } \
2676 } else { \
2677 if (alpha == 1.0) { \
2678 RDR_GEMV_LOOP(type, pfx, \
2679 AXPBY_DOW(beta, yvec[dof], 1.0, sum, yvec[dof])); \
2680 } else { \
2681 RDR_GEMV_LOOP(type, pfx, \
2682 AXPBY_DOW(beta, yvec[dof], alpha, sum, yvec[dof])); \
2683 } \
2684 }
2685
2686 if (a->type == MATENT_REAL_D) {
2687 RDR_GEMV_CASES(REAL_D, DM);
2688 } else {
2689 RDR_GEMV_CASES(REAL, SCM);
2690 }
2691 #undef RDR_GEMV_LOOP
2692 #undef RDR_GEMV_CASES
2693 } else { /* Transposed operation */
2694 TEST_EXIT(m_admin == x_admin,
2695 "matrix- and x-admins do not match: %p %p (Transpose).\n",
2696 m_admin, x_admin);
2697
2698 #define RDR_GEMV_LOOP(mat_type, mat_pfx, add_insn) \
2699 for (dof = 0; dof < m_admin->size_used; dof++) { \
2700 REAL xval = xvec[dof]; \
2701 FOR_ALL_MAT_COLS(mat_type, a->matrix_row[dof], \
2702 if (mask == NULL || mask->vec[col_dof] < DIRICHLET) { \
2703 DM##mat_pfx##add_insn; \
2704 }); \
2705 }
2706 #define RDR_GEMV_CASES(type, pfx) \
2707 if (beta == 1.0) { \
2708 /* y is unchanged, just add to it */ \
2709 } else if (beta == 0.0) { \
2710 FOR_ALL_DOFS(y_admin, SET_DOW(0.0, yvec[dof])); \
2711 } else { \
2712 FOR_ALL_DOFS(y_admin, SCAL_DOW(beta, yvec[dof])); \
2713 } \
2714 if (alpha == 1.0) { \
2715 RDR_GEMV_LOOP(type, pfx, \
2716 AXPY_DOW(xval, row->entry[col_idx], yvec[col_dof])); \
2717 } else { \
2718 RDR_GEMV_LOOP(type, pfx, \
2719 AXPY_DOW(alpha*xval, row->entry[col_idx], \
2720 yvec[col_dof])); \
2721 }
2722
2723 if (a->type == MATENT_REAL_D) {
2724 RDR_GEMV_CASES(REAL_D, DM);
2725 } else {
2726 RDR_GEMV_CASES(REAL, SCM);
2727 }
2728 #undef RDR_GEMV_LOOP
2729 #undef RDR_GEMV_CASES
2730 }
2731 }
2732
2733 FLATTEN_ATTR
dof_gemv_rdr(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,REAL beta,DOF_REAL_D_VEC * y)2734 void dof_gemv_rdr(MatrixTranspose transpose,
2735 REAL alpha,
2736 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
2737 const DOF_REAL_VEC *x, REAL beta, DOF_REAL_D_VEC *y)
2738 {
2739 const DOF_MATRIX *A_chain;
2740
2741 if (transpose == NoTranspose) {
2742 COL_CHAIN_DO(A, const DOF_MATRIX) {
2743 if (mask) {
2744 __dof_gemv_rdr(transpose, alpha, A, mask, x, beta, y);
2745 } else {
2746 __dof_gemv_rdr(transpose, alpha, A, NULL, x, beta, y);
2747 }
2748 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2749 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2750 if (mask) {
2751 __dof_gemv_rdr(transpose, alpha, A_chain, mask, x, 1.0, y);
2752 } else {
2753 __dof_gemv_rdr(transpose, alpha, A_chain, NULL, x, 1.0, y);
2754 }
2755 }
2756 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2757 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
2758 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2759 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
2760 } else {
2761 ROW_CHAIN_DO(A, const DOF_MATRIX) {
2762 if (mask) {
2763 __dof_gemv_rdr(transpose, alpha, A, mask, x, beta, y);
2764 } else {
2765 __dof_gemv_rdr(transpose, alpha, A, NULL, x, beta, y);
2766 }
2767 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2768 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2769 if (mask) {
2770 __dof_gemv_rdr(transpose, alpha, A_chain, mask, x, 1.0, y);
2771 } else {
2772 __dof_gemv_rdr(transpose, alpha, A_chain, NULL, x, 1.0, y);
2773 }
2774 }
2775 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2776 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
2777 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2778 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
2779 }
2780 }
2781
2782 FORCE_INLINE_ATTR
2783 static inline
__dof_mv_rdr(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,DOF_REAL_D_VEC * y)2784 void __dof_mv_rdr(MatrixTranspose transpose,
2785 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2786 const DOF_REAL_VEC *x,
2787 DOF_REAL_D_VEC *y)
2788 {
2789 /* The compiler should be clever enough now ... (see abbove) */
2790 __dof_gemv_rdr(transpose, 1.0, a, mask, x, 0.0, y);
2791 }
2792
2793 FLATTEN_ATTR
dof_mv_rdr(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,DOF_REAL_D_VEC * y)2794 void dof_mv_rdr(MatrixTranspose transpose,
2795 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2796 const DOF_REAL_VEC *x, DOF_REAL_D_VEC *y)
2797 {
2798 /* The compiler should be clever enough now ... (see abbove) */
2799 dof_gemv_rdr(transpose, 1.0, a, mask, x, 0.0, y);
2800 }
2801
2802 /* >>> */
2803
2804 /* <<< variable stride */
2805
2806 FORCE_INLINE_ATTR
2807 static inline
__dof_gemv_dow_scl(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,REAL beta,DOF_REAL_VEC_D * y)2808 void __dof_gemv_dow_scl(MatrixTranspose transpose, REAL alpha,
2809 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2810 const DOF_REAL_VEC *x, REAL beta, DOF_REAL_VEC_D *y)
2811 {
2812 if (y->stride != 1) {
2813 __dof_gemv_rdr(transpose, alpha, a, mask, x, beta, (DOF_REAL_D_VEC *)y);
2814 } else {
2815 __dof_gemv(transpose, alpha, a, mask, x, beta, (DOF_REAL_VEC *)y);
2816 }
2817 }
2818
2819 FLATTEN_ATTR
dof_gemv_dow_scl(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,REAL beta,DOF_REAL_VEC_D * y)2820 void dof_gemv_dow_scl(MatrixTranspose transpose, REAL alpha,
2821 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
2822 const DOF_REAL_VEC *x,
2823 REAL beta, DOF_REAL_VEC_D *y)
2824 {
2825 const DOF_MATRIX *A_chain;
2826
2827 if (transpose == NoTranspose) {
2828 COL_CHAIN_DO(A, const DOF_MATRIX) {
2829 if (mask) {
2830 __dof_gemv_dow_scl(transpose, alpha, A, mask, x, beta, y);
2831 } else {
2832 __dof_gemv_dow_scl(transpose, alpha, A, NULL, x, beta, y);
2833 }
2834 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2835 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2836 if (mask) {
2837 __dof_gemv_dow_scl(transpose, alpha, A_chain, mask, x, 1.0, y);
2838 } else {
2839 __dof_gemv_dow_scl(transpose, alpha, A_chain, NULL, x, 1.0, y);
2840 }
2841 }
2842 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2843 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2844 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2845 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
2846 } else {
2847 ROW_CHAIN_DO(A, const DOF_MATRIX) {
2848 if (mask) {
2849 __dof_gemv_dow_scl(transpose, alpha, A, mask, x, beta, y);
2850 } else {
2851 __dof_gemv_dow_scl(transpose, alpha, A, NULL, x, beta, y);
2852 }
2853 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2854 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2855 if (mask) {
2856 __dof_gemv_dow_scl(transpose, alpha, A_chain, mask, x, 1.0, y);
2857 } else {
2858 __dof_gemv_dow_scl(transpose, alpha, A_chain, NULL, x, 1.0, y);
2859 }
2860 }
2861 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2862 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2863 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2864 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
2865 }
2866 }
2867
2868 FORCE_INLINE_ATTR
2869 static inline
__dof_mv_dow_scl(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,DOF_REAL_VEC_D * y)2870 void __dof_mv_dow_scl(MatrixTranspose transpose,
2871 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2872 const DOF_REAL_VEC *x,
2873 DOF_REAL_VEC_D *y)
2874 {
2875 if (y->stride != 1) {
2876 __dof_mv_rdr(transpose, a, mask, x, (DOF_REAL_D_VEC *)y);
2877 } else {
2878 __dof_mv(transpose, a, mask, x, (DOF_REAL_VEC *)y);
2879 }
2880 }
2881
2882 FLATTEN_ATTR
dof_mv_dow_scl(MatrixTranspose transpose,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC * x,DOF_REAL_VEC_D * y)2883 void dof_mv_dow_scl(MatrixTranspose transpose,
2884 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
2885 const DOF_REAL_VEC *x, DOF_REAL_VEC_D *y)
2886 {
2887 const DOF_MATRIX *A_chain;
2888
2889 if (transpose == NoTranspose) {
2890 COL_CHAIN_DO(A, const DOF_MATRIX) {
2891 if (mask) {
2892 __dof_mv_dow_scl(transpose, A, mask, x, y);
2893 } else {
2894 __dof_mv_dow_scl(transpose, A, NULL, x, y);
2895 }
2896 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2897 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2898 if (mask) {
2899 __dof_gemv_dow_scl(transpose, 1.0, A_chain, mask, x, 1.0, y);
2900 } else {
2901 __dof_gemv_dow_scl(transpose, 1.0, A_chain, NULL, x, 1.0, y);
2902 }
2903 }
2904 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2905 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2906 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2907 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
2908 } else {
2909 ROW_CHAIN_DO(A, const DOF_MATRIX) {
2910 if (mask) {
2911 __dof_mv_dow_scl(transpose, A, mask, x, y);
2912 } else {
2913 __dof_mv_dow_scl(transpose, A, NULL, x, y);
2914 }
2915 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
2916 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2917 if (mask) {
2918 __dof_gemv_dow_scl(transpose, 1.0, A_chain, mask, x, 1.0, y);
2919 } else {
2920 __dof_gemv_dow_scl(transpose, 1.0, A_chain, NULL, x, 1.0, y);
2921 }
2922 }
2923 x = CHAIN_NEXT(x, const DOF_REAL_VEC);
2924 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
2925 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
2926 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
2927 }
2928 }
2929
2930 /* >>> */
2931
2932 /* >>> */
2933
2934 /* <<< REAL x REAL_D */
2935
2936 /* <<< fixed stride */
2937
2938 FORCE_INLINE_ATTR
2939 static inline
__dof_gemv_rrd(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,REAL beta,DOF_REAL_VEC * y)2940 void __dof_gemv_rrd(MatrixTranspose transpose, REAL alpha,
2941 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
2942 const DOF_REAL_D_VEC *x, REAL beta, DOF_REAL_VEC *y)
2943 {
2944 FUNCNAME("dof_gemv_rrd");
2945 int ysize, dof;
2946 REAL sum;
2947 REAL_D *xvec;
2948 REAL *yvec;
2949 const DOF_ADMIN *m_admin, *y_admin, *x_admin;
2950
2951 MV_TEST_AND_GET_ADMINS(a, x, y, m_admin, x_admin, y_admin);
2952
2953 DEBUG_TEST_EXIT(a->type == MATENT_REAL || a->type == MATENT_REAL_D,
2954 "incompatible block-matrix type");
2955
2956 xvec = x->vec;
2957 yvec = y->vec;
2958
2959 ysize = y->size;
2960
2961 FOR_ALL_FREE_DOFS(y_admin, if (dof < ysize) yvec[dof] = 0.0);
2962
2963 if (a->is_diagonal) {
2964 if (a->type == MATENT_REAL_D) {
2965 if (x_admin == y_admin) {
2966 FOR_ALL_DOFS(m_admin, {
2967 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
2968 yvec[dof] =
2969 beta * yvec[dof] + alpha * SCP_DOW(a->diagonal.real_d->vec[dof], xvec[dof]);
2970 } else {
2971 yvec[dof] *= beta;
2972 }
2973 });
2974 } else if (transpose == NoTranspose) {
2975 DOF *col_dofs = a->diag_cols->vec;
2976 FOR_ALL_DOFS(m_admin, {
2977 DOF col_dof = col_dofs[dof];
2978 if (ENTRY_USED(col_dof) &&
2979 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
2980 yvec[dof] =
2981 beta * yvec[dof]
2982 +
2983 alpha * SCP_DOW(a->diagonal.real_d->vec[dof], xvec[col_dof]);
2984 } else {
2985 yvec[dof] *= beta;
2986 }
2987 });
2988 } else {
2989 DOF *col_dofs = a->diag_cols->vec;
2990 FOR_ALL_DOFS(m_admin, {
2991 DOF col_dof = col_dofs[dof];
2992 if (ENTRY_USED(col_dof) &&
2993 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
2994 yvec[col_dof] =
2995 beta * yvec[col_dof]
2996 +
2997 alpha * SCP_DOW(a->diagonal.real_d->vec[dof], xvec[dof]);
2998 } else {
2999 yvec[dof] *= beta;
3000 }
3001 });
3002 }
3003 } else {
3004 // Matrix is block-wise scalar
3005 if (x_admin == y_admin) {
3006 FOR_ALL_DOFS(m_admin, {
3007 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3008 yvec[dof] =
3009 beta * yvec[dof] + alpha * a->diagonal.real->vec[dof] * SUM_DOW(xvec[dof]);
3010 } else {
3011 yvec[dof] *= beta;
3012 }
3013 });
3014 } else if (transpose == NoTranspose) {
3015 DOF *col_dofs = a->diag_cols->vec;
3016 FOR_ALL_DOFS(m_admin, {
3017 DOF col_dof = col_dofs[dof];
3018 if (ENTRY_USED(col_dof) &&
3019 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3020 yvec[dof] =
3021 beta * yvec[dof]
3022 +
3023 alpha * a->diagonal.real->vec[dof] * SUM_DOW(xvec[col_dof]);
3024 } else {
3025 yvec[dof] *= beta;
3026 }
3027 });
3028 } else {
3029 DOF *col_dofs = a->diag_cols->vec;
3030 FOR_ALL_DOFS(m_admin, {
3031 DOF col_dof = col_dofs[dof];
3032 if (ENTRY_USED(col_dof) &&
3033 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3034 yvec[col_dof] =
3035 beta * yvec[col_dof]
3036 +
3037 alpha * a->diagonal.real->vec[dof] * SUM_DOW(xvec[dof]);
3038 } else {
3039 yvec[dof] *= beta;
3040 }
3041 });
3042 }
3043 }
3044 } else if (transpose == NoTranspose) { // not diagonal
3045
3046 TEST_EXIT(m_admin == y_admin,
3047 "matrix- and y-admins do not match: %p %p.\n",
3048 m_admin, y_admin);
3049
3050 /* Handle some very special cases efficiently (s.t. the compiler can
3051 * optimize them).
3052 */
3053 #define RRD_GEMV_LOOP(mat_type, MV_INSN, ADD_INSN) \
3054 for (dof = 0; dof < m_admin->size_used; dof++) { \
3055 sum = 0.0; \
3056 if (mask == NULL || mask->vec[dof] < DIRICHLET) { \
3057 FOR_ALL_MAT_COLS(mat_type, a->matrix_row[dof], \
3058 sum += MV_INSN(row->entry[col_idx], \
3059 xvec[col_dof])); \
3060 } \
3061 ADD_INSN; \
3062 }
3063 #define RRD_GEMV_CASES(mat_type, MV_INSN) \
3064 if (beta == 0.0) { \
3065 if (alpha == 1.0) { \
3066 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[dof] = sum); \
3067 } else { \
3068 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[dof] = alpha*sum); \
3069 } \
3070 } else if (beta == 1.0) { \
3071 if (alpha == 1.0) { \
3072 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[dof] += sum); \
3073 } else { \
3074 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[dof] += alpha*sum); \
3075 } \
3076 } else { \
3077 if (alpha == 1.0) { \
3078 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[dof] = beta * yvec[dof] + sum); \
3079 } else { \
3080 RRD_GEMV_LOOP(mat_type, MV_INSN, \
3081 yvec[dof] = beta * yvec[dof] + alpha * sum); \
3082 } \
3083 }
3084
3085 if (a->type == MATENT_REAL_D) {
3086 RRD_GEMV_CASES(REAL_D, SCP_DOW);
3087 } else {
3088 #define MV_INSN(a, b) a * SUM_DOW(b)
3089 RRD_GEMV_CASES(REAL, MV_INSN);
3090 #undef MV_INSN
3091 }
3092 #undef RRD_GEMV_LOOP
3093 #undef RRD_GEMV_CASES
3094 } else { /* transposed matrix operation */
3095 TEST_EXIT(m_admin == x_admin,
3096 "matrix- and x-admins do not match: %p %p (Transpose).\n",
3097 m_admin, x_admin);
3098
3099 #define RRD_GEMV_LOOP(mat_type, MV_INSN, ADD_INSN) \
3100 for (dof = 0; dof < m_admin->size_used; dof++) { \
3101 REAL *xval = xvec[dof]; \
3102 REAL Mx; \
3103 FOR_ALL_MAT_COLS(mat_type, a->matrix_row[dof], \
3104 if (mask == NULL || mask->vec[col_dof] < DIRICHLET) { \
3105 Mx = MV_INSN(row->entry[col_idx], xval); \
3106 ADD_INSN; \
3107 }); \
3108 }
3109 #define RRD_GEMV_CASES(mat_type, MV_INSN) \
3110 if (beta == 1.0) { \
3111 /* y is unchanged */ \
3112 } else if (beta == 0.0) { \
3113 FOR_ALL_DOFS(y_admin, yvec[dof] = 0.0); \
3114 } else { \
3115 FOR_ALL_DOFS(y_admin, yvec[dof] *= beta); \
3116 } \
3117 if (alpha == 1.0) { \
3118 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[col_dof] += Mx); \
3119 } else { \
3120 RRD_GEMV_LOOP(mat_type, MV_INSN, yvec[col_dof] += alpha*Mx); \
3121 }
3122
3123 if (a->type == MATENT_REAL_D) {
3124 RRD_GEMV_CASES(REAL_D, SCP_DOW);
3125 } else {
3126 #define MV_INSN(a, b) a * SUM_DOW(b)
3127 RRD_GEMV_CASES(REAL, MV_INSN);
3128 #undef MV_INSN
3129 }
3130 #undef RRD_GEMV_LOOP
3131 #undef RRD_GEMV_CASES
3132 }
3133 }
3134
3135 FLATTEN_ATTR
dof_gemv_rrd(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,REAL beta,DOF_REAL_VEC * y)3136 void dof_gemv_rrd(MatrixTranspose transpose, REAL alpha,
3137 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3138 const DOF_REAL_D_VEC *x, REAL beta, DOF_REAL_VEC *y)
3139 {
3140 const DOF_MATRIX *A_chain;
3141
3142 if (transpose == NoTranspose) {
3143 COL_CHAIN_DO(A, const DOF_MATRIX) {
3144 if (mask) {
3145 __dof_gemv_rrd(transpose, alpha, A, mask, x, beta, y);
3146 } else {
3147 __dof_gemv_rrd(transpose, alpha, A, NULL, x, beta, y);
3148 }
3149 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3150 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3151 if (mask) {
3152 __dof_gemv_rrd(transpose, alpha, A_chain, mask, x, 1.0, y);
3153 } else {
3154 __dof_gemv_rrd(transpose, alpha, A_chain, NULL, x, 1.0, y);
3155 }
3156 }
3157 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3158 y = CHAIN_NEXT(y, DOF_REAL_VEC);
3159 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3160 } else {
3161 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3162 if (mask) {
3163 __dof_gemv_rrd(transpose, alpha, A, mask, x, beta, y);
3164 } else {
3165 __dof_gemv_rrd(transpose, alpha, A, NULL, x, beta, y);
3166 }
3167 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3168 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3169 if (mask) {
3170 __dof_gemv_rrd(transpose, alpha, A_chain, mask, x, 1.0, y);
3171 } else {
3172 __dof_gemv_rrd(transpose, alpha, A_chain, NULL, x, 1.0, y);
3173 }
3174 }
3175 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3176 y = CHAIN_NEXT(y, DOF_REAL_VEC);
3177 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3178 }
3179 }
3180
3181 FORCE_INLINE_ATTR
3182 static inline
__dof_mv_rrd(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,DOF_REAL_VEC * y)3183 void __dof_mv_rrd(MatrixTranspose transpose,
3184 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3185 const DOF_REAL_D_VEC *x,
3186 DOF_REAL_VEC *y)
3187 {
3188 /* The compiler should be clever enough now ... (see abbove) */
3189 __dof_gemv_rrd(transpose, 1.0, a, mask, x, 0.0, y);
3190 }
3191
3192 FLATTEN_ATTR
dof_mv_rrd(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,DOF_REAL_VEC * y)3193 void dof_mv_rrd(MatrixTranspose transpose,
3194 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3195 const DOF_REAL_D_VEC *x, DOF_REAL_VEC *y)
3196 {
3197 /* The compiler should be clever enough now ... (see abbove) */
3198 dof_gemv_rrd(transpose, 1.0, a, mask, x, 0.0, y);
3199 }
3200
3201 /* >>> */
3202
3203 /* <<< variable stride */
3204
3205 FORCE_INLINE_ATTR
3206 static inline
__dof_gemv_scl_dow(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,REAL beta,DOF_REAL_VEC * y)3207 void __dof_gemv_scl_dow(MatrixTranspose transpose, REAL alpha,
3208 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3209 const DOF_REAL_VEC_D *x,
3210 REAL beta, DOF_REAL_VEC *y)
3211 {
3212 if (x->stride != 1) {
3213 __dof_gemv_rrd(transpose,
3214 alpha, a, mask, (const DOF_REAL_D_VEC *)x, beta, y);
3215 } else {
3216 __dof_gemv(transpose,
3217 alpha, a, mask, (const DOF_REAL_VEC *)x, beta, y);
3218 }
3219 }
3220
3221 FLATTEN_ATTR
dof_gemv_scl_dow(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,REAL beta,DOF_REAL_VEC * y)3222 void dof_gemv_scl_dow(MatrixTranspose transpose, REAL alpha,
3223 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3224 const DOF_REAL_VEC_D *x, REAL beta, DOF_REAL_VEC *y)
3225 {
3226 const DOF_MATRIX *A_chain;
3227
3228 if (transpose == NoTranspose) {
3229 COL_CHAIN_DO(A, const DOF_MATRIX) {
3230 if (mask) {
3231 __dof_gemv_scl_dow(transpose, alpha, A, mask, x, beta, y);
3232 } else {
3233 __dof_gemv_scl_dow(transpose, alpha, A, NULL, x, beta, y);
3234 }
3235 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3236 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3237 if (mask) {
3238 __dof_gemv_scl_dow(transpose, alpha, A_chain, mask, x, 1.0, y);
3239 } else {
3240 __dof_gemv_scl_dow(transpose, alpha, A_chain, NULL, x, 1.0, y);
3241 }
3242 }
3243 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3244 y = CHAIN_NEXT(y, DOF_REAL_VEC);
3245 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3246 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3247 } else {
3248 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3249 if (mask) {
3250 __dof_gemv_scl_dow(transpose, alpha, A, mask, x, beta, y);
3251 } else {
3252 __dof_gemv_scl_dow(transpose, alpha, A, NULL, x, beta, y);
3253 }
3254 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3255 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3256 if (mask) {
3257 __dof_gemv_scl_dow(transpose, alpha, A_chain, mask, x, 1.0, y);
3258 } else {
3259 __dof_gemv_scl_dow(transpose, alpha, A_chain, NULL, x, 1.0, y);
3260 }
3261 }
3262 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3263 y = CHAIN_NEXT(y, DOF_REAL_VEC);
3264 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3265 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3266 }
3267 }
3268
3269 FORCE_INLINE_ATTR
3270 static inline
__dof_mv_scl_dow(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,DOF_REAL_VEC * y)3271 void __dof_mv_scl_dow(MatrixTranspose transpose,
3272 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3273 const DOF_REAL_VEC_D *x, DOF_REAL_VEC *y)
3274 {
3275 if (x->stride != 1) {
3276 __dof_mv_rrd(transpose, a, mask, (const DOF_REAL_D_VEC *)x, y);
3277 } else {
3278 __dof_mv(transpose, a, mask, (const DOF_REAL_VEC *)x, y);
3279 }
3280 }
3281
3282 FLATTEN_ATTR
dof_mv_scl_dow(MatrixTranspose transpose,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,DOF_REAL_VEC * y)3283 void dof_mv_scl_dow(MatrixTranspose transpose,
3284 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3285 const DOF_REAL_VEC_D *x, DOF_REAL_VEC *y)
3286 {
3287 const DOF_MATRIX *A_chain;
3288
3289 if (transpose == NoTranspose) {
3290 COL_CHAIN_DO(A, const DOF_MATRIX) {
3291 if (mask) {
3292 __dof_mv_scl_dow(transpose, A, mask, x, y);
3293 } else {
3294 __dof_mv_scl_dow(transpose, A, NULL, x, y);
3295 }
3296 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3297 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3298 if (mask) {
3299 __dof_gemv_scl_dow(transpose, 1.0, A_chain, mask, x, 1.0, y);
3300 } else {
3301 __dof_gemv_scl_dow(transpose, 1.0, A_chain, NULL, x, 1.0, y);
3302 }
3303 }
3304 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3305 y = CHAIN_NEXT(y, DOF_REAL_VEC);
3306 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3307 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3308 } else {
3309 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3310 if (mask) {
3311 __dof_mv_scl_dow(transpose, A, mask, x, y);
3312 } else {
3313 __dof_mv_scl_dow(transpose, A, NULL, x, y);
3314 }
3315 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3316 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3317 if (mask) {
3318 __dof_gemv_scl_dow(transpose, 1.0, A_chain, mask, x, 1.0, y);
3319 } else {
3320 __dof_gemv_scl_dow(transpose, 1.0, A_chain, NULL, x, 1.0, y);
3321 }
3322 }
3323 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3324 y = CHAIN_NEXT(y, DOF_REAL_VEC);
3325 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3326 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3327 }
3328 }
3329
3330 /* >>> */
3331
3332 /* >>> */
3333
3334 /* <<< REAL_D x REAL_D */
3335
3336 /* <<< fixed stride */
3337
3338 FORCE_INLINE_ATTR
3339 static inline
__dof_gemv_d(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,REAL beta,DOF_REAL_D_VEC * y)3340 void __dof_gemv_d(MatrixTranspose transpose, REAL alpha,
3341 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3342 const DOF_REAL_D_VEC *x, REAL beta, DOF_REAL_D_VEC *y)
3343 {
3344 FUNCNAME("dof_gemv_d");
3345 int dof, ysize;
3346 REAL_D sum, ax;
3347 REAL_D *xvec, *yvec;
3348 const DOF_ADMIN *m_admin, *y_admin, *x_admin;
3349
3350 MV_TEST_AND_GET_ADMINS(a, x, y, m_admin, x_admin, y_admin);
3351
3352 xvec = x->vec;
3353 yvec = y->vec;
3354
3355 ysize = y->size;
3356 FOR_ALL_FREE_DOFS(y_admin, if (dof < ysize) SET_DOW(0.0, yvec[dof]));
3357
3358 if (a->is_diagonal) {
3359 switch (a->type) {
3360 case MATENT_REAL:
3361 if (x_admin == y_admin) {
3362 FOR_ALL_DOFS(m_admin, {
3363 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3364 SCMGEMV_DOW(alpha, a->diagonal.real->vec[dof], xvec[dof],
3365 beta, yvec[dof]);
3366 } else {
3367 SCAL_DOW(beta, yvec[dof]);
3368 }
3369 });
3370 } else if (transpose == NoTranspose) {
3371 DOF *col_dofs = a->diag_cols->vec;
3372 FOR_ALL_DOFS(m_admin, {
3373 DOF col_dof = col_dofs[dof];
3374 if (ENTRY_USED(col_dof) &&
3375 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3376 SCMGEMV_DOW(alpha, a->diagonal.real->vec[dof], xvec[col_dof],
3377 beta, yvec[dof]);
3378 } else {
3379 SCAL_DOW(beta, yvec[dof]);
3380 }
3381 });
3382 } else {
3383 DOF *col_dofs = a->diag_cols->vec;
3384 FOR_ALL_DOFS(m_admin, {
3385 DOF col_dof = col_dofs[dof];
3386 if (ENTRY_USED(col_dof) &&
3387 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3388 SCMGEMV_DOW(alpha, a->diagonal.real->vec[dof], xvec[dof],
3389 beta, yvec[col_dof]);
3390 } else {
3391 SCAL_DOW(beta, yvec[col_dof]);
3392 }
3393 });
3394 }
3395 break;
3396 case MATENT_REAL_D:
3397 if (x_admin == y_admin) {
3398 FOR_ALL_DOFS(m_admin, {
3399 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3400 DMGEMV_DOW(alpha, a->diagonal.real_d->vec[dof], xvec[dof],
3401 beta, yvec[dof]);
3402 } else {
3403 SCAL_DOW(beta, yvec[dof]);
3404 }
3405 });
3406 } else if (transpose == NoTranspose) {
3407 DOF *col_dofs = a->diag_cols->vec;
3408 FOR_ALL_DOFS(m_admin, {
3409 DOF col_dof = col_dofs[dof];
3410 if (ENTRY_USED(col_dof) &&
3411 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3412 DMGEMV_DOW(alpha,
3413 a->diagonal.real_d->vec[dof], xvec[col_dof],
3414 beta, yvec[dof]);
3415 } else {
3416 SCAL_DOW(beta, yvec[dof]);
3417 }
3418 });
3419 } else {
3420 DOF *col_dofs = a->diag_cols->vec;
3421 FOR_ALL_DOFS(m_admin, {
3422 DOF col_dof = col_dofs[dof];
3423 if (ENTRY_USED(col_dof) &&
3424 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3425 DMGEMV_DOW(alpha,
3426 a->diagonal.real_d->vec[dof], xvec[dof],
3427 beta, yvec[col_dof]);
3428 } else {
3429 SCAL_DOW(beta, yvec[col_dof]);
3430 }
3431 });
3432 }
3433 break;
3434 case MATENT_REAL_DD:
3435 if (x_admin == y_admin) {
3436 if (transpose == NoTranspose) {
3437 FOR_ALL_DOFS(m_admin, {
3438 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3439 MGEMV_DOW(alpha, (const REAL_D *)a->diagonal.real_dd->vec[dof],
3440 xvec[dof], beta, yvec[dof]);
3441 } else {
3442 SCAL_DOW(beta, yvec[dof]);
3443 }
3444 });
3445 } else {
3446 FOR_ALL_DOFS(m_admin, {
3447 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3448 MGEMTV_DOW(alpha, (const REAL_D *)a->diagonal.real_dd->vec[dof],
3449 xvec[dof], beta, yvec[dof]);
3450 } else {
3451 SCAL_DOW(beta, yvec[dof]);
3452 }
3453 });
3454 }
3455 } else { /* x_admin != y_admin */
3456 DOF *col_dofs = a->diag_cols->vec;
3457 if (transpose == NoTranspose) {
3458 FOR_ALL_DOFS(m_admin, {
3459 DOF col_dof = col_dofs[dof];
3460 if (ENTRY_USED(col_dof) &&
3461 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3462 MGEMV_DOW(
3463 alpha, (const REAL_D *)a->diagonal.real_dd->vec[dof],
3464 xvec[col_dof], beta, yvec[dof]);
3465 } else {
3466 SCAL_DOW(beta, yvec[dof]);
3467 }
3468 });
3469 } else {
3470 FOR_ALL_DOFS(m_admin, {
3471 DOF col_dof = col_dofs[dof];
3472 if (ENTRY_USED(col_dof) &&
3473 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3474 MGEMTV_DOW(alpha, (const REAL_D *)a->diagonal.real_dd->vec[dof],
3475 xvec[dof], beta, yvec[col_dof]);
3476 } else {
3477 SCAL_DOW(beta, yvec[col_dof]);
3478 }
3479 });
3480 }
3481 }
3482 break;
3483 default:
3484 break;
3485 }
3486 } else if (transpose == NoTranspose) {
3487
3488 TEST_EXIT(m_admin == y_admin,
3489 "matrix- and y-admins do not match: %p %p.\n",
3490 m_admin, y_admin);
3491
3492 #undef MAT_BODY
3493 #define MAT_BODY(F, CC, C, S, TYPE) \
3494 for (dof = 0; dof < m_admin->size_used; dof++) { \
3495 SET_DOW(0.0, sum); \
3496 if (mask == NULL || mask->vec[dof] < DIRICHLET) { \
3497 FOR_ALL_MAT_COLS(TYPE, a->matrix_row[dof], \
3498 F##V_DOW(CC row->entry[col_idx], \
3499 xvec[col_dof], sum)); \
3500 } \
3501 AXPBY_DOW(alpha, sum, beta, yvec[dof], yvec[dof]); \
3502 }
3503 MAT_EMIT_BODY_SWITCH(a->type);
3504
3505 } else if (transpose == Transpose) {
3506
3507 TEST_EXIT(m_admin == x_admin,
3508 "matrix- and x-admins do not match: %p %p.\n",
3509 m_admin, x_admin);
3510
3511 FOR_ALL_DOFS(y_admin, AX_DOW(beta, yvec[dof]));
3512
3513 #undef MAT_BODY
3514 #define MAT_BODY(F, CC, C, S, TYPE) \
3515 for (dof = 0; dof < m_admin->size_used; dof++) { \
3516 COPY_DOW(xvec[dof], ax); AX_DOW(alpha, ax); \
3517 FOR_ALL_MAT_COLS(TYPE, a->matrix_row[dof], \
3518 if (mask == NULL || mask->vec[col_dof] < DIRICHLET) { \
3519 F##TV_DOW(CC row->entry[col_idx], ax, yvec[col_dof]); \
3520 }); \
3521 }
3522 MAT_EMIT_BODY_SWITCH(a->type);
3523 }
3524 else {
3525 ERROR_EXIT("transpose=%d\n", transpose);
3526 }
3527 }
3528
3529 FLATTEN_ATTR
dof_gemv_d(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,REAL beta,DOF_REAL_D_VEC * y)3530 void dof_gemv_d(MatrixTranspose transpose, REAL alpha,
3531 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3532 const DOF_REAL_D_VEC *x, REAL beta, DOF_REAL_D_VEC *y)
3533 {
3534 const DOF_MATRIX *A_chain;
3535
3536 if (transpose == NoTranspose) {
3537 COL_CHAIN_DO(A, const DOF_MATRIX) {
3538 if (mask) {
3539 __dof_gemv_d(transpose, alpha, A, mask, x, beta, y);
3540 } else {
3541 __dof_gemv_d(transpose, alpha, A, NULL, x, beta, y);
3542 }
3543 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3544 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3545 if (mask) {
3546 __dof_gemv_d(transpose, alpha, A_chain, mask, x, 1.0, y);
3547 } else {
3548 __dof_gemv_d(transpose, alpha, A_chain, NULL, x, 1.0, y);
3549 }
3550 }
3551 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3552 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
3553 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3554 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3555 } else {
3556 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3557 if (mask) {
3558 __dof_gemv_d(transpose, alpha, A, mask, x, beta, y);
3559 } else {
3560 __dof_gemv_d(transpose, alpha, A, NULL, x, beta, y);
3561 }
3562 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3563 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3564 if (mask) {
3565 __dof_gemv_d(transpose, alpha, A_chain, mask, x, 1.0, y);
3566 } else {
3567 __dof_gemv_d(transpose, alpha, A_chain, NULL, x, 1.0, y);
3568 }
3569 }
3570 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3571 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
3572 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3573 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3574 }
3575 }
3576
3577 /* Multiply a DOF_REAL_D_VEC with a DOF_MATRIX, which may consist of
3578 * scalar, diagonal or DOWxDOW blocks.
3579 */
3580 FORCE_INLINE_ATTR
3581 static inline
__dof_mv_d(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)3582 void __dof_mv_d(MatrixTranspose transpose,
3583 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3584 const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
3585 {
3586 FUNCNAME("dof_mv_d");
3587 int j, dof;
3588 REAL_D sum;
3589 REAL_D *xvec, *yvec;
3590 const DOF_ADMIN *m_admin, *y_admin, *x_admin;
3591
3592 MV_TEST_AND_GET_ADMINS(a, x, y, m_admin, x_admin, y_admin);
3593
3594 xvec = x->vec;
3595 yvec = y->vec;
3596
3597 if (a->is_diagonal) {
3598 switch (a->type) {
3599 case MATENT_REAL:
3600 if (x_admin == y_admin) {
3601 FOR_ALL_DOFS(m_admin, {
3602 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3603 SCMVEQ_DOW(a->diagonal.real->vec[dof], xvec[dof], yvec[dof]);
3604 }
3605 });
3606 } else if (transpose == NoTranspose) {
3607 DOF *col_dofs = a->diag_cols->vec;
3608 FOR_ALL_DOFS(m_admin, {
3609 DOF col_dof = col_dofs[dof];
3610 if (ENTRY_USED(col_dof) &&
3611 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3612 SCMVEQ_DOW(a->diagonal.real->vec[dof], xvec[col_dof], yvec[dof]);
3613 }
3614 });
3615 } else {
3616 DOF *col_dofs = a->diag_cols->vec;
3617 FOR_ALL_DOFS(m_admin, {
3618 DOF col_dof = col_dofs[dof];
3619 if (ENTRY_USED(col_dof) &&
3620 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3621 SCMVEQ_DOW(a->diagonal.real->vec[dof], xvec[dof], yvec[col_dof]);
3622 }
3623 });
3624 }
3625 break;
3626 case MATENT_REAL_D:
3627 if (x_admin == y_admin) {
3628 FOR_ALL_DOFS(m_admin, {
3629 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3630 DMVEQ_DOW(a->diagonal.real_d->vec[dof], xvec[dof], yvec[dof]);
3631 }
3632 });
3633 } else if (transpose == NoTranspose) {
3634 DOF *col_dofs = a->diag_cols->vec;
3635 FOR_ALL_DOFS(m_admin, {
3636 DOF col_dof = col_dofs[dof];
3637 if (ENTRY_USED(col_dof) &&
3638 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3639 DMVEQ_DOW(a->diagonal.real_d->vec[dof], xvec[col_dof], yvec[dof]);
3640 }
3641 });
3642 } else {
3643 DOF *col_dofs = a->diag_cols->vec;
3644 FOR_ALL_DOFS(m_admin, {
3645 DOF col_dof = col_dofs[dof];
3646 if (ENTRY_USED(col_dof) &&
3647 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3648 DMVEQ_DOW(a->diagonal.real_d->vec[dof], xvec[dof], yvec[col_dof]);
3649 }
3650 });
3651 }
3652 break;
3653 case MATENT_REAL_DD:
3654 if (x_admin == y_admin) {
3655 if (transpose == NoTranspose) {
3656 FOR_ALL_DOFS(m_admin, {
3657 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3658 MVEQ_DOW((const REAL_D *)a->diagonal.real_dd->vec[dof],
3659 xvec[dof], yvec[dof]);
3660 }
3661 });
3662 } else {
3663 FOR_ALL_DOFS(m_admin, {
3664 if (mask == NULL || mask->vec[dof] < DIRICHLET) {
3665 MTVEQ_DOW((const REAL_D *)a->diagonal.real_dd->vec[dof],
3666 xvec[dof], yvec[dof]);
3667 }
3668 });
3669 }
3670 } else {
3671 DOF *col_dofs = a->diag_cols->vec;
3672 if (transpose == NoTranspose) {
3673 FOR_ALL_DOFS(m_admin, {
3674 DOF col_dof = col_dofs[dof];
3675 if (ENTRY_USED(col_dof) &&
3676 (mask == NULL || mask->vec[dof] < DIRICHLET)) {
3677 MVEQ_DOW((const REAL_D *)a->diagonal.real_dd->vec[dof],
3678 xvec[col_dof], yvec[dof]);
3679 }
3680 });
3681 } else {
3682 FOR_ALL_DOFS(m_admin, {
3683 DOF col_dof = col_dofs[dof];
3684 if (ENTRY_USED(col_dof) &&
3685 (mask == NULL || mask->vec[col_dof] < DIRICHLET)) {
3686 MTVEQ_DOW((const REAL_D *)a->diagonal.real_dd->vec[dof],
3687 xvec[dof], yvec[col_dof]);
3688 }
3689 });
3690 }
3691 }
3692 break;
3693 default:
3694 break;
3695 }
3696 } else if (transpose == NoTranspose) {
3697 int ysize = y->size;
3698
3699 TEST_EXIT(m_admin == y_admin,
3700 "matrix- and y-admins do not match: %p %p.\n",
3701 m_admin, y_admin);
3702
3703 FOR_ALL_FREE_DOFS(y_admin, if (dof < ysize) SET_DOW(0.0, yvec[dof]));
3704
3705 #undef MAT_BODY
3706 #define MAT_BODY(F, CC, C, S, TYPE) \
3707 for (dof = 0; dof < m_admin->size_used; dof++) { \
3708 SET_DOW(0.0, sum); \
3709 if (mask == NULL || mask->vec[dof] < DIRICHLET) { \
3710 FOR_ALL_MAT_COLS(TYPE, a->matrix_row[dof], \
3711 F##V_DOW(CC row->entry[col_idx], \
3712 xvec[col_dof], sum)); \
3713 } \
3714 COPY_DOW(sum, yvec[dof]); \
3715 }
3716 MAT_EMIT_BODY_SWITCH(a->type);
3717 } else if (transpose == Transpose) {
3718
3719 TEST_EXIT(m_admin == x_admin,
3720 "matrix- and x-admins do not match: %p %p.\n",
3721 m_admin, x_admin);
3722
3723 for (j = 0; j < y_admin->size_used; j++) {
3724 SET_DOW(0.0, yvec[j]);
3725 }
3726
3727 #undef MAT_BODY
3728 #define MAT_BODY(F, CC, C, S, TYPE) \
3729 for (dof = 0; dof < m_admin->size_used; dof++) { \
3730 FOR_ALL_MAT_COLS(TYPE, a->matrix_row[dof], \
3731 if (mask == NULL || mask->vec[col_dof] < DIRICHLET) { \
3732 F##TV_DOW(CC row->entry[col_idx], \
3733 xvec[dof], \
3734 yvec[col_dof]); \
3735 }); \
3736 }
3737 MAT_EMIT_BODY_SWITCH(a->type);
3738 } else {
3739 ERROR_EXIT("transpose=%d\n", transpose);
3740 }
3741 }
3742
3743 FLATTEN_ATTR
dof_mv_d(MatrixTranspose transpose,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_D_VEC * x,DOF_REAL_D_VEC * y)3744 void dof_mv_d(MatrixTranspose transpose,
3745 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3746 const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y)
3747 {
3748 const DOF_MATRIX *A_chain;
3749
3750 if (transpose == NoTranspose) {
3751 COL_CHAIN_DO(A, const DOF_MATRIX) {
3752 if (mask) {
3753 __dof_mv_d(transpose, A, mask, x, y);
3754 } else {
3755 __dof_mv_d(transpose, A, NULL, x, y);
3756 }
3757 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3758 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3759 if (mask) {
3760 __dof_gemv_d(transpose, 1.0, A, mask, x, 1.0, y);
3761 } else {
3762 __dof_gemv_d(transpose, 1.0, A, NULL, x, 1.0, y);
3763 }
3764 }
3765 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3766 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
3767 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3768 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3769 } else {
3770 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3771 if (mask) {
3772 __dof_mv_d(transpose, A, mask, x, y);
3773 } else {
3774 __dof_mv_d(transpose, A, NULL, x, y);
3775 }
3776 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3777 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3778 if (mask) {
3779 __dof_gemv_d(transpose, 1.0, A_chain, mask, x, 1.0, y);
3780 } else {
3781 __dof_gemv_d(transpose, 1.0, A_chain, NULL, x, 1.0, y);
3782 }
3783 }
3784 x = CHAIN_NEXT(x, const DOF_REAL_D_VEC);
3785 y = CHAIN_NEXT(y, DOF_REAL_D_VEC);
3786 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3787 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3788 }
3789 }
3790
3791 /* >>> */
3792
3793 /* <<< variable stride */
3794
3795 FORCE_INLINE_ATTR
3796 static inline
__dof_gemv_dow(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,REAL beta,DOF_REAL_VEC_D * y)3797 void __dof_gemv_dow(MatrixTranspose transpose, REAL alpha,
3798 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3799 const DOF_REAL_VEC_D *x, REAL beta, DOF_REAL_VEC_D *y)
3800 {
3801 if (y->stride == 1) {
3802 __dof_gemv_scl_dow(transpose, alpha, a, mask, x, beta, (DOF_REAL_VEC *)y);
3803 } else if (x->stride == 1) {
3804 __dof_gemv_dow_scl(transpose, alpha, a, mask, (const DOF_REAL_VEC *)x,
3805 beta, y);
3806 } else {
3807 __dof_gemv_d(transpose, alpha, a, mask, (const DOF_REAL_D_VEC *)x,
3808 beta, (DOF_REAL_D_VEC *)y);
3809 }
3810 }
3811
3812 FLATTEN_ATTR
dof_gemv_dow(MatrixTranspose transpose,REAL alpha,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,REAL beta,DOF_REAL_VEC_D * y)3813 void dof_gemv_dow(MatrixTranspose transpose, REAL alpha,
3814 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3815 const DOF_REAL_VEC_D *x, REAL beta, DOF_REAL_VEC_D *y)
3816 {
3817 const DOF_MATRIX *A_chain;
3818
3819 if (transpose == NoTranspose) {
3820 COL_CHAIN_DO(A, const DOF_MATRIX) {
3821 if (mask) {
3822 __dof_gemv_dow(transpose, alpha, A, mask, x, beta, y);
3823 } else {
3824 __dof_gemv_dow(transpose, alpha, A, NULL, x, beta, y);
3825 }
3826 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3827 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3828 if (mask) {
3829 __dof_gemv_dow(transpose, alpha, A_chain, mask, x, 1.0, y);
3830 } else {
3831 __dof_gemv_dow(transpose, alpha, A_chain, NULL, x, 1.0, y);
3832 }
3833 }
3834 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3835 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
3836 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3837 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3838 } else {
3839 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3840 if (mask) {
3841 __dof_gemv_dow(transpose, alpha, A, mask, x, beta, y);
3842 } else {
3843 __dof_gemv_dow(transpose, alpha, A, NULL, x, beta, y);
3844 }
3845 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3846 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3847 if (mask) {
3848 __dof_gemv_dow(transpose, alpha, A_chain, mask, x, 1.0, y);
3849 } else {
3850 __dof_gemv_dow(transpose, alpha, A_chain, NULL, x, 1.0, y);
3851 }
3852 }
3853 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3854 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
3855 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3856 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3857 }
3858 }
3859
3860 FORCE_INLINE_ATTR
3861 static inline
__dof_mv_dow(MatrixTranspose transpose,const DOF_MATRIX * a,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,DOF_REAL_VEC_D * y)3862 void __dof_mv_dow(MatrixTranspose transpose,
3863 const DOF_MATRIX *a, const DOF_SCHAR_VEC *mask,
3864 const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y)
3865 {
3866 if (y->stride == 1) {
3867 __dof_mv_scl_dow(transpose, a, mask, x, (DOF_REAL_VEC *)y);
3868 } else if (x->stride == 1) {
3869 __dof_mv_dow_scl(transpose, a, mask, (const DOF_REAL_VEC *)x, y);
3870 } else {
3871 __dof_mv_d(transpose, a, mask,
3872 (const DOF_REAL_D_VEC *)x, (DOF_REAL_D_VEC *)y);
3873 }
3874 }
3875
3876 FLATTEN_ATTR
dof_mv_dow(MatrixTranspose transpose,const DOF_MATRIX * A,const DOF_SCHAR_VEC * mask,const DOF_REAL_VEC_D * x,DOF_REAL_VEC_D * y)3877 void dof_mv_dow(MatrixTranspose transpose,
3878 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
3879 const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y)
3880 {
3881 const DOF_MATRIX *A_chain;
3882
3883 if (transpose == NoTranspose) {
3884 COL_CHAIN_DO(A, const DOF_MATRIX) {
3885 if (mask) {
3886 __dof_mv_dow(transpose, A, mask, x, y);
3887 } else {
3888 __dof_mv_dow(transpose, A, NULL, x, y);
3889 }
3890 ROW_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3891 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3892 if (mask) {
3893 __dof_gemv_dow(transpose, 1.0, A_chain, mask, x, 1.0, y);
3894 } else {
3895 __dof_gemv_dow(transpose, 1.0, A_chain, NULL, x, 1.0, y);
3896 }
3897 }
3898 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3899 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
3900 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3901 } COL_CHAIN_WHILE(A, const DOF_MATRIX);
3902 } else {
3903 ROW_CHAIN_DO(A, const DOF_MATRIX) {
3904 if (mask) {
3905 __dof_mv_dow(transpose, A, mask, x, y);
3906 } else {
3907 __dof_mv_dow(transpose, A, NULL, x, y);
3908 }
3909 COL_CHAIN_FOREACH(A_chain, A, const DOF_MATRIX) {
3910 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3911 if (mask) {
3912 __dof_gemv_dow(transpose, 1.0, A_chain, mask, x, 1.0, y);
3913 } else {
3914 __dof_gemv_dow(transpose, 1.0, A_chain, NULL, x, 1.0, y);
3915 }
3916 }
3917 x = CHAIN_NEXT(x, const DOF_REAL_VEC_D);
3918 y = CHAIN_NEXT(y, DOF_REAL_VEC_D);
3919 mask = mask ? CHAIN_NEXT(mask, DOF_SCHAR_VEC) : NULL;
3920 } ROW_CHAIN_WHILE(A, const DOF_MATRIX);
3921 }
3922 }
3923
3924 /* >>> */
3925
3926 /* >>> */
3927
3928 /*--------------------------------------------------------------------------*/
3929
3930 /*--------------------------------------------------------------------------*/
3931 /* print_dof_matrix: print matrix in compressed format */
3932 /*--------------------------------------------------------------------------*/
3933
3934 static inline
__print_dof_matrix_row_real(const DOF_MATRIX * matrix,int i)3935 void __print_dof_matrix_row_real(const DOF_MATRIX *matrix, int i)
3936 {
3937 FUNCNAME("print_dof_matrix");
3938 int j, jcol;
3939 MATRIX_ROW *row;
3940
3941 if (matrix->is_diagonal) {
3942 if (i < matrix->diagonal.real->size) {
3943 MSG("row %3d: (%d, %.8e)", i, i, matrix->diagonal.real->vec[i]);
3944 }
3945 return;
3946 }
3947 for (row = matrix->matrix_row[i]; row; row = row->next) {
3948 MSG("row %3d:", i);
3949 for (j=0; j<ROW_LENGTH; j++) {
3950 jcol = row->col[j];
3951 if (ENTRY_USED(jcol)) {
3952 print_msg(" (%3d, %.8e)", jcol, row->entry.real[j]);
3953 } else if (jcol == NO_MORE_ENTRIES) {
3954 break;
3955 }
3956 }
3957 print_msg("\n");
3958 if (jcol == NO_MORE_ENTRIES) {
3959 break;
3960 }
3961 }
3962 }
3963
3964 static inline
__print_dof_matrix_real(const DOF_MATRIX * matrix)3965 void __print_dof_matrix_real(const DOF_MATRIX *matrix)
3966 {
3967 FUNCNAME("print_dof_matrix");
3968 int i;
3969
3970 if (matrix->is_diagonal) {
3971 print_dof_real_vec(matrix->diagonal.real);
3972 } else for (i = 0; i < matrix->size; i++) {
3973 __print_dof_matrix_row_real(matrix, i);
3974 }
3975 }
3976
3977 static inline
__print_dof_matrix_row_real_d(const DOF_MATRIX * matrix,int i)3978 void __print_dof_matrix_row_real_d(const DOF_MATRIX *matrix, int i)
3979 {
3980 FUNCNAME("print_dof_rdr_matrix");
3981 int j, jcol;
3982 MATRIX_ROW *row;
3983
3984 if (matrix->is_diagonal) {
3985 if (i < matrix->diagonal.real_d->size) {
3986 MSG("row %3d: (%d, "FORMAT_DOW")\n",
3987 i, i, EXPAND_DOW(matrix->diagonal.real_d->vec[i]));
3988 }
3989 } else for (row = matrix->matrix_row[i]; row; row = row->next) {
3990 MSG("row %3d:",i);
3991 for (j=0; j<ROW_LENGTH; j++) {
3992 jcol = row->col[j];
3993 if (ENTRY_USED(jcol)) {
3994 print_msg(" (%3d, "FORMAT_DOW")",
3995 jcol, EXPAND_DOW(row->entry.real_d[j]));
3996 } else if (jcol == NO_MORE_ENTRIES) {
3997 break;
3998 }
3999 }
4000 print_msg("\n");
4001 if (jcol == NO_MORE_ENTRIES) {
4002 break;
4003 }
4004 }
4005 }
4006
4007 static inline
__print_dof_matrix_real_d(const DOF_MATRIX * matrix)4008 void __print_dof_matrix_real_d(const DOF_MATRIX *matrix)
4009 {
4010 FUNCNAME("print_dof_rdr_matrix");
4011 int i;
4012
4013 if (matrix->is_diagonal) {
4014 print_dof_real_d_vec(matrix->diagonal.real_d);
4015 } else for (i= 0; i < matrix->size; i++) {
4016 __print_dof_matrix_row_real_d(matrix, i);
4017 }
4018 }
4019
4020 static inline
__print_dof_matrix_row_real_dd(const DOF_MATRIX * matrix,int i)4021 void __print_dof_matrix_row_real_dd(const DOF_MATRIX *matrix, int i)
4022 {
4023 FUNCNAME("print_dof_matrix");
4024 int j, jcol, n, m;
4025 MATRIX_ROW *row;
4026
4027 if (matrix->is_diagonal) {
4028 if (i < matrix->diagonal.real_d->size) {
4029 MSG("row %3d: (%d, "MFORMAT_DOW")\n",
4030 i, i, MEXPAND_DOW(matrix->diagonal.real_dd->vec[i]));
4031 }
4032 } else if (matrix->matrix_row[i]) {
4033 for (n = 0; n < DIM_OF_WORLD; n++) {
4034 if (n == 0) {
4035 MSG("row %3d:",i);
4036 } else {
4037 MSG(" ");
4038 }
4039 for (row = matrix->matrix_row[i]; row; row = row->next) {
4040 for (j=0; j<ROW_LENGTH; j++) {
4041 jcol = row->col[j];
4042 if (ENTRY_USED(jcol)) {
4043 if (n == 0) {
4044 print_msg(" |%3d", jcol);
4045 } else {
4046 print_msg(" | ");
4047 }
4048 for (m = 0; m < DIM_OF_WORLD; m++) {
4049 print_msg(" % .2e", row->entry.real_dd[j][n][m]);
4050 }
4051 } else if (jcol == NO_MORE_ENTRIES) {
4052 break;
4053 }
4054 }
4055 if (jcol == NO_MORE_ENTRIES) {
4056 break;
4057 }
4058 }
4059 print_msg("\n");
4060 }
4061 }
4062 }
4063
4064 static inline
__print_dof_matrix_real_dd(const DOF_MATRIX * matrix)4065 void __print_dof_matrix_real_dd(const DOF_MATRIX *matrix)
4066 {
4067 FUNCNAME("print_dof_matrix");
4068 int i;
4069
4070 if (matrix->is_diagonal) {
4071 print_dof_real_dd_vec(matrix->diagonal.real_dd);
4072 } else for (i= 0 ; i < matrix->size; i++) {
4073 __print_dof_matrix_row_real_dd(matrix, i);
4074 }
4075 }
4076
4077 static inline
__print_dof_matrix_row(const DOF_MATRIX * matrix,int row)4078 void __print_dof_matrix_row(const DOF_MATRIX *matrix, int row)
4079 {
4080 FUNCNAME("print_dof_matrix");
4081
4082 switch (matrix->type) {
4083 case MATENT_REAL:
4084 __print_dof_matrix_row_real(matrix, row); break;
4085 case MATENT_REAL_D:
4086 __print_dof_matrix_row_real_d(matrix, row); break;
4087 case MATENT_REAL_DD:
4088 __print_dof_matrix_row_real_dd(matrix, row); break;
4089 case MATENT_NONE:
4090 MSG("Attempt to print uninitialized dof-matrix.");
4091 break;
4092 default:
4093 ERROR_EXIT("Unknown MATENT_TYPE: %d\n", matrix->type);
4094 break;
4095 }
4096 }
4097
4098 static inline
__print_dof_matrix(const DOF_MATRIX * matrix)4099 void __print_dof_matrix(const DOF_MATRIX *matrix)
4100 {
4101 FUNCNAME("print_dof_matrix");
4102
4103 switch (matrix->type) {
4104 case MATENT_REAL:
4105 __print_dof_matrix_real(matrix); break;
4106 case MATENT_REAL_D:
4107 __print_dof_matrix_real_d(matrix); break;
4108 case MATENT_REAL_DD:
4109 __print_dof_matrix_real_dd(matrix); break;
4110 case MATENT_NONE:
4111 MSG("Attempt to print uninitialized dof-matrix.");
4112 break;
4113 default:
4114 ERROR_EXIT("Unknown MATENT_TYPE: %d\n", matrix->type);
4115 break;
4116 }
4117 }
4118
print_dof_matrix(const DOF_MATRIX * matrix)4119 void print_dof_matrix(const DOF_MATRIX *matrix)
4120 {
4121 FUNCNAME("print_dof_matrix");
4122 int i, j;
4123
4124 i = 0;
4125 COL_CHAIN_DO(matrix, const DOF_MATRIX) {
4126 j = 0;
4127 ROW_CHAIN_DO(matrix, const DOF_MATRIX) {
4128 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix)) {
4129 MSG("BLOCK(%d,%d):\n", i , j);
4130 }
4131 __print_dof_matrix(matrix);
4132 ++j;
4133 } ROW_CHAIN_WHILE(matrix, const DOF_MATRIX);
4134 ++i;
4135 } COL_CHAIN_WHILE(matrix, const DOF_MATRIX);
4136 }
4137
print_dof_matrix_row(const DOF_MATRIX * matrix,int i)4138 void print_dof_matrix_row(const DOF_MATRIX *matrix, int i)
4139 {
4140 FUNCNAME("print_dof_matrix_row");
4141 int j;
4142
4143 j = 0;
4144 ROW_CHAIN_DO(matrix, const DOF_MATRIX) {
4145 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix)) {
4146 MSG("BLOCK(%d):\n", j);
4147 }
4148 __print_dof_matrix_row(matrix, i);
4149 ++j;
4150 } ROW_CHAIN_WHILE(matrix, const DOF_MATRIX);
4151 }
4152
4153 static inline
__print_dof_real_vec(const DOF_REAL_VEC * drv)4154 void __print_dof_real_vec(const DOF_REAL_VEC *drv)
4155 {
4156 FUNCNAME("print_dof_real_vec");
4157 int i, j;
4158 const DOF_ADMIN *admin = NULL;
4159 const char *format;
4160
4161 if (drv->fe_space) admin = drv->fe_space->admin;
4162
4163 MSG("Vec `%s':\n", drv->name);
4164 j = 0;
4165 if (admin)
4166 {
4167 if (admin->size_used > 100)
4168 format = "%s(%3d,%10.5le)";
4169 else if (admin->size_used > 10)
4170 format = "%s(%2d,%10.5le)";
4171 else
4172 format = "%s(%1d,%10.5le)";
4173
4174 FOR_ALL_DOFS(admin,
4175 if ((j % 3) == 0) {
4176 if (j) print_msg("\n");
4177 MSG(format, "", dof, drv->vec[dof]);
4178 }
4179 else
4180 print_msg(format, " ", dof, drv->vec[dof]);
4181 j++;
4182 );
4183 print_msg("\n");
4184 }
4185 else
4186 {
4187 MSG("no DOF_ADMIN, print whole vector.\n");
4188
4189 for (i = 0; i < drv->size; i++) {
4190 if ((j % 3) == 0)
4191 {
4192 if (j) print_msg("\n");
4193 MSG("(%d,%10.5le)",i,drv->vec[i]);
4194 }
4195 else
4196 print_msg(" (%d,%10.5le)",i,drv->vec[i]);
4197 j++;
4198 }
4199 print_msg("\n");
4200 }
4201 return;
4202 }
4203
print_dof_real_vec(const DOF_REAL_VEC * vec)4204 void print_dof_real_vec(const DOF_REAL_VEC *vec)
4205 {
4206 FUNCNAME("print_dof_real_vec");
4207 int i;
4208
4209 i = 0;
4210 CHAIN_DO(vec, const DOF_REAL_VEC) {
4211 if (!CHAIN_SINGLE(vec)) {
4212 MSG("BLOCK(%d):\n", i);
4213 }
4214 __print_dof_real_vec(vec);
4215 ++i;
4216 } CHAIN_WHILE(vec, const DOF_REAL_VEC);
4217 }
4218
4219 static inline
__print_dof_real_d_vec(const DOF_REAL_D_VEC * drdv)4220 void __print_dof_real_d_vec(const DOF_REAL_D_VEC *drdv)
4221 {
4222 FUNCNAME("print_dof_real_d_vec");
4223 int i, j, k;
4224 const DOF_ADMIN *admin = NULL;
4225 #if DIM_OF_WORLD < 2
4226 static int per_line = 4;
4227 #elif DIM_OF_WORLD < 4
4228 static int per_line = 2;
4229 #else
4230 static int per_line = 1;
4231 #endif
4232
4233 if (drdv->fe_space) admin = drdv->fe_space->admin;
4234
4235 MSG("Vec `%s':\n", drdv->name);
4236 j = 0;
4237 if (admin)
4238 {
4239 FOR_ALL_DOFS(admin,
4240 if ((j % per_line) == 0) {
4241 if (j) print_msg("\n");
4242 MSG("(%3d:",dof);
4243 }
4244 else
4245 print_msg(" (%3d:", dof);
4246 for (k=0; k<DIM_OF_WORLD; k++)
4247 print_msg("%c%10.5le", (k > 0 ? ',':' '), drdv->vec[dof][k]);
4248 print_msg(")");
4249 j++;
4250 );
4251 print_msg("\n");
4252 }
4253 else
4254 {
4255 MSG("no DOF_ADMIN, print whole vector.\n");
4256
4257 for (i = 0; i < drdv->size; i++) {
4258 if ((j % per_line) == 0)
4259 {
4260 if (j) print_msg("\n");
4261 MSG("(%3d:",i);
4262 }
4263 else
4264 print_msg(" (%3d:", i);
4265 for (k=0; k<DIM_OF_WORLD; k++)
4266 print_msg("%c%10.5le", (k > 0 ? ',':' '), drdv->vec[i][k]);
4267 print_msg(")");
4268 }
4269 print_msg("\n");
4270 }
4271 return;
4272 }
4273
print_dof_real_d_vec(const DOF_REAL_D_VEC * vec)4274 void print_dof_real_d_vec(const DOF_REAL_D_VEC *vec)
4275 {
4276
4277 FUNCNAME("print_dof_real_d_vec");
4278 int i;
4279
4280 i = 0;
4281 CHAIN_DO(vec, const DOF_REAL_D_VEC) {
4282 if (!CHAIN_SINGLE(vec)) {
4283 MSG("BLOCK(%d):\n", i);
4284 }
4285 __print_dof_real_d_vec(vec);
4286 ++i;
4287 } CHAIN_WHILE(vec, const DOF_REAL_D_VEC);
4288 }
4289
print_dof_real_vec_dow(const DOF_REAL_VEC_D * vec)4290 void print_dof_real_vec_dow(const DOF_REAL_VEC_D *vec)
4291 {
4292 FUNCNAME("print_dof_real_vec_dow");
4293 int i;
4294
4295 i = 0;
4296 CHAIN_DO(vec, const DOF_REAL_VEC_D) {
4297 if (!CHAIN_SINGLE(vec)) {
4298 MSG("BLOCK(%d):\n", i);
4299 }
4300 if (vec->stride != 1) {
4301 __print_dof_real_d_vec((const DOF_REAL_D_VEC *)vec);
4302 } else{
4303 __print_dof_real_vec((const DOF_REAL_VEC *)vec);
4304 }
4305 ++i;
4306 } CHAIN_WHILE(vec, const DOF_REAL_VEC_D);
4307 }
4308
4309 static inline
__print_dof_real_dd_vec(const DOF_REAL_DD_VEC * drdv)4310 void __print_dof_real_dd_vec(const DOF_REAL_DD_VEC *drdv)
4311 {
4312 FUNCNAME("print_dof_real_dd_vec");
4313 int i;
4314 const DOF_ADMIN *admin = NULL;
4315
4316 if (drdv->fe_space) admin = drdv->fe_space->admin;
4317
4318 MSG("Vec `%s':\n", drdv->name);
4319 if (admin) {
4320 FOR_ALL_DOFS(admin, {
4321 MSG("(%3d: "MFORMAT_DOW")\n", dof, MEXPAND_DOW(drdv->vec[dof]));
4322 });
4323 } else {
4324 MSG("no DOF_ADMIN, print whole vector.\n");
4325
4326 for (i = 0; i < drdv->size; i++) {
4327 MSG("(%3d: "MFORMAT_DOW")\n", i, MEXPAND_DOW(drdv->vec[i]));
4328 }
4329 }
4330 }
4331
print_dof_real_dd_vec(const DOF_REAL_DD_VEC * vec)4332 void print_dof_real_dd_vec(const DOF_REAL_DD_VEC *vec)
4333 {
4334 FUNCNAME("print_dof_real_dd_vec");
4335 int i;
4336
4337 i = 0;
4338 CHAIN_DO(vec, const DOF_REAL_DD_VEC) {
4339 if (!CHAIN_SINGLE(vec)) {
4340 MSG("BLOCK(%d):\n", i);
4341 }
4342 __print_dof_real_dd_vec(vec);
4343 ++i;
4344 } CHAIN_WHILE(vec, const DOF_REAL_DD_VEC);
4345 }
4346
print_real_vec_maple(REAL * vector,int size,const char * vec_name)4347 void print_real_vec_maple(REAL *vector, int size, const char *vec_name)
4348 {
4349 fprint_real_vec_maple(stdout, vector, size, vec_name);
4350 }
4351
print_dof_real_vec_maple(const DOF_REAL_VEC * vec,const char * vec_name)4352 void print_dof_real_vec_maple(const DOF_REAL_VEC *vec, const char *vec_name)
4353 {
4354 if (vec_name == NULL)
4355 vec_name = vec->name;
4356
4357 fprint_dof_real_vec_maple(stdout, vec, vec_name);
4358 }
4359
print_dof_real_d_vec_maple(const DOF_REAL_D_VEC * vec,const char * vec_name)4360 void print_dof_real_d_vec_maple(const DOF_REAL_D_VEC *vec, const char *vec_name)
4361 {
4362 if (vec_name == NULL)
4363 vec_name = vec->name;
4364
4365 fprint_dof_real_d_vec_maple(stdout, vec, vec_name);
4366 }
4367
print_dof_real_vec_dow_maple(const DOF_REAL_VEC_D * vec,const char * vec_name)4368 void print_dof_real_vec_dow_maple(const DOF_REAL_VEC_D *vec, const char *vec_name)
4369 {
4370 if (vec_name == NULL)
4371 vec_name = vec->name;
4372
4373 fprint_dof_real_vec_dow_maple(stdout, vec, vec_name);
4374 }
4375
print_dof_matrix_maple(const DOF_MATRIX * matrix,const char * matrix_name)4376 void print_dof_matrix_maple(const DOF_MATRIX *matrix, const char *matrix_name)
4377 {
4378 if (matrix_name == NULL)
4379 matrix_name = matrix->name;
4380
4381 fprint_dof_matrix_maple(stdout, matrix, matrix_name);
4382 }
4383
4384
fprint_real_vec_maple(FILE * fp,REAL * vector,int size,const char * vec_name)4385 void fprint_real_vec_maple(FILE *fp, REAL *vector, int size, const char *vec_name)
4386 {
4387 int i;
4388 char _vec_name[24];
4389
4390
4391 if (vec_name == NULL) {
4392 sprintf(_vec_name, "REAL_VEC");
4393 vec_name = _vec_name;
4394 }
4395
4396 fprintf(fp, "\n#REAL_VEC \"%s\" in maple-format:\n\n", vec_name);
4397 fflush(fp);
4398
4399 fprintf(fp, "%s:=Vector(%d,proc(i) 0 end):\n\n", vec_name, size);
4400 fflush(fp);
4401
4402 for (i = 0; i < size; i++)
4403 {
4404 fprintf(fp, " %s[%d]:=%.17e:\n", vec_name, i+1, vector[i]);
4405 fflush(fp);
4406 }
4407 fprintf(fp, "\n%s:=Vector([%s]);\n\n\n\n\n", vec_name, vec_name);
4408 fflush(fp);
4409 }
4410
fprint_dof_real_vec_maple(FILE * fp,const DOF_REAL_VEC * vec,const char * vec_name)4411 void fprint_dof_real_vec_maple(FILE *fp, const DOF_REAL_VEC *vec, const char *vec_name)
4412 {
4413 fprint_dof_real_vec_dow_maple(fp, (const DOF_REAL_VEC_D *)vec, vec_name);
4414 }
4415
fprint_dof_real_d_vec_maple(FILE * fp,const DOF_REAL_D_VEC * vec,const char * vec_name)4416 void fprint_dof_real_d_vec_maple(FILE *fp, const DOF_REAL_D_VEC *vec, const char *vec_name)
4417 {
4418 fprint_dof_real_vec_dow_maple(fp, (const DOF_REAL_VEC_D *)vec, vec_name);
4419 }
4420
fprint_dof_real_vec_dow_maple(FILE * fp,const DOF_REAL_VEC_D * _vec,const char * vec_name)4421 void fprint_dof_real_vec_dow_maple(FILE *fp,
4422 const DOF_REAL_VEC_D *_vec,
4423 const char *vec_name)
4424 {
4425 char Chain_name[24];
4426
4427 int i, k, bdim = 0;
4428 if (vec_name == NULL)
4429 vec_name = _vec->name;
4430
4431 fprintf(fp, "\n#DOF_REAL_VEC_D %s in maple-format:\n\n", vec_name);
4432 fflush(fp);
4433
4434
4435 i = 0;
4436 CHAIN_DO(_vec, const DOF_REAL_VEC_D) {
4437 int count = 0;
4438 int dim = 0;
4439 fprintf(fp, "%s", vec_name);
4440 fflush(fp);
4441 *Chain_name = '\0';
4442 if (!CHAIN_SINGLE(_vec))
4443 sprintf(Chain_name,"_Chain%d", i);
4444 if (_vec->stride != 1) {
4445 const DOF_REAL_D_VEC *vec = (const DOF_REAL_D_VEC *)_vec;
4446 dim = DIM_OF_WORLD*vec->fe_space->admin->size_used;
4447 fprintf(fp, "%s", Chain_name);
4448 fprintf(fp, ":=Vector(%d,proc(i) 0 end):\n\n", dim);
4449 fflush(fp);
4450
4451 FOR_ALL_DOFS(vec->fe_space->admin,
4452 for (k = 0; k < DIM_OF_WORLD; k++) {
4453 fprintf(fp, " ");
4454 fprintf(fp, "%s", vec_name);
4455 fprintf(fp, "%s", Chain_name);
4456 fprintf(fp, "[%d]:=%.17e:\n", count+1, vec->vec[dof][k]);
4457
4458 count++;
4459 }
4460 fflush(fp);
4461 );
4462 fprintf(fp, "\n\n\n\n");
4463 fflush(fp);
4464 } else {
4465 const DOF_REAL_VEC *vec = (const DOF_REAL_VEC *)_vec;
4466 dim = vec->fe_space->admin->size_used;
4467 fprintf(fp, "%s", Chain_name);
4468 fprintf(fp, ":=Vector(%d,proc(i) 0 end):\n\n", dim);
4469 fflush(fp);
4470
4471 FOR_ALL_DOFS(vec->fe_space->admin,
4472 fprintf(fp, " ");
4473 fprintf(fp, "%s", vec_name);
4474 fprintf(fp, "%s", Chain_name);
4475 fprintf(fp, "[%d]:=%.17e:\n", dof+1, vec->vec[dof]);
4476 fflush(fp);
4477 );
4478
4479 fprintf(fp, "\n\n\n\n");
4480 fflush(fp);
4481 }
4482 ++i;
4483 bdim += dim;
4484 } CHAIN_WHILE(_vec, const DOF_REAL_VEC_D);
4485
4486 /* output as blockvector */
4487 fprintf(fp, "%s", vec_name);
4488 fprintf(fp, ":=Vector([");
4489 for (k = 0; k < i; k++) {
4490 if (k != 0)
4491 fprintf(fp, ",");
4492 fprintf(fp, "%s", vec_name);
4493 if (i >= 2)
4494 fprintf(fp, "_Chain%d", k);
4495 }
4496 fprintf(fp, "]);\n");
4497 fprintf(fp, "\n\n\n\n\n");
4498 fflush(fp);
4499
4500 }
4501
4502
fprint_dof_matrix_maple(FILE * fp,const DOF_MATRIX * matrix,const char * matrix_name)4503 void fprint_dof_matrix_maple(FILE *fp, const DOF_MATRIX *matrix, const char *matrix_name)
4504 {
4505 FUNCNAME("fprint_dof_matrix_maple");
4506
4507
4508 if (matrix_name == NULL)
4509 matrix_name = matrix->name;
4510
4511
4512
4513 fprintf(fp, "\n");
4514 fprintf(fp, "#DOF_MATRIX ");
4515 fprintf(fp, "%s", matrix_name);
4516 fprintf(fp, " in maple-format:\n\n");
4517 fflush(fp);
4518
4519
4520
4521
4522 int i, j, k, l, n, m, row_shift, col_shift;
4523 n = 0;
4524 COL_CHAIN_DO(matrix, const DOF_MATRIX) {
4525 m = 0;
4526 ROW_CHAIN_DO(matrix, const DOF_MATRIX) {
4527
4528 int row_dim = matrix->row_fe_space->admin->size_used;
4529 int col_dim = matrix->col_fe_space->admin->size_used;
4530
4531 switch (matrix->type) {
4532 case MATENT_REAL:
4533 if ((matrix->row_fe_space->rdim == DIM_OF_WORLD) && (matrix->col_fe_space->rdim == DIM_OF_WORLD) &&
4534 (matrix->row_fe_space->bas_fcts->rdim == 1) && (matrix->col_fe_space->bas_fcts->rdim == 1))
4535 {
4536 if (matrix->row_fe_space == matrix->col_fe_space) {
4537 fprintf(fp, "%s", matrix_name);
4538
4539 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4540 fprintf(fp, "_Chain%d%d", n, m);
4541 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) if i<>j then 0; else 1; end; end):\n\n",
4542 DIM_OF_WORLD*row_dim, DIM_OF_WORLD*col_dim);
4543 fflush(fp);
4544 } else {
4545 fprintf(fp, "%s", matrix_name);
4546
4547 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix)) {
4548 fprintf(fp, "_Chain%d%d", n, m);
4549 }
4550 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) 0 end):\n\n",
4551 DIM_OF_WORLD*row_dim, DIM_OF_WORLD*col_dim);
4552 fflush(fp);
4553 }
4554 row_shift = 0;
4555
4556 if (matrix->is_diagonal) {
4557 for (i = 0; i < matrix->row_fe_space->admin->size_used; i++) {
4558 row_shift = i*DIM_OF_WORLD;
4559 col_shift = row_shift;
4560 fprintf(fp, " ");
4561 fprintf(fp, "%s", matrix_name);
4562 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4563 fprintf(fp, "_Chain%d%d", n, m);
4564 fprintf(fp, "[%d,%d]:=%.17e:\n", j + row_shift + 1, j + col_shift + 1, matrix->diagonal.real->vec[i]);
4565 fflush(fp);
4566 fprintf(fp, "\n");
4567 fflush(fp);
4568 }
4569 } else {
4570 for (i = 0; i < matrix->size; i++) {
4571 row_shift = i*DIM_OF_WORLD;
4572 FOR_ALL_MAT_COLS(REAL, matrix->matrix_row[i], {
4573 col_shift = col_dof*DIM_OF_WORLD;
4574 for (j = 0; j < DIM_OF_WORLD; j++){
4575 fprintf(fp, " ");
4576 fprintf(fp, "%s", matrix_name);
4577 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4578 fprintf(fp, "_Chain%d%d", n, m);
4579 fprintf(fp, "[%d,%d]:=%.17e:\n", j + row_shift + 1, j + col_shift + 1, row->entry[col_idx]);
4580 }
4581 fflush(fp);
4582 });
4583 if (matrix->matrix_row[i] != NULL) {
4584 fprintf(fp, "\n");
4585 fflush(fp);
4586 }
4587 }
4588 }
4589 } else {
4590 if (matrix->row_fe_space == matrix->col_fe_space) {
4591 fprintf(fp, "%s", matrix_name);
4592 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4593 fprintf(fp, "_Chain%d%d", n, m);
4594 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) if i<>j then 0; else 1; end; end):\n\n",
4595 row_dim, col_dim);
4596 fflush(fp);
4597 } else {
4598 fprintf(fp, "%s", matrix_name);
4599 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4600 fprintf(fp, "_Chain%d%d", n, m);
4601 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) 0 end):\n\n",
4602 row_dim, col_dim);
4603 fflush(fp);
4604 }
4605 if (matrix->is_diagonal) {
4606 for (i = 0; i < matrix->row_fe_space->admin->size_used; i++) {
4607 fprintf(fp, " ");
4608 fprintf(fp, "%s", matrix_name);
4609 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4610 fprintf(fp, "_Chain%d%d", n, m);
4611 fprintf(fp, "[%d,%d]:=%.17e:\n", i+1, i+1, matrix->diagonal.real->vec[i]);
4612 fflush(fp);
4613 fprintf(fp, "\n");
4614 fflush(fp);
4615 }
4616 } else {
4617 for (i = 0; i < matrix->size; i++) {
4618 FOR_ALL_MAT_COLS(REAL, matrix->matrix_row[i], {
4619 fprintf(fp, " ");
4620 fprintf(fp, "%s", matrix_name);
4621 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4622 fprintf(fp, "_Chain%d%d", n, m);
4623 fprintf(fp, "[%d,%d]:=%.17e:\n", i+1, col_dof+1, row->entry[col_idx]);
4624 fflush(fp);
4625 });
4626 if (matrix->matrix_row[i] != NULL) {
4627 fprintf(fp, "\n");
4628 fflush(fp);
4629 }
4630 }
4631 }
4632 }
4633 break;
4634 case MATENT_REAL_D:
4635 if ((matrix->row_fe_space->rdim == DIM_OF_WORLD &&
4636 matrix->col_fe_space->rdim == 1)
4637 ||
4638 (matrix->row_fe_space->rdim == DIM_OF_WORLD &&
4639 matrix->col_fe_space->rdim == DIM_OF_WORLD &&
4640 matrix->col_fe_space->bas_fcts->rdim == DIM_OF_WORLD)) {
4641 fprintf(fp, "%s", matrix_name);
4642 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4643 fprintf(fp, "_Chain%d%d", n, m);
4644 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) 0 end):\n\n",
4645 DIM_OF_WORLD*row_dim, col_dim);
4646 fflush(fp);
4647 if (matrix->is_diagonal) {
4648 for (i = 0; i < matrix->row_fe_space->admin->size_used; i++) {
4649 for (k = 0; k < DIM_OF_WORLD; k++) {
4650 fprintf(fp, " ");
4651 fprintf(fp, "%s", matrix_name);
4652 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4653 fprintf(fp, "_Chain%d%d", n, m);
4654 fprintf(fp, "[%d,%d]:=%.17e:\n",
4655 i*DIM_OF_WORLD+1+k, i+1, matrix->diagonal.real_d->vec[i][k]);
4656 }
4657 fprintf(fp, "\n");
4658 fflush(fp);
4659 }
4660 } else {
4661 for (i = 0; i < matrix->size; i++) {
4662 FOR_ALL_MAT_COLS(REAL_D, matrix->matrix_row[i], {
4663 for (k = 0; k < DIM_OF_WORLD; k++) {
4664 fprintf(fp, " ");
4665 fprintf(fp, "%s", matrix_name);
4666 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4667 fprintf(fp, "_Chain%d%d", n, m);
4668 fprintf(fp, "[%d,%d]:=%.17e:\n",
4669 i*DIM_OF_WORLD+1+k, col_dof+1, row->entry[col_idx][k]);
4670 }
4671 fflush(fp);
4672 });
4673 if (matrix->matrix_row[i] != NULL) {
4674 fprintf(fp, "\n");
4675 fflush(fp);
4676 }
4677 }
4678 }
4679 } else if ((matrix->row_fe_space->rdim == 1 &&
4680 matrix->col_fe_space->rdim == DIM_OF_WORLD)
4681 ||
4682 (matrix->row_fe_space->rdim == DIM_OF_WORLD &&
4683 matrix->row_fe_space->bas_fcts->rdim == DIM_OF_WORLD &&
4684 matrix->col_fe_space->rdim == DIM_OF_WORLD)) {
4685 fprintf(fp, "%s", matrix_name);
4686 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4687 fprintf(fp, "_Chain%d%d", n, m);
4688 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) 0 end):\n\n",
4689 row_dim, col_dim*DIM_OF_WORLD);
4690 fflush(fp);
4691 if (matrix->is_diagonal) {
4692 for (i = 0; i < matrix->row_fe_space->admin->size_used; i++) {
4693 for (k = 0; k < DIM_OF_WORLD; k++) {
4694 fprintf(fp, " ");
4695 fprintf(fp, "%s", matrix_name);
4696 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4697 fprintf(fp, "_Chain%d%d", n, m);
4698 fprintf(fp, "[%d,%d]:=%.17e:\n",
4699 i+1, i*DIM_OF_WORLD+1+k, matrix->diagonal.real_d->vec[i][k]);
4700 }
4701 fprintf(fp, "\n");
4702 fflush(fp);
4703 }
4704 } else {
4705 for (i = 0; i < matrix->size; i++) {
4706 FOR_ALL_MAT_COLS(REAL_D, matrix->matrix_row[i], {
4707 for (k = 0; k < DIM_OF_WORLD; k++) {
4708 fprintf(fp, " ");
4709 fprintf(fp, "%s", matrix_name);
4710 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4711 fprintf(fp, "_Chain%d%d", n, m);
4712 fprintf(fp, "[%d,%d]:=%.17e:\n",
4713 i+1, col_dof*DIM_OF_WORLD+1+k, row->entry[col_idx][k]);
4714 }
4715 fflush(fp);
4716 });
4717 if (matrix->matrix_row[i] != NULL) {
4718 fprintf(fp, "\n");
4719 fflush(fp);
4720 }
4721 }
4722 }
4723 }
4724 break;
4725 case MATENT_REAL_DD:
4726 fprintf(fp, "%s", matrix_name);
4727 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4728 fprintf(fp, "_Chain%d%d", n, m);
4729 fprintf(fp, ":=Matrix(%d,%d,proc(i,j) 0 end):\n\n",
4730 row_dim*DIM_OF_WORLD, col_dim*DIM_OF_WORLD);
4731 fflush(fp);
4732 if (matrix->is_diagonal) {
4733 for (i = 0; i < matrix->row_fe_space->admin->size_used; i++) {
4734 for (k = 0; k < DIM_OF_WORLD; k++) {
4735 for (l = 0; l < DIM_OF_WORLD; l++) {
4736 fprintf(fp, " ");
4737 fprintf(fp, "%s", matrix_name);
4738 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4739 fprintf(fp, "_Chain%d%d", n, m);
4740 fprintf(fp, "[%d,%d]:=%.17e:\n",
4741 i*DIM_OF_WORLD+1+k, i*DIM_OF_WORLD+1+l, matrix->diagonal.real_dd->vec[i][k][l]);
4742 }
4743 }
4744 fprintf(fp, "\n");
4745 fflush(fp);
4746 }
4747 } else {
4748 for (i = 0; i < matrix->size; i++) {
4749 FOR_ALL_MAT_COLS(REAL_DD, matrix->matrix_row[i], {
4750 for (k = 0; k < DIM_OF_WORLD; k++) {
4751 for (l = 0; l < DIM_OF_WORLD; l++) {
4752 fprintf(fp, " ");
4753 fprintf(fp, "%s", matrix_name);
4754 if (!COL_CHAIN_SINGLE(matrix) || !ROW_CHAIN_SINGLE(matrix))
4755 fprintf(fp, "_Chain%d%d", n, m);
4756 fprintf(fp, "[%d,%d]:=%.17e:\n",
4757 i*DIM_OF_WORLD+1+k, col_dof*DIM_OF_WORLD+1+l,
4758 row->entry[col_idx][k][l]);
4759 }
4760 }
4761 fflush(fp);
4762 });
4763 if (matrix->matrix_row[i] != NULL) {
4764 fprintf(fp, "\n");
4765 fflush(fp);
4766 }
4767 }
4768 }
4769 break;
4770 default:
4771 ERROR("Unknown matrix type: %d\n", matrix->type);
4772 }
4773 fprintf(fp, "\n");
4774 fflush(fp);
4775
4776 ++m;
4777 } ROW_CHAIN_WHILE(matrix, const DOF_MATRIX);
4778 ++n;
4779 } COL_CHAIN_WHILE(matrix, const DOF_MATRIX);
4780
4781 /* output as blockmatrix */
4782 fprintf(fp, "%s", matrix_name);
4783 fprintf(fp, ":=Matrix([");
4784 for (i = 0; i < n; i++) {
4785 if (i != 0)
4786 fprintf(fp, ",");
4787 fprintf(fp, "[");
4788 for (j = 0; j < m; j++) {
4789 if (j != 0)
4790 fprintf(fp, ",");
4791 fprintf(fp, "evalm(");
4792 fprintf(fp, "%s", matrix_name);
4793 if ((n >= 2)||(m >= 2))
4794 fprintf(fp, "_Chain%d%d", i, j);
4795 fprintf(fp, ")");
4796 }
4797 fprintf(fp, "]");
4798 }
4799 fprintf(fp, "]);\n");
4800 fprintf(fp, "\n\n\n\n\n");
4801 fflush(fp);
4802
4803 }
4804
4805
file_print_real_vec_maple(const char * file_name,const char * mode,REAL * vector,int size,const char * vec_name)4806 void file_print_real_vec_maple(const char *file_name,
4807 const char *mode,
4808 REAL *vector, int size, const char *vec_name)
4809 {
4810 FILE *fp;
4811
4812 fp = fopen(file_name, mode);
4813 fprint_real_vec_maple(fp, vector, size, vec_name);
4814 fclose(fp);
4815 }
4816
file_print_dof_real_vec_maple(const char * file_name,const char * mode,const DOF_REAL_VEC * vec,const char * vec_name)4817 void file_print_dof_real_vec_maple(const char *file_name,
4818 const char *mode,
4819 const DOF_REAL_VEC *vec,
4820 const char *vec_name)
4821 {
4822 if (vec_name == NULL)
4823 vec_name = vec->name;
4824
4825 FILE *fp;
4826
4827 fp = fopen(file_name, mode);
4828 fprint_dof_real_vec_maple(fp, vec, vec_name);
4829 fclose(fp);
4830 }
4831
file_print_dof_real_d_vec_maple(const char * file_name,const char * mode,const DOF_REAL_D_VEC * vec,const char * vec_name)4832 void file_print_dof_real_d_vec_maple(const char *file_name,
4833 const char *mode,
4834 const DOF_REAL_D_VEC *vec,
4835 const char *vec_name)
4836 {
4837 if (vec_name == NULL)
4838 vec_name = vec->name;
4839
4840 FILE *fp;
4841
4842 fp = fopen(file_name, mode);
4843 fprint_dof_real_d_vec_maple(fp, vec, vec_name);
4844 fclose(fp);
4845 }
4846
file_print_dof_real_vec_dow_maple(const char * file_name,const char * mode,const DOF_REAL_VEC_D * vec,const char * vec_name)4847 void file_print_dof_real_vec_dow_maple(const char *file_name,
4848 const char *mode,
4849 const DOF_REAL_VEC_D *vec,
4850 const char *vec_name)
4851 {
4852 if (vec_name == NULL)
4853 vec_name = vec->name;
4854
4855 FILE *fp;
4856
4857 fp = fopen(file_name, mode);
4858 fprint_dof_real_vec_dow_maple(fp, vec, vec_name);
4859 fclose(fp);
4860 }
4861
file_print_dof_matrix_maple(const char * file_name,const char * mode,const DOF_MATRIX * matrix,const char * matrix_name)4862 void file_print_dof_matrix_maple(const char *file_name,
4863 const char *mode,
4864 const DOF_MATRIX *matrix,
4865 const char *matrix_name)
4866 {
4867 if (matrix_name == NULL)
4868 matrix_name = matrix->name;
4869
4870 FILE *fp;
4871
4872 fp = fopen(file_name, mode);
4873 fprint_dof_matrix_maple(fp, matrix, matrix_name);
4874 fclose(fp);
4875 }
4876
4877
4878 static inline
__print_dof_ptr_vec(const DOF_PTR_VEC * dpv)4879 void __print_dof_ptr_vec(const DOF_PTR_VEC *dpv)
4880 {
4881 FUNCNAME("print_dof_ptr_vec");
4882 int i, j;
4883 const DOF_ADMIN *admin = NULL;
4884 const char *format;
4885
4886 if (dpv->fe_space) admin = dpv->fe_space->admin;
4887
4888 MSG("Vector `%s':\n", dpv->name);
4889 j = 0;
4890
4891
4892 if (admin) {
4893 if (admin->size_used > 100)
4894 format = "%s(%3d,%p)";
4895 else if (admin->size_used > 10)
4896 format = "%s(%2d,%p)";
4897 else
4898 format = "%s(%1d,%p)";
4899
4900 FOR_ALL_DOFS(admin,
4901 if ((j % 5) == 0)
4902 {
4903 if (j) print_msg("\n");
4904 MSG(format, "", dof,dpv->vec[dof]);
4905 }
4906 else
4907 print_msg(format, " ",dof,dpv->vec[dof]);
4908 j++;
4909 );
4910 print_msg("\n");
4911 }
4912 else {
4913 if (dpv->size > 100)
4914 format = "%s(%3d,%p)";
4915 else if (dpv->size > 10)
4916 format = "%s(%2d,%p)";
4917 else
4918 format = "%s(%1d,%p)";
4919
4920 for (i = 0; i < dpv->size; i++) {
4921 if ((j % 5) == 0)
4922 {
4923 if (j) print_msg("\n");
4924 MSG(format, "", i, dpv->vec[i]);
4925 }
4926 else
4927 print_msg(format, " ", i, dpv->vec[i]);
4928 j++;
4929 }
4930 print_msg("\n");
4931 }
4932 return;
4933 }
4934
print_dof_ptr_vec(const DOF_PTR_VEC * vec)4935 void print_dof_ptr_vec(const DOF_PTR_VEC *vec)
4936 {
4937
4938 FUNCNAME("print_dof_ptr_vec");
4939 int i;
4940
4941 i = 0;
4942 CHAIN_DO(vec, const DOF_PTR_VEC) {
4943 if (!CHAIN_SINGLE(vec)) {
4944 MSG("BLOCK(%d):\n", i);
4945 }
4946 __print_dof_ptr_vec(vec);
4947 ++i;
4948 } CHAIN_WHILE(vec, const DOF_PTR_VEC);
4949 }
4950
4951 static inline
__print_dof_int_vec(const DOF_INT_VEC * div)4952 void __print_dof_int_vec(const DOF_INT_VEC *div)
4953 {
4954 FUNCNAME("print_dof_int_vec");
4955 int i, j;
4956 const DOF_ADMIN *admin = NULL;
4957 const char *format;
4958
4959 if (div->fe_space) admin = div->fe_space->admin;
4960
4961 MSG("Vector `%s':\n", div->name);
4962 j = 0;
4963 if (admin) {
4964 if (admin->size_used > 100)
4965 format = "%s(%3d,%3d)";
4966 else if (admin->size_used > 10)
4967 format = "%s(%2d,%3d)";
4968 else
4969 format = "%s(%1d,%3d)";
4970
4971 FOR_ALL_DOFS(admin,
4972 if ((j % 5) == 0)
4973 {
4974 if (j) print_msg("\n");
4975 MSG(format, "", dof,div->vec[dof]);
4976 }
4977 else
4978 print_msg(format, " ", dof, div->vec[dof]);
4979 j++;
4980 );
4981 print_msg("\n");
4982 }
4983 else {
4984 if (div->size > 100)
4985 format = "%s(%3d,%3d)";
4986 else if (div->size > 10)
4987 format = "%s(%2d,%3d)";
4988 else
4989 format = "%s(%1d,%3d)";
4990
4991 for (i = 0; i < div->size; i++) {
4992 if ((j % 5) == 0)
4993 {
4994 if (j) print_msg("\n");
4995 MSG(format, "", i, div->vec[i]);
4996 }
4997 else
4998 print_msg(format, " ", i, div->vec[i]);
4999 j++;
5000 }
5001 print_msg("\n");
5002 }
5003 return;
5004 }
5005
print_dof_int_vec(const DOF_INT_VEC * vec)5006 void print_dof_int_vec(const DOF_INT_VEC *vec)
5007 {
5008
5009 FUNCNAME("print_dof_int_vec");
5010 int i;
5011
5012 i = 0;
5013 CHAIN_DO(vec, const DOF_INT_VEC) {
5014 if (!CHAIN_SINGLE(vec)) {
5015 MSG("BLOCK(%d):\n", i);
5016 }
5017 __print_dof_int_vec(vec);
5018 ++i;
5019 } CHAIN_WHILE(vec, const DOF_INT_VEC);
5020 }
5021
5022 static inline
__print_dof_uchar_vec(const DOF_UCHAR_VEC * duv)5023 void __print_dof_uchar_vec(const DOF_UCHAR_VEC *duv)
5024 {
5025 FUNCNAME("print_dof_uchar_vec");
5026 int i, j;
5027 const DOF_ADMIN *admin = NULL;
5028 const char *format;
5029
5030 if (duv->fe_space) admin = duv->fe_space->admin;
5031
5032 MSG("Vector `%s':\n", duv->name);
5033 j = 0;
5034
5035 if (admin) {
5036 if (admin->size_used > 100)
5037 format = "%s(%3d,0x%02X)";
5038 else if (admin->size_used > 10)
5039 format = "%s(%2d,0x%02X)";
5040 else
5041 format = "%s(%1d,0x%02X)";
5042
5043 FOR_ALL_DOFS(admin,
5044 if ((j % 5) == 0) {
5045 if (j) print_msg("\n");
5046 MSG(format, "", dof, duv->vec[dof]);
5047 }
5048 else
5049 print_msg(format, " ", dof, duv->vec[dof]);
5050 j++;
5051 );
5052 print_msg("\n");
5053 }
5054 else {
5055 if (duv->size > 100)
5056 format = "%s(%3d,0x%20X)";
5057 else if (duv->size > 10)
5058 format = "%s(%2d,0x%02X)";
5059 else
5060 format = "%s(%1d,0x%02X)";
5061
5062 for (i = 0; i < duv->size; i++) {
5063 if ((j % 5) == 0)
5064 {
5065 if (j) print_msg("\n");
5066 MSG(format, "", i, duv->vec[i]);
5067 }
5068 else
5069 print_msg(format, " ", i, duv->vec[i]);
5070 j++;
5071 }
5072 print_msg("\n");
5073 }
5074 return;
5075 }
5076
print_dof_uchar_vec(const DOF_UCHAR_VEC * vec)5077 void print_dof_uchar_vec(const DOF_UCHAR_VEC *vec)
5078 {
5079
5080 FUNCNAME("print_dof_uchar_vec");
5081 int i;
5082
5083 i = 0;
5084 CHAIN_DO(vec, const DOF_UCHAR_VEC) {
5085 if (!CHAIN_SINGLE(vec)) {
5086 MSG("BLOCK(%d):\n", i);
5087 }
5088 __print_dof_uchar_vec(vec);
5089 ++i;
5090 } CHAIN_WHILE(vec, const DOF_UCHAR_VEC);
5091 }
5092
5093 static inline
__print_dof_schar_vec(const DOF_SCHAR_VEC * dsv)5094 void __print_dof_schar_vec(const DOF_SCHAR_VEC *dsv)
5095 {
5096 FUNCNAME("print_dof_schar_vec");
5097 int i, j;
5098 const DOF_ADMIN *admin = NULL;
5099 const char *format;
5100
5101 if (dsv->fe_space) admin = dsv->fe_space->admin;
5102
5103 MSG("Vector `%s':\n", dsv->name);
5104 j = 0;
5105
5106
5107 if (admin) {
5108 if (admin->size_used > 100)
5109 format = "%s(%3d,0x%02X)";
5110 else if (admin->size_used > 10)
5111 format = "%s(%2d,0x%02X)";
5112 else
5113 format = "%s(%1d,0x%02X)";
5114
5115 FOR_ALL_DOFS(admin,
5116 if ((j % 5) == 0)
5117 {
5118 if (j) print_msg("\n");
5119 MSG(format, "", dof, dsv->vec[dof] & 0xFF);
5120 }
5121 else
5122 print_msg(format, " ", dof, dsv->vec[dof] & 0xFF);
5123 j++;
5124 );
5125 print_msg("\n");
5126 }
5127 else {
5128 if (dsv->size > 100)
5129 format = "%s(%3d,0x%02X)";
5130 else if (dsv->size > 10)
5131 format = "%s(%2d,0x%02X)";
5132 else
5133 format = "%s(%1d,0x%02X)";
5134
5135 for (i = 0; i < dsv->size; i++) {
5136 if ((j % 5) == 0)
5137 {
5138 if (j) print_msg("\n");
5139 MSG(format, "", i, dsv->vec[i] & 0xFF);
5140 }
5141 else
5142 print_msg(format, " ", i, dsv->vec[i] & 0xFF);
5143 j++;
5144 }
5145 print_msg("\n");
5146 }
5147 return;
5148 }
5149
print_dof_schar_vec(const DOF_SCHAR_VEC * vec)5150 void print_dof_schar_vec(const DOF_SCHAR_VEC *vec)
5151 {
5152
5153 FUNCNAME("print_dof_schar_vec");
5154 int i;
5155
5156 i = 0;
5157 CHAIN_DO(vec, const DOF_SCHAR_VEC) {
5158 if (!CHAIN_SINGLE(vec)) {
5159 MSG("BLOCK(%d):\n", i);
5160 }
5161 __print_dof_schar_vec(vec);
5162 ++i;
5163 } CHAIN_WHILE(vec, const DOF_SCHAR_VEC);
5164 }
5165
5166 /*--------------------------------------------------------------------------*/
5167 /* clear_dof_matrix: remove all entries from dof_matrix */
5168 /*--------------------------------------------------------------------------*/
5169
_AI_clear_dof_matrix_single(DOF_MATRIX * matrix)5170 static inline void _AI_clear_dof_matrix_single(DOF_MATRIX *matrix)
5171 {
5172 int i;
5173 MATRIX_ROW *row, *next;
5174
5175 if (!matrix->is_diagonal) {
5176 if (matrix->matrix_row) {
5177 for (i=0; i < matrix->size; i++) {
5178 for (row = matrix->matrix_row[i]; row; row = next) {
5179 next = row->next;
5180 free_matrix_row(matrix->row_fe_space, row);
5181 }
5182 matrix->matrix_row[i] = NULL;
5183 }
5184 }
5185 } else {
5186 if (matrix->diagonal.real) {
5187 MAT_SWITCH_TYPE(matrix->type,
5188 free_dof_real_dd_vec(matrix->diagonal.real_dd),
5189 free_dof_real_d_vec(matrix->diagonal.real_d),
5190 free_dof_real_vec(matrix->diagonal.real));
5191 matrix->diagonal.real = NULL;
5192 if (matrix->unchained) {
5193 ((DOF_MATRIX *)matrix->unchained)->diagonal.real = NULL;
5194 }
5195 }
5196 if (matrix->inv_diag.real) {
5197 MAT_SWITCH_TYPE(matrix->type,
5198 free_dof_real_dd_vec(matrix->inv_diag.real_dd),
5199 free_dof_real_d_vec(matrix->inv_diag.real_d),
5200 free_dof_real_vec(matrix->inv_diag.real));
5201 matrix->inv_diag.real = NULL;
5202 if (matrix->unchained) {
5203 ((DOF_MATRIX *)matrix->unchained)->inv_diag.real = NULL;
5204 }
5205 }
5206 FOR_ALL_DOFS(matrix->row_fe_space->admin,
5207 matrix->diag_cols->vec[dof] = UNUSED_ENTRY);
5208 }
5209 matrix->type = MATENT_NONE; /* set by add_element_matrix() as required */
5210 matrix->n_entries = 0;
5211 }
5212
clear_dof_matrix(DOF_MATRIX * matrix)5213 void clear_dof_matrix(DOF_MATRIX *matrix)
5214 {
5215 COL_CHAIN_DO(matrix, DOF_MATRIX) {
5216 ROW_CHAIN_DO(matrix, DOF_MATRIX) {
5217 _AI_clear_dof_matrix_single(matrix);
5218 } ROW_CHAIN_WHILE(matrix, DOF_MATRIX);
5219 } COL_CHAIN_WHILE(matrix, DOF_MATRIX);
5220 }
5221
5222 /* Initialize new DOFs to UNUSED_ENTRY, and also mark the (now unused)
5223 * parental DOF as UNUSED_ENTRY.
5224 */
refine_diag_cols(DOF_INT_VEC * vec,RC_LIST_EL * rclist,int n)5225 static void refine_diag_cols(DOF_INT_VEC *vec, RC_LIST_EL *rclist, int n)
5226 {
5227 const DOF_ADMIN *admin = vec->fe_space->admin;
5228 int n0, node, n_dofs;
5229 int i;
5230
5231 node = admin->mesh->node[CENTER];
5232 n0 = admin->n0_dof[CENTER];
5233 n_dofs = admin->n_dof[CENTER];
5234
5235 for (i = 0; i < n; i++) {
5236 int j, k;
5237 EL *ch, *el;
5238 el = rclist[i].el_info.el;
5239 for (j = 0; j < 2; j++) {
5240 ch = el->child[j];
5241 for (k = 0; k < n_dofs; k++) {
5242 DOF dof = ch->dof[node][n0+k];
5243 vec->vec[dof] = UNUSED_ENTRY;
5244 }
5245 }
5246 for (k = 0; k < n_dofs; k++) {
5247 DOF dof = el->dof[node][n0+k];
5248 vec->vec[dof] = UNUSED_ENTRY;
5249 }
5250 }
5251 }
5252
dof_matrix_set_diagonal(DOF_MATRIX * matrix,bool diag)5253 void dof_matrix_set_diagonal(DOF_MATRIX *matrix, bool diag)
5254 {
5255 matrix->is_diagonal = diag;
5256 if (matrix->is_diagonal) {
5257 if (matrix->matrix_row) {
5258 MEM_FREE(matrix->matrix_row, matrix->size, MATRIX_ROW *);
5259 matrix->matrix_row = NULL;
5260 }
5261 if (matrix->diag_cols == NULL) {
5262 matrix->diag_cols = get_dof_int_vec("diag cols", matrix->row_fe_space);
5263 matrix->diag_cols->refine_interpol = refine_diag_cols;
5264 FOR_ALL_DOFS(matrix->row_fe_space->admin,
5265 matrix->diag_cols->vec[dof] = UNUSED_ENTRY);
5266 }
5267 } else {
5268 if (matrix->matrix_row == NULL) {
5269 matrix->matrix_row = MEM_CALLOC(matrix->size, MATRIX_ROW *);
5270 }
5271 if (matrix->diag_cols) {
5272 free_dof_int_vec(matrix->diag_cols);
5273 matrix->diag_cols = NULL;
5274 }
5275 }
5276 }
5277
5278 /* Set parts of a DOF_MATRIX chain to "diagonal", if the undlerlying
5279 * spaces have exactly one DOF per element which belongs to the
5280 * interior.
5281 */
dof_matrix_try_diagonal(DOF_MATRIX * matrix)5282 void dof_matrix_try_diagonal(DOF_MATRIX *matrix)
5283 {
5284 COL_CHAIN_DO(matrix, DOF_MATRIX) {
5285 ROW_CHAIN_DO(matrix, DOF_MATRIX) {
5286 const FE_SPACE *row_fe_space = matrix->row_fe_space;
5287 const FE_SPACE *col_fe_space = matrix->col_fe_space;
5288 if (row_fe_space->admin->n_dof[CENTER] == 1 &&
5289 row_fe_space->admin->n_dof[VERTEX] == 0 &&
5290 row_fe_space->admin->n_dof[EDGE] == 0 &&
5291 row_fe_space->admin->n_dof[FACE] == 0 &&
5292 (col_fe_space == NULL ||
5293 (col_fe_space->admin->n_dof[CENTER] == 1 &&
5294 col_fe_space->admin->n_dof[VERTEX] == 0 &&
5295 col_fe_space->admin->n_dof[EDGE] == 0 &&
5296 col_fe_space->admin->n_dof[FACE] == 0))) {
5297 dof_matrix_set_diagonal(matrix, true);
5298 }
5299 } ROW_CHAIN_WHILE(matrix, DOF_MATRIX);
5300 } COL_CHAIN_WHILE(matrix, DOF_MATRIX);
5301 }
5302
5303 /*--------------------------------------------------------------------------*/
5304 /* test_dof_matrix */
5305 /*--------------------------------------------------------------------------*/
5306
test_dof_matrix(DOF_MATRIX * matrix)5307 void test_dof_matrix(DOF_MATRIX *matrix)
5308 {
5309 FUNCNAME("test_dof_matrix");
5310 int i, j, jcol, k,kcol;
5311 MATRIX_ROW *row, *row2;
5312 int non_symmetric, found;
5313 MATENT_TYPE type = matrix->type;
5314 REAL sum;
5315
5316 /* test symmetry */
5317 non_symmetric = 0;
5318 for (i=0; i<matrix->size; i++) {
5319 sum = 0.0;
5320 for (row = matrix->matrix_row[i]; row; row = row->next) {
5321 for (j=0; j<ROW_LENGTH; j++) {
5322 jcol = row->col[j];
5323 if (ENTRY_USED(jcol)) {
5324 found = 0;
5325 for (row2 = matrix->matrix_row[jcol]; row2;
5326 row2 = row2 ? row2->next : NULL) {
5327 for (k=0; k<ROW_LENGTH; k++) {
5328 kcol = row2->col[k];
5329 if (ENTRY_USED(kcol)) {
5330 if (kcol==i) {
5331 found = 1;
5332 switch (type) {
5333 case MATENT_REAL_DD:
5334 if (MDST2_DOW((const REAL_D *)row2->entry.real_dd[k],
5335 (const REAL_D *)row->entry.real_dd[j])
5336 > 1.E-10) {
5337 non_symmetric = 1;
5338 MSG("mat[%d,%d]="MFORMAT_DOW \
5339 " != mat[%d,%d]="MFORMAT_DOW"\n",
5340 i, jcol, MEXPAND_DOW(row->entry.real_dd[j]),
5341 jcol, i, MEXPAND_DOW(row2->entry.real_dd[k]));
5342 }
5343 row2 = NULL;
5344 break;
5345 case MATENT_REAL_D:
5346 if (DST2_DOW(row2->entry.real_d[k],
5347 row->entry.real_d[j]) > 1.E-10) {
5348 non_symmetric = 1;
5349 MSG("mat[%d,%d]="DMFORMAT_DOW \
5350 " != mat[%d,%d]="DMFORMAT_DOW"\n",
5351 i, jcol, DMEXPAND_DOW(row->entry.real_d[j]),
5352 jcol, i, DMEXPAND_DOW(row2->entry.real_d[k]));
5353 }
5354 row2 = NULL;
5355 break;
5356 case MATENT_REAL:
5357 if (fabs(row2->entry.real[k]
5358 -
5359 row->entry.real[j]) > 1.E-10) {
5360 non_symmetric = 1;
5361 MSG("mat[%d,%d]="SCMFORMAT_DOW \
5362 " != mat[%d,%d]="SCMFORMAT_DOW"\n",
5363 i, jcol, SCMEXPAND_DOW(row->entry.real[j]),
5364 jcol, i, SCMEXPAND_DOW(row2->entry.real[k]));
5365 }
5366 row2 = NULL;
5367 break;
5368 default:
5369 ERROR_EXIT("Unknown or invalid MATENT_TYPE: %d\n",
5370 matrix->type);
5371 }
5372 if (row2 == NULL) {
5373 break;
5374 }
5375 }
5376 }
5377 }
5378 }
5379 if (!found) {
5380 non_symmetric = 1;
5381 MSG("mat[%d,%d] not found\n",jcol,i);
5382 }
5383 }
5384 }
5385 if (ABS(sum) > 1.E-5) {
5386 MSG("Row sum[%d] = %10.5le\n", i, sum);
5387 }
5388 }
5389 }
5390 if (non_symmetric) {
5391 MSG("matrix `%s' not symmetric.\n",matrix->name);
5392 WAIT;
5393 }
5394 else {
5395 MSG("matrix `%s' is symmetric.\n",matrix->name);
5396 }
5397 return;
5398 }
5399
5400 /*--------------------------------------------------------------------------*/
5401 /* the new update routines: */
5402 /*--------------------------------------------------------------------------*/
5403
5404 FORCE_INLINE_ATTR
5405 static inline void _AI_add_element_matrix(DOF_MATRIX *matrix,
5406 MATENT_TYPE mat_type,
5407 REAL factor,
5408 const EL_MATRIX *el_matrix,
5409 MATENT_TYPE elm_type,
5410 bool transpose,
5411 const EL_DOF_VEC *row_dof,
5412 const EL_DOF_VEC *col_dof,
5413 const EL_SCHAR_VEC *bound);
5414
5415 FORCE_INLINE_ATTR
__AI_check_matrix_types(MATENT_TYPE mat_type,MATENT_TYPE elm_type)5416 static inline bool __AI_check_matrix_types(MATENT_TYPE mat_type,
5417 MATENT_TYPE elm_type)
5418 {
5419 FUNCNAME("add_element_matrix");
5420 bool result = false;
5421
5422 if (mat_type == MATENT_NONE) {
5423 return true;
5424 }
5425
5426 switch (mat_type) {
5427 case MATENT_REAL:
5428 result = (elm_type == MATENT_REAL);
5429 DEBUG_TEST_EXIT(result,
5430 "Trying to add non-scalar element matrix "
5431 "to scalar DOF_MATRIX\n");
5432 break;
5433 case MATENT_REAL_D:
5434 result = (elm_type == MATENT_REAL_D || elm_type == MATENT_REAL);
5435 DEBUG_TEST_EXIT(result,
5436 "Trying to add REAL_DD element matrix to non-REAL_DD "
5437 "DOF_MATRIX\n");
5438 break;
5439 case MATENT_REAL_DD:
5440 result = (elm_type == MATENT_REAL ||
5441 elm_type == MATENT_REAL_D ||
5442 elm_type == MATENT_REAL_DD);
5443 if (!result) {
5444 ERROR_EXIT("Unsupported MATENT-type %d in element matrix\n",
5445 elm_type);
5446 }
5447 break;
5448 default:
5449 ERROR_EXIT("Unsupported MATENT-type %d in DOF_MATRIX\n", mat_type);
5450 break;
5451 }
5452 return result;
5453 }
5454
_AI_check_matrix_types(MATENT_TYPE mat_type,MATENT_TYPE elm_type)5455 bool _AI_check_matrix_types(MATENT_TYPE mat_type, MATENT_TYPE elm_type)
5456 {
5457 return __AI_check_matrix_types(mat_type, elm_type);
5458 }
5459
5460 FORCE_INLINE_ATTR
5461 static inline void
_AI_add_element_matrix_single(DOF_MATRIX * matrix,REAL factor,const EL_MATRIX * el_matrix,MatrixTranspose transpose,const EL_DOF_VEC * row_dof,const EL_DOF_VEC * col_dof,const EL_SCHAR_VEC * bound)5462 _AI_add_element_matrix_single(DOF_MATRIX *matrix,
5463 REAL factor,
5464 const EL_MATRIX *el_matrix,
5465 MatrixTranspose transpose,
5466 const EL_DOF_VEC *row_dof,
5467 const EL_DOF_VEC *col_dof,
5468 const EL_SCHAR_VEC *bound)
5469 {
5470 FUNCNAME("add_element_matrix_single");
5471 bool tr = (transpose == Transpose);
5472
5473 if (matrix->type == MATENT_NONE) {
5474 matrix->type = el_matrix->type;
5475 }
5476
5477 if (!__AI_check_matrix_types(matrix->type, el_matrix->type)) {
5478 ERROR_EXIT("Non-matching matrix/element-matrix type");
5479 }
5480
5481 switch (matrix->type) {
5482 case MATENT_REAL:
5483 _AI_add_element_matrix(matrix, MATENT_REAL,
5484 factor, el_matrix, MATENT_REAL, tr,
5485 row_dof, col_dof, bound);
5486 break;
5487 case MATENT_REAL_D:
5488 if (el_matrix->type == MATENT_REAL) {
5489 _AI_add_element_matrix(matrix, MATENT_REAL_D,
5490 factor, el_matrix, MATENT_REAL, tr,
5491 row_dof, col_dof, bound);
5492 } else {
5493 _AI_add_element_matrix(matrix, MATENT_REAL_D,
5494 factor, el_matrix, MATENT_REAL_D, tr,
5495 row_dof, col_dof, bound);
5496 }
5497 break;
5498 case MATENT_REAL_DD:
5499 switch (el_matrix->type) {
5500 case MATENT_REAL:
5501 _AI_add_element_matrix(matrix, MATENT_REAL_DD,
5502 factor, el_matrix, MATENT_REAL, tr,
5503 row_dof, col_dof, bound);
5504 break;
5505 case MATENT_REAL_D:
5506 _AI_add_element_matrix(matrix, MATENT_REAL_DD,
5507 factor, el_matrix, MATENT_REAL_D, tr,
5508 row_dof, col_dof, bound);
5509 break;
5510 case MATENT_REAL_DD:
5511 _AI_add_element_matrix(matrix, MATENT_REAL_DD,
5512 factor, el_matrix, MATENT_REAL_DD, tr,
5513 row_dof, col_dof, bound);
5514 break;
5515 default:
5516 break;
5517 }
5518 default:
5519 break;
5520 }
5521 }
5522
5523 FORCE_INLINE_ATTR
_AI_add_element_matrix(DOF_MATRIX * matrix,MATENT_TYPE mat_type,REAL factor,const EL_MATRIX * el_mat,MATENT_TYPE elm_type,bool transpose,const EL_DOF_VEC * row_dof_vec,const EL_DOF_VEC * col_dof_vec,const EL_SCHAR_VEC * bound_vec)5524 static inline void _AI_add_element_matrix(DOF_MATRIX *matrix,
5525 MATENT_TYPE mat_type,
5526 REAL factor,
5527 const EL_MATRIX *el_mat,
5528 MATENT_TYPE elm_type,
5529 bool transpose,
5530 const EL_DOF_VEC *row_dof_vec,
5531 const EL_DOF_VEC *col_dof_vec,
5532 const EL_SCHAR_VEC *bound_vec)
5533 {
5534 FUNCNAME("_AI_add_element_matrix");
5535 DOF i, j, I, J, k, irow, jcol;
5536 int free_col = 0, n_row, n_col;
5537 REAL_DD *const*el_mat_dd = el_mat->data.real_dd;
5538 REAL_D *const*el_mat_d = el_mat->data.real_d;
5539 REAL *const*el_mat_s = el_mat->data.real;
5540 const DOF *row_dof;
5541 const DOF *col_dof;
5542 const S_CHAR *bound;
5543 MATRIX_ROW *row, *free_row;
5544 bool same_space = false;
5545
5546 DEBUG_TEST_EXIT(matrix, "no matrix\n");
5547
5548 if (transpose) {
5549 n_row = el_mat->n_col;
5550 n_col = el_mat->n_row;
5551 } else {
5552 n_row = el_mat->n_row;
5553 n_col = el_mat->n_col;
5554 }
5555
5556 if (n_col == 0 || n_row == 0) {
5557 return;
5558 }
5559
5560 row_dof = row_dof_vec->vec;
5561 bound = bound_vec ? bound_vec->vec : NULL;
5562
5563 if (!col_dof_vec || n_col < 0) {
5564 col_dof = row_dof;
5565 n_col = n_row;
5566 same_space = true;
5567 } else {
5568 col_dof = col_dof_vec->vec;
5569 same_space = (row_dof_vec == col_dof_vec);
5570 }
5571
5572 if (matrix->is_diagonal) {
5573 #if 0
5574 DEBUG_TEST_EXIT(matrix->col_fe_space == NULL ||
5575 fe_space_is_eq(matrix->row_fe_space, matrix->col_fe_space),
5576 "Diagonal matrices are only supported if "
5577 "row_fe_space == col_fe_space.\n");
5578 #else
5579 DEBUG_TEST_EXIT(matrix->row_fe_space->admin->n_dof[CENTER] == 1 &&
5580 matrix->row_fe_space->admin->n_dof[VERTEX] == 0 &&
5581 matrix->row_fe_space->admin->n_dof[EDGE] == 0 &&
5582 matrix->row_fe_space->admin->n_dof[FACE] == 0 &&
5583 (matrix->col_fe_space == NULL ||
5584 (matrix->col_fe_space->admin->n_dof[CENTER] == 1 &&
5585 matrix->col_fe_space->admin->n_dof[VERTEX] == 0 &&
5586 matrix->col_fe_space->admin->n_dof[EDGE] == 0 &&
5587 matrix->col_fe_space->admin->n_dof[FACE] == 0)),
5588 "Incompatible FE-space for diagonal matrix.\n");
5589 #endif
5590 DEBUG_TEST_EXIT(n_col == 1 && n_row == 1,
5591 "This matrix \"%s\" is structurally non-diagonal\n",
5592 matrix->name);
5593
5594 if (matrix->diagonal.real == NULL) {
5595 switch (mat_type) {
5596 case MATENT_REAL:
5597 matrix->diagonal.real = get_dof_real_vec(
5598 "matrix diagonal", matrix->row_fe_space->unchained);
5599 dof_set(0.0, matrix->diagonal.real);
5600 break;
5601 case MATENT_REAL_D:
5602 matrix->diagonal.real_d = get_dof_real_d_vec(
5603 "matrix diagonal", matrix->row_fe_space->unchained);
5604 dof_set_d(0.0, matrix->diagonal.real_d);
5605 break;
5606 case MATENT_REAL_DD:
5607 matrix->diagonal.real_dd = get_dof_real_dd_vec(
5608 "matrix diagonal", matrix->row_fe_space->unchained);
5609 dof_set_dd(0.0, matrix->diagonal.real_dd);
5610 break;
5611 default:
5612 break;
5613 }
5614 if (matrix->unchained) {
5615 ((DOF_MATRIX *)matrix->unchained)->diagonal.real =
5616 matrix->diagonal.real;
5617 }
5618 matrix->n_entries = matrix->row_fe_space->admin->used_count;
5619 }
5620
5621 irow = row_dof[0];
5622 jcol = col_dof[0];
5623
5624 matrix->diag_cols->vec[irow] = jcol;
5625
5626 if (bound && bound[0] >= DIRICHLET) {
5627 switch (mat_type) {
5628 case MATENT_REAL:
5629 matrix->diagonal.real->vec[irow] = 1.0;
5630 break;
5631 case MATENT_REAL_D:
5632 DMSET_DOW(1.0, matrix->diagonal.real_d->vec[irow]);
5633 break;
5634 case MATENT_REAL_DD:
5635 DMSET_DOW(1.0, matrix->diagonal.real_d->vec[irow]);
5636 break;
5637 default:
5638 break;
5639 }
5640 return;
5641 }
5642
5643 if (factor == 0.0) {
5644 return;
5645 }
5646
5647 switch (mat_type) {
5648 case MATENT_REAL:
5649 switch (elm_type) {
5650 case MATENT_REAL:
5651 matrix->diagonal.real->vec[irow] += factor * el_mat->data.real[0][0];
5652 break;
5653 default:
5654 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5655 break;
5656 }
5657 break;
5658 case MATENT_REAL_D:
5659 switch (elm_type) {
5660 case MATENT_REAL:
5661 DMSCMAXPY_DOW(factor, el_mat->data.real[0][0],
5662 matrix->diagonal.real_d->vec[irow]);
5663 break;
5664 case MATENT_REAL_D:
5665 AXPY_DOW(factor, el_mat->data.real_d[0][0],
5666 matrix->diagonal.real_d->vec[irow]);
5667 break;
5668 default:
5669 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5670 break;
5671 }
5672 break;
5673 case MATENT_REAL_DD:
5674 switch (elm_type) {
5675 case MATENT_REAL_DD:
5676 if (transpose) {
5677 MMAXTPY_DOW(factor, (const REAL_D *)el_mat->data.real_dd[0][0],
5678 matrix->diagonal.real_dd->vec[irow]);
5679 } else {
5680 MMAXPY_DOW(factor, (const REAL_D *)el_mat->data.real_dd[0][0],
5681 matrix->diagonal.real_dd->vec[irow]);
5682 }
5683 break;
5684 case MATENT_REAL_D:
5685 MDMAXPY_DOW(factor, (const REAL *)el_mat->data.real_d[0][0],
5686 matrix->diagonal.real_dd->vec[irow]);
5687 break;
5688 case MATENT_REAL:
5689 MSCMAXPY_DOW(factor, el_mat->data.real[0][0],
5690 matrix->diagonal.real_dd->vec[irow]);
5691 break;
5692 default:
5693 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5694 break;
5695 }
5696 break;
5697 default:
5698 break;
5699 }
5700
5701 return;
5702 }
5703
5704 /* NOTE: if row_dof != col_dof DIRICHLET rows are simply skipped.
5705 */
5706
5707 for (i = 0; i < n_row; i++) {
5708 irow = row_dof[i];
5709
5710 if (matrix->matrix_row[irow] == NULL) {
5711 if (same_space) {
5712 row = matrix->matrix_row[irow] =
5713 get_matrix_row(matrix->row_fe_space, mat_type);
5714 row->col[0] = irow; /* first entry is diagonal element */
5715 ++matrix->n_entries;
5716 if (bound && bound[i] >= DIRICHLET) {
5717 #undef MAT_BODY
5718 #define MAT_BODY(F, CC, C, S, TYPE) F##SET_DOW(1.0, row->entry.S[0])
5719 MAT_EMIT_BODY_SWITCH(mat_type);
5720 continue; /* done with this row */
5721 } else {
5722 #undef MAT_BODY
5723 #define MAT_BODY(F, CC, C, S, TYPE) F##SET_DOW(0.0, row->entry.S[0])
5724 MAT_EMIT_BODY_SWITCH(mat_type);
5725 }
5726 } else {
5727 if (bound && bound[i] >= DIRICHLET) {
5728 continue;
5729 }
5730 row = matrix->matrix_row[irow] =
5731 get_matrix_row(matrix->row_fe_space, mat_type);
5732 row->col[0] = UNUSED_ENTRY;
5733 }
5734 } else if (bound && bound[i] >= DIRICHLET) {
5735 continue;
5736 }
5737
5738 if (factor == 0.0) {
5739 continue;
5740 }
5741
5742 for (j = 0; j < n_col; j++) {
5743
5744 if (transpose) {
5745 I = j; J = i;
5746 } else {
5747 I = i; J = j;
5748 }
5749
5750 #if 1
5751 /* Omit zero entries (???) Should we? */
5752 # undef MAT_BODY
5753 # define MAT_BODY(F, CC, C, S, TYPE) \
5754 if (F##CMP_DOW(0.0, CC el_mat->data.S[I][J])) { \
5755 continue; \
5756 }
5757 MAT_EMIT_BODY_SWITCH(elm_type);
5758 #endif
5759
5760 jcol = col_dof[j];
5761 row = matrix->matrix_row[irow];
5762 free_row = NULL;
5763 do {
5764 /* <<< add entries */
5765 for (k=0; k < ROW_LENGTH; k++) {
5766 if (row->col[k] == jcol) {
5767 switch(mat_type) {
5768 case MATENT_REAL_DD:
5769 switch (elm_type) {
5770 case MATENT_REAL_DD:
5771 if (transpose) { /* don't forget to transpose the blocks */
5772 MMAXTPY_DOW(factor, (const REAL_D *)el_mat_dd[I][J],
5773 row->entry.real_dd[k]);
5774 } else {
5775 MMAXPY_DOW(factor, (const REAL_D *)el_mat_dd[I][J],
5776 row->entry.real_dd[k]);
5777 }
5778 break;
5779 case MATENT_REAL_D:
5780 MDMAXPY_DOW(factor, el_mat_d[I][J], row->entry.real_dd[k]);
5781 break;
5782 case MATENT_REAL:
5783 MSCMAXPY_DOW(factor, el_mat_s[I][J], row->entry.real_dd[k]);
5784 break;
5785 default:
5786 ERROR_EXIT("Unknown MATENT_TYPE (%d)\n", elm_type);
5787 break;
5788 }
5789 break;
5790 case MATENT_REAL_D:
5791 switch (elm_type) {
5792 case MATENT_REAL_D:
5793 DMDMAXPY_DOW(factor, el_mat_d[I][J], row->entry.real_d[k]);
5794 break;
5795 case MATENT_REAL:
5796 DMSCMAXPY_DOW(factor, el_mat_s[I][J], row->entry.real_d[k]);
5797 break;
5798 default:
5799 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5800 break;
5801 }
5802 break;
5803 case MATENT_REAL:
5804 SCMSCMAXPY_DOW(factor, el_mat_s[I][J], row->entry.real[k]);
5805 break;
5806 default:
5807 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5808 break;
5809 }
5810 break; /* for-loop break */
5811 }
5812 if (ENTRY_NOT_USED(row->col[k])) {
5813 free_col = k;
5814 free_row = row;
5815 if (row->col[k] == NO_MORE_ENTRIES) {
5816 k = ROW_LENGTH;
5817 break;
5818 }
5819 }
5820 }
5821 /* >>> */
5822 if (k < ROW_LENGTH) { /* done? */
5823 break;
5824 }
5825 /* <<< allocate memory if necessary */
5826 if (row->next || free_row) {
5827 row = row->next;
5828 } else {
5829 free_row =
5830 row->next = get_matrix_row(matrix->row_fe_space, mat_type);
5831 free_col = 0;
5832 row = NULL;
5833 }
5834 /* >>> */
5835 } while (row);
5836
5837 if (k >= ROW_LENGTH) { /* not done? */
5838 /* <<< occupy new column slot */
5839 free_row->col[free_col] = jcol;
5840 ++matrix->n_entries;
5841 switch (mat_type) {
5842 case MATENT_REAL_DD:
5843 switch (elm_type) {
5844 case MATENT_REAL_DD:
5845 MMAXEY_DOW(factor, (const REAL_D *)el_mat_dd[I][J],
5846 free_row->entry.real_dd[free_col]);
5847 break;
5848 case MATENT_REAL_D:
5849 MDMAXEY_DOW(factor, el_mat_d[I][J],
5850 free_row->entry.real_dd[free_col]);
5851 break;
5852 case MATENT_REAL:
5853 MSCMAXEY_DOW(factor, el_mat_s[I][J],
5854 free_row->entry.real_dd[free_col]);
5855 break;
5856 default:
5857 ERROR_EXIT("Unknown MATENT_TYPE (%d)\n", elm_type);
5858 break;
5859 }
5860 break;
5861 case MATENT_REAL_D:
5862 switch (elm_type) {
5863 case MATENT_REAL_D:
5864 DMDMAXEY_DOW(factor, el_mat_d[I][J],
5865 free_row->entry.real_d[free_col]);
5866 break;
5867 case MATENT_REAL:
5868 DMSCMAXEY_DOW(factor, el_mat_s[I][J],
5869 free_row->entry.real_d[free_col]);
5870 break;
5871 default:
5872 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5873 break;
5874 }
5875 break;
5876 case MATENT_REAL:
5877 SCMSCMAXEY_DOW(factor, el_mat_s[I][J],
5878 free_row->entry.real[free_col]);
5879 break;
5880 default:
5881 ERROR_EXIT("Unknown or invalid MATENT_TYPE (%d)\n", elm_type);
5882 break;
5883 }
5884 /* >>> */
5885 }
5886 }
5887 }
5888 }
5889
5890 /* Careful: when adding element-matrices in transpose-mode, we must
5891 * also transpose the block-structure.
5892 *
5893 * TRANSPOSE only refers to the element matrix; row_dof and col_dof
5894 * are relative to the row- and column fe-spaces of the matrix (row ==
5895 * range fe-space, col == domain fe-space).
5896 */
5897 FLATTEN_ATTR
add_element_matrix(DOF_MATRIX * matrix,REAL factor,const EL_MATRIX * el_matrix,MatrixTranspose transpose,const EL_DOF_VEC * row_dof,const EL_DOF_VEC * col_dof,const EL_SCHAR_VEC * bound)5898 void add_element_matrix(DOF_MATRIX *matrix,
5899 REAL factor,
5900 const EL_MATRIX *el_matrix,
5901 MatrixTranspose transpose,
5902 const EL_DOF_VEC *row_dof,
5903 const EL_DOF_VEC *col_dof,
5904 const EL_SCHAR_VEC *bound)
5905 {
5906 if (transpose == NoTranspose) {
5907 COL_CHAIN_DO(matrix, DOF_MATRIX) {
5908 ROW_CHAIN_DO(matrix, DOF_MATRIX) {
5909 _AI_add_element_matrix_single(matrix, factor, el_matrix, NoTranspose,
5910 row_dof, col_dof, bound);
5911 col_dof = CHAIN_NEXT(col_dof, const EL_DOF_VEC);
5912 el_matrix = ROW_CHAIN_NEXT(el_matrix, const EL_MATRIX);
5913 } ROW_CHAIN_WHILE(matrix, DOF_MATRIX);
5914 row_dof = CHAIN_NEXT(row_dof, const EL_DOF_VEC);
5915 el_matrix = COL_CHAIN_NEXT(el_matrix, const EL_MATRIX);
5916 bound = bound != NULL ? CHAIN_NEXT(bound, EL_SCHAR_VEC) : NULL;
5917 } COL_CHAIN_WHILE(matrix, DOF_MATRIX);
5918 } else {
5919 COL_CHAIN_DO(matrix, DOF_MATRIX) {
5920 ROW_CHAIN_DO(matrix, DOF_MATRIX) {
5921 _AI_add_element_matrix_single(matrix, factor, el_matrix, Transpose,
5922 row_dof, col_dof, bound);
5923 col_dof = CHAIN_NEXT(col_dof, const EL_DOF_VEC);
5924 el_matrix = COL_CHAIN_NEXT(el_matrix, const EL_MATRIX);
5925 } ROW_CHAIN_WHILE(matrix, DOF_MATRIX);
5926 row_dof = CHAIN_NEXT(row_dof, const EL_DOF_VEC);
5927 el_matrix = ROW_CHAIN_NEXT(el_matrix, const EL_MATRIX);
5928 bound = bound != NULL ? CHAIN_NEXT(bound, EL_SCHAR_VEC) : NULL;
5929 } COL_CHAIN_WHILE(matrix, DOF_MATRIX);
5930 }
5931 }
5932
5933 static inline void
_AI_add_element_vec_single(DOF_REAL_VEC * drv,REAL factor,const EL_REAL_VEC * el_vec,const EL_DOF_VEC * dof,const EL_SCHAR_VEC * bound)5934 _AI_add_element_vec_single(DOF_REAL_VEC *drv, REAL factor,
5935 const EL_REAL_VEC *el_vec,
5936 const EL_DOF_VEC *dof,
5937 const EL_SCHAR_VEC *bound)
5938 {
5939 int i;
5940
5941 if (bound) {
5942 for (i = 0; i < el_vec->n_components; i++) {
5943 if (bound->vec[i] < DIRICHLET) {
5944 drv->vec[dof->vec[i]] += factor*el_vec->vec[i];
5945 }
5946 }
5947 } else {
5948 for (i = 0; i < el_vec->n_components; i++) {
5949 drv->vec[dof->vec[i]] += factor*el_vec->vec[i];
5950 }
5951 }
5952 }
5953
5954 static inline void
_AI_add_element_d_vec_single(DOF_REAL_D_VEC * drdv,REAL factor,const EL_REAL_D_VEC * el_vec,const EL_DOF_VEC * dof,const EL_SCHAR_VEC * bound)5955 _AI_add_element_d_vec_single(DOF_REAL_D_VEC *drdv, REAL factor,
5956 const EL_REAL_D_VEC *el_vec,
5957 const EL_DOF_VEC *dof,
5958 const EL_SCHAR_VEC *bound)
5959 {
5960 int i;
5961
5962 if (bound) {
5963 for (i = 0; i < el_vec->n_components; i++) {
5964 if (bound->vec[i] < DIRICHLET) {
5965 AXPY_DOW(factor, el_vec->vec[i], drdv->vec[dof->vec[i]]);
5966 }
5967 }
5968 } else {
5969 for (i = 0; i < el_vec->n_components; i++) {
5970 AXPY_DOW(factor, el_vec->vec[i], drdv->vec[dof->vec[i]]);
5971 }
5972 }
5973 }
5974
5975 static inline void
_AI_add_element_vec_dow_single(DOF_REAL_VEC_D * drv,REAL factor,const EL_REAL_VEC_D * el_vec,const EL_DOF_VEC * dof,const EL_SCHAR_VEC * bound)5976 _AI_add_element_vec_dow_single(DOF_REAL_VEC_D *drv, REAL factor,
5977 const EL_REAL_VEC_D *el_vec,
5978 const EL_DOF_VEC *dof,
5979 const EL_SCHAR_VEC *bound)
5980 {
5981 if (drv->stride != 1) {
5982 _AI_add_element_d_vec_single(
5983 (DOF_REAL_D_VEC *)drv, factor, (const EL_REAL_D_VEC *)el_vec, dof, bound);
5984 } else {
5985 _AI_add_element_vec_single(
5986 (DOF_REAL_VEC *)drv, factor, (const EL_REAL_VEC *)el_vec, dof, bound);
5987 }
5988 }
5989
add_element_vec(DOF_REAL_VEC * drv,REAL factor,const EL_REAL_VEC * el_vec,const EL_DOF_VEC * dof,const EL_SCHAR_VEC * bound)5990 void add_element_vec(DOF_REAL_VEC *drv, REAL factor,
5991 const EL_REAL_VEC *el_vec,
5992 const EL_DOF_VEC *dof,
5993 const EL_SCHAR_VEC *bound)
5994 {
5995 CHAIN_DO(el_vec, const EL_REAL_VEC) {
5996 _AI_add_element_vec_single(drv, factor, el_vec, dof, bound);
5997 drv = CHAIN_NEXT(drv, DOF_REAL_VEC);
5998 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
5999 bound = bound ? CHAIN_NEXT(bound, EL_SCHAR_VEC) : NULL;
6000 } CHAIN_WHILE(el_vec, const EL_REAL_VEC);
6001 }
6002
add_element_d_vec(DOF_REAL_D_VEC * drv,REAL factor,const EL_REAL_D_VEC * el_vec,const EL_DOF_VEC * dof,const EL_SCHAR_VEC * bound)6003 void add_element_d_vec(DOF_REAL_D_VEC *drv, REAL factor,
6004 const EL_REAL_D_VEC *el_vec,
6005 const EL_DOF_VEC *dof,
6006 const EL_SCHAR_VEC *bound)
6007 {
6008 CHAIN_DO(el_vec, const EL_REAL_D_VEC) {
6009 _AI_add_element_d_vec_single(drv, factor, el_vec, dof, bound);
6010 drv = CHAIN_NEXT(drv, DOF_REAL_D_VEC);
6011 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
6012 bound = bound ? CHAIN_NEXT(bound, EL_SCHAR_VEC) : NULL;
6013 } CHAIN_WHILE(el_vec, const EL_REAL_D_VEC);
6014 }
6015
add_element_vec_dow(DOF_REAL_VEC_D * drv,REAL factor,const EL_REAL_VEC_D * el_vec,const EL_DOF_VEC * dof,const EL_SCHAR_VEC * bound)6016 void add_element_vec_dow(DOF_REAL_VEC_D *drv, REAL factor,
6017 const EL_REAL_VEC_D *el_vec,
6018 const EL_DOF_VEC *dof,
6019 const EL_SCHAR_VEC *bound)
6020 {
6021 CHAIN_DO(el_vec, cosnt EL_REAL_VEC_D) {
6022 _AI_add_element_vec_dow_single(drv, factor, el_vec, dof, bound);
6023 drv = CHAIN_NEXT(drv, DOF_REAL_VEC_D);
6024 dof = CHAIN_NEXT(dof, EL_DOF_VEC);
6025 bound = bound ? CHAIN_NEXT(bound, EL_SCHAR_VEC) : NULL;
6026 } CHAIN_WHILE(el_vec, const EL_REAL_VEC_D);
6027 }
6028
update_matrix(DOF_MATRIX * dof_matrix,const EL_MATRIX_INFO * minfo,MatrixTranspose transpose)6029 void update_matrix(DOF_MATRIX *dof_matrix,
6030 const EL_MATRIX_INFO *minfo,
6031 MatrixTranspose transpose)
6032 {
6033 FUNCNAME("update_matrix");
6034 MESH *mesh;
6035 const FE_SPACE *row_fe_space, *col_fe_space, *neigh_fe_space = NULL;
6036 const BAS_FCTS *row_bfcts, *col_bfcts, *neigh_bfcts;
6037 const DOF_ADMIN *row_admin;
6038 bool use_get_bound;
6039 int dim;
6040 EL_DOF_VEC *row_dof, *col_dof, *neigh_dof = NULL;
6041 EL_SCHAR_VEC *bound = NULL;
6042 EL_BNDRY_VEC *bndry_bits = NULL;
6043 FLAGS fill_flag;
6044
6045 TEST_EXIT(minfo, "no EL_MATRIX_INFO\n");
6046 TEST_EXIT(minfo->el_matrix_fct, "no el_matrix_fct in EL_MATRIX_INFO\n");
6047 TEST_EXIT(dof_matrix, "no DOF_MATRIX\n");
6048
6049 mesh = minfo->row_fe_space->mesh;
6050
6051 COL_CHAIN_DO(dof_matrix, DOF_MATRIX) {
6052 ROW_CHAIN_DO(dof_matrix, DOF_MATRIX) {
6053 BNDRY_FLAGS_CPY(dof_matrix->dirichlet_bndry, minfo->dirichlet_bndry);
6054 } ROW_CHAIN_WHILE(dof_matrix, DOF_MATRIX);
6055 } COL_CHAIN_WHILE(dof_matrix, DOF_MATRIX);
6056
6057 col_fe_space = NULL;
6058 if (transpose == NoTranspose) {
6059 row_fe_space = minfo->row_fe_space;
6060 if (minfo->col_fe_space && minfo->col_fe_space != row_fe_space) {
6061 col_fe_space = minfo->col_fe_space;
6062 }
6063 } else {
6064 if (minfo->col_fe_space && minfo->col_fe_space != minfo->row_fe_space) {
6065 row_fe_space = minfo->col_fe_space;
6066 col_fe_space = minfo->row_fe_space;
6067 } else {
6068 row_fe_space = minfo->col_fe_space;
6069 }
6070 }
6071
6072 row_bfcts = row_fe_space->bas_fcts;
6073 row_admin = row_fe_space->admin;
6074
6075 if (col_fe_space) {
6076 col_bfcts = col_fe_space->bas_fcts;
6077 } else {
6078 col_bfcts = NULL;
6079 }
6080
6081 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(dof_matrix->dirichlet_bndry);
6082 if (use_get_bound) {
6083 fill_flag = minfo->fill_flag|FILL_BOUND;
6084 if (mesh->is_periodic && !(row_admin->flags & ADM_PERIODIC)) {
6085 fill_flag |= FILL_NON_PERIODIC;
6086 }
6087 } else {
6088 fill_flag = minfo->fill_flag;
6089 }
6090
6091 minfo->el_matrix_fct(NULL, minfo->fill_info);
6092
6093 row_dof = get_el_dof_vec(row_bfcts);
6094 if (use_get_bound) {
6095 bound = get_el_schar_vec(row_bfcts);
6096 bndry_bits = get_el_bndry_vec(row_bfcts);
6097 }
6098
6099 if (col_bfcts) {
6100 col_dof = get_el_dof_vec(col_bfcts);
6101 } else {
6102 col_dof = row_dof;
6103 }
6104
6105 if (minfo->neigh_el_mat_fcts) {
6106 neigh_fe_space = col_fe_space ? col_fe_space : row_fe_space;
6107 neigh_bfcts = neigh_fe_space->bas_fcts;
6108 neigh_dof = get_el_dof_vec(neigh_bfcts);
6109 }
6110
6111 dim = mesh->dim;
6112 TRAVERSE_FIRST(mesh, -1, fill_flag) {
6113 const EL_MATRIX *mat;
6114 const EL *el = el_info->el;
6115
6116 if ((mat = minfo->el_matrix_fct(el_info, minfo->fill_info)) == NULL) {
6117 continue;
6118 }
6119
6120 get_dof_indices(row_dof, row_fe_space, el);
6121 if (col_bfcts) {
6122 get_dof_indices(col_dof, col_fe_space, el);
6123 }
6124 if (use_get_bound) {
6125 get_bound(bndry_bits, row_bfcts, el_info);
6126 dirichlet_map(bound, bndry_bits, dof_matrix->dirichlet_bndry);
6127 }
6128 add_element_matrix(dof_matrix, minfo->factor, mat, transpose,
6129 row_dof, col_dof, use_get_bound ? bound : NULL);
6130
6131 if (minfo->neigh_el_mat_fcts) {
6132 /* Possibly a DG-discretisation, add jump contributions */
6133 int wall;
6134 for (wall = 0; wall < N_WALLS(dim); wall++) {
6135 mat = minfo->neigh_el_mat_fcts[wall](el_info, minfo->neigh_fill_info);
6136 if (mat == NULL) {
6137 continue;
6138 }
6139 /* It is always the space of _ANSATZ_ functions which belongs
6140 * to the neighbour element, so we need the col-DOFs on the
6141 * respective neighbour elements.
6142 */
6143 DEBUG_TEST_EXIT(el_info->neigh[wall] != NULL,
6144 "Jump contribution, but no neighbour????\n");
6145
6146 get_dof_indices(neigh_dof, neigh_fe_space, el_info->neigh[wall]);
6147 add_element_matrix(dof_matrix, minfo->factor, mat, transpose,
6148 row_dof, neigh_dof, use_get_bound ? bound : NULL);
6149 }
6150 }
6151
6152 } TRAVERSE_NEXT();
6153
6154 free_el_dof_vec(row_dof);
6155 if (col_bfcts) {
6156 free_el_dof_vec(col_dof);
6157 }
6158 if (minfo->neigh_el_mat_fcts) {
6159 free_el_dof_vec(neigh_dof);
6160 }
6161
6162 if (use_get_bound) {
6163 free_el_schar_vec(bound);
6164 free_el_bndry_vec(bndry_bits);
6165 }
6166
6167 }
6168
update_real_vec(DOF_REAL_VEC * drv,const EL_VEC_INFO * vec_info)6169 void update_real_vec(DOF_REAL_VEC *drv, const EL_VEC_INFO *vec_info)
6170 {
6171 FUNCNAME("update_real_vec");
6172 MESH *mesh;
6173 const BAS_FCTS *bfcts;
6174 const FE_SPACE *fe_space;
6175 const DOF_ADMIN *admin;
6176 bool use_get_bound;
6177 EL_DOF_VEC *dof;
6178 EL_SCHAR_VEC *bound = NULL;
6179 EL_BNDRY_VEC *bndry_bits = NULL;
6180 FLAGS fill_flag;
6181
6182 TEST_EXIT(vec_info,"no EL_VEC_INFO\n");
6183 TEST_EXIT(vec_info->el_vec_fct,"no el_vec_fct in EL_VEC_INFO\n");
6184 TEST_EXIT(drv,"no DOF_REAL_VEC\n");
6185
6186 fe_space = vec_info->fe_space;
6187 mesh = fe_space->mesh;
6188 bfcts = fe_space->bas_fcts;
6189 admin = fe_space->admin;
6190
6191 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(vec_info->dirichlet_bndry);
6192 if (use_get_bound) {
6193 fill_flag = vec_info->fill_flag|FILL_BOUND;
6194 if (mesh->is_periodic && !(admin->flags & ADM_PERIODIC))
6195 fill_flag |= FILL_NON_PERIODIC;
6196 } else {
6197 fill_flag = vec_info->fill_flag;
6198 }
6199
6200 vec_info->el_vec_fct(NULL, vec_info->fill_info);
6201
6202 dof = get_el_dof_vec(bfcts);
6203 if (use_get_bound) {
6204 bound = get_el_schar_vec(bfcts);
6205 bndry_bits = get_el_bndry_vec(bfcts);
6206 }
6207
6208 TRAVERSE_FIRST(mesh, -1, fill_flag) {
6209 const EL_REAL_VEC *vec;
6210
6211 vec = vec_info->el_vec_fct(el_info, vec_info->fill_info);
6212 if (vec == NULL) {
6213 continue;
6214 }
6215
6216 get_dof_indices(dof, fe_space, el_info->el);
6217 if (use_get_bound) {
6218 get_bound(bndry_bits, bfcts, el_info);
6219 dirichlet_map(bound, bndry_bits, vec_info->dirichlet_bndry);
6220 }
6221
6222 add_element_vec(drv, vec_info->factor, vec, dof,
6223 use_get_bound ? bound : NULL);
6224 } TRAVERSE_NEXT();
6225
6226 free_el_dof_vec(dof);
6227 if (use_get_bound) {
6228 free_el_schar_vec(bound);
6229 free_el_bndry_vec(bndry_bits);
6230 }
6231
6232 }
6233
update_real_d_vec(DOF_REAL_D_VEC * drdv,const EL_VEC_D_INFO * vec_info)6234 void update_real_d_vec(DOF_REAL_D_VEC *drdv, const EL_VEC_D_INFO *vec_info)
6235 {
6236 FUNCNAME("update_real_d_vec");
6237 MESH *mesh;
6238 const BAS_FCTS *bfcts;
6239 const FE_SPACE *fe_space;
6240 const DOF_ADMIN *admin;
6241 bool use_get_bound;
6242 EL_DOF_VEC *dof;
6243 EL_SCHAR_VEC *bound = NULL;
6244 EL_BNDRY_VEC *bndry_bits = NULL;
6245 FLAGS fill_flag;
6246
6247 TEST_EXIT(vec_info,"no EL_VEC_D_INFO\n");
6248 TEST_EXIT(vec_info->el_vec_fct,"no el_vec_fct in EL_VEC_D_INFO\n");
6249 TEST_EXIT(drdv,"no DOF_REAL_D_VEC\n");
6250
6251
6252 fe_space = vec_info->fe_space;
6253 mesh = fe_space->mesh;
6254 bfcts = fe_space->bas_fcts;
6255 admin = fe_space->admin;
6256
6257 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(vec_info->dirichlet_bndry);
6258 if (use_get_bound) {
6259 fill_flag = vec_info->fill_flag|FILL_BOUND;
6260 if (mesh->is_periodic && !(admin->flags & ADM_PERIODIC)) {
6261 fill_flag |= FILL_NON_PERIODIC;
6262 }
6263 } else {
6264 fill_flag = vec_info->fill_flag;
6265 }
6266
6267 vec_info->el_vec_fct(NULL, vec_info->fill_info);
6268
6269 dof = get_el_dof_vec(bfcts);
6270 if (use_get_bound) {
6271 bound = get_el_schar_vec(bfcts);
6272 bndry_bits = get_el_bndry_vec(bfcts);
6273 }
6274
6275 TRAVERSE_FIRST(mesh, -1, fill_flag) {
6276 const EL_REAL_D_VEC *vec;
6277
6278 vec = vec_info->el_vec_fct(el_info, vec_info->fill_info);
6279 if (vec == NULL) {
6280 continue;
6281 }
6282
6283 get_dof_indices(dof, fe_space, el_info->el);
6284 if (use_get_bound) {
6285 get_bound(bndry_bits, bfcts, el_info);
6286 dirichlet_map(bound, bndry_bits, vec_info->dirichlet_bndry);
6287 }
6288
6289 add_element_d_vec(drdv, vec_info->factor, vec, dof,
6290 use_get_bound ? bound : NULL);
6291
6292 } TRAVERSE_NEXT();
6293
6294 free_el_dof_vec(dof);
6295 if (use_get_bound) {
6296 free_el_schar_vec(bound);
6297 free_el_bndry_vec(bndry_bits);
6298 }
6299 }
6300
update_real_vec_dow(DOF_REAL_VEC_D * drdv,const EL_VEC_INFO_D * vec_info)6301 void update_real_vec_dow(DOF_REAL_VEC_D *drdv, const EL_VEC_INFO_D *vec_info)
6302 {
6303 FUNCNAME("update_real_vec_dow");
6304 MESH *mesh;
6305 const BAS_FCTS *bfcts;
6306 const FE_SPACE *fe_space;
6307 const DOF_ADMIN *admin;
6308 bool use_get_bound;
6309 EL_DOF_VEC *dof;
6310 EL_SCHAR_VEC *bound = NULL;
6311 EL_BNDRY_VEC *bndry_bits = NULL;
6312 FLAGS fill_flag;
6313
6314 TEST_EXIT(vec_info, "no EL_VEC_INFO_D\n");
6315 TEST_EXIT(vec_info->el_vec_fct, "no el_vec_fct in EL_VEC_INFO_D\n");
6316 TEST_EXIT(drdv, "no DOF_REAL_VEC_D\n");
6317
6318
6319 fe_space = vec_info->fe_space;
6320 mesh = fe_space->mesh;
6321 bfcts = fe_space->bas_fcts;
6322 admin = fe_space->admin;
6323
6324 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(vec_info->dirichlet_bndry);
6325 if (use_get_bound) {
6326 fill_flag = vec_info->fill_flag|FILL_BOUND;
6327 if (mesh->is_periodic && !(admin->flags & ADM_PERIODIC)) {
6328 fill_flag |= FILL_NON_PERIODIC;
6329 }
6330 } else {
6331 fill_flag = vec_info->fill_flag;
6332 }
6333
6334 vec_info->el_vec_fct(NULL, vec_info->fill_info);
6335
6336 dof = get_el_dof_vec(bfcts);
6337 if (use_get_bound) {
6338 bound = get_el_schar_vec(bfcts);
6339 bndry_bits = get_el_bndry_vec(bfcts);
6340 }
6341
6342 TRAVERSE_FIRST(mesh, -1, fill_flag) {
6343 const EL_REAL_VEC_D *vec;
6344
6345 vec = vec_info->el_vec_fct(el_info, vec_info->fill_info);
6346 if (vec == NULL) {
6347 continue;
6348 }
6349
6350 get_dof_indices(dof, fe_space, el_info->el);
6351 if (use_get_bound) {
6352 get_bound(bndry_bits, bfcts, el_info);
6353 dirichlet_map(bound, bndry_bits, vec_info->dirichlet_bndry);
6354 }
6355
6356 add_element_vec_dow(drdv, vec_info->factor, vec, dof,
6357 use_get_bound ? bound : NULL);
6358
6359 } TRAVERSE_NEXT();
6360
6361 free_el_dof_vec(dof);
6362 if (use_get_bound) {
6363 free_el_schar_vec(bound);
6364 free_el_bndry_vec(bndry_bits);
6365 }
6366 }
6367
6368 /*--------------------------------------------------------------------------*/
6369 /* attempt to find a DOF_ADMIN structure which stores vertex DOFs */
6370 /*--------------------------------------------------------------------------*/
6371
get_vertex_admin(MESH * mesh,FLAGS flags)6372 const DOF_ADMIN *get_vertex_admin(MESH *mesh, FLAGS flags)
6373 {
6374 int i, n_admin = mesh->n_dof_admin;
6375 DOF_ADMIN **admins = mesh->dof_admin;
6376 const DOF_ADMIN *admin = NULL;
6377
6378 if (!mesh->is_periodic) {
6379 flags &= ~ADM_PERIODIC;
6380 }
6381
6382 for (i = 0; i < n_admin; i++) {
6383 if (admins[i]->n_dof[VERTEX] && (admins[i]->flags == flags)) {
6384 if (!admin)
6385 admin = admins[i];
6386 else if (admins[i]->used_count < admin->used_count)
6387 admin = admins[i];
6388 }
6389 }
6390
6391 if(!admin) {
6392 const int n_dof[N_NODE_TYPES] = {1, 0, 0, 0};
6393 const FE_SPACE *fe_space;
6394
6395 fe_space = get_dof_space(mesh, "Vertex DOF admin", n_dof, flags);
6396
6397 admin = fe_space->admin;
6398
6399 free_fe_space(fe_space);
6400 }
6401
6402 return admin;
6403 }
6404
6405 /* attempt to find a minimal admin with admin->n_dof[type] >= n_dof[type] */
6406 const DOF_ADMIN *
get_minimal_admin(MESH * mesh,const int n_dof[N_NODE_TYPES],FLAGS flags)6407 get_minimal_admin(MESH *mesh, const int n_dof[N_NODE_TYPES], FLAGS flags)
6408 {
6409 int i, j, n_admin = mesh->n_dof_admin;
6410 DOF_ADMIN **admins = mesh->dof_admin;
6411 const DOF_ADMIN *admin = NULL;
6412
6413 if (!mesh->is_periodic) {
6414 flags &= ~ADM_PERIODIC;
6415 }
6416
6417 for (i = 0; i < n_admin; i++) {
6418 if (admins[i]->flags != flags) {
6419 continue;
6420 }
6421 for (j = 0; j < N_NODE_TYPES && admins[i]->n_dof[j] >= n_dof[j]; j++);
6422 if (j < N_NODE_TYPES) {
6423 continue;
6424 }
6425 if (admin == NULL) {
6426 admin = admins[i];
6427 } else if (admins[i]->used_count < admin->used_count) {
6428 admin = admins[i];
6429 }
6430 }
6431
6432 if (admin == NULL) {
6433 const FE_SPACE *fe_space;
6434
6435 fe_space = get_dof_space(mesh, "minimal admin", n_dof, flags);
6436
6437 admin = fe_space->admin;
6438
6439 free_fe_space(fe_space);
6440 }
6441
6442 return admin;
6443 }
6444
6445 /* Give a summary about the number of currently attached objects
6446 * (matrices and vectors).
6447 */
summarize_dof_admin(const DOF_ADMIN * admin)6448 void summarize_dof_admin(const DOF_ADMIN *admin)
6449 {
6450 FUNCNAME("summarize_dof_admin");
6451 unsigned cnt;
6452 DOF_INT_VEC *dof_int_vec; /* linked list of int vectors */
6453 DOF_DOF_VEC *dof_dof_vec; /* linked list of dof vectors */
6454 DOF_DOF_VEC *int_dof_vec; /* linked list of dof vectors */
6455 DOF_UCHAR_VEC *dof_uchar_vec; /* linked list of u_char vectors */
6456 DOF_SCHAR_VEC *dof_schar_vec; /* linked list of s_char vectors */
6457 DOF_REAL_VEC *dof_real_vec; /* linked list of real vectors */
6458 DOF_REAL_D_VEC *dof_real_d_vec; /* linked list of real_d vectors */
6459 DOF_PTR_VEC *dof_ptr_vec; /* linked list of void * vectors */
6460 DOF_MATRIX *dof_matrix; /* linked list of matrices */
6461
6462 MSG("DOF_ADMIN \"%s@%s\"\n", admin->name, admin->mesh->name);
6463 MSG("size : %d\n", admin->size);
6464 MSG("used_count: %d\n", admin->used_count);
6465 MSG("size_used : %d\n", admin->size_used);
6466 MSG("hole_count: %d\n", admin->hole_count);
6467 #ifdef COUNT_OBJECTS
6468 # undef COUNT_OBJECTS
6469 #endif
6470 #define COUNT_OBJECTS(name) \
6471 for (cnt = 0, name = admin->name; name != NULL; name = name->next, ++cnt); \
6472 if (cnt > 0) { \
6473 MSG(#name": %d\n", cnt); \
6474 }
6475
6476 COUNT_OBJECTS(dof_int_vec);
6477 COUNT_OBJECTS(dof_dof_vec);
6478 COUNT_OBJECTS(int_dof_vec);
6479 COUNT_OBJECTS(dof_uchar_vec);
6480 COUNT_OBJECTS(dof_schar_vec);
6481 COUNT_OBJECTS(dof_real_vec);
6482 COUNT_OBJECTS(dof_real_d_vec);
6483 COUNT_OBJECTS(dof_ptr_vec);
6484 COUNT_OBJECTS(dof_matrix);
6485 }
6486
summarize_all_admins(MESH * mesh)6487 void summarize_all_admins(MESH *mesh)
6488 {
6489 FUNCNAME("summarize_all_admins");
6490 int i;
6491
6492 MSG("******************** Admins@%s ************** \n", mesh->name);
6493 for (i = 0; i < mesh->n_dof_admin; ++i) {
6494 summarize_dof_admin(mesh->dof_admin[i]);
6495 MSG("\n");
6496 }
6497 }
6498