Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
BlaSchwarzSetup.c
Go to the documentation of this file.
1
15#include <math.h>
16#include <time.h>
17#include "fasp.h"
18#include "fasp_functs.h"
19
20/*---------------------------------*/
21/*-- Declare Private Functions --*/
22/*---------------------------------*/
23
24static void SWZ_level(const INT, dCSRmat*, INT*, INT*, INT*, INT*, const INT);
25static void SWZ_block(SWZ_data*, const INT, const INT*, const INT*, INT*);
26
27/*---------------------------------*/
28/*-- Public Functions --*/
29/*---------------------------------*/
30
46INT fasp_swz_dcsr_setup(SWZ_data* swzdata, SWZ_param* swzparam)
47{
48 // information about A
49 dCSRmat A = swzdata->A;
50 INT n = A.row;
51
52 INT blksolver = swzparam->SWZ_blksolver;
53 INT maxlev = swzparam->SWZ_maxlvl;
54
55 // local variables
56 INT i;
57 INT inroot = -10, nsizei = -10, nsizeall = -10, nlvl = 0;
58 INT* jb = NULL;
59 ivector MIS;
60
61 // data for Schwarz method
62 INT nblk;
63 INT *iblock = NULL, *jblock = NULL, *mask = NULL, *maxa = NULL;
64
65 // return
66 INT flag = FASP_SUCCESS;
67
68 swzdata->swzparam = swzparam;
69
70#if DEBUG_MODE > 0
71 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
72#endif
73
74 // allocate memory
75 maxa = (INT*)fasp_mem_calloc(n, sizeof(INT));
76 mask = (INT*)fasp_mem_calloc(n, sizeof(INT));
77 iblock = (INT*)fasp_mem_calloc(n, sizeof(INT));
78 jblock = (INT*)fasp_mem_calloc(n, sizeof(INT));
79
80 nsizeall = 0;
81 memset(mask, 0, sizeof(INT) * n);
82 memset(iblock, 0, sizeof(INT) * n);
83 memset(maxa, 0, sizeof(INT) * n);
84
85 maxa[0] = 0;
86
87 // select root nodes
88 MIS = fasp_sparse_mis(&A);
89
90 /*-------------------------------------------*/
91 // find the blocks
92 /*-------------------------------------------*/
93
94 // first pass: do a maxlev level sets out for each node
95 for (i = 0; i < MIS.row; i++) {
96 inroot = MIS.val[i];
97 SWZ_level(inroot, &A, mask, &nlvl, maxa, jblock, maxlev);
98 nsizei = maxa[nlvl];
99 nsizeall += nsizei;
100 }
101
102#if DEBUG_MODE > 1
103 printf("### DEBUG: nsizeall = %d\n", nsizeall);
104#endif
105
106 // calculated the size of jblock up to here
107 jblock = (INT*)fasp_mem_realloc(jblock, (nsizeall + n) * sizeof(INT));
108
109 // second pass: redo the same again, but this time we store in jblock
110 maxa[0] = 0;
111 iblock[0] = 0;
112 nsizeall = 0;
113 jb = jblock;
114 for (i = 0; i < MIS.row; i++) {
115 inroot = MIS.val[i];
116 SWZ_level(inroot, &A, mask, &nlvl, maxa, jb, maxlev);
117 nsizei = maxa[nlvl];
118 iblock[i + 1] = iblock[i] + nsizei;
119 nsizeall += nsizei;
120 jb += nsizei;
121 }
122 nblk = MIS.row;
123
124#if DEBUG_MODE > 1
125 printf("### DEBUG: nsizeall = %d, %d\n", nsizeall, iblock[nblk]);
126#endif
127
128 /*-------------------------------------------*/
129 // LU decomposition of blocks
130 /*-------------------------------------------*/
131
132 memset(mask, 0, sizeof(INT) * n);
133
134 swzdata->blk_data = (dCSRmat*)fasp_mem_calloc(nblk, sizeof(dCSRmat));
135
136 SWZ_block(swzdata, nblk, iblock, jblock, mask);
137
138 // Setup for each block solver
139 switch (blksolver) {
140
141#if WITH_MUMPS
142 case SOLVER_MUMPS:
143 {
144 /* use MUMPS direct solver on each block */
145 dCSRmat* blk = swzdata->blk_data;
146 Mumps_data* mumps =
147 (Mumps_data*)fasp_mem_calloc(nblk, sizeof(Mumps_data));
148 for (i = 0; i < nblk; ++i)
149 mumps[i] = fasp_mumps_factorize(&blk[i], NULL, NULL, PRINT_NONE);
150 swzdata->mumps = mumps;
151
152 break;
153 }
154#endif
155
156#if WITH_UMFPACK
157 case SOLVER_UMFPACK:
158 {
159 /* use UMFPACK direct solver on each block */
160 dCSRmat* blk = swzdata->blk_data;
161 void** numeric = (void**)fasp_mem_calloc(nblk, sizeof(void*));
162 dCSRmat Ac_tran;
163 for (i = 0; i < nblk; ++i) {
164 Ac_tran = fasp_dcsr_create(blk[i].row, blk[i].col, blk[i].nnz);
165 fasp_dcsr_transz(&blk[i], NULL, &Ac_tran);
166 fasp_dcsr_cp(&Ac_tran, &blk[i]);
167 numeric[i] = fasp_umfpack_factorize(&blk[i], 0);
168 }
169 swzdata->numeric = numeric;
170 fasp_dcsr_free(&Ac_tran);
171
172 break;
173 }
174#endif
175
176 default:
177 {
178 /* do nothing for iterative methods */
179 }
180 }
181
182#if DEBUG_MODE > 1
183 printf("### DEBUG: n = %d, #blocks = %d, max block size = %d\n", n, nblk,
184 swzdata->maxbs);
185#endif
186
187 /*-------------------------------------------*/
188 // return
189 /*-------------------------------------------*/
190 swzdata->nblk = nblk;
191 swzdata->iblock = iblock;
192 swzdata->jblock = jblock;
193 swzdata->mask = mask;
194 swzdata->maxa = maxa;
195 swzdata->SWZ_type = swzparam->SWZ_type;
196
197#if DEBUG_MODE > 0
198 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
199#endif
200
201 return flag;
202}
203
218void fasp_dcsr_swz_forward(SWZ_data* swzdata, SWZ_param* swzparam, dvector* x,
219 dvector* b)
220{
221 INT i, j, iblk, ki, kj, kij, is, ibl0, ibl1, nloc, iaa, iab;
222
223 // Schwarz partition
224 INT nblk = swzdata->nblk;
225 dCSRmat* blk = swzdata->blk_data;
226 INT* iblock = swzdata->iblock;
227 INT* jblock = swzdata->jblock;
228 INT* mask = swzdata->mask;
229 INT blksolver = swzparam->SWZ_blksolver;
230
231 // Schwarz data
232 dCSRmat A = swzdata->A;
233 INT* ia = A.IA;
234 INT* ja = A.JA;
235 REAL* val = A.val;
236
237 // Local solution and right hand vectors
238 dvector rhs = swzdata->rhsloc1;
239 dvector u = swzdata->xloc1;
240
241#if WITH_UMFPACK
242 void** numeric = swzdata->numeric;
243#endif
244
245#if WITH_MUMPS
246 Mumps_data* mumps = swzdata->mumps;
247#endif
248
249 for (is = 0; is < nblk; ++is) {
250 // Form the right hand of eack block
251 ibl0 = iblock[is];
252 ibl1 = iblock[is + 1];
253 nloc = ibl1 - ibl0;
254 for (i = 0; i < nloc; ++i) {
255 iblk = ibl0 + i;
256 ki = jblock[iblk];
257 mask[ki] = i + 1;
258 }
259
260 for (i = 0; i < nloc; ++i) {
261 iblk = ibl0 + i;
262 ki = jblock[iblk];
263 rhs.val[i] = b->val[ki];
264 iaa = ia[ki] - 1;
265 iab = ia[ki + 1] - 1;
266 for (kij = iaa; kij < iab; ++kij) {
267 kj = ja[kij] - 1;
268 j = mask[kj];
269 if (j == 0) {
270 rhs.val[i] -= val[kij] * x->val[kj];
271 }
272 }
273 }
274
275 // Solve each block
276 switch (blksolver) {
277
278#if WITH_MUMPS
279 case SOLVER_MUMPS:
280 {
281 /* use MUMPS direct solver on each block */
282 fasp_mumps_solve(&blk[is], &rhs, &u, mumps[is], 0);
283 break;
284 }
285#endif
286
287#if WITH_UMFPACK
288 case SOLVER_UMFPACK:
289 {
290 /* use UMFPACK direct solver on each block */
291 fasp_umfpack_solve(&blk[is], &rhs, &u, numeric[is], 0);
292 break;
293 }
294#endif
295 default:
296 /* use iterative solver on each block */
297 u.row = blk[is].row;
298 rhs.row = blk[is].row;
299 fasp_dvec_set(u.row, &u, 0);
300 fasp_solver_dcsr_pvgmres(&blk[is], &rhs, &u, NULL, 1e-8, 1e-20, 100, 20,
301 1, 0);
302 }
303
304 // zero the mask so that everyting is as it was
305 for (i = 0; i < nloc; ++i) {
306 iblk = ibl0 + i;
307 ki = jblock[iblk];
308 mask[ki] = 0;
309 x->val[ki] = u.val[i];
310 }
311 }
312}
313
328void fasp_dcsr_swz_backward(SWZ_data* swzdata, SWZ_param* swzparam, dvector* x,
329 dvector* b)
330{
331 INT i, j, iblk, ki, kj, kij, is, ibl0, ibl1, nloc, iaa, iab;
332
333 // Schwarz partition
334 INT nblk = swzdata->nblk;
335 dCSRmat* blk = swzdata->blk_data;
336 INT* iblock = swzdata->iblock;
337 INT* jblock = swzdata->jblock;
338 INT* mask = swzdata->mask;
339 INT blksolver = swzparam->SWZ_blksolver;
340
341 // Schwarz data
342 dCSRmat A = swzdata->A;
343 INT* ia = A.IA;
344 INT* ja = A.JA;
345 REAL* val = A.val;
346
347 // Local solution and right hand vectors
348 dvector rhs = swzdata->rhsloc1;
349 dvector u = swzdata->xloc1;
350
351#if WITH_UMFPACK
352 void** numeric = swzdata->numeric;
353#endif
354
355#if WITH_MUMPS
356 Mumps_data* mumps = swzdata->mumps;
357#endif
358
359 for (is = nblk - 1; is >= 0; --is) {
360 // Form the right hand of eack block
361 ibl0 = iblock[is];
362 ibl1 = iblock[is + 1];
363 nloc = ibl1 - ibl0;
364 for (i = 0; i < nloc; ++i) {
365 iblk = ibl0 + i;
366 ki = jblock[iblk];
367 mask[ki] = i + 1;
368 }
369
370 for (i = 0; i < nloc; ++i) {
371 iblk = ibl0 + i;
372 ki = jblock[iblk];
373 rhs.val[i] = b->val[ki];
374 iaa = ia[ki] - 1;
375 iab = ia[ki + 1] - 1;
376 for (kij = iaa; kij < iab; ++kij) {
377 kj = ja[kij] - 1;
378 j = mask[kj];
379 if (j == 0) {
380 rhs.val[i] -= val[kij] * x->val[kj];
381 }
382 }
383 }
384
385 // Solve each block
386 switch (blksolver) {
387
388#if WITH_MUMPS
389 case SOLVER_MUMPS:
390 {
391 /* use MUMPS direct solver on each block */
392 fasp_mumps_solve(&blk[is], &rhs, &u, mumps[is], 0);
393 break;
394 }
395#endif
396
397#if WITH_UMFPACK
398 case SOLVER_UMFPACK:
399 {
400 /* use UMFPACK direct solver on each block */
401 fasp_umfpack_solve(&blk[is], &rhs, &u, numeric[is], 0);
402 break;
403 }
404#endif
405 default:
406 /* use iterative solver on each block */
407 rhs.row = blk[is].row;
408 u.row = blk[is].row;
409 fasp_dvec_set(u.row, &u, 0);
410 fasp_solver_dcsr_pvgmres(&blk[is], &rhs, &u, NULL, 1e-8, 1e-20, 100, 20,
411 1, 0);
412 }
413
414 // zero the mask so that everyting is as it was
415 for (i = 0; i < nloc; ++i) {
416 iblk = ibl0 + i;
417 ki = jblock[iblk];
418 mask[ki] = 0;
419 x->val[ki] = u.val[i];
420 }
421 }
422}
423
424/*---------------------------------*/
425/*-- Private Functions --*/
426/*---------------------------------*/
427
445static void SWZ_level(const INT inroot, dCSRmat* A, INT* mask, INT* nlvl, INT* iblock,
446 INT* jblock, const INT maxlev)
447{
448 INT* ia = A->IA;
449 INT* ja = A->JA;
450 INT nnz = A->nnz;
451 INT i, j, lvl, lbegin, lvlend, nsize, node;
452 INT jstrt, jstop, nbr, lvsize;
453
454 // This is diagonal
455 if (ia[inroot + 1] - ia[inroot] <= 1) {
456 lvl = 0;
457 iblock[lvl] = 0;
458 jblock[iblock[lvl]] = inroot;
459 lvl++;
460 iblock[lvl] = 1;
461 } else {
462 // input node as root node (level 0)
463 lvl = 0;
464 jblock[0] = inroot;
465 lvlend = 0;
466 nsize = 1;
467 // mark root node
468 mask[inroot] = 1;
469
470 lvsize = nnz;
471
472 // form the level hierarchy for root node(level1, level2, ... maxlev)
473 while (lvsize > 0 && lvl < maxlev) {
474 lbegin = lvlend;
475 lvlend = nsize;
476 iblock[lvl] = lbegin;
477 lvl++;
478 for (i = lbegin; i < lvlend; ++i) {
479 node = jblock[i];
480 jstrt = ia[node] - 1;
481 jstop = ia[node + 1] - 1;
482 for (j = jstrt; j < jstop; ++j) {
483 nbr = ja[j] - 1;
484 if (mask[nbr] == 0) {
485 jblock[nsize] = nbr;
486 mask[nbr] = lvl;
487 nsize++;
488 }
489 }
490 }
491 lvsize = nsize - lvlend;
492 }
493
494 iblock[lvl] = nsize;
495
496 // reset mask array
497 for (i = 0; i < nsize; ++i) {
498 node = jblock[i];
499 mask[node] = 0;
500 }
501 }
502
503 *nlvl = lvl;
504}
505
521static void SWZ_block(SWZ_data* swzdata, const INT nblk, const INT* iblock,
522 const INT* jblock, INT* mask)
523{
524 INT i, j, iblk, ki, kj, kij, is, ibl0, ibl1, nloc, iaa, iab;
525 INT maxbs = 0, count, nnz;
526
527 dCSRmat A = swzdata->A;
528 dCSRmat* blk = swzdata->blk_data;
529
530 INT* ia = A.IA;
531 INT* ja = A.JA;
532 REAL* val = A.val;
533
534 // get maximal block size
535 for (is = 0; is < nblk; ++is) {
536 ibl0 = iblock[is];
537 ibl1 = iblock[is + 1];
538 nloc = ibl1 - ibl0;
539 maxbs = MAX(maxbs, nloc);
540 }
541
542 swzdata->maxbs = maxbs;
543
544 // allocate memory for each sub_block's right hand
545 swzdata->xloc1 = fasp_dvec_create(maxbs);
546 swzdata->rhsloc1 = fasp_dvec_create(maxbs);
547
548 for (is = 0; is < nblk; ++is) {
549 ibl0 = iblock[is];
550 ibl1 = iblock[is + 1];
551 nloc = ibl1 - ibl0;
552 count = 0;
553 for (i = 0; i < nloc; ++i) {
554 iblk = ibl0 + i;
555 ki = jblock[iblk];
556 iaa = ia[ki] - 1;
557 iab = ia[ki + 1] - 1;
558 count += iab - iaa;
559 mask[ki] = i + 1;
560 }
561
562 blk[is] = fasp_dcsr_create(nloc, nloc, count);
563 blk[is].IA[0] = 0;
564 nnz = 0;
565
566 for (i = 0; i < nloc; ++i) {
567 iblk = ibl0 + i;
568 ki = jblock[iblk];
569 iaa = ia[ki] - 1;
570 iab = ia[ki + 1] - 1;
571 for (kij = iaa; kij < iab; ++kij) {
572 kj = ja[kij] - 1;
573 j = mask[kj];
574 if (j != 0) {
575 blk[is].JA[nnz] = j - 1;
576 blk[is].val[nnz] = val[kij];
577 nnz++;
578 }
579 }
580 blk[is].IA[i + 1] = nnz;
581 }
582
583 blk[is].nnz = nnz;
584
585 // zero the mask so that everyting is as it was
586 for (i = 0; i < nloc; ++i) {
587 iblk = ibl0 + i;
588 ki = jblock[iblk];
589 mask[ki] = 0;
590 }
591 }
592}
593
594/*---------------------------------*/
595/*-- End of File --*/
596/*---------------------------------*/
fasp_mem_realloc
void * fasp_mem_realloc(void *oldmem, const LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: AuxMemory.c:113
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_dvec_create
dvector fasp_dvec_create(const INT m)
Create dvector data space of REAL type.
Definition: AuxVector.c:62
fasp_dvec_set
void fasp_dvec_set(INT n, dvector *x, const REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: AuxVector.c:222
fasp_swz_dcsr_setup
INT fasp_swz_dcsr_setup(SWZ_data *swzdata, SWZ_param *swzparam)
Setup phase for the Schwarz methods.
Definition: BlaSchwarzSetup.c:46
fasp_dcsr_swz_forward
void fasp_dcsr_swz_forward(SWZ_data *swzdata, SWZ_param *swzparam, dvector *x, dvector *b)
Schwarz smoother: forward sweep.
Definition: BlaSchwarzSetup.c:218
fasp_dcsr_swz_backward
void fasp_dcsr_swz_backward(SWZ_data *swzdata, SWZ_param *swzparam, dvector *x, dvector *b)
Schwarz smoother: backward sweep.
Definition: BlaSchwarzSetup.c:328
fasp_dcsr_create
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:47
fasp_dcsr_free
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:184
fasp_dcsr_transz
void fasp_dcsr_transz(dCSRmat *A, INT *p, dCSRmat *AT)
Generalized transpose of A: (n x m) matrix given in dCSRmat format.
Definition: BlaSparseCSR.c:1416
fasp_dcsr_cp
void fasp_dcsr_cp(const dCSRmat *A, dCSRmat *B)
copy a dCSRmat to a new one B=A
Definition: BlaSparseCSR.c:851
fasp_sparse_mis
ivector fasp_sparse_mis(dCSRmat *A)
Get the maximal independet set of a CSR matrix.
Definition: BlaSparseUtil.c:907
fasp_solver_dcsr_pvgmres
INT fasp_solver_dcsr_pvgmres(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during it...
Definition: KryPvgmres.c:66
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
SOLVER_MUMPS
#define SOLVER_MUMPS
Definition: fasp_const.h:125
SOLVER_UMFPACK
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
Mumps_data
Data for MUMPS interface.
Definition: fasp.h:607
SWZ_data
Data for Schwarz methods.
Definition: fasp.h:726
SWZ_data::mask
INT * mask
mask
Definition: fasp.h:769
SWZ_data::maxbs
INT maxbs
maximal block size
Definition: fasp.h:772
SWZ_data::A
dCSRmat A
pointer to the original coefficient matrix
Definition: fasp.h:731
SWZ_data::rhsloc1
dvector rhsloc1
local right hand side
Definition: fasp.h:748
SWZ_data::xloc1
dvector xloc1
local solution
Definition: fasp.h:751
SWZ_data::swzparam
SWZ_param * swzparam
param for Schwarz
Definition: fasp.h:794
SWZ_data::nblk
INT nblk
number of blocks
Definition: fasp.h:736
SWZ_data::blk_data
dCSRmat * blk_data
matrix for each partition
Definition: fasp.h:778
SWZ_data::SWZ_type
INT SWZ_type
Schwarz method type.
Definition: fasp.h:760
SWZ_data::jblock
INT * jblock
column index of blocks
Definition: fasp.h:742
SWZ_data::mumps
Mumps_data * mumps
param for MUMPS
Definition: fasp.h:791
SWZ_data::maxa
INT * maxa
maxa
Definition: fasp.h:775
SWZ_data::iblock
INT * iblock
row index of blocks
Definition: fasp.h:739
SWZ_param
Parameters for Schwarz method.
Definition: fasp.h:430
SWZ_param::SWZ_maxlvl
INT SWZ_maxlvl
maximal level for constructing the blocks
Definition: fasp.h:439
SWZ_param::SWZ_blksolver
INT SWZ_blksolver
type of Schwarz block solver
Definition: fasp.h:445
SWZ_param::SWZ_type
SHORT SWZ_type
type for Schwarz method
Definition: fasp.h:436
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
dvector::row
INT row
number of rows
Definition: fasp.h:357
ivector
Vector with n entries of INT type.
Definition: fasp.h:368
ivector::row
INT row
number of rows
Definition: fasp.h:371
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374