Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
precond_bcsr.c
Go to the documentation of this file.
1 
6 #include "fasp.h"
7 #include "fasp_block.h"
8 #include "fasp_functs.h"
9 
10 /*---------------------------------*/
11 /*-- Public Functions --*/
12 /*---------------------------------*/
13 
27  REAL *z,
28  void *data)
29 {
30 
31  precond_block_data *precdata=(precond_block_data *)data;
32  dCSRmat *A_diag = precdata->A_diag;
33  dvector *tempr = &(precdata->r);
34 
35  const INT N0 = A_diag[0].row;
36  const INT N1 = A_diag[1].row;
37  const INT N2 = A_diag[2].row;
38  const INT N = N0 + N1 + N2;
39 
40  // back up r, setup z;
41  fasp_array_cp(N, r, tempr->val);
42  fasp_array_set(N, z, 0.0);
43 
44  // prepare
45 #if WITH_UMFPACK
46  void **LU_diag = precdata->LU_diag;
47  dvector r0, r1, r2, z0, z1, z2;
48 
49  r0.row = N0; z0.row = N0;
50  r1.row = N1; z1.row = N1;
51  r2.row = N2; z2.row = N2;
52 
53  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
54  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
55 #elif WITH_SuperLU
56  dvector r0, r1, r2, z0, z1, z2;
57 
58  r0.row = N0; z0.row = N0;
59  r1.row = N1; z1.row = N1;
60  r2.row = N2; z2.row = N2;
61 
62  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
63  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
64 #endif
65 
66  // Preconditioning A00 block
67 #if WITH_UMFPACK
68  /* use UMFPACK direct solver */
69  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
70 #elif WITH_SuperLU
71  /* use SuperLU direct solver on the coarsest level */
72  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
73 #endif
74 
75  // Preconditioning A11 block
76 #if WITH_UMFPACK
77  /* use UMFPACK direct solver */
78  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
79 #elif WITH_SuperLU
80  /* use SuperLU direct solver on the coarsest level */
81  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
82 #endif
83 
84  // Preconditioning A22 block
85 #if WITH_UMFPACK
86  /* use UMFPACK direct solver */
87  fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
88 #elif WITH_SuperLU
89  /* use SuperLU direct solver on the coarsest level */
90  fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
91 #endif
92 
93  // restore r
94  fasp_array_cp(N, tempr->val, r);
95 
96 }
97 
111  REAL *z,
112  void *data)
113 {
114 
115  precond_block_data *precdata=(precond_block_data *)data;
116  block_dCSRmat *A = precdata->Abcsr;
117  dvector *tempr = &(precdata->r);
118 
119  AMG_param *amgparam = precdata->amgparam;
120  AMG_data **mgl = precdata->mgl;
121 
122  const INT N0 = A->blocks[0]->row;
123  const INT N1 = A->blocks[4]->row;
124  const INT N2 = A->blocks[8]->row;
125  const INT N = N0 + N1 + N2;
126 
127  // back up r, setup z;
128  fasp_array_cp(N, r, tempr->val);
129  fasp_array_set(N, z, 0.0);
130 
131  // prepare
132  dvector r0, r1, r2, z0, z1, z2;
133  r0.row = N0; z0.row = N0; r1.row = N1; z1.row = N1; r2.row = N2; z2.row = N2;
134  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); z0.val = z;
135  z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
136 
137  // Preconditioning A00 block
138  mgl[0]->b.row=N0; fasp_array_cp(N0, r0.val, mgl[0]->b.val); // residual is an input
139  mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
140 
141  fasp_solver_mgcycle(mgl[0], amgparam);
142  fasp_array_cp(N0, mgl[0]->x.val, z0.val);
143 
144  // Preconditioning A11 block
145  mgl[1]->b.row=N1; fasp_array_cp(N1, r1.val, mgl[1]->b.val); // residual is an input
146  mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
147 
148  fasp_solver_mgcycle(mgl[1], amgparam);
149  fasp_array_cp(N1, mgl[1]->x.val, z1.val);
150 
151  // Preconditioning A22 block
152  mgl[2]->b.row=N2; fasp_array_cp(N2, r2.val, mgl[2]->b.val); // residual is an input
153  mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
154 
155  fasp_solver_mgcycle(mgl[2], amgparam);
156  fasp_array_cp(N2, mgl[2]->x.val, z2.val);
157 
158  // restore r
159  fasp_array_cp(N, tempr->val, r);
160 
161 }
162 
176  REAL *z,
177  void *data)
178 {
179 
180  precond_block_data *precdata=(precond_block_data *)data;
181  dCSRmat *A_diag = precdata->A_diag;
182  dvector *tempr = &(precdata->r);
183 
184  const INT N0 = A_diag[0].row;
185  const INT N1 = A_diag[1].row;
186  const INT N2 = A_diag[2].row;
187  const INT N3 = A_diag[3].row;
188  const INT N = N0 + N1 + N2 + N3;
189 
190  // back up r, setup z;
191  fasp_array_cp(N, r, tempr->val);
192  fasp_array_set(N, z, 0.0);
193 
194  // prepare
195 #if WITH_UMFPACK
196  void **LU_diag = precdata->LU_diag;
197  dvector r0, r1, r2, r3, z0, z1, z2, z3;
198 
199  r0.row = N0; z0.row = N0;
200  r1.row = N1; z1.row = N1;
201  r2.row = N2; z2.row = N2;
202  r3.row = N3; z3.row = N3;
203 
204  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); r3.val = &(r[N0+N1+N2]);
205  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]); z3.val = &(z[N0+N1+N2]);
206 #elif WITH_SuperLU
207  dvector r0, r1, r2, r3, z0, z1, z2, z3;
208 
209  r0.row = N0; z0.row = N0;
210  r1.row = N1; z1.row = N1;
211  r2.row = N2; z2.row = N2;
212  r3.row = N3; z3.row = N3;
213 
214  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); r3.val = &(r[N0+N1+N2]);
215  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]); z3.val = &(z[N0+N1+N2]);
216 #endif
217 
218  // Preconditioning A00 block
219 #if WITH_UMFPACK
220  /* use UMFPACK direct solver */
221  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
222 #elif WITH_SuperLU
223  /* use SuperLU direct solver on the coarsest level */
224  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
225 #endif
226 
227  // Preconditioning A11 block
228 #if WITH_UMFPACK
229  /* use UMFPACK direct solver */
230  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
231 #elif WITH_SuperLU
232  /* use SuperLU direct solver on the coarsest level */
233  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
234 #endif
235 
236  // Preconditioning A22 block
237 #if WITH_UMFPACK
238  /* use UMFPACK direct solver */
239  fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
240 #elif WITH_SuperLU
241  /* use SuperLU direct solver on the coarsest level */
242  fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
243 #endif
244 
245  // Preconditioning A33 block
246 #if WITH_UMFPACK
247  /* use UMFPACK direct solver */
248  fasp_umfpack_solve(&A_diag[3], &r3, &z3, LU_diag[3], 0);
249 #elif WITH_SuperLU
250  /* use SuperLU direct solver on the coarsest level */
251  fasp_solver_superlu(&A_diag[3], &r3, &z3, 0);
252 #endif
253 
254  // restore r
255  fasp_array_cp(N, tempr->val, r);
256 
257 }
258 
272  REAL *z,
273  void *data)
274 {
275 
276  precond_block_data *precdata=(precond_block_data *)data;
277  block_dCSRmat *A = precdata->Abcsr;
278  dCSRmat *A_diag = precdata->A_diag;
279  void **LU_diag = precdata->LU_diag;
280 
281  dvector *tempr = &(precdata->r);
282 
283  const INT N0 = A_diag[0].row;
284  const INT N1 = A_diag[1].row;
285  const INT N2 = A_diag[2].row;
286  const INT N = N0 + N1 + N2;
287 
288  // back up r, setup z;
289  fasp_array_cp(N, r, tempr->val);
290  fasp_array_set(N, z, 0.0);
291 
292  // prepare
293  dvector r0, r1, r2, z0, z1, z2;
294 
295  r0.row = N0; z0.row = N0;
296  r1.row = N1; z1.row = N1;
297  r2.row = N2; z2.row = N2;
298 
299  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
300  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
301 
302  // Preconditioning A00 block
303 #if WITH_UMFPACK
304  /* use UMFPACK direct solver */
305  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
306 #elif WITH_SuperLU
307  /* use SuperLU direct solver on the coarsest level */
308  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
309 #endif
310 
311  // r1 = r1 - A3*z0
312  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
313 
314  // Preconditioning A11 block
315 #if WITH_UMFPACK
316  /* use UMFPACK direct solver */
317  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
318 #elif WITH_SuperLU
319  /* use SuperLU direct solver on the coarsest level */
320  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
321 #endif
322 
323  // r2 = r2 - A6*z0 - A7*z1
324  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
325  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
326 
327  // Preconditioning A22 block
328 #if WITH_UMFPACK
329  /* use UMFPACK direct solver */
330  fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
331 #elif WITH_SuperLU
332  /* use SuperLU direct solver on the coarsest level */
333  fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
334 #endif
335 
336  // restore r
337  fasp_array_cp(N, tempr->val, r);
338 
339 }
340 
354  REAL *z,
355  void *data)
356 {
357 
358  precond_block_data *precdata=(precond_block_data *)data;
359  block_dCSRmat *A = precdata->Abcsr;
360  dvector *tempr = &(precdata->r);
361 
362  AMG_param *amgparam = precdata->amgparam;
363  AMG_data **mgl = precdata->mgl;
364 
365  const INT N0 = A->blocks[0]->row;
366  const INT N1 = A->blocks[4]->row;
367  const INT N2 = A->blocks[8]->row;
368  const INT N = N0 + N1 + N2;
369 
370  INT i;
371 
372  // back up r, setup z;
373  fasp_array_cp(N, r, tempr->val);
374  fasp_array_set(N, z, 0.0);
375 
376  // prepare
377  dvector r0, r1, r2, z0, z1, z2;
378  r0.row = N0; z0.row = N0; r1.row = N1; z1.row = N1; r2.row = N2; z2.row = N2;
379  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); z0.val = z;
380  z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
381 
382  // Preconditioning A00 block
383  mgl[0]->b.row=N0; fasp_array_cp(N0, r0.val, mgl[0]->b.val); // residual is an input
384  mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
385 
386  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
387  fasp_array_cp(N0, mgl[0]->x.val, z0.val);
388 
389  // r1 = r1 - A10*z0
390  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
391 
392  // Preconditioning A11 block
393  mgl[1]->b.row=N1; fasp_array_cp(N1, r1.val, mgl[1]->b.val); // residual is an input
394  mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
395 
396  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
397  fasp_array_cp(N1, mgl[1]->x.val, z1.val);
398 
399  // r2 = r2 - A20*z0 - A21*z1
400  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
401  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
402 
403  // Preconditioning A22 block
404  mgl[2]->b.row=N2; fasp_array_cp(N2, r2.val, mgl[2]->b.val); // residual is an input
405  mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
406 
407  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
408  fasp_array_cp(N2, mgl[2]->x.val, z2.val);
409 
410  // restore r
411  fasp_array_cp(N, tempr->val, r);
412 
413 }
414 
428  REAL *z,
429  void *data)
430 {
431 
432  precond_block_data *precdata=(precond_block_data *)data;
433  block_dCSRmat *A = precdata->Abcsr;
434  dCSRmat *A_diag = precdata->A_diag;
435  void **LU_diag = precdata->LU_diag;
436 
437  dvector *tempr = &(precdata->r);
438 
439  const INT N0 = A_diag[0].row;
440  const INT N1 = A_diag[1].row;
441  const INT N2 = A_diag[2].row;
442  const INT N3 = A_diag[3].row;
443  const INT N = N0 + N1 + N2 + N3;
444 
445  // back up r, setup z;
446  fasp_array_cp(N, r, tempr->val);
447  fasp_array_set(N, z, 0.0);
448 
449  // prepare
450  dvector r0, r1, r2, r3, z0, z1, z2, z3;
451 
452  r0.row = N0; z0.row = N0;
453  r1.row = N1; z1.row = N1;
454  r2.row = N2; z2.row = N2;
455  r3.row = N3; z3.row = N3;
456 
457  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); r3.val = &(r[N0+N1+N2]);
458  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]); z3.val = &(z[N0+N1+N2]);
459 
460  // Preconditioning A00 block
461 #if WITH_UMFPACK
462  /* use UMFPACK direct solver */
463  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
464 #elif WITH_SuperLU
465  /* use SuperLU direct solver on the coarsest level */
466  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
467 #endif
468 
469  // r1 = r1 - A4*z0
470  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[4], z0.val, r1.val);
471 
472  // Preconditioning A11 block
473 #if WITH_UMFPACK
474  /* use UMFPACK direct solver */
475  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
476 #elif WITH_SuperLU
477  /* use SuperLU direct solver on the coarsest level */
478  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
479 #endif
480 
481  // r2 = r2 - A8*z0 - A9*z1
482  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[8], z0.val, r2.val);
483  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[9], z1.val, r2.val);
484 
485  // Preconditioning A22 block
486 #if WITH_UMFPACK
487  /* use UMFPACK direct solver */
488  fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
489 #elif WITH_SuperLU
490  /* use SuperLU direct solver on the coarsest level */
491  fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
492 #endif
493 
494  // r3 = r3 - A12*z0 - A13*z1-A14*z2
495  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[12], z0.val, r3.val);
496  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[13], z1.val, r3.val);
497  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[14], z2.val, r3.val);
498 
499  // Preconditioning A33 block
500 #if WITH_UMFPACK
501  /* use UMFPACK direct solver */
502  fasp_umfpack_solve(&A_diag[3], &r3, &z3, LU_diag[3], 0);
503 #elif WITH_SuperLU
504  /* use SuperLU direct solver on the coarsest level */
505  fasp_solver_superlu(&A_diag[3], &r3, &z3, 0);
506 #endif
507 
508  // restore r
509  fasp_array_cp(N, tempr->val, r);
510 
511 }
512 
526  REAL *z,
527  void *data)
528 {
529 
530  precond_block_data *precdata=(precond_block_data *)data;
531  block_dCSRmat *A = precdata->Abcsr;
532  dCSRmat *A_diag = precdata->A_diag;
533  void **LU_diag = precdata->LU_diag;
534 
535  dvector *tempr = &(precdata->r);
536 
537  const INT N0 = A_diag[0].row;
538  const INT N1 = A_diag[1].row;
539  const INT N2 = A_diag[2].row;
540  const INT N = N0 + N1 + N2;
541 
542  // back up r, setup z;
543  fasp_array_cp(N, r, tempr->val);
544  fasp_array_set(N, z, 0.0);
545 
546  // prepare
547  dvector r0, r1, r2, z0, z1, z2;
548 
549  r0.row = N0; z0.row = N0;
550  r1.row = N1; z1.row = N1;
551  r2.row = N2; z2.row = N2;
552 
553  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
554  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
555 
556  // Preconditioning A22 block
557 #if WITH_UMFPACK
558  /* use UMFPACK direct solver */
559  fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
560 #elif WITH_SuperLU
561  /* use SuperLU direct solver on the coarsest level */
562  fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
563 #endif
564 
565  // r1 = r1 - A5*z2
566  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
567 
568  // Preconditioning A11 block
569 #if WITH_UMFPACK
570  /* use UMFPACK direct solver */
571  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
572 #elif WITH_SuperLU
573  /* use SuperLU direct solver on the coarsest level */
574  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
575 #endif
576 
577  // r0 = r0 - A1*z1 - A2*z2
578  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
579  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
580 
581  // Preconditioning A00 block
582 #if WITH_UMFPACK
583  /* use UMFPACK direct solver */
584  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
585 #elif WITH_SuperLU
586  /* use SuperLU direct solver on the coarsest level */
587  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
588 #endif
589 
590  // restore r
591  fasp_array_cp(N, tempr->val, r);
592 
593 }
594 
608  REAL *z,
609  void *data)
610 {
611 
612  precond_block_data *precdata=(precond_block_data *)data;
613  block_dCSRmat *A = precdata->Abcsr;
614  dCSRmat *A_diag = precdata->A_diag;
615  //void **LU_diag = precdata->LU_diag;
616 
617  AMG_param *amgparam = precdata->amgparam;
618  AMG_data **mgl = precdata->mgl;
619 
620  dvector *tempr = &(precdata->r);
621 
622  const INT N0 = A_diag[0].row;
623  const INT N1 = A_diag[1].row;
624  const INT N2 = A_diag[2].row;
625  const INT N = N0 + N1 + N2;
626 
627  INT i;
628 
629  // back up r, setup z;
630  fasp_array_cp(N, r, tempr->val);
631  fasp_array_set(N, z, 0.0);
632 
633  // prepare
634  dvector r0, r1, r2, z0, z1, z2;
635 
636  r0.row = N0; z0.row = N0;
637  r1.row = N1; z1.row = N1;
638  r2.row = N2; z2.row = N2;
639 
640  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
641  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
642 
643  // Preconditioning A22 block
644  mgl[2]->b.row=N2; fasp_array_cp(N2, r2.val, mgl[2]->b.val); // residual is an input
645  mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
646 
647  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
648  fasp_array_cp(N2, mgl[2]->x.val, z2.val);
649 
650  // r1 = r1 - A5*z2
651  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
652 
653  // Preconditioning A11 block
654  mgl[1]->b.row=N1; fasp_array_cp(N1, r1.val, mgl[1]->b.val); // residual is an input
655  mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
656 
657  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
658  fasp_array_cp(N1, mgl[1]->x.val, z1.val);
659 
660  // r0 = r0 - A1*z1 - A2*z2
661  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
662  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
663 
664  // Preconditioning A00 block
665  mgl[0]->b.row=N0; fasp_array_cp(N0, r0.val, mgl[0]->b.val); // residual is an input
666  mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
667 
668  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
669  fasp_array_cp(N0, mgl[0]->x.val, z0.val);
670 
671  // restore r
672  fasp_array_cp(N, tempr->val, r);
673 
674 }
675 
689  REAL *z,
690  void *data)
691 {
692 
693  precond_block_data *precdata=(precond_block_data *)data;
694  block_dCSRmat *A = precdata->Abcsr;
695  dCSRmat *A_diag = precdata->A_diag;
696  void **LU_diag = precdata->LU_diag;
697 
698  dvector *tempr = &(precdata->r);
699 
700  const INT N0 = A_diag[0].row;
701  const INT N1 = A_diag[1].row;
702  const INT N2 = A_diag[2].row;
703  const INT N = N0 + N1 + N2;
704 
705  // back up r, setup z;
706  fasp_array_cp(N, r, tempr->val);
707  fasp_array_set(N, z, 0.0);
708 
709  // prepare
710  dvector r0, r1, r2, z0, z1, z2;
711 
712  r0.row = N0; z0.row = N0;
713  r1.row = N1; z1.row = N1;
714  r2.row = N2; z2.row = N2;
715 
716  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
717  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
718 
719  // Preconditioning A00 block
720 #if WITH_UMFPACK
721  /* use UMFPACK direct solver */
722  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
723 #elif WITH_SuperLU
724  /* use SuperLU direct solver on the coarsest level */
725  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
726 #endif
727 
728  // r1 = r1 - A3*z0
729  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
730 
731  // Preconditioning A11 block
732 #if WITH_UMFPACK
733  /* use UMFPACK direct solver */
734  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
735 #elif WITH_SuperLU
736  /* use SuperLU direct solver on the coarsest level */
737  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
738 #endif
739 
740  // r2 = r2 - A6*z0 - A7*z1
741  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
742  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
743 
744  // Preconditioning A22 block
745 #if WITH_UMFPACK
746  /* use UMFPACK direct solver */
747  fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
748 #elif WITH_SuperLU
749  /* use SuperLU direct solver on the coarsest level */
750  fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
751 #endif
752 
753 
754  // Preconditioning A22 block
755 //#if WITH_UMFPACK
756 // /* use UMFPACK direct solver */
757 // fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
758 //#elif WITH_SuperLU
759 // /* use SuperLU direct solver on the coarsest level */
760 // fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
761 //#endif
762 
763 
764  // r1 = r1 - A5*z2
765  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
766 
767  // Preconditioning A11 block
768 #if WITH_UMFPACK
769  /* use UMFPACK direct solver */
770  fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
771 #elif WITH_SuperLU
772  /* use SuperLU direct solver on the coarsest level */
773  fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
774 #endif
775 
776  // r0 = r0 - A1*z1 - A2*z2
777  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
778  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
779 
780  // Preconditioning A00 block
781 #if WITH_UMFPACK
782  /* use UMFPACK direct solver */
783  fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
784 #elif WITH_SuperLU
785  /* use SuperLU direct solver on the coarsest level */
786  fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
787 #endif
788 
789  // restore r
790  fasp_array_cp(N, tempr->val, r);
791 
792 }
793 
807  REAL *z,
808  void *data)
809 {
810 
811  precond_block_data *precdata=(precond_block_data *)data;
812  block_dCSRmat *A = precdata->Abcsr;
813  dCSRmat *A_diag = precdata->A_diag;
814  //void **LU_diag = precdata->LU_diag;
815 
816  AMG_param *amgparam = precdata->amgparam;
817  AMG_data **mgl = precdata->mgl;
818 
819  INT i;
820 
821  dvector *tempr = &(precdata->r);
822 
823  const INT N0 = A_diag[0].row;
824  const INT N1 = A_diag[1].row;
825  const INT N2 = A_diag[2].row;
826  const INT N = N0 + N1 + N2;
827 
828  // back up r, setup z;
829  fasp_array_cp(N, r, tempr->val);
830  fasp_array_set(N, z, 0.0);
831 
832  // prepare
833  dvector r0, r1, r2, z0, z1, z2;
834 
835  r0.row = N0; z0.row = N0;
836  r1.row = N1; z1.row = N1;
837  r2.row = N2; z2.row = N2;
838 
839  r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
840  z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
841 
842  // Preconditioning A00 block
843  mgl[0]->b.row=N0; fasp_array_cp(N0, r0.val, mgl[0]->b.val); // residual is an input
844  mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
845 
846  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
847  fasp_array_cp(N0, mgl[0]->x.val, z0.val);
848 
849  // r1 = r1 - A3*z0
850  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
851 
852  // Preconditioning A11 block
853  mgl[1]->b.row=N1; fasp_array_cp(N1, r1.val, mgl[1]->b.val); // residual is an input
854  mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
855 
856  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
857  fasp_array_cp(N1, mgl[1]->x.val, z1.val);
858 
859  // r2 = r2 - A6*z0 - A7*z1
860  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
861  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
862 
863  // Preconditioning A22 block
864  mgl[2]->b.row=N2; fasp_array_cp(N2, r2.val, mgl[2]->b.val); // residual is an input
865  mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
866 
867  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
868  fasp_array_cp(N2, mgl[2]->x.val, z2.val);
869 
870  // Preconditioning A22 block
871  /*
872  mgl[2]->b.row=N2; fasp_array_cp(N2, r2.val, mgl[2]->b.val); // residual is an input
873  mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
874 
875  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
876  fasp_array_cp(N2, mgl[2]->x.val, z2.val);
877  */
878 
879  // r1 = r1 - A5*z2
880  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
881 
882  // Preconditioning A11 block
883  mgl[1]->b.row=N1; fasp_array_cp(N1, r1.val, mgl[1]->b.val); // residual is an input
884  mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
885 
886  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
887  fasp_array_cp(N1, mgl[1]->x.val, z1.val);
888 
889  // r0 = r0 - A1*z1 - A2*z2
890  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
891  fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
892 
893  // Preconditioning A00 block
894  mgl[0]->b.row=N0; fasp_array_cp(N0, r0.val, mgl[0]->b.val); // residual is an input
895  mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
896 
897  for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
898  fasp_array_cp(N0, mgl[0]->x.val, z0.val);
899 
900  // restore r
901  fasp_array_cp(N, tempr->val, r);
902 
903 }
904 
917  REAL *z,
918  void *data)
919 {
920 
922 
923  INT NumLayers = precdata->NumLayers;
924  block_dCSRmat *A = precdata->A;
925  block_dCSRmat *Ai = precdata->Ai;
926  dCSRmat *local_A = precdata->local_A;
927  ivector *local_index = precdata->local_index;
928  void **local_LU = precdata->local_LU;
929 
930  dvector *r_backup = &(precdata->r);
931  REAL *w = precdata->w;
932 
933  // local veriables
934  INT i,l;
935  dvector temp_r;
936  dvector temp_e;
937 
938  dvector *local_r = (dvector *)fasp_mem_calloc(NumLayers, sizeof(dvector));
939  dvector *local_e = (dvector *)fasp_mem_calloc(NumLayers, sizeof(dvector));
940 
941  // calculate the size and generate block local_r and local_z
942  INT N=0;
943 
944  for (l=0;l<NumLayers; l++) {
945 
946  local_r[l].row = A->blocks[l*NumLayers+l]->row;
947  local_r[l].val = r+N;
948 
949  local_e[l].row = A->blocks[l*NumLayers+l]->col;
950  local_e[l].val = z+N;
951 
952  N = N+A->blocks[l*NumLayers+l]->col;
953 
954  }
955 
956  temp_r.val = w;
957  temp_e.val = w+N;
958 
959  // back up r, setup z;
960  fasp_array_cp(N, r, r_backup->val);
961  fasp_array_cp(N, r, z);
962 
963  // L^{-1}r
964  for (l=0; l<NumLayers-1; l++){
965 
966  temp_r.row = local_A[l].row;
967  temp_e.row = local_A[l].row;
968 
969  fasp_dvec_set(local_A[l].row, &temp_r, 0.0);
970 
971  for (i=0; i<local_e[l].row; i++){
972  temp_r.val[local_index[l].val[i]] = local_e[l].val[i];
973  }
974 
975 #if WITH_UMFPACK
976  /* use UMFPACK direct solver */
977  fasp_umfpack_solve(&local_A[l], &temp_r, &temp_e, local_LU[l], 0);
978 #elif WITH_SuperLU
979  /* use SuperLU direct solver on the coarsest level */
980  fasp_solver_superlu(&local_A[l], &temp_r, &temp_e, 0);
981 #endif
982 
983  for (i=0; i<local_r[l].row; i++){
984  local_r[l].val[i] = temp_e.val[local_index[l].val[i]];
985  }
986 
987  fasp_blas_dcsr_aAxpy(-1.0, Ai->blocks[(l+1)*NumLayers+l], local_r[l].val,
988  local_e[l+1].val);
989 
990  }
991 
992  // D^{-1}L^{-1}r
993  for (l=0; l<NumLayers; l++){
994 
995  temp_r.row = local_A[l].row;
996  temp_e.row = local_A[l].row;
997 
998  fasp_dvec_set(local_A[l].row, &temp_r, 0.0);
999 
1000  for (i=0; i<local_e[l].row; i++){
1001  temp_r.val[local_index[l].val[i]] = local_e[l].val[i];
1002  }
1003 
1004 #if WITH_UMFPACK
1005  /* use UMFPACK direct solver */
1006  fasp_umfpack_solve(&local_A[l], &temp_r, &temp_e, local_LU[l], 0);
1007 #elif WITH_SuperLU
1008  /* use SuperLU direct solver on the coarsest level */
1009  fasp_solver_superlu(&local_A[l], &temp_r, &temp_e, 0);
1010 #endif
1011 
1012  for (i=0; i<local_e[l].row; i++) {
1013  local_e[l].val[i] = temp_e.val[local_index[l].val[i]];
1014  }
1015 
1016  }
1017 
1018  // L^{-t}D^{-1}L^{-1}u
1019  for (l=NumLayers-2; l>=0; l--){
1020 
1021  temp_r.row = local_A[l].row;
1022  temp_e.row = local_A[l].row;
1023 
1024  fasp_dvec_set(local_A[l].row, &temp_r, 0.0);
1025 
1026  fasp_blas_dcsr_mxv (Ai->blocks[l*NumLayers+l+1], local_e[l+1].val, local_r[l].val);
1027 
1028  for (i=0; i<local_r[l].row; i++){
1029  temp_r.val[local_index[l].val[i]] = local_r[l].val[i];
1030  }
1031 
1032 #if WITH_UMFPACK
1033  /* use UMFPACK direct solver */
1034  fasp_umfpack_solve(&local_A[l], &temp_r, &temp_e, local_LU[l], 0);
1035 #elif WITH_SuperLU
1036  /* use SuperLU direct solver on the coarsest level */
1037  fasp_solver_superlu(&local_A[l], &temp_r, &temp_e, 0);
1038 #endif
1039 
1040  for (i=0; i<local_e[l].row; i++){
1041  local_e[l].val[i] = local_e[l].val[i] - temp_e.val[local_index[l].val[i]];
1042  }
1043 
1044  }
1045 
1046  // restore r
1047  fasp_array_cp(N, r_backup->val, r);
1048 
1049 }
1050 
1051 /*---------------------------------*/
1052 /*-- End of File --*/
1053 /*---------------------------------*/
int fasp_solver_superlu(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by SuperLU.
block_dCSRmat * A
Definition: fasp_block.h:652
void fasp_precond_block_upper_3(REAL *r, REAL *z, void *data)
block upper triangular preconditioning (3x3 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:525
dCSRmat ** blocks
blocks of dCSRmat, point to blocks[brow][bcol]
Definition: fasp_block.h:93
Parameters for AMG solver.
Definition: fasp.h:583
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
dCSRmat * A_diag
Definition: fasp_block.h:509
dvector b
pointer to the right-hand side at level level_num
Definition: fasp.h:744
AMG_data ** mgl
Definition: fasp_block.h:520
void fasp_precond_block_upper_3_amg(REAL *r, REAL *z, void *data)
block upper triangular preconditioning (3x3 block matrix, each diagonal block is solved AMG) ...
Definition: precond_bcsr.c:607
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
block_dCSRmat * Ai
Definition: fasp_block.h:653
void fasp_precond_block_SGS_3(REAL *r, REAL *z, void *data)
block symmetric GS preconditioning (3x3 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:688
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:747
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
INT * val
actual vector entries
Definition: fasp.h:362
INT row
row number of matrix A, m
Definition: fasp.h:151
void fasp_precond_block_lower_3_amg(REAL *r, REAL *z, void *data)
block lower triangular preconditioning (3x3 block matrix, each diagonal block is solved by AMG) ...
Definition: precond_bcsr.c:353
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
void fasp_precond_block_diag_3(REAL *r, REAL *z, void *data)
block diagonal preconditioning (3x3 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:26
void fasp_precond_block_diag_3_amg(REAL *r, REAL *z, void *data)
block diagonal preconditioning (3x3 block matrix, each diagonal block is solved by AMG) ...
Definition: precond_bcsr.c:110
INT row
number of rows
Definition: fasp.h:345
Block REAL CSR matrix format.
Definition: fasp_block.h:84
Data passed to the preconditioner for sweeping preconditioning.
Definition: fasp_block.h:648
void fasp_solver_mgcycle(AMG_data *mgl, AMG_param *param)
#include "forts_ns.h"
Definition: mgcycle.c:41
Data for AMG solvers.
Definition: fasp.h:722
void fasp_precond_block_SGS_3_amg(REAL *r, REAL *z, void *data)
block symmetric GS preconditioning (3x3 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:806
void fasp_blas_dcsr_aAxpy(const REAL alpha, dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_csr.c:479
void fasp_dvec_set(INT n, dvector *x, REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: vec.c:235
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
Header file for FASP block matrices.
void fasp_precond_block_diag_4(REAL *r, REAL *z, void *data)
block diagonal preconditioning (4x4 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:175
AMG_param * amgparam
Definition: fasp_block.h:521
block_dCSRmat * Abcsr
Definition: fasp_block.h:507
void fasp_precond_block_lower_3(REAL *r, REAL *z, void *data)
block lower triangular preconditioning (3x3 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:271
Data passed to the preconditioner for block preconditioning for block_dCSRmat format.
Definition: fasp_block.h:502
void fasp_precond_block_lower_4(REAL *r, REAL *z, void *data)
block lower triangular preconditioning (4x4 block matrix, each diagonal block is solved exactly) ...
Definition: precond_bcsr.c:427
void fasp_precond_sweeping(REAL *r, REAL *z, void *data)
sweeping preconditioner for Maxwell equations
Definition: precond_bcsr.c:916
void fasp_blas_dcsr_mxv(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_csr.c:225