Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreBLC.c
Go to the documentation of this file.
1
16#include "fasp.h"
17#include "fasp_block.h"
18#include "fasp_functs.h"
19
20/*---------------------------------*/
21/*-- Public Functions --*/
22/*---------------------------------*/
23
38void fasp_precond_dblc_diag_3 (REAL *r,
39 REAL *z,
40 void *data)
41{
42#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
43
44 precond_data_blc *precdata = (precond_data_blc *)data;
45 dCSRmat *A_diag = precdata->A_diag;
46 dvector *tempr = &(precdata->r);
47
48 const INT N0 = A_diag[0].row;
49 const INT N1 = A_diag[1].row;
50 const INT N2 = A_diag[2].row;
51 const INT N = N0 + N1 + N2;
52
53 // back up r, setup z;
54 fasp_darray_cp(N, r, tempr->val);
55 fasp_darray_set(N, z, 0.0);
56
57 // prepare
58#if WITH_UMFPACK
59 void **LU_diag = precdata->LU_diag;
60 dvector r0, r1, r2, z0, z1, z2;
61
62 r0.row = N0; z0.row = N0;
63 r1.row = N1; z1.row = N1;
64 r2.row = N2; z2.row = N2;
65
66 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
67 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
68#elif WITH_SuperLU
69 dvector r0, r1, r2, z0, z1, z2;
70
71 r0.row = N0; z0.row = N0;
72 r1.row = N1; z1.row = N1;
73 r2.row = N2; z2.row = N2;
74
75 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
76 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
77#endif
78
79 // Preconditioning A00 block
80#if WITH_UMFPACK
81 /* use UMFPACK direct solver */
82 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
83#elif WITH_SuperLU
84 /* use SuperLU direct solver */
85 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
86#endif
87
88 // Preconditioning A11 block
89#if WITH_UMFPACK
90 /* use UMFPACK direct solver */
91 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
92#elif WITH_SuperLU
93 /* use SuperLU direct solver */
94 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
95#endif
96
97 // Preconditioning A22 block
98#if WITH_UMFPACK
99 /* use UMFPACK direct solver */
100 fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
101#elif WITH_SuperLU
102 /* use SuperLU direct solver */
103 fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
104#endif
105
106 // restore r
107 fasp_darray_cp(N, tempr->val, r);
108
109#endif
110}
111
126void fasp_precond_dblc_diag_3_amg (REAL *r,
127 REAL *z,
128 void *data)
129{
130 precond_data_blc *precdata = (precond_data_blc *)data;
131 dBLCmat *A = precdata->Ablc;
132 dvector *tempr = &(precdata->r);
133
134 AMG_param *amgparam = precdata->amgparam;
135 AMG_data **mgl = precdata->mgl;
136
137 const INT N0 = A->blocks[0]->row;
138 const INT N1 = A->blocks[4]->row;
139 const INT N2 = A->blocks[8]->row;
140 const INT N = N0 + N1 + N2;
141
142 // back up r, setup z;
143 fasp_darray_cp(N, r, tempr->val);
144 fasp_darray_set(N, z, 0.0);
145
146 // prepare
147 dvector r0, r1, r2, z0, z1, z2;
148 r0.row = N0; z0.row = N0; r1.row = N1; z1.row = N1; r2.row = N2; z2.row = N2;
149 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); z0.val = z;
150 z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
151
152 // Preconditioning A00 block
153 mgl[0]->b.row=N0; fasp_darray_cp(N0, r0.val, mgl[0]->b.val);
154 mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
155
156 fasp_solver_mgcycle(mgl[0], amgparam);
157 fasp_darray_cp(N0, mgl[0]->x.val, z0.val);
158
159 // Preconditioning A11 block
160 mgl[1]->b.row=N1; fasp_darray_cp(N1, r1.val, mgl[1]->b.val);
161 mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
162
163 fasp_solver_mgcycle(mgl[1], amgparam);
164 fasp_darray_cp(N1, mgl[1]->x.val, z1.val);
165
166 // Preconditioning A22 block
167 mgl[2]->b.row=N2; fasp_darray_cp(N2, r2.val, mgl[2]->b.val);
168 mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
169
170 fasp_solver_mgcycle(mgl[2], amgparam);
171 fasp_darray_cp(N2, mgl[2]->x.val, z2.val);
172
173 // restore r
174 fasp_darray_cp(N, tempr->val, r);
175}
176
191void fasp_precond_dblc_diag_4 (REAL *r,
192 REAL *z,
193 void *data)
194{
195#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
196
197 precond_data_blc *precdata=(precond_data_blc *)data;
198 dCSRmat *A_diag = precdata->A_diag;
199 dvector *tempr = &(precdata->r);
200
201 const INT N0 = A_diag[0].row;
202 const INT N1 = A_diag[1].row;
203 const INT N2 = A_diag[2].row;
204 const INT N3 = A_diag[3].row;
205 const INT N = N0 + N1 + N2 + N3;
206
207 // back up r, setup z;
208 fasp_darray_cp(N, r, tempr->val);
209 fasp_darray_set(N, z, 0.0);
210
211 // prepare
212#if WITH_UMFPACK
213 void **LU_diag = precdata->LU_diag;
214 dvector r0, r1, r2, r3, z0, z1, z2, z3;
215
216 r0.row = N0; z0.row = N0;
217 r1.row = N1; z1.row = N1;
218 r2.row = N2; z2.row = N2;
219 r3.row = N3; z3.row = N3;
220
221 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); r3.val = &(r[N0+N1+N2]);
222 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]); z3.val = &(z[N0+N1+N2]);
223#elif WITH_SuperLU
224 dvector r0, r1, r2, r3, z0, z1, z2, z3;
225
226 r0.row = N0; z0.row = N0;
227 r1.row = N1; z1.row = N1;
228 r2.row = N2; z2.row = N2;
229 r3.row = N3; z3.row = N3;
230
231 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); r3.val = &(r[N0+N1+N2]);
232 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]); z3.val = &(z[N0+N1+N2]);
233#endif
234
235 // Preconditioning A00 block
236#if WITH_UMFPACK
237 /* use UMFPACK direct solver */
238 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
239#elif WITH_SuperLU
240 /* use SuperLU direct solver */
241 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
242#endif
243
244 // Preconditioning A11 block
245#if WITH_UMFPACK
246 /* use UMFPACK direct solver */
247 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
248#elif WITH_SuperLU
249 /* use SuperLU direct solver */
250 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
251#endif
252
253 // Preconditioning A22 block
254#if WITH_UMFPACK
255 /* use UMFPACK direct solver */
256 fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
257#elif WITH_SuperLU
258 /* use SuperLU direct solver */
259 fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
260#endif
261
262 // Preconditioning A33 block
263#if WITH_UMFPACK
264 /* use UMFPACK direct solver */
265 fasp_umfpack_solve(&A_diag[3], &r3, &z3, LU_diag[3], 0);
266#elif WITH_SuperLU
267 /* use SuperLU direct solver */
268 fasp_solver_superlu(&A_diag[3], &r3, &z3, 0);
269#endif
270
271 // restore r
272 fasp_darray_cp(N, tempr->val, r);
273
274#endif
275}
276
291void fasp_precond_dblc_lower_3 (REAL *r,
292 REAL *z,
293 void *data)
294{
295#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
296
297 precond_data_blc *precdata = (precond_data_blc *)data;
298 dBLCmat *A = precdata->Ablc;
299 dCSRmat *A_diag = precdata->A_diag;
300 dvector *tempr = &(precdata->r);
301
302#if WITH_UMFPACK
303 void **LU_diag = precdata->LU_diag;
304#endif
305
306 const INT N0 = A_diag[0].row;
307 const INT N1 = A_diag[1].row;
308 const INT N2 = A_diag[2].row;
309 const INT N = N0 + N1 + N2;
310
311 // back up r, setup z;
312 fasp_darray_cp(N, r, tempr->val);
313 fasp_darray_set(N, z, 0.0);
314
315 // prepare
316 dvector r0, r1, r2, z0, z1, z2;
317
318 r0.row = N0; z0.row = N0;
319 r1.row = N1; z1.row = N1;
320 r2.row = N2; z2.row = N2;
321
322 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
323 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
324
325 // Preconditioning A00 block
326#if WITH_UMFPACK
327 /* use UMFPACK direct solver */
328 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
329#elif WITH_SuperLU
330 /* use SuperLU direct solver */
331 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
332#endif
333
334 // r1 = r1 - A3*z0
335 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
336
337 // Preconditioning A11 block
338#if WITH_UMFPACK
339 /* use UMFPACK direct solver */
340 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
341#elif WITH_SuperLU
342 /* use SuperLU direct solver */
343 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
344#endif
345
346 // r2 = r2 - A6*z0 - A7*z1
347 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
348 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
349
350 // Preconditioning A22 block
351#if WITH_UMFPACK
352 /* use UMFPACK direct solver */
353 fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
354#elif WITH_SuperLU
355 /* use SuperLU direct solver */
356 fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
357#endif
358
359 // restore r
360 fasp_darray_cp(N, tempr->val, r);
361
362#endif
363}
364
379void fasp_precond_dblc_lower_3_amg (REAL *r,
380 REAL *z,
381 void *data)
382{
383 precond_data_blc *precdata = (precond_data_blc *)data;
384 dBLCmat *A = precdata->Ablc;
385 dvector *tempr = &(precdata->r);
386
387 AMG_param *amgparam = precdata->amgparam;
388 AMG_data **mgl = precdata->mgl;
389
390 const INT N0 = A->blocks[0]->row;
391 const INT N1 = A->blocks[4]->row;
392 const INT N2 = A->blocks[8]->row;
393 const INT N = N0 + N1 + N2;
394
395 INT i;
396
397 // back up r, setup z;
398 fasp_darray_cp(N, r, tempr->val);
399 fasp_darray_set(N, z, 0.0);
400
401 // prepare
402 dvector r0, r1, r2, z0, z1, z2;
403 r0.row = N0; z0.row = N0; r1.row = N1; z1.row = N1; r2.row = N2; z2.row = N2;
404 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); z0.val = z;
405 z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
406
407 // Preconditioning A00 block
408 mgl[0]->b.row=N0; fasp_darray_cp(N0, r0.val, mgl[0]->b.val);
409 mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
410
411 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
412 fasp_darray_cp(N0, mgl[0]->x.val, z0.val);
413
414 // r1 = r1 - A10*z0
415 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
416
417 // Preconditioning A11 block
418 mgl[1]->b.row=N1; fasp_darray_cp(N1, r1.val, mgl[1]->b.val);
419 mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
420
421 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
422 fasp_darray_cp(N1, mgl[1]->x.val, z1.val);
423
424 // r2 = r2 - A20*z0 - A21*z1
425 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
426 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
427
428 // Preconditioning A22 block
429 mgl[2]->b.row=N2; fasp_darray_cp(N2, r2.val, mgl[2]->b.val);
430 mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
431
432 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
433 fasp_darray_cp(N2, mgl[2]->x.val, z2.val);
434
435 // restore r
436 fasp_darray_cp(N, tempr->val, r);
437}
438
453void fasp_precond_dblc_lower_4 (REAL *r,
454 REAL *z,
455 void *data)
456{
457#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
458
459 precond_data_blc *precdata = (precond_data_blc *)data;
460 dBLCmat *A = precdata->Ablc;
461 dCSRmat *A_diag = precdata->A_diag;
462 dvector *tempr = &(precdata->r);
463
464#if WITH_UMFPACK
465 void **LU_diag = precdata->LU_diag;
466#endif
467
468 const INT N0 = A_diag[0].row;
469 const INT N1 = A_diag[1].row;
470 const INT N2 = A_diag[2].row;
471 const INT N3 = A_diag[3].row;
472 const INT N = N0 + N1 + N2 + N3;
473
474 // back up r, setup z;
475 fasp_darray_cp(N, r, tempr->val);
476 fasp_darray_set(N, z, 0.0);
477
478 // prepare
479 dvector r0, r1, r2, r3, z0, z1, z2, z3;
480
481 r0.row = N0; z0.row = N0;
482 r1.row = N1; z1.row = N1;
483 r2.row = N2; z2.row = N2;
484 r3.row = N3; z3.row = N3;
485
486 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]); r3.val = &(r[N0+N1+N2]);
487 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]); z3.val = &(z[N0+N1+N2]);
488
489 // Preconditioning A00 block
490#if WITH_UMFPACK
491 /* use UMFPACK direct solver */
492 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
493#elif WITH_SuperLU
494 /* use SuperLU direct solver */
495 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
496#endif
497
498 // r1 = r1 - A4*z0
499 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[4], z0.val, r1.val);
500
501 // Preconditioning A11 block
502#if WITH_UMFPACK
503 /* use UMFPACK direct solver */
504 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
505#elif WITH_SuperLU
506 /* use SuperLU direct solver */
507 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
508#endif
509
510 // r2 = r2 - A8*z0 - A9*z1
511 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[8], z0.val, r2.val);
512 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[9], z1.val, r2.val);
513
514 // Preconditioning A22 block
515#if WITH_UMFPACK
516 /* use UMFPACK direct solver */
517 fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
518#elif WITH_SuperLU
519 /* use SuperLU direct solver */
520 fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
521#endif
522
523 // r3 = r3 - A12*z0 - A13*z1-A14*z2
524 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[12], z0.val, r3.val);
525 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[13], z1.val, r3.val);
526 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[14], z2.val, r3.val);
527
528 // Preconditioning A33 block
529#if WITH_UMFPACK
530 /* use UMFPACK direct solver */
531 fasp_umfpack_solve(&A_diag[3], &r3, &z3, LU_diag[3], 0);
532#elif WITH_SuperLU
533 /* use SuperLU direct solver */
534 fasp_solver_superlu(&A_diag[3], &r3, &z3, 0);
535#endif
536
537 // restore r
538 fasp_darray_cp(N, tempr->val, r);
539
540#endif
541}
542
557void fasp_precond_dblc_upper_3 (REAL *r,
558 REAL *z,
559 void *data)
560{
561#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
562
563 precond_data_blc *precdata = (precond_data_blc *)data;
564 dBLCmat *A = precdata->Ablc;
565 dCSRmat *A_diag = precdata->A_diag;
566 dvector *tempr = &(precdata->r);
567
568#if WITH_UMFPACK
569 void **LU_diag = precdata->LU_diag;
570#endif
571
572 const INT N0 = A_diag[0].row;
573 const INT N1 = A_diag[1].row;
574 const INT N2 = A_diag[2].row;
575 const INT N = N0 + N1 + N2;
576
577 // back up r, setup z;
578 fasp_darray_cp(N, r, tempr->val);
579 fasp_darray_set(N, z, 0.0);
580
581 // prepare
582 dvector r0, r1, r2, z0, z1, z2;
583
584 r0.row = N0; z0.row = N0;
585 r1.row = N1; z1.row = N1;
586 r2.row = N2; z2.row = N2;
587
588 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
589 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
590
591 // Preconditioning A22 block
592#if WITH_UMFPACK
593 /* use UMFPACK direct solver */
594 fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
595#elif WITH_SuperLU
596 /* use SuperLU direct solver */
597 fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
598#endif
599
600 // r1 = r1 - A5*z2
601 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
602
603 // Preconditioning A11 block
604#if WITH_UMFPACK
605 /* use UMFPACK direct solver */
606 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
607#elif WITH_SuperLU
608 /* use SuperLU direct solver */
609 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
610#endif
611
612 // r0 = r0 - A1*z1 - A2*z2
613 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
614 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
615
616 // Preconditioning A00 block
617#if WITH_UMFPACK
618 /* use UMFPACK direct solver */
619 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
620#elif WITH_SuperLU
621 /* use SuperLU direct solver */
622 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
623#endif
624
625 // restore r
626 fasp_darray_cp(N, tempr->val, r);
627
628#endif
629}
630
645void fasp_precond_dblc_upper_3_amg (REAL *r,
646 REAL *z,
647 void *data)
648{
649 precond_data_blc *precdata = (precond_data_blc *)data;
650 dBLCmat *A = precdata->Ablc;
651 dCSRmat *A_diag = precdata->A_diag;
652
653 AMG_param *amgparam = precdata->amgparam;
654 AMG_data **mgl = precdata->mgl;
655
656 dvector *tempr = &(precdata->r);
657
658 const INT N0 = A_diag[0].row;
659 const INT N1 = A_diag[1].row;
660 const INT N2 = A_diag[2].row;
661 const INT N = N0 + N1 + N2;
662
663 INT i;
664
665 // back up r, setup z;
666 fasp_darray_cp(N, r, tempr->val);
667 fasp_darray_set(N, z, 0.0);
668
669 // prepare
670 dvector r0, r1, r2, z0, z1, z2;
671
672 r0.row = N0; z0.row = N0;
673 r1.row = N1; z1.row = N1;
674 r2.row = N2; z2.row = N2;
675
676 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
677 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
678
679 // Preconditioning A22 block
680 mgl[2]->b.row=N2; fasp_darray_cp(N2, r2.val, mgl[2]->b.val);
681 mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
682
683 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
684 fasp_darray_cp(N2, mgl[2]->x.val, z2.val);
685
686 // r1 = r1 - A5*z2
687 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
688
689 // Preconditioning A11 block
690 mgl[1]->b.row=N1; fasp_darray_cp(N1, r1.val, mgl[1]->b.val);
691 mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
692
693 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
694 fasp_darray_cp(N1, mgl[1]->x.val, z1.val);
695
696 // r0 = r0 - A1*z1 - A2*z2
697 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
698 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
699
700 // Preconditioning A00 block
701 mgl[0]->b.row=N0; fasp_darray_cp(N0, r0.val, mgl[0]->b.val);
702 mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
703
704 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
705 fasp_darray_cp(N0, mgl[0]->x.val, z0.val);
706
707 // restore r
708 fasp_darray_cp(N, tempr->val, r);
709}
710
725void fasp_precond_dblc_SGS_3 (REAL *r,
726 REAL *z,
727 void *data)
728{
729#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
730
731 precond_data_blc *precdata = (precond_data_blc *)data;
732 dBLCmat *A = precdata->Ablc;
733 dCSRmat *A_diag = precdata->A_diag;
734 dvector *tempr = &(precdata->r);
735
736#if WITH_UMFPACK
737 void **LU_diag = precdata->LU_diag;
738#endif
739
740 const INT N0 = A_diag[0].row;
741 const INT N1 = A_diag[1].row;
742 const INT N2 = A_diag[2].row;
743 const INT N = N0 + N1 + N2;
744
745 // back up r, setup z;
746 fasp_darray_cp(N, r, tempr->val);
747 fasp_darray_set(N, z, 0.0);
748
749 // prepare
750 dvector r0, r1, r2, z0, z1, z2;
751
752 r0.row = N0; z0.row = N0;
753 r1.row = N1; z1.row = N1;
754 r2.row = N2; z2.row = N2;
755
756 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
757 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
758
759 // Preconditioning A00 block
760#if WITH_UMFPACK
761 /* use UMFPACK direct solver */
762 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
763#elif WITH_SuperLU
764 /* use SuperLU direct solver */
765 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
766#endif
767
768 // r1 = r1 - A3*z0
769 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
770
771 // Preconditioning A11 block
772#if WITH_UMFPACK
773 /* use UMFPACK direct solver */
774 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
775#elif WITH_SuperLU
776 /* use SuperLU direct solver */
777 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
778#endif
779
780 // r2 = r2 - A6*z0 - A7*z1
781 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
782 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
783
784 // Preconditioning A22 block
785#if WITH_UMFPACK
786 /* use UMFPACK direct solver */
787 fasp_umfpack_solve(&A_diag[2], &r2, &z2, LU_diag[2], 0);
788#elif WITH_SuperLU
789 /* use SuperLU direct solver */
790 fasp_solver_superlu(&A_diag[2], &r2, &z2, 0);
791#endif
792
793 // r1 = r1 - A5*z2
794 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
795
796 // Preconditioning A11 block
797#if WITH_UMFPACK
798 /* use UMFPACK direct solver */
799 fasp_umfpack_solve(&A_diag[1], &r1, &z1, LU_diag[1], 0);
800#elif WITH_SuperLU
801 /* use SuperLU direct solver */
802 fasp_solver_superlu(&A_diag[1], &r1, &z1, 0);
803#endif
804
805 // r0 = r0 - A1*z1 - A2*z2
806 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
807 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
808
809 // Preconditioning A00 block
810#if WITH_UMFPACK
811 /* use UMFPACK direct solver */
812 fasp_umfpack_solve(&A_diag[0], &r0, &z0, LU_diag[0], 0);
813#elif WITH_SuperLU
814 /* use SuperLU direct solver */
815 fasp_solver_superlu(&A_diag[0], &r0, &z0, 0);
816#endif
817
818 // restore r
819 fasp_darray_cp(N, tempr->val, r);
820
821#endif
822}
823
838void fasp_precond_dblc_SGS_3_amg (REAL *r,
839 REAL *z,
840 void *data)
841{
842 precond_data_blc *precdata = (precond_data_blc *)data;
843 dBLCmat *A = precdata->Ablc;
844 dCSRmat *A_diag = precdata->A_diag;
845
846 AMG_param *amgparam = precdata->amgparam;
847 AMG_data **mgl = precdata->mgl;
848
849 INT i;
850
851 dvector *tempr = &(precdata->r);
852
853 const INT N0 = A_diag[0].row;
854 const INT N1 = A_diag[1].row;
855 const INT N2 = A_diag[2].row;
856 const INT N = N0 + N1 + N2;
857
858 // back up r, setup z;
859 fasp_darray_cp(N, r, tempr->val);
860 fasp_darray_set(N, z, 0.0);
861
862 // prepare
863 dvector r0, r1, r2, z0, z1, z2;
864
865 r0.row = N0; z0.row = N0;
866 r1.row = N1; z1.row = N1;
867 r2.row = N2; z2.row = N2;
868
869 r0.val = r; r1.val = &(r[N0]); r2.val = &(r[N0+N1]);
870 z0.val = z; z1.val = &(z[N0]); z2.val = &(z[N0+N1]);
871
872 // Preconditioning A00 block
873 mgl[0]->b.row=N0; fasp_darray_cp(N0, r0.val, mgl[0]->b.val);
874 mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
875
876 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
877 fasp_darray_cp(N0, mgl[0]->x.val, z0.val);
878
879 // r1 = r1 - A3*z0
880 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[3], z0.val, r1.val);
881
882 // Preconditioning A11 block
883 mgl[1]->b.row=N1; fasp_darray_cp(N1, r1.val, mgl[1]->b.val);
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_darray_cp(N1, mgl[1]->x.val, z1.val);
888
889 // r2 = r2 - A6*z0 - A7*z1
890 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[6], z0.val, r2.val);
891 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[7], z1.val, r2.val);
892
893 // Preconditioning A22 block
894 mgl[2]->b.row=N2; fasp_darray_cp(N2, r2.val, mgl[2]->b.val);
895 mgl[2]->x.row=N2; fasp_dvec_set(N2, &mgl[2]->x,0.0);
896
897 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[2], amgparam);
898 fasp_darray_cp(N2, mgl[2]->x.val, z2.val);
899
900 // r1 = r1 - A5*z2
901 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[5], z2.val, r1.val);
902
903 // Preconditioning A11 block
904 mgl[1]->b.row=N1; fasp_darray_cp(N1, r1.val, mgl[1]->b.val);
905 mgl[1]->x.row=N1; fasp_dvec_set(N1, &mgl[1]->x,0.0);
906
907 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[1], amgparam);
908 fasp_darray_cp(N1, mgl[1]->x.val, z1.val);
909
910 // r0 = r0 - A1*z1 - A2*z2
911 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[1], z1.val, r0.val);
912 fasp_blas_dcsr_aAxpy(-1.0, A->blocks[2], z2.val, r0.val);
913
914 // Preconditioning A00 block
915 mgl[0]->b.row=N0; fasp_darray_cp(N0, r0.val, mgl[0]->b.val);
916 mgl[0]->x.row=N0; fasp_dvec_set(N0, &mgl[0]->x, 0.0);
917
918 for(i=0;i<1;++i) fasp_solver_mgcycle(mgl[0], amgparam);
919 fasp_darray_cp(N0, mgl[0]->x.val, z0.val);
920
921 // restore r
922 fasp_darray_cp(N, tempr->val, r);
923}
924
939void fasp_precond_dblc_sweeping (REAL *r,
940 REAL *z,
941 void *data)
942{
943#if WITH_UMFPACK || WITH_SuperLU // Must use direct solvers for this method!
944
945 precond_data_sweeping *precdata = (precond_data_sweeping *)data;
946
947 INT NumLayers = precdata->NumLayers;
948 dBLCmat *A = precdata->A;
949 dBLCmat *Ai = precdata->Ai;
950 dCSRmat *local_A = precdata->local_A;
951 ivector *local_index = precdata->local_index;
952
953#if WITH_UMFPACK
954 void **local_LU = precdata->local_LU;
955#endif
956
957 dvector *r_backup = &(precdata->r);
958 REAL *w = precdata->w;
959
960 // local veriables
961 INT i,l;
962 dvector temp_r;
963 dvector temp_e;
964
965 dvector *local_r = (dvector *)fasp_mem_calloc(NumLayers, sizeof(dvector));
966 dvector *local_e = (dvector *)fasp_mem_calloc(NumLayers, sizeof(dvector));
967
968 // calculate the size and generate block local_r and local_z
969 INT N=0;
970
971 for (l=0;l<NumLayers; l++) {
972
973 local_r[l].row = A->blocks[l*NumLayers+l]->row;
974 local_r[l].val = r+N;
975
976 local_e[l].row = A->blocks[l*NumLayers+l]->col;
977 local_e[l].val = z+N;
978
979 N = N+A->blocks[l*NumLayers+l]->col;
980
981 }
982
983 temp_r.val = w;
984 temp_e.val = w+N;
985
986 // back up r, setup z;
987 fasp_darray_cp(N, r, r_backup->val);
988 fasp_darray_cp(N, r, z);
989
990 // L^{-1}r
991 for (l=0; l<NumLayers-1; l++){
992
993 temp_r.row = local_A[l].row;
994 temp_e.row = local_A[l].row;
995
996 fasp_dvec_set(local_A[l].row, &temp_r, 0.0);
997
998 for (i=0; i<local_e[l].row; i++){
999 temp_r.val[local_index[l].val[i]] = local_e[l].val[i];
1000 }
1001
1002#if WITH_UMFPACK
1003 /* use UMFPACK direct solver */
1004 fasp_umfpack_solve(&local_A[l], &temp_r, &temp_e, local_LU[l], 0);
1005#elif WITH_SuperLU
1006 /* use SuperLU direct solver */
1007 fasp_solver_superlu(&local_A[l], &temp_r, &temp_e, 0);
1008#endif
1009
1010 for (i=0; i<local_r[l].row; i++){
1011 local_r[l].val[i] = temp_e.val[local_index[l].val[i]];
1012 }
1013
1014 fasp_blas_dcsr_aAxpy(-1.0, Ai->blocks[(l+1)*NumLayers+l], local_r[l].val,
1015 local_e[l+1].val);
1016
1017 }
1018
1019 // D^{-1}L^{-1}r
1020 for (l=0; l<NumLayers; l++){
1021
1022 temp_r.row = local_A[l].row;
1023 temp_e.row = local_A[l].row;
1024
1025 fasp_dvec_set(local_A[l].row, &temp_r, 0.0);
1026
1027 for (i=0; i<local_e[l].row; i++){
1028 temp_r.val[local_index[l].val[i]] = local_e[l].val[i];
1029 }
1030
1031#if WITH_UMFPACK
1032 /* use UMFPACK direct solver */
1033 fasp_umfpack_solve(&local_A[l], &temp_r, &temp_e, local_LU[l], 0);
1034#elif WITH_SuperLU
1035 /* use SuperLU direct solver */
1036 fasp_solver_superlu(&local_A[l], &temp_r, &temp_e, 0);
1037#endif
1038
1039 for (i=0; i<local_e[l].row; i++) {
1040 local_e[l].val[i] = temp_e.val[local_index[l].val[i]];
1041 }
1042
1043 }
1044
1045 // L^{-t}D^{-1}L^{-1}u
1046 for (l=NumLayers-2; l>=0; l--){
1047
1048 temp_r.row = local_A[l].row;
1049 temp_e.row = local_A[l].row;
1050
1051 fasp_dvec_set(local_A[l].row, &temp_r, 0.0);
1052
1053 fasp_blas_dcsr_mxv (Ai->blocks[l*NumLayers+l+1], local_e[l+1].val, local_r[l].val);
1054
1055 for (i=0; i<local_r[l].row; i++){
1056 temp_r.val[local_index[l].val[i]] = local_r[l].val[i];
1057 }
1058
1059#if WITH_UMFPACK
1060 /* use UMFPACK direct solver */
1061 fasp_umfpack_solve(&local_A[l], &temp_r, &temp_e, local_LU[l], 0);
1062#elif WITH_SuperLU
1063 /* use SuperLU direct solver */
1064 fasp_solver_superlu(&local_A[l], &temp_r, &temp_e, 0);
1065#endif
1066
1067 for (i=0; i<local_e[l].row; i++){
1068 local_e[l].val[i] = local_e[l].val[i] - temp_e.val[local_index[l].val[i]];
1069 }
1070
1071 }
1072
1073 // restore r
1074 fasp_darray_cp(N, r_backup->val, r);
1075
1076#endif
1077}
1078
1079/*---------------------------------*/
1080/*-- End of File --*/
1081/*---------------------------------*/
fasp_darray_set
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:41
fasp_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
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_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_blas_dcsr_mxv
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:242
fasp_blas_dcsr_aAxpy
void fasp_blas_dcsr_aAxpy(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvCSR.c:494
fasp_precond_dblc_upper_3
void fasp_precond_dblc_upper_3(REAL *r, REAL *z, void *data)
block upper triangular preconditioning (3x3 blocks)
Definition: PreBLC.c:557
fasp_precond_dblc_SGS_3
void fasp_precond_dblc_SGS_3(REAL *r, REAL *z, void *data)
Block symmetric GS preconditioning (3x3 blocks)
Definition: PreBLC.c:725
fasp_precond_dblc_diag_4
void fasp_precond_dblc_diag_4(REAL *r, REAL *z, void *data)
Block diagonal preconditioning (4x4 blocks)
Definition: PreBLC.c:191
fasp_precond_dblc_lower_3
void fasp_precond_dblc_lower_3(REAL *r, REAL *z, void *data)
block lower triangular preconditioning (3x3 blocks)
Definition: PreBLC.c:291
fasp_precond_dblc_SGS_3_amg
void fasp_precond_dblc_SGS_3_amg(REAL *r, REAL *z, void *data)
Block symmetric GS preconditioning (3x3 blocks)
Definition: PreBLC.c:838
fasp_precond_dblc_diag_3_amg
void fasp_precond_dblc_diag_3_amg(REAL *r, REAL *z, void *data)
Block diagonal preconditioning (3x3 blocks)
Definition: PreBLC.c:126
fasp_precond_dblc_lower_4
void fasp_precond_dblc_lower_4(REAL *r, REAL *z, void *data)
block lower triangular preconditioning (4x4 blocks)
Definition: PreBLC.c:453
fasp_precond_dblc_upper_3_amg
void fasp_precond_dblc_upper_3_amg(REAL *r, REAL *z, void *data)
block upper triangular preconditioning (3x3 blocks)
Definition: PreBLC.c:645
fasp_precond_dblc_lower_3_amg
void fasp_precond_dblc_lower_3_amg(REAL *r, REAL *z, void *data)
block lower triangular preconditioning (3x3 blocks)
Definition: PreBLC.c:379
fasp_precond_dblc_diag_3
void fasp_precond_dblc_diag_3(REAL *r, REAL *z, void *data)
Block diagonal preconditioner (3x3 blocks)
Definition: PreBLC.c:38
fasp_precond_dblc_sweeping
void fasp_precond_dblc_sweeping(REAL *r, REAL *z, void *data)
Sweeping preconditioner for Maxwell equations.
Definition: PreBLC.c:939
fasp_solver_mgcycle
void fasp_solver_mgcycle(AMG_data *mgl, AMG_param *param)
Solve Ax=b with non-recursive multigrid cycle.
Definition: PreMGCycle.c:48
fasp_solver_superlu
int fasp_solver_superlu(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by SuperLU.
Definition: XtrSuperlu.c:47
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
INT
#define INT
Definition: fasp.h:72
fasp_block.h
Header file for FASP block matrices.
AMG_data
Data for AMG methods.
Definition: fasp.h:804
AMG_data::b
dvector b
pointer to the right-hand side at level level_num
Definition: fasp.h:826
AMG_data::x
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:829
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
dBLCmat
Block REAL CSR matrix format.
Definition: fasp_block.h:74
dBLCmat::blocks
dCSRmat ** blocks
blocks of dCSRmat, point to blocks[brow][bcol]
Definition: fasp_block.h:83
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
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::val
INT * val
actual vector entries
Definition: fasp.h:374
precond_data_blc
Data for block preconditioners in dBLCmat format.
Definition: fasp_block.h:363
precond_data_blc::Ablc
dBLCmat * Ablc
Definition: fasp_block.h:368
precond_data_blc::A_diag
dCSRmat * A_diag
Definition: fasp_block.h:370
precond_data_blc::mgl
AMG_data ** mgl
Definition: fasp_block.h:382
precond_data_blc::LU_diag
void ** LU_diag
Definition: fasp_block.h:379
precond_data_blc::amgparam
AMG_param * amgparam
Definition: fasp_block.h:384
precond_data_blc::r
dvector r
Definition: fasp_block.h:372
precond_data_sweeping
Data for sweeping preconditioner.
Definition: fasp_block.h:398
precond_data_sweeping::local_index
ivector * local_index
Definition: fasp_block.h:408
precond_data_sweeping::Ai
dBLCmat * Ai
Definition: fasp_block.h:403
precond_data_sweeping::w
REAL * w
Definition: fasp_block.h:412
precond_data_sweeping::A
dBLCmat * A
Definition: fasp_block.h:402
precond_data_sweeping::NumLayers
INT NumLayers
Definition: fasp_block.h:400
precond_data_sweeping::local_LU
void ** local_LU
Definition: fasp_block.h:406
precond_data_sweeping::local_A
dCSRmat * local_A
Definition: fasp_block.h:405
precond_data_sweeping::r
dvector r
Definition: fasp_block.h:411