Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
schwarz_setup.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <time.h>
8 
9 #include "fasp.h"
10 #include "fasp_functs.h"
12 #include "mg_util.inl"
13 
14 static void Schwarz_levels (INT, dCSRmat *, INT *, INT *, INT *, INT *, INT);
15 
16 /*---------------------------------*/
17 /*-- Public Functions --*/
18 /*---------------------------------*/
19 
36  INT nblk,
37  INT *iblock,
38  INT *jblock,
39  INT *mask)
40 {
41  INT i, j, iblk, ki, kj, kij, is, ibl0, ibl1, nloc, iaa, iab;
42  INT maxbs = 0, count, nnz;
43 
44  dCSRmat A = Schwarz->A;
45  dCSRmat *blk = Schwarz->blk_data;
46 
47  INT *ia = A.IA;
48  INT *ja = A.JA;
49  REAL *val = A.val;
50 
51  // get maximal block size
52  for (is=0; is<nblk; ++is) {
53  ibl0 = iblock[is];
54  ibl1 = iblock[is+1];
55  nloc = ibl1-ibl0;
56  maxbs = MAX(maxbs, nloc);
57  }
58 
59  Schwarz->maxbs = maxbs;
60 
61  // allocate memory for each sub_block's right hand
62  Schwarz->xloc1 = fasp_dvec_create(maxbs);
63  Schwarz->rhsloc1 = fasp_dvec_create(maxbs);
64 
65  for (is=0; is<nblk; ++is) {
66  ibl0 = iblock[is];
67  ibl1 = iblock[is+1];
68  nloc = ibl1-ibl0;
69  count = 0;
70  for (i=0; i<nloc; ++i ) {
71  iblk = ibl0 + i;
72  ki = jblock[iblk];
73  iaa = ia[ki]-1;
74  iab = ia[ki+1]-1;
75  count += iab - iaa;
76  mask[ki] = i+1;
77  }
78 
79  blk[is] = fasp_dcsr_create(nloc, nloc, count);
80  blk[is].IA[0] = 0;
81  nnz = 0;
82 
83  for (i=0; i<nloc; ++i) {
84  iblk = ibl0 + i;
85  ki = jblock[iblk];
86  iaa = ia[ki]-1;
87  iab = ia[ki+1]-1;
88  for (kij = iaa; kij<iab; ++kij) {
89  kj = ja[kij]-1;
90  j = mask[kj];
91  if(j != 0) {
92  blk[is].JA[nnz] = j-1;
93  blk[is].val[nnz] = val[kij];
94  nnz ++;
95  }
96  }
97  blk[is].IA[i+1] = nnz;
98  }
99 
100  blk[is].nnz = nnz;
101 
102  // zero the mask so that everyting is as it was
103  for (i=0; i<nloc; ++i) {
104  iblk = ibl0 + i;
105  ki = jblock[iblk];
106  mask[ki] = 0;
107  }
108  }
109 }
110 
127  Schwarz_param *param)
128 {
129  // information about A
130  dCSRmat A = Schwarz->A;
131  INT n = A.row;
132 
133  INT block_solver = param->Schwarz_blksolver;
134  INT maxlev = param->Schwarz_maxlvl;
135  Schwarz->swzparam = param;
136 
137  // local variables
138  INT i;
139  INT inroot = -10, nsizei = -10, nsizeall = -10, nlvl = 0;
140  INT *jb=NULL;
141  ivector MIS;
142 
143  // data for Schwarz method
144  INT nblk;
145  INT *iblock = NULL, *jblock = NULL, *mask = NULL, *maxa = NULL;
146 
147  // return
148  INT flag = FASP_SUCCESS;
149 
150 #if DEBUG_MODE > 0
151  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
152 #endif
153 
154  // allocate memory
155  maxa = (INT *)fasp_mem_calloc(n,sizeof(INT));
156  mask = (INT *)fasp_mem_calloc(n,sizeof(INT));
157  iblock = (INT *)fasp_mem_calloc(n,sizeof(INT));
158  jblock = (INT *)fasp_mem_calloc(n,sizeof(INT));
159 
160  nsizeall=0;
161  memset(mask, 0, sizeof(INT)*n);
162  memset(iblock, 0, sizeof(INT)*n);
163  memset(maxa, 0, sizeof(INT)*n);
164 
165  maxa[0]=0;
166 
167  // select root nodes.
168  MIS = fasp_sparse_MIS(&A);
169 
170  /*-------------------------------------------*/
171  // find the blocks
172  /*-------------------------------------------*/
173 
174  // first pass: do a maxlev level sets out for each node
175  for ( i = 0; i < MIS.row; i++ ) {
176  inroot = MIS.val[i];
177  Schwarz_levels(inroot,&A,mask,&nlvl,maxa,jblock,maxlev);
178  nsizei=maxa[nlvl];
179  nsizeall+=nsizei;
180  }
181 
182 #if DEBUG_MODE > 1
183  fprintf(stdout,"### DEBUG: nsizeall = %d\n",nsizeall);
184 #endif
185 
186  /* We only calculated the size of this up to here. So we can reallocate jblock */
187  jblock = (INT *)fasp_mem_realloc(jblock,(nsizeall+n)*sizeof(INT));
188 
189  // second pass: redo the same again, but this time we store in jblock
190  maxa[0]=0;
191  iblock[0]=0;
192  nsizeall=0;
193  jb=jblock;
194  for (i=0;i<MIS.row;i++) {
195  inroot = MIS.val[i];
196  Schwarz_levels(inroot,&A,mask,&nlvl,maxa,jb,maxlev);
197  nsizei=maxa[nlvl];
198  iblock[i+1]=iblock[i]+nsizei;
199  nsizeall+=nsizei;
200  jb+=nsizei;
201  }
202  nblk = MIS.row;
203 
204 #if DEBUG_MODE > 1
205  fprintf(stdout,"### DEBUG: nsizeall = %d, %d\n",nsizeall,iblock[nblk]);
206 #endif
207 
208  /*-------------------------------------------*/
209  // LU decomposition of blocks
210  /*-------------------------------------------*/
211 
212  memset(mask, 0, sizeof(INT)*n);
213 
214  Schwarz->blk_data = (dCSRmat*)fasp_mem_calloc(nblk, sizeof(dCSRmat));
215 
216  fasp_Schwarz_get_block_matrix(Schwarz, nblk, iblock, jblock, mask);
217 
218  // Setup for each block solver
219  switch (block_solver) {
220 
221 #if WITH_MUMPS
222  case SOLVER_MUMPS: {
223  /* use MUMPS direct solver on each block */
224  dCSRmat *blk = Schwarz->blk_data;
225  Mumps_data *mumps = (Mumps_data*)fasp_mem_calloc(nblk, sizeof(Mumps_data));
226  for (i=0; i<nblk; ++i)
227  mumps[i] = fasp_mumps_factorize(&blk[i], NULL, NULL, PRINT_NONE);
228  Schwarz->mumps = mumps;
229 
230  break;
231  }
232 #endif
233 
234 #if WITH_UMFPACK
235  case SOLVER_UMFPACK: {
236  /* use UMFPACK direct solver on each block */
237  dCSRmat *blk = Schwarz->blk_data;
238  void **numeric = (void**)fasp_mem_calloc(nblk, sizeof(void*));
239  dCSRmat Ac_tran;
240  for (i=0; i<nblk; ++i) {
241  Ac_tran = fasp_dcsr_create(blk[i].row, blk[i].col, blk[i].nnz);
242  fasp_dcsr_transz(&blk[i], NULL, &Ac_tran);
243  fasp_dcsr_cp(&Ac_tran, &blk[i]);
244  numeric[i] = fasp_umfpack_factorize(&blk[i], 0);
245  }
246  Schwarz->numeric = numeric;
247  fasp_dcsr_free(&Ac_tran);
248 
249  break;
250  }
251 #endif
252 
253  default: {
254  /* do nothing for iterative methods */
255  }
256  }
257 
258 #if DEBUG_MODE > 1
259  fprintf(stdout, "### DEBUG: n = %d, #blocks = %d, max block size = %d\n",
260  n, nblk, Schwarz->maxbs);
261 #endif
262 
263  /*-------------------------------------------*/
264  // return
265  /*-------------------------------------------*/
266  Schwarz->nblk = nblk;
267  Schwarz->iblock = iblock;
268  Schwarz->jblock = jblock;
269  Schwarz->mask = mask;
270  Schwarz->maxa = maxa;
271  Schwarz->Schwarz_type = param->Schwarz_type;
272 
273 #if DEBUG_MODE > 0
274  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
275 #endif
276 
277  return flag;
278 }
279 
296  Schwarz_param *param,
297  dvector *x,
298  dvector *b)
299 {
300  INT i, j, iblk, ki, kj, kij, is, ibl0, ibl1, nloc, iaa, iab;
301 
302  // Schwarz partition
303  INT nblk = Schwarz->nblk;
304  dCSRmat *blk = Schwarz->blk_data;
305  INT *iblock = Schwarz->iblock;
306  INT *jblock = Schwarz->jblock;
307  INT *mask = Schwarz->mask;
308  INT block_solver = param->Schwarz_blksolver;
309 
310  // Schwarz data
311  dCSRmat A = Schwarz->A;
312  INT *ia = A.IA;
313  INT *ja = A.JA;
314  REAL *val = A.val;
315 
316  // Local solution and right hand vectors
317  dvector rhs = Schwarz->rhsloc1;
318  dvector u = Schwarz->xloc1;
319 
320 #if WITH_UMFPACK
321  void **numeric = Schwarz->numeric;
322 #endif
323 
324 #if WITH_MUMPS
325  Mumps_data *mumps = Schwarz->mumps;
326 #endif
327 
328  for (is=0; is<nblk; ++is) {
329  // Form the right hand of eack block
330  ibl0 = iblock[is];
331  ibl1 = iblock[is+1];
332  nloc = ibl1-ibl0;
333  for (i=0; i<nloc; ++i ) {
334  iblk = ibl0 + i;
335  ki = jblock[iblk];
336  mask[ki] = i+1;
337  }
338 
339  for (i=0; i<nloc; ++i) {
340  iblk = ibl0 + i;
341  ki = jblock[iblk];
342  rhs.val[i] = b->val[ki];
343  iaa = ia[ki]-1;
344  iab = ia[ki+1]-1;
345  for (kij = iaa; kij<iab; ++kij) {
346  kj = ja[kij]-1;
347  j = mask[kj];
348  if(j == 0) {
349  rhs.val[i] -= val[kij]*x->val[kj];
350  }
351  }
352  }
353 
354  // Solve each block
355  switch (block_solver) {
356 
357 #if WITH_MUMPS
358  case SOLVER_MUMPS: {
359  /* use MUMPS direct solver on each block */
360  fasp_mumps_solve(&blk[is], &rhs, &u, mumps[is], 0);
361  break;
362  }
363 #endif
364 
365 #if WITH_UMFPACK
366  case SOLVER_UMFPACK: {
367  /* use UMFPACK direct solver on each block */
368  fasp_umfpack_solve(&blk[is], &rhs, &u, numeric[is], 0);
369  break;
370  }
371 #endif
372  default:
373  /* use iterative solver on each block */
374  u.row = blk[is].row;
375  rhs.row = blk[is].row;
376  fasp_dvec_set(u.row, &u, 0);
377  fasp_solver_dcsr_pvgmres(&blk[is], &rhs, &u, NULL, 1e-8, 100, 20, 1, 0);
378  }
379 
380  //zero the mask so that everyting is as it was
381  for (i=0; i<nloc; ++i) {
382  iblk = ibl0 + i;
383  ki = jblock[iblk];
384  mask[ki] = 0;
385  x->val[ki] = u.val[i];
386  }
387  }
388 }
389 
406  Schwarz_param *param,
407  dvector *x,
408  dvector *b)
409 {
410  INT i, j, iblk, ki, kj, kij, is, ibl0, ibl1, nloc, iaa, iab;
411 
412  // Schwarz partition
413  INT nblk = Schwarz->nblk;
414  dCSRmat *blk = Schwarz->blk_data;
415  INT *iblock = Schwarz->iblock;
416  INT *jblock = Schwarz->jblock;
417  INT *mask = Schwarz->mask;
418  INT block_solver = param->Schwarz_blksolver;
419 
420  // Schwarz data
421  dCSRmat A = Schwarz->A;
422  INT *ia = A.IA;
423  INT *ja = A.JA;
424  REAL *val = A.val;
425 
426  // Local solution and right hand vectors
427  dvector rhs = Schwarz->rhsloc1;
428  dvector u = Schwarz->xloc1;
429 
430 #if WITH_UMFPACK
431  void **numeric = Schwarz->numeric;
432 #endif
433 
434 #if WITH_MUMPS
435  Mumps_data *mumps = Schwarz->mumps;
436 #endif
437 
438  for (is=nblk-1; is>=0; --is) {
439  // Form the right hand of eack block
440  ibl0 = iblock[is];
441  ibl1 = iblock[is+1];
442  nloc = ibl1-ibl0;
443  for (i=0; i<nloc; ++i ) {
444  iblk = ibl0 + i;
445  ki = jblock[iblk];
446  mask[ki] = i+1;
447  }
448 
449  for (i=0; i<nloc; ++i) {
450  iblk = ibl0 + i;
451  ki = jblock[iblk];
452  rhs.val[i] = b->val[ki];
453  iaa = ia[ki]-1;
454  iab = ia[ki+1]-1;
455  for (kij = iaa; kij<iab; ++kij) {
456  kj = ja[kij]-1;
457  j = mask[kj];
458  if(j == 0) {
459  rhs.val[i] -= val[kij]*x->val[kj];
460  }
461  }
462  }
463 
464  // Solve each block
465  switch (block_solver) {
466 
467 #if WITH_MUMPS
468  case SOLVER_MUMPS: {
469  /* use MUMPS direct solver on each block */
470  fasp_mumps_solve(&blk[is], &rhs, &u, mumps[is], 0);
471  break;
472  }
473 #endif
474 
475 #if WITH_UMFPACK
476  case SOLVER_UMFPACK: {
477  /* use UMFPACK direct solver on each block */
478  fasp_umfpack_solve(&blk[is], &rhs, &u, numeric[is], 0);
479  break;
480  }
481 #endif
482  default:
483  /* use iterative solver on each block */
484  rhs.row = blk[is].row;
485  u.row = blk[is].row;
486  fasp_dvec_set(u.row, &u, 0);
487  fasp_solver_dcsr_pvgmres (&blk[is], &rhs, &u, NULL, 1e-8, 100, 20, 1, 0);
488  }
489 
490  //zero the mask so that everyting is as it was
491  for (i=0; i<nloc; ++i) {
492  iblk = ibl0 + i;
493  ki = jblock[iblk];
494  mask[ki] = 0;
495  x->val[ki] = u.val[i];
496  }
497  }
498 }
499 
500 /*---------------------------------*/
501 /*-- Private Functions --*/
502 /*---------------------------------*/
503 
521 static void Schwarz_levels (INT inroot,
522  dCSRmat *A,
523  INT *mask,
524  INT *nlvl,
525  INT *iblock,
526  INT *jblock,
527  INT maxlev)
528 {
529  INT *ia = A->IA;
530  INT *ja = A->JA;
531  INT nnz = A->nnz;
532  INT i, j, lvl, lbegin, lvlend, nsize, node;
533  INT jstrt, jstop, nbr, lvsize;
534 
535  // This is diagonal
536  if (ia[inroot+1]-ia[inroot] <= 1) {
537  lvl = 0;
538  iblock[lvl] = 0;
539  jblock[iblock[lvl]] = inroot;
540  lvl ++;
541  iblock[lvl] = 1;
542  }
543  else {
544  // input node as root node (level 0)
545  lvl = 0;
546  jblock[0] = inroot;
547  lvlend = 0;
548  nsize = 1;
549  // mark root node
550  mask[inroot] = 1;
551 
552  lvsize = nnz;
553 
554  // start to form the level hierarchy for root node(level1, level2, ... maxlev)
555  while (lvsize > 0 && lvl < maxlev) {
556  lbegin = lvlend;
557  lvlend = nsize;
558  iblock[lvl] = lbegin;
559  lvl ++;
560  for(i=lbegin; i<lvlend; ++i) {
561  node = jblock[i];
562  jstrt = ia[node]-1;
563  jstop = ia[node+1]-1;
564  for (j = jstrt; j<jstop; ++j) {
565  nbr = ja[j]-1;
566  if (mask[nbr] == 0) {
567  jblock[nsize] = nbr;
568  mask[nbr] = lvl;
569  nsize ++;
570  }
571  }
572  }
573  lvsize = nsize - lvlend;
574  }
575 
576  iblock[lvl] = nsize;
577 
578  // reset mask array
579  for (i = 0; i< nsize; ++i) {
580  node = jblock[i];
581  mask[node] = 0;
582  }
583  }
584 
585  *nlvl = lvl;
586 }
587 
588 /*---------------------------------*/
589 /*-- End of File --*/
590 /*---------------------------------*/
void fasp_Schwarz_get_block_matrix(Schwarz_data *Schwarz, INT nblk, INT *iblock, INT *jblock, INT *mask)
Form Schwarz partition data.
Definition: schwarz_setup.c:35
dvector fasp_dvec_create(const INT m)
Create dvector data space of REAL type.
Definition: vec.c:56
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: sparse_csr.c:34
void fasp_dcsr_Schwarz_forward_smoother(Schwarz_data *Schwarz, Schwarz_param *param, dvector *x, dvector *b)
Schwarz smoother: forward sweep.
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
INT maxbs
maximal block size
Definition: fasp.h:551
INT * iblock
row index of blocks
Definition: fasp.h:518
INT fasp_Schwarz_setup(Schwarz_data *Schwarz, Schwarz_param *param)
Setup phase for the Schwarz methods.
Mumps_data * mumps
param for MUMPS
Definition: fasp.h:570
void fasp_dcsr_Schwarz_backward_smoother(Schwarz_data *Schwarz, Schwarz_param *param, dvector *x, dvector *b)
Schwarz smoother: backward sweep.
INT nblk
number of blocks
Definition: fasp.h:515
INT * maxa
maxa
Definition: fasp.h:554
REAL * val
nonzero entries of A
Definition: fasp.h:166
INT Schwarz_maxlvl
maximal level for constructing the blocks
Definition: fasp.h:443
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
INT Schwarz_type
Schwarz method type.
Definition: fasp.h:539
Parameters for MUMPS interface.
Definition: fasp.h:459
void fasp_dcsr_cp(dCSRmat *A, dCSRmat *B)
copy a dCSRmat to a new one B=A
Definition: sparse_csr.c:723
void fasp_dcsr_transz(dCSRmat *A, INT *p, dCSRmat *AT)
Generalized transpose of A: (n x m) matrix given in dCSRmat format.
Definition: sparse_csr.c:1366
Schwarz_param * swzparam
param for Schwarz
Definition: fasp.h:573
INT nnz
number of nonzero entries
Definition: fasp.h:157
SHORT Schwarz_type
type for Schwarz method
Definition: fasp.h:440
Main header file for FASP.
void * fasp_mem_realloc(void *oldmem, LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: memory.c:110
#define SOLVER_MUMPS
Definition: fasp_const.h:125
ivector fasp_sparse_MIS(dCSRmat *A)
get the maximal independet set of a CSR matrix
Definition: sparse_util.c:909
INT * val
actual vector entries
Definition: fasp.h:362
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
INT * mask
mask
Definition: fasp.h:548
dvector rhsloc1
local right hand side
Definition: fasp.h:527
INT row
row number of matrix A, m
Definition: fasp.h:151
INT Schwarz_blksolver
type of Schwarz block solver
Definition: fasp.h:449
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT row
number of rows
Definition: fasp.h:345
dCSRmat A
pointer to the matrix
Definition: fasp.h:510
Parameters for Schwarz method.
Definition: fasp.h:434
Data for Schwarz methods.
Definition: fasp.h:505
INT count
INT row
number of rows
Definition: fasp.h:359
dCSRmat * blk_data
matrix for each partition
Definition: fasp.h:557
void fasp_dvec_set(INT n, dvector *x, REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: vec.c:235
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
INT fasp_solver_dcsr_pvgmres(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during th...
Definition: pvgmres.c:51
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:166
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
dvector xloc1
local solution
Definition: fasp.h:530
INT * jblock
column index of blocks
Definition: fasp.h:521