Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreBSR.c
Go to the documentation of this file.
1
16#ifdef _OPENMP
17#include <omp.h>
18#endif
19
20#include "fasp.h"
21#include "fasp_functs.h"
22
23/*---------------------------------*/
24/*-- Declare Private Functions --*/
25/*---------------------------------*/
26
27#include "PreMGUtil.inl"
28
29/*---------------------------------*/
30/*-- Public Functions --*/
31/*---------------------------------*/
32
49void fasp_precond_dbsr_diag(REAL* r, REAL* z, void* data)
50{
51 precond_diag_bsr* diag = (precond_diag_bsr*)data;
52 const INT nb = diag->nb;
53
54 switch (nb) {
55
56 case 2:
57 fasp_precond_dbsr_diag_nc2(r, z, diag);
58 break;
59 case 3:
60 fasp_precond_dbsr_diag_nc3(r, z, diag);
61 break;
62
63 case 4:
64 fasp_precond_dbsr_diag_nc4(r, z, diag);
65 break;
66
67 case 5:
68 fasp_precond_dbsr_diag_nc5(r, z, diag);
69 break;
70
71 case 7:
72 fasp_precond_dbsr_diag_nc7(r, z, diag);
73 break;
74
75 default:
76 {
77 REAL* diagptr = diag->diag.val;
78 const INT nb2 = nb * nb;
79 const INT m = diag->diag.row / nb2;
80 INT i;
81
82#ifdef _OPENMP
83 if (m > OPENMP_HOLDS) {
84 INT myid, mybegin, myend;
85 INT nthreads = fasp_get_num_threads();
86#pragma omp parallel for private(myid, mybegin, myend, i)
87 for (myid = 0; myid < nthreads; myid++) {
88 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
89 for (i = mybegin; i < myend; ++i) {
90 fasp_blas_smat_mxv(&(diagptr[i * nb2]), &(r[i * nb]),
91 &(z[i * nb]), nb);
92 }
93 }
94 } else {
95#endif
96 for (i = 0; i < m; ++i) {
97 fasp_blas_smat_mxv(&(diagptr[i * nb2]), &(r[i * nb]),
98 &(z[i * nb]), nb);
99 }
100#ifdef _OPENMP
101 }
102#endif
103 break;
104 }
105 }
106}
107
124void fasp_precond_dbsr_diag_nc2(REAL* r, REAL* z, void* data)
125{
126 precond_diag_bsr* diag = (precond_diag_bsr*)data;
127 REAL* diagptr = diag->diag.val;
128
129 INT i;
130 const INT m = diag->diag.row / 4;
131
132#ifdef _OPENMP
133 if (m > OPENMP_HOLDS) {
134 INT myid, mybegin, myend;
135 INT nthreads = fasp_get_num_threads();
136#pragma omp parallel for private(myid, mybegin, myend, i)
137 for (myid = 0; myid < nthreads; myid++) {
138 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
139 for (i = mybegin; i < myend; ++i) {
140 fasp_blas_smat_mxv_nc2(&(diagptr[i * 4]), &(r[i * 2]), &(z[i * 2]));
141 }
142 }
143 } else {
144#endif
145 for (i = 0; i < m; ++i) {
146 fasp_blas_smat_mxv_nc2(&(diagptr[i * 4]), &(r[i * 2]), &(z[i * 2]));
147 }
148#ifdef _OPENMP
149 }
150#endif
151}
152
169void fasp_precond_dbsr_diag_nc3(REAL* r, REAL* z, void* data)
170{
171 precond_diag_bsr* diag = (precond_diag_bsr*)data;
172 REAL* diagptr = diag->diag.val;
173
174 const INT m = diag->diag.row / 9;
175 INT i;
176
177#ifdef _OPENMP
178 if (m > OPENMP_HOLDS) {
179 INT myid, mybegin, myend;
180 INT nthreads = fasp_get_num_threads();
181#pragma omp parallel for private(myid, mybegin, myend, i)
182 for (myid = 0; myid < nthreads; myid++) {
183 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
184 for (i = mybegin; i < myend; ++i) {
185 fasp_blas_smat_mxv_nc3(&(diagptr[i * 9]), &(r[i * 3]), &(z[i * 3]));
186 }
187 }
188 } else {
189#endif
190 for (i = 0; i < m; ++i) {
191 fasp_blas_smat_mxv_nc3(&(diagptr[i * 9]), &(r[i * 3]), &(z[i * 3]));
192 }
193#ifdef _OPENMP
194 }
195#endif
196}
197
214void fasp_precond_dbsr_diag_nc4(REAL* r, REAL* z, void* data)
215{
216 precond_diag_bsr* diag = (precond_diag_bsr*)data;
217 REAL* diagptr = diag->diag.val;
218
219 const INT m = diag->diag.row / 16;
220 INT i;
221
222#ifdef _OPENMP
223 if (m > OPENMP_HOLDS) {
224 INT myid, mybegin, myend;
225 INT nthreads = fasp_get_num_threads();
226#pragma omp parallel for private(myid, mybegin, myend, i)
227 for (myid = 0; myid < nthreads; myid++) {
228 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
229 for (i = mybegin; i < myend; ++i) {
230 fasp_blas_smat_mxv_nc4(&(diagptr[i * 16]), &(r[i * 4]), &(z[i * 4]));
231 }
232 }
233 } else {
234#endif
235 for (i = 0; i < m; ++i) {
236 fasp_blas_smat_mxv_nc4(&(diagptr[i * 16]), &(r[i * 4]), &(z[i * 4]));
237 }
238#ifdef _OPENMP
239 }
240#endif
241}
242
259void fasp_precond_dbsr_diag_nc5(REAL* r, REAL* z, void* data)
260{
261 precond_diag_bsr* diag = (precond_diag_bsr*)data;
262 REAL* diagptr = diag->diag.val;
263
264 const INT m = diag->diag.row / 25;
265 INT i;
266
267#ifdef _OPENMP
268 if (m > OPENMP_HOLDS) {
269 INT myid, mybegin, myend;
270 INT nthreads = fasp_get_num_threads();
271#pragma omp parallel for private(myid, mybegin, myend, i)
272 for (myid = 0; myid < nthreads; myid++) {
273 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
274 for (i = mybegin; i < myend; ++i) {
275 fasp_blas_smat_mxv_nc5(&(diagptr[i * 25]), &(r[i * 5]), &(z[i * 5]));
276 }
277 }
278 } else {
279#endif
280 for (i = 0; i < m; ++i) {
281 fasp_blas_smat_mxv_nc5(&(diagptr[i * 25]), &(r[i * 5]), &(z[i * 5]));
282 }
283#ifdef _OPENMP
284 }
285#endif
286}
287
304void fasp_precond_dbsr_diag_nc7(REAL* r, REAL* z, void* data)
305{
306 precond_diag_bsr* diag = (precond_diag_bsr*)data;
307 REAL* diagptr = diag->diag.val;
308
309 const INT m = diag->diag.row / 49;
310 INT i;
311
312#ifdef _OPENMP
313 if (m > OPENMP_HOLDS) {
314 INT myid, mybegin, myend;
315 INT nthreads = fasp_get_num_threads();
316#pragma omp parallel for private(myid, mybegin, myend, i)
317 for (myid = 0; myid < nthreads; myid++) {
318 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
319 for (i = mybegin; i < myend; ++i) {
320 fasp_blas_smat_mxv_nc7(&(diagptr[i * 49]), &(r[i * 7]), &(z[i * 7]));
321 }
322 }
323 } else {
324#endif
325 for (i = 0; i < m; ++i) {
326 fasp_blas_smat_mxv_nc7(&(diagptr[i * 49]), &(r[i * 7]), &(z[i * 7]));
327 }
328#ifdef _OPENMP
329 }
330#endif
331}
332
347void fasp_precond_dbsr_ilu(REAL* r, REAL* z, void* data)
348{
349 const ILU_data* iludata = (ILU_data*)data;
350 const INT m = iludata->row, mm1 = m - 1, mm2 = m - 2, memneed = 2 * m;
351 const INT nb = iludata->nb, nb2 = nb * nb, size = m * nb;
352
353 INT* ijlu = iludata->ijlu;
354 REAL* lu = iludata->luval;
355
356 INT ib, ibstart, ibstart1;
357 INT i, j, jj, begin_row, end_row;
358 REAL *zz, *zr, *mult;
359
360 if (iludata->nwork < memneed) {
361 printf("### ERROR: Need %d memory, only %d available!\n", memneed,
362 iludata->nwork);
363 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
364 }
365
366 zz = iludata->work;
367 zr = zz + size;
368 mult = zr + size;
369
370 memcpy(zr, r, size * sizeof(REAL));
371
372 switch (nb) {
373
374 case 1:
375
376 // forward sweep: solve unit lower matrix equation L*zz=zr
377 zz[0] = zr[0];
378 for (i = 1; i <= mm1; ++i) {
379 begin_row = ijlu[i];
380 end_row = ijlu[i + 1] - 1;
381 for (j = begin_row; j <= end_row; ++j) {
382 jj = ijlu[j];
383 if (jj < i)
384 zr[i] -= lu[j] * zz[jj];
385 else
386 break;
387 }
388 zz[i] = zr[i];
389 }
390
391 // backward sweep: solve upper matrix equation U*z=zz
392 z[mm1] = zz[mm1] * lu[mm1];
393 for (i = mm2; i >= 0; i--) {
394 begin_row = ijlu[i];
395 end_row = ijlu[i + 1] - 1;
396 for (j = end_row; j >= begin_row; j--) {
397 jj = ijlu[j];
398 if (jj > i)
399 zz[i] -= lu[j] * z[jj];
400 else
401 break;
402 }
403 z[i] = zz[i] * lu[i];
404 }
405
406 break; // end (if nb==1)
407
408 case 3:
409
410 // forward sweep: solve unit lower matrix equation L*zz=zr
411 zz[0] = zr[0];
412 zz[1] = zr[1];
413 zz[2] = zr[2];
414
415 for (i = 1; i <= mm1; ++i) {
416 begin_row = ijlu[i];
417 end_row = ijlu[i + 1] - 1;
418 ibstart = i * nb;
419 for (j = begin_row; j <= end_row; ++j) {
420 jj = ijlu[j];
421 if (jj < i) {
422 fasp_blas_smat_mxv_nc3(&(lu[j * nb2]), &(zz[jj * nb]), mult);
423 for (ib = 0; ib < nb; ++ib) zr[ibstart + ib] -= mult[ib];
424 } else
425 break;
426 }
427
428 zz[ibstart] = zr[ibstart];
429 zz[ibstart + 1] = zr[ibstart + 1];
430 zz[ibstart + 2] = zr[ibstart + 2];
431 }
432
433 // backward sweep: solve upper matrix equation U*z=zz
434 ibstart = mm1 * nb2;
435 ibstart1 = mm1 * nb;
436 fasp_blas_smat_mxv_nc3(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]));
437
438 for (i = mm2; i >= 0; i--) {
439 begin_row = ijlu[i];
440 end_row = ijlu[i + 1] - 1;
441 ibstart = i * nb2;
442 ibstart1 = i * nb;
443 for (j = end_row; j >= begin_row; j--) {
444 jj = ijlu[j];
445 if (jj > i) {
446 fasp_blas_smat_mxv_nc3(&(lu[j * nb2]), &(z[jj * nb]), mult);
447 for (ib = 0; ib < nb; ++ib) zz[ibstart1 + ib] -= mult[ib];
448 }
449
450 else
451 break;
452 }
453
454 fasp_blas_smat_mxv_nc3(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]));
455 }
456
457 break; // end (if nb=3)
458
459 case 5:
460
461 // forward sweep: solve unit lower matrix equation L*zz=zr
462 fasp_darray_cp(nb, &(zr[0]), &(zz[0]));
463
464 for (i = 1; i <= mm1; ++i) {
465 begin_row = ijlu[i];
466 end_row = ijlu[i + 1] - 1;
467 ibstart = i * nb;
468 for (j = begin_row; j <= end_row; ++j) {
469 jj = ijlu[j];
470 if (jj < i) {
471 fasp_blas_smat_mxv_nc5(&(lu[j * nb2]), &(zz[jj * nb]), mult);
472 for (ib = 0; ib < nb; ++ib) zr[ibstart + ib] -= mult[ib];
473 } else
474 break;
475 }
476
477 fasp_darray_cp(nb, &(zr[ibstart]), &(zz[ibstart]));
478 }
479
480 // backward sweep: solve upper matrix equation U*z=zz
481 ibstart = mm1 * nb2;
482 ibstart1 = mm1 * nb;
483 fasp_blas_smat_mxv_nc5(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]));
484
485 for (i = mm2; i >= 0; i--) {
486 begin_row = ijlu[i];
487 end_row = ijlu[i + 1] - 1;
488 ibstart = i * nb2;
489 ibstart1 = i * nb;
490 for (j = end_row; j >= begin_row; j--) {
491 jj = ijlu[j];
492 if (jj > i) {
493 fasp_blas_smat_mxv_nc5(&(lu[j * nb2]), &(z[jj * nb]), mult);
494 for (ib = 0; ib < nb; ++ib) zz[ibstart1 + ib] -= mult[ib];
495 }
496
497 else
498 break;
499 }
500
501 fasp_blas_smat_mxv_nc5(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]));
502 }
503
504 break; // end (if nb==5)
505
506 case 7:
507
508 // forward sweep: solve unit lower matrix equation L*zz=zr
509 fasp_darray_cp(nb, &(zr[0]), &(zz[0]));
510
511 for (i = 1; i <= mm1; ++i) {
512 begin_row = ijlu[i];
513 end_row = ijlu[i + 1] - 1;
514 ibstart = i * nb;
515 for (j = begin_row; j <= end_row; ++j) {
516 jj = ijlu[j];
517 if (jj < i) {
518 fasp_blas_smat_mxv_nc7(&(lu[j * nb2]), &(zz[jj * nb]), mult);
519 for (ib = 0; ib < nb; ++ib) zr[ibstart + ib] -= mult[ib];
520 } else
521 break;
522 }
523
524 fasp_darray_cp(nb, &(zr[ibstart]), &(zz[ibstart]));
525 }
526
527 // backward sweep: solve upper matrix equation U*z=zz
528 ibstart = mm1 * nb2;
529 ibstart1 = mm1 * nb;
530 fasp_blas_smat_mxv_nc7(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]));
531
532 for (i = mm2; i >= 0; i--) {
533 begin_row = ijlu[i];
534 end_row = ijlu[i + 1] - 1;
535 ibstart = i * nb2;
536 ibstart1 = i * nb;
537 for (j = end_row; j >= begin_row; j--) {
538 jj = ijlu[j];
539 if (jj > i) {
540 fasp_blas_smat_mxv_nc7(&(lu[j * nb2]), &(z[jj * nb]), mult);
541 for (ib = 0; ib < nb; ++ib) zz[ibstart1 + ib] -= mult[ib];
542 }
543
544 else
545 break;
546 }
547
548 fasp_blas_smat_mxv_nc7(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]));
549 }
550
551 break; // end (if nb==7)
552
553 default:
554
555 // forward sweep: solve unit lower matrix equation L*zz=zr
556 fasp_darray_cp(nb, &(zr[0]), &(zz[0]));
557
558 for (i = 1; i <= mm1; ++i) {
559 begin_row = ijlu[i];
560 end_row = ijlu[i + 1] - 1;
561 ibstart = i * nb;
562 for (j = begin_row; j <= end_row; ++j) {
563 jj = ijlu[j];
564 if (jj < i) {
565 fasp_blas_smat_mxv(&(lu[j * nb2]), &(zz[jj * nb]), mult, nb);
566 for (ib = 0; ib < nb; ++ib) zr[ibstart + ib] -= mult[ib];
567 } else
568 break;
569 }
570
571 fasp_darray_cp(nb, &(zr[ibstart]), &(zz[ibstart]));
572 }
573
574 // backward sweep: solve upper matrix equation U*z=zz
575 ibstart = mm1 * nb2;
576 ibstart1 = mm1 * nb;
577 fasp_blas_smat_mxv(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]), nb);
578
579 for (i = mm2; i >= 0; i--) {
580 begin_row = ijlu[i];
581 end_row = ijlu[i + 1] - 1;
582 ibstart = i * nb2;
583 ibstart1 = i * nb;
584 for (j = end_row; j >= begin_row; j--) {
585 jj = ijlu[j];
586 if (jj > i) {
587 fasp_blas_smat_mxv(&(lu[j * nb2]), &(z[jj * nb]), mult, nb);
588 for (ib = 0; ib < nb; ++ib) zz[ibstart1 + ib] -= mult[ib];
589 }
590
591 else
592 break;
593 }
594
595 fasp_blas_smat_mxv(&(lu[ibstart]), &(zz[ibstart1]), &(z[ibstart1]), nb);
596 }
597
598 break; // end everything else
599 }
600
601 return;
602}
603
618void fasp_precond_dbsr_ilu_mc_omp(REAL* r, REAL* z, void* data)
619{
620#ifdef _OPENMP
621 const ILU_data* iludata = (ILU_data*)data;
622 const INT m = iludata->row, memneed = 2 * m;
623 const INT nb = iludata->nb, nb2 = nb * nb, size = m * nb;
624
625 INT* ijlu = iludata->ijlu;
626 REAL* lu = iludata->luval;
627 INT ncolors = iludata->nlevL;
628 INT* ic = iludata->ilevL;
629
630 INT ib, ibstart, ibstart1;
631 INT i, j, jj, k, begin_row, end_row;
632 REAL *zz, *zr, *mult;
633
634 if (iludata->nwork < memneed) {
635 printf("### ERROR: Need %d memory, only %d available!\n", memneed,
636 iludata->nwork);
637 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
638 }
639
640 zz = iludata->work;
641 zr = zz + size;
642
643 memcpy(zr, r, size * sizeof(REAL));
644
645 switch (nb) {
646
647 case 1:
648 // forward sweep: solve unit lower matrix equation L*zz=zr
649 for (k = 0; k < ncolors; ++k) {
650#pragma omp parallel for private(i, begin_row, end_row, j, jj)
651 for (i = ic[k]; i < ic[k + 1]; ++i) {
652 begin_row = ijlu[i];
653 end_row = ijlu[i + 1] - 1;
654 for (j = begin_row; j <= end_row; ++j) {
655 jj = ijlu[j];
656 if (jj < i)
657 zr[i] -= lu[j] * zz[jj];
658 else
659 break;
660 }
661 zz[i] = zr[i];
662 }
663 }
664 // backward sweep: solve upper matrix equation U*z=zz
665 for (k = ncolors - 1; k >= 0; k--) {
666#pragma omp parallel for private(i, begin_row, end_row, j, jj)
667 for (i = ic[k + 1] - 1; i >= ic[k]; i--) {
668 begin_row = ijlu[i];
669 end_row = ijlu[i + 1] - 1;
670 for (j = end_row; j >= begin_row; j--) {
671 jj = ijlu[j];
672 if (jj > i)
673 zz[i] -= lu[j] * z[jj];
674 else
675 break;
676 }
677 z[i] = zz[i] * lu[i];
678 }
679 }
680
681 break; // end (if nb==1)
682
683 case 2:
684
685 for (k = 0; k < ncolors; ++k) {
686#pragma omp parallel private(i, begin_row, end_row, ibstart, j, jj, ib, mult)
687 {
688 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
689#pragma omp for
690 for (i = ic[k]; i < ic[k + 1]; ++i) {
691 begin_row = ijlu[i];
692 end_row = ijlu[i + 1] - 1;
693 ibstart = i * nb;
694 for (j = begin_row; j <= end_row; ++j) {
695 jj = ijlu[j];
696 if (jj < i) {
697 fasp_blas_smat_mxv_nc2(&(lu[j * nb2]), &(zz[jj * nb]),
698 mult);
699 for (ib = 0; ib < nb; ++ib)
700 zr[ibstart + ib] -= mult[ib];
701 } else
702 break;
703 }
704
705 zz[ibstart] = zr[ibstart];
706 zz[ibstart + 1] = zr[ibstart + 1];
707 }
708
709 fasp_mem_free(mult);
710 mult = NULL;
711 }
712 }
713
714 for (k = ncolors - 1; k >= 0; k--) {
715#pragma omp parallel private(i, begin_row, end_row, ibstart, ibstart1, j, jj, ib, mult)
716 {
717 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
718#pragma omp for
719 for (i = ic[k + 1] - 1; i >= ic[k]; i--) {
720 begin_row = ijlu[i];
721 end_row = ijlu[i + 1] - 1;
722 ibstart = i * nb2;
723 ibstart1 = i * nb;
724 for (j = end_row; j >= begin_row; j--) {
725 jj = ijlu[j];
726 if (jj > i) {
727 fasp_blas_smat_mxv_nc2(&(lu[j * nb2]), &(z[jj * nb]),
728 mult);
729 for (ib = 0; ib < nb; ++ib)
730 zz[ibstart1 + ib] -= mult[ib];
731 }
732
733 else
734 break;
735 }
736
737 fasp_blas_smat_mxv_nc2(&(lu[ibstart]), &(zz[ibstart1]),
738 &(z[ibstart1]));
739 }
740
741 fasp_mem_free(mult);
742 mult = NULL;
743 }
744 }
745
746 break; // end (if nb=2)
747 case 3:
748
749 for (k = 0; k < ncolors; ++k) {
750#pragma omp parallel private(i, begin_row, end_row, ibstart, j, jj, ib, mult)
751 {
752 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
753#pragma omp for
754 for (i = ic[k]; i < ic[k + 1]; ++i) {
755 begin_row = ijlu[i];
756 end_row = ijlu[i + 1] - 1;
757 ibstart = i * nb;
758 for (j = begin_row; j <= end_row; ++j) {
759 jj = ijlu[j];
760 if (jj < i) {
761 fasp_blas_smat_mxv_nc3(&(lu[j * nb2]), &(zz[jj * nb]),
762 mult);
763 for (ib = 0; ib < nb; ++ib)
764 zr[ibstart + ib] -= mult[ib];
765 } else
766 break;
767 }
768
769 zz[ibstart] = zr[ibstart];
770 zz[ibstart + 1] = zr[ibstart + 1];
771 zz[ibstart + 2] = zr[ibstart + 2];
772 }
773
774 fasp_mem_free(mult);
775 mult = NULL;
776 }
777 }
778
779 for (k = ncolors - 1; k >= 0; k--) {
780#pragma omp parallel private(i, begin_row, end_row, ibstart, ibstart1, j, jj, ib, mult)
781 {
782 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
783#pragma omp for
784 for (i = ic[k + 1] - 1; i >= ic[k]; i--) {
785 begin_row = ijlu[i];
786 end_row = ijlu[i + 1] - 1;
787 ibstart = i * nb2;
788 ibstart1 = i * nb;
789 for (j = end_row; j >= begin_row; j--) {
790 jj = ijlu[j];
791 if (jj > i) {
792 fasp_blas_smat_mxv_nc3(&(lu[j * nb2]), &(z[jj * nb]),
793 mult);
794 for (ib = 0; ib < nb; ++ib)
795 zz[ibstart1 + ib] -= mult[ib];
796 }
797
798 else
799 break;
800 }
801
802 fasp_blas_smat_mxv_nc3(&(lu[ibstart]), &(zz[ibstart1]),
803 &(z[ibstart1]));
804 }
805
806 fasp_mem_free(mult);
807 mult = NULL;
808 }
809 }
810
811 break; // end (if nb=3)
812
813 default:
814 {
815 if (nb > 3) {
816 printf("### ERROR: Multi-thread Parallel ILU for %d components \
817 has not yet been implemented!!!",
818 nb);
819 fasp_chkerr(ERROR_UNKNOWN, __FUNCTION__);
820 }
821 break;
822 }
823 }
824
825 return;
826#endif
827}
828
844void fasp_precond_dbsr_ilu_ls_omp(REAL* r, REAL* z, void* data)
845{
846#ifdef _OPENMP
847 const ILU_data* iludata = (ILU_data*)data;
848 const INT m = iludata->row, memneed = 2 * m;
849 const INT nb = iludata->nb, nb2 = nb * nb, size = m * nb;
850
851 INT* ijlu = iludata->ijlu;
852 REAL* lu = iludata->luval;
853 INT nlevL = iludata->nlevL;
854 INT* ilevL = iludata->ilevL;
855 INT* jlevL = iludata->jlevL;
856 INT nlevU = iludata->nlevU;
857 INT* ilevU = iludata->ilevU;
858 INT* jlevU = iludata->jlevU;
859
860 INT ib, ibstart, ibstart1;
861 INT i, ii, j, jj, k, begin_row, end_row;
862 REAL *zz, *zr, *mult;
863
864 if (iludata->nwork < memneed) {
865 printf("### ERROR: Need %d memory, only %d available!\n", memneed,
866 iludata->nwork);
867 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
868 }
869
870 zz = iludata->work;
871 zr = zz + size;
872 // mult = zr + size;
873
874 memcpy(zr, r, size * sizeof(REAL));
875
876 switch (nb) {
877
878 case 1:
879 // forward sweep: solve unit lower matrix equation L*zz=zr
880 for (k = 0; k < nlevL; ++k) {
881#pragma omp parallel for private(i, ii, begin_row, end_row, j, jj)
882 for (ii = ilevL[k]; ii < ilevL[k + 1]; ++ii) {
883 i = jlevL[ii];
884 begin_row = ijlu[i];
885 end_row = ijlu[i + 1] - 1;
886 for (j = begin_row; j <= end_row; ++j) {
887 jj = ijlu[j];
888 if (jj < i)
889 zr[i] -= lu[j] * zz[jj];
890 else
891 break;
892 }
893 zz[i] = zr[i];
894 }
895 }
896 // backward sweep: solve upper matrix equation U*z=zz
897 for (k = 0; k < nlevU; k++) {
898#pragma omp parallel for private(i, ii, begin_row, end_row, j, jj)
899 for (ii = ilevU[k + 1] - 1; ii >= ilevU[k]; ii--) {
900 i = jlevU[ii];
901 begin_row = ijlu[i];
902 end_row = ijlu[i + 1] - 1;
903 for (j = end_row; j >= begin_row; j--) {
904 jj = ijlu[j];
905 if (jj > i)
906 zz[i] -= lu[j] * z[jj];
907 else
908 break;
909 }
910 z[i] = zz[i] * lu[i];
911 }
912 }
913
914 break; // end (if nb==1)
915
916 case 2:
917
918 for (k = 0; k < nlevL; ++k) {
919#pragma omp parallel private(i, ii, begin_row, end_row, ibstart, j, jj, ib, mult)
920 {
921 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
922#pragma omp for
923 for (ii = ilevL[k]; ii < ilevL[k + 1]; ++ii) {
924 i = jlevL[ii];
925 begin_row = ijlu[i];
926 end_row = ijlu[i + 1] - 1;
927 ibstart = i * nb;
928 for (j = begin_row; j <= end_row; ++j) {
929 jj = ijlu[j];
930 if (jj < i) {
931 fasp_blas_smat_mxv_nc2(&(lu[j * nb2]), &(zz[jj * nb]),
932 mult);
933 for (ib = 0; ib < nb; ++ib)
934 zr[ibstart + ib] -= mult[ib];
935 } else
936 break;
937 }
938
939 zz[ibstart] = zr[ibstart];
940 zz[ibstart + 1] = zr[ibstart + 1];
941 }
942
943 fasp_mem_free(mult);
944 mult = NULL;
945 }
946 }
947
948 for (k = 0; k < nlevU; k++) {
949#pragma omp parallel private(i, ii, begin_row, end_row, ibstart, ibstart1, j, jj, ib, \
950 mult)
951 {
952 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
953#pragma omp for
954 for (ii = ilevU[k + 1] - 1; ii >= ilevU[k]; ii--) {
955 i = jlevU[ii];
956 begin_row = ijlu[i];
957 end_row = ijlu[i + 1] - 1;
958 ibstart = i * nb2;
959 ibstart1 = i * nb;
960 for (j = end_row; j >= begin_row; j--) {
961 jj = ijlu[j];
962 if (jj > i) {
963 fasp_blas_smat_mxv_nc2(&(lu[j * nb2]), &(z[jj * nb]),
964 mult);
965 for (ib = 0; ib < nb; ++ib)
966 zz[ibstart1 + ib] -= mult[ib];
967 }
968
969 else
970 break;
971 }
972
973 fasp_blas_smat_mxv_nc2(&(lu[ibstart]), &(zz[ibstart1]),
974 &(z[ibstart1]));
975 }
976
977 fasp_mem_free(mult);
978 mult = NULL;
979 }
980 }
981
982 break; // end (if nb=2)
983 case 3:
984
985 for (k = 0; k < nlevL; ++k) {
986#pragma omp parallel private(i, ii, begin_row, end_row, ibstart, j, jj, ib, mult)
987 {
988 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
989#pragma omp for
990 for (ii = ilevL[k]; ii < ilevL[k + 1]; ++ii) {
991 i = jlevL[ii];
992 begin_row = ijlu[i];
993 end_row = ijlu[i + 1] - 1;
994 ibstart = i * nb;
995 for (j = begin_row; j <= end_row; ++j) {
996 jj = ijlu[j];
997 if (jj < i) {
998 fasp_blas_smat_mxv_nc3(&(lu[j * nb2]), &(zz[jj * nb]),
999 mult);
1000 for (ib = 0; ib < nb; ++ib)
1001 zr[ibstart + ib] -= mult[ib];
1002 } else
1003 break;
1004 }
1005
1006 zz[ibstart] = zr[ibstart];
1007 zz[ibstart + 1] = zr[ibstart + 1];
1008 zz[ibstart + 2] = zr[ibstart + 2];
1009 }
1010
1011 fasp_mem_free(mult);
1012 mult = NULL;
1013 }
1014 }
1015
1016 for (k = 0; k < nlevU; k++) {
1017#pragma omp parallel private(i, ii, begin_row, end_row, ibstart, ibstart1, j, jj, ib, \
1018 mult)
1019 {
1020 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
1021#pragma omp for
1022 for (ii = ilevU[k + 1] - 1; ii >= ilevU[k]; ii--) {
1023 i = jlevU[ii];
1024 begin_row = ijlu[i];
1025 end_row = ijlu[i + 1] - 1;
1026 ibstart = i * nb2;
1027 ibstart1 = i * nb;
1028 for (j = end_row; j >= begin_row; j--) {
1029 jj = ijlu[j];
1030 if (jj > i) {
1031 fasp_blas_smat_mxv_nc3(&(lu[j * nb2]), &(z[jj * nb]),
1032 mult);
1033 for (ib = 0; ib < nb; ++ib)
1034 zz[ibstart1 + ib] -= mult[ib];
1035 }
1036
1037 else
1038 break;
1039 }
1040
1041 fasp_blas_smat_mxv_nc3(&(lu[ibstart]), &(zz[ibstart1]),
1042 &(z[ibstart1]));
1043 }
1044
1045 fasp_mem_free(mult);
1046 mult = NULL;
1047 }
1048 }
1049
1050 break; // end (if nb=3)
1051
1052 default:
1053 {
1054
1055 for (k = 0; k < nlevL; ++k) {
1056#pragma omp parallel private(i, ii, begin_row, end_row, ibstart, j, jj, ib, mult)
1057 {
1058 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
1059#pragma omp for
1060 for (ii = ilevL[k]; ii < ilevL[k + 1]; ++ii) {
1061 i = jlevL[ii];
1062 begin_row = ijlu[i];
1063 end_row = ijlu[i + 1] - 1;
1064 ibstart = i * nb;
1065 for (j = begin_row; j <= end_row; ++j) {
1066 jj = ijlu[j];
1067 if (jj < i) {
1068 fasp_blas_smat_mxv(&(lu[j * nb2]), &(zz[jj * nb]),
1069 mult, nb);
1070 for (ib = 0; ib < nb; ++ib)
1071 zr[ibstart + ib] -= mult[ib];
1072 } else
1073 break;
1074 }
1075
1076 for (j = 0; j < nb; j++)
1077 zz[ibstart + j] =
1078 zr[ibstart + j]; // Li Zhao, 09/19/2022
1079 }
1080
1081 fasp_mem_free(mult);
1082 mult = NULL;
1083 }
1084 }
1085
1086 for (k = 0; k < nlevU; k++) {
1087#pragma omp parallel private(i, ii, begin_row, end_row, ibstart, ibstart1, j, jj, ib, \
1088 mult)
1089 {
1090 mult = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
1091#pragma omp for
1092 for (ii = ilevU[k + 1] - 1; ii >= ilevU[k]; ii--) {
1093 i = jlevU[ii];
1094 begin_row = ijlu[i];
1095 end_row = ijlu[i + 1] - 1;
1096 ibstart = i * nb2;
1097 ibstart1 = i * nb;
1098 for (j = end_row; j >= begin_row; j--) {
1099 jj = ijlu[j];
1100 if (jj > i) {
1101 fasp_blas_smat_mxv(&(lu[j * nb2]), &(z[jj * nb]),
1102 mult, nb);
1103 for (ib = 0; ib < nb; ++ib)
1104 zz[ibstart1 + ib] -= mult[ib];
1105 }
1106
1107 else
1108 break;
1109 }
1110
1111 fasp_blas_smat_mxv(&(lu[ibstart]), &(zz[ibstart1]),
1112 &(z[ibstart1]), nb);
1113 }
1114
1115 fasp_mem_free(mult);
1116 mult = NULL;
1117 }
1118 }
1119
1120 break;
1121
1122 /*
1123 if (nb > 3) {
1124 printf("### ERROR: Multi-thread Parallel ILU for %d components \
1125 has not yet been implemented!!!", nb);
1126 fasp_chkerr(ERROR_UNKNOWN, __FUNCTION__);
1127 }
1128 break;
1129 */
1130 }
1131 }
1132
1133 return;
1134#endif
1135}
1136
1149void fasp_precond_dbsr_amg(REAL* r, REAL* z, void* data)
1150{
1151 precond_data_bsr* predata = (precond_data_bsr*)data;
1152 const INT row = predata->mgl_data[0].A.ROW;
1153 const INT nb = predata->mgl_data[0].A.nb;
1154 const INT maxit = predata->maxit;
1155 const INT m = row * nb;
1156
1157 INT i;
1158
1159 AMG_param amgparam;
1160 fasp_param_amg_init(&amgparam);
1161 amgparam.cycle_type = predata->cycle_type;
1162 amgparam.smoother = predata->smoother;
1163 amgparam.smooth_order = predata->smooth_order;
1164 amgparam.presmooth_iter = predata->presmooth_iter;
1165 amgparam.postsmooth_iter = predata->postsmooth_iter;
1166 amgparam.relaxation = predata->relaxation;
1167 amgparam.coarse_scaling = predata->coarse_scaling;
1168 amgparam.tentative_smooth = predata->tentative_smooth;
1169 amgparam.ILU_levels = predata->mgl_data->ILU_levels;
1170
1171 AMG_data_bsr* mgl = predata->mgl_data;
1172 mgl->b.row = m;
1173 fasp_darray_cp(m, r, mgl->b.val); // residual is an input
1174 mgl->x.row = m;
1175 fasp_dvec_set(m, &mgl->x, 0.0);
1176
1177 for (i = maxit; i--;) fasp_solver_mgcycle_bsr(mgl, &amgparam);
1178
1179 fasp_darray_cp(m, mgl->x.val, z);
1180}
1181
1194void fasp_precond_dbsr_amg_nk(REAL* r, REAL* z, void* data)
1195{
1196 precond_data_bsr* predata = (precond_data_bsr*)data;
1197 const INT row = predata->mgl_data[0].A.ROW;
1198 const INT nb = predata->mgl_data[0].A.nb;
1199 const INT maxit = predata->maxit;
1200 const INT m = row * nb;
1201
1202 INT i;
1203
1204 dCSRmat* A_nk = predata->A_nk;
1205 dCSRmat* P_nk = predata->P_nk;
1206 dCSRmat* R_nk = predata->R_nk;
1207
1208 fasp_darray_set(m, z, 0.0);
1209
1210 // local variables
1211 dvector r_nk, z_nk;
1212 fasp_dvec_alloc(A_nk->row, &r_nk);
1213 fasp_dvec_alloc(A_nk->row, &z_nk);
1214
1215 //----------------------
1216 // extra kernel solve
1217 //----------------------
1218 // r_nk = R_nk*r
1219 fasp_blas_dcsr_mxv(R_nk, r, r_nk.val);
1220
1221 // z_nk = A_nk^{-1}*r_nk
1222#if WITH_UMFPACK // use UMFPACK directly
1223 fasp_solver_umfpack(A_nk, &r_nk, &z_nk, 0);
1224#else
1225 fasp_coarse_itsolver(A_nk, &r_nk, &z_nk, 1e-12, 0);
1226#endif
1227
1228 // z = z + P_nk*z_nk;
1229 fasp_blas_dcsr_aAxpy(1.0, P_nk, z_nk.val, z);
1230
1231 //----------------------
1232 // AMG solve
1233 //----------------------
1234 AMG_param amgparam;
1235 fasp_param_amg_init(&amgparam);
1236 amgparam.cycle_type = predata->cycle_type;
1237 amgparam.smoother = predata->smoother;
1238 amgparam.smooth_order = predata->smooth_order;
1239 amgparam.presmooth_iter = predata->presmooth_iter;
1240 amgparam.postsmooth_iter = predata->postsmooth_iter;
1241 amgparam.relaxation = predata->relaxation;
1242 amgparam.coarse_scaling = predata->coarse_scaling;
1243 amgparam.tentative_smooth = predata->tentative_smooth;
1244 amgparam.ILU_levels = predata->mgl_data->ILU_levels;
1245
1246 AMG_data_bsr* mgl = predata->mgl_data;
1247 mgl->b.row = m;
1248 fasp_darray_cp(m, r, mgl->b.val); // residual is an input
1249 mgl->x.row = m; // fasp_dvec_set(m,&mgl->x,0.0);
1250 fasp_darray_cp(m, z, mgl->x.val);
1251
1252 for (i = maxit; i--;) fasp_solver_mgcycle_bsr(mgl, &amgparam);
1253
1254 fasp_darray_cp(m, mgl->x.val, z);
1255
1256 //----------------------
1257 // extra kernel solve
1258 //----------------------
1259 // r = r - A*z
1260 fasp_blas_dbsr_aAxpy(-1.0, &(predata->mgl_data[0].A), z, mgl->b.val);
1261
1262 // r_nk = R_nk*r
1263 fasp_blas_dcsr_mxv(R_nk, mgl->b.val, r_nk.val);
1264
1265 // z_nk = A_nk^{-1}*r_nk
1266#if WITH_UMFPACK // use UMFPACK directly
1267 fasp_solver_umfpack(A_nk, &r_nk, &z_nk, 0);
1268#else
1269 fasp_coarse_itsolver(A_nk, &r_nk, &z_nk, 1e-12, 0);
1270#endif
1271
1272 // z = z + P_nk*z_nk;
1273 fasp_blas_dcsr_aAxpy(1.0, P_nk, z_nk.val, z);
1274}
1275
1276double PreSmoother_time_zl = 0.0;
1277double PostSmoother_time_zl = 0.0;
1278double Krylov_time_zl = 0.0;
1279double Coarsen_time_zl = 0.0;
1280double AMLI_cycle_time_zl = 0.0;
1281
1294void fasp_precond_dbsr_namli(REAL* r, REAL* z, void* data)
1295{
1296 precond_data_bsr* pcdata = (precond_data_bsr*)data;
1297 const INT row = pcdata->mgl_data[0].A.ROW;
1298 const INT nb = pcdata->mgl_data[0].A.nb;
1299 const INT maxit = pcdata->maxit;
1300 const SHORT num_levels = pcdata->max_levels;
1301 const INT m = row * nb;
1302
1303 INT i;
1304
1305 AMG_param amgparam;
1306 fasp_param_amg_init(&amgparam);
1307 fasp_param_precbsr_to_amg(&amgparam, pcdata);
1308
1309 AMG_data_bsr* mgl = pcdata->mgl_data;
1310 mgl->b.row = m;
1311 fasp_darray_cp(m, r, mgl->b.val); // residual is an input
1312 mgl->x.row = m;
1313 fasp_dvec_set(m, &mgl->x, 0.0);
1314
1315 // REAL start_time, end_time; //! zhaoli
1316 // fasp_gettime(&start_time); //! zhaoli
1317
1318 for (i = maxit; i--;) fasp_solver_namli_bsr(mgl, &amgparam, 0, num_levels);
1319
1320 // fasp_gettime(&end_time); //! zhaoli
1321 // AMLI_cycle_time_zl += end_time - start_time;
1322 // printf("nonlinear AMLI-cycle time: %.4f\n", AMLI_cycle_time_zl); //! zhaoli
1323 // printf("PreSmoother_time_zl: %.4f\n", PreSmoother_time_zl); //! zhaoli
1324 // printf("PostSmoother_time_zl: %.4f\n", PostSmoother_time_zl); //! zhaoli
1325 // printf("Krylov_time_zl: %.4f\n", Krylov_time_zl); //! zhaoli
1326 // printf("Coarsen_time_zl: %.4f\n", Coarsen_time_zl); //! zhaoli
1327
1328 fasp_darray_cp(m, mgl->x.val, z);
1329}
1330
1331/*---------------------------------*/
1332/*-- End of File --*/
1333/*---------------------------------*/
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_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
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_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_param_precbsr_to_amg
void fasp_param_precbsr_to_amg(AMG_param *amgparam, const precond_data_bsr *pcdata)
Set AMG_param with precond_data.
Definition: AuxParam.c:882
fasp_param_amg_init
void fasp_param_amg_init(AMG_param *amgparam)
Initialize AMG parameters.
Definition: AuxParam.c:431
fasp_get_start_end
void fasp_get_start_end(const INT procid, const INT nprocs, const INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: AuxThreads.c:93
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_dvec_alloc
void fasp_dvec_alloc(const INT m, dvector *u)
Create dvector data space of REAL type.
Definition: AuxVector.c:105
fasp_blas_smat_mxv_nc7
void fasp_blas_smat_mxv_nc7(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 7*7 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:200
fasp_blas_smat_mxv_nc5
void fasp_blas_smat_mxv_nc5(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 5*5 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:176
fasp_blas_smat_mxv_nc4
void fasp_blas_smat_mxv_nc4(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 4*4 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:154
fasp_blas_smat_mxv
void fasp_blas_smat_mxv(const REAL *a, const REAL *b, REAL *c, const INT n)
Compute the product of a small full matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:238
fasp_blas_smat_mxv_nc3
void fasp_blas_smat_mxv_nc3(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 3*3 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:133
fasp_blas_smat_mxv_nc2
void fasp_blas_smat_mxv_nc2(const REAL *a, const REAL *b, REAL *c)
Compute the product of a 2*2 matrix a and a array b, stored in c.
Definition: BlaSmallMat.c:113
fasp_blas_dbsr_aAxpy
void fasp_blas_dbsr_aAxpy(const REAL alpha, const dBSRmat *A, const REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: BlaSpmvBSR.c:514
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_dbsr_diag_nc7
void fasp_precond_dbsr_diag_nc7(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreBSR.c:304
fasp_precond_dbsr_amg
void fasp_precond_dbsr_amg(REAL *r, REAL *z, void *data)
AMG preconditioner.
Definition: PreBSR.c:1149
fasp_precond_dbsr_ilu_mc_omp
void fasp_precond_dbsr_ilu_mc_omp(REAL *r, REAL *z, void *data)
Multi-thread Parallel ILU preconditioner based on graph coloring.
Definition: PreBSR.c:618
fasp_precond_dbsr_diag_nc5
void fasp_precond_dbsr_diag_nc5(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreBSR.c:259
fasp_precond_dbsr_diag_nc4
void fasp_precond_dbsr_diag_nc4(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreBSR.c:214
fasp_precond_dbsr_ilu_ls_omp
void fasp_precond_dbsr_ilu_ls_omp(REAL *r, REAL *z, void *data)
Multi-thread Parallel ILU preconditioner based on level schedule strategy.
Definition: PreBSR.c:844
fasp_precond_dbsr_amg_nk
void fasp_precond_dbsr_amg_nk(REAL *r, REAL *z, void *data)
AMG with extra near kernel solve preconditioner.
Definition: PreBSR.c:1194
fasp_precond_dbsr_diag
void fasp_precond_dbsr_diag(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreBSR.c:49
fasp_precond_dbsr_ilu
void fasp_precond_dbsr_ilu(REAL *r, REAL *z, void *data)
ILU preconditioner.
Definition: PreBSR.c:347
fasp_precond_dbsr_diag_nc2
void fasp_precond_dbsr_diag_nc2(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreBSR.c:124
fasp_precond_dbsr_diag_nc3
void fasp_precond_dbsr_diag_nc3(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreBSR.c:169
fasp_precond_dbsr_namli
void fasp_precond_dbsr_namli(REAL *r, REAL *z, void *data)
Nonlinear AMLI-cycle AMG preconditioner.
Definition: PreBSR.c:1294
fasp_solver_mgcycle_bsr
void fasp_solver_mgcycle_bsr(AMG_data_bsr *mgl, AMG_param *param)
Solve Ax=b with non-recursive multigrid cycle.
Definition: PreMGCycle.c:287
fasp_solver_namli_bsr
void fasp_solver_namli_bsr(AMG_data_bsr *mgl, AMG_param *param, INT l, INT num_levels)
Solve Ax=b with recursive nonlinear AMLI-cycle.
Definition: PreMGRecurAMLI.c:528
fasp_solver_umfpack
INT fasp_solver_umfpack(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by UMFpack.
Definition: XtrUmfpack.c:44
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
INT
#define INT
Definition: fasp.h:72
OPENMP_HOLDS
#define OPENMP_HOLDS
Definition: fasp_const.h:269
ERROR_ALLOC_MEM
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
ERROR_UNKNOWN
#define ERROR_UNKNOWN
Definition: fasp_const.h:56
AMG_data_bsr
Data for multigrid levels in dBSRmat format.
Definition: fasp_block.h:146
AMG_data_bsr::A
dBSRmat A
pointer to the matrix at level level_num
Definition: fasp_block.h:155
AMG_data_bsr::b
dvector b
pointer to the right-hand side at level level_num
Definition: fasp_block.h:164
AMG_data_bsr::ILU_levels
INT ILU_levels
number of levels use ILU smoother
Definition: fasp_block.h:217
AMG_data_bsr::x
dvector x
pointer to the iterative solution at level level_num
Definition: fasp_block.h:167
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
AMG_param::coarse_scaling
SHORT coarse_scaling
switch of scaling of the coarse grid correction
Definition: fasp.h:503
AMG_param::ILU_levels
SHORT ILU_levels
number of levels use ILU smoother
Definition: fasp.h:560
AMG_param::relaxation
REAL relaxation
relaxation parameter for Jacobi and SOR smoother
Definition: fasp.h:494
AMG_param::smoother
SHORT smoother
smoother type
Definition: fasp.h:482
AMG_param::cycle_type
SHORT cycle_type
type of AMG cycle
Definition: fasp.h:476
AMG_param::tentative_smooth
REAL tentative_smooth
relaxation parameter for smoothing the tentative prolongation
Definition: fasp.h:551
AMG_param::postsmooth_iter
SHORT postsmooth_iter
number of postsmoothers
Definition: fasp.h:491
AMG_param::presmooth_iter
SHORT presmooth_iter
number of presmoothers
Definition: fasp.h:488
AMG_param::smooth_order
SHORT smooth_order
smoother order
Definition: fasp.h:485
ILU_data
Data for ILU setup.
Definition: fasp.h:651
ILU_data::ijlu
INT * ijlu
integer array of row pointers and column indexes, the size is nzlu
Definition: fasp.h:669
ILU_data::luval
REAL * luval
nonzero entries of LU
Definition: fasp.h:672
ILU_data::jlevU
INT * jlevU
mapping from row to color for upper triangle
Definition: fasp.h:716
ILU_data::nlevU
INT nlevU
number of colors for upper triangle
Definition: fasp.h:704
ILU_data::nwork
INT nwork
work space size
Definition: fasp.h:678
ILU_data::nb
INT nb
block size for BSR type only
Definition: fasp.h:675
ILU_data::row
INT row
row number of matrix LU, m
Definition: fasp.h:660
ILU_data::nlevL
INT nlevL
number of colors for lower triangle
Definition: fasp.h:701
ILU_data::jlevL
INT * jlevL
mapping from row to color for lower triangle
Definition: fasp.h:713
ILU_data::work
REAL * work
work space
Definition: fasp.h:681
ILU_data::ilevL
INT * ilevL
number of vertices in each color for lower triangle
Definition: fasp.h:707
ILU_data::ilevU
INT * ilevU
number of vertices in each color for upper triangle
Definition: fasp.h:710
dBSRmat::nb
INT nb
dimension of each sub-block
Definition: fasp_block.h:46
dBSRmat::ROW
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:37
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
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
precond_data_bsr
Data for preconditioners in dBSRmat format.
Definition: fasp_block.h:271
precond_data_bsr::coarse_scaling
SHORT coarse_scaling
switch of scaling of the coarse grid correction
Definition: fasp_block.h:313
precond_data_bsr::A_nk
dCSRmat * A_nk
Matrix data for near kernal.
Definition: fasp_block.h:342
precond_data_bsr::relaxation
REAL relaxation
relaxation parameter for SOR smoother
Definition: fasp_block.h:307
precond_data_bsr::smoother
SHORT smoother
AMG smoother type.
Definition: fasp_block.h:292
precond_data_bsr::cycle_type
SHORT cycle_type
AMG cycle type.
Definition: fasp_block.h:289
precond_data_bsr::tentative_smooth
REAL tentative_smooth
smooth factor for smoothing the tentative prolongation
Definition: fasp_block.h:322
precond_data_bsr::postsmooth_iter
SHORT postsmooth_iter
number of postsmoothing
Definition: fasp_block.h:301
precond_data_bsr::R_nk
dCSRmat * R_nk
Resriction for near kernal.
Definition: fasp_block.h:348
precond_data_bsr::mgl_data
AMG_data_bsr * mgl_data
AMG preconditioner data.
Definition: fasp_block.h:328
precond_data_bsr::max_levels
INT max_levels
max number of AMG levels
Definition: fasp_block.h:283
precond_data_bsr::P_nk
dCSRmat * P_nk
Prolongation for near kernal.
Definition: fasp_block.h:345
precond_data_bsr::presmooth_iter
SHORT presmooth_iter
number of presmoothing
Definition: fasp_block.h:298
precond_data_bsr::maxit
INT maxit
max number of iterations of AMG preconditioner
Definition: fasp_block.h:280
precond_data_bsr::smooth_order
SHORT smooth_order
AMG smoother ordering.
Definition: fasp_block.h:295
precond_diag_bsr
Data for diagnal preconditioners in dBSRmat format.
Definition: fasp_block.h:255
precond_diag_bsr::nb
INT nb
dimension of each sub-block
Definition: fasp_block.h:258
precond_diag_bsr::diag
dvector diag
diagnal elements
Definition: fasp_block.h:261