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