Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
ItrSmootherBSR.c
Go to the documentation of this file.
1
17#include <math.h>
18
19#ifdef _OPENMP
20#include <omp.h>
21#endif
22
23#include "fasp.h"
24#include "fasp_functs.h"
25
26/*---------------------------------*/
27/*-- Declare Private Functions --*/
28/*---------------------------------*/
29
30#ifdef _OPENMP
31
32#if ILU_MC_OMP
33static inline void perm(const INT, const INT, const REAL*, const INT*, REAL*);
34static inline void invperm(const INT, const INT, const REAL*, const INT*, REAL*);
35#endif
36
37#endif
38
39REAL ilu_solve_time = 0.0;
41/*---------------------------------*/
42/*-- Public Functions --*/
43/*---------------------------------*/
44
59void fasp_smoother_dbsr_jacobi(dBSRmat* A, dvector* b, dvector* u)
60{
61 // members of A
62 const INT ROW = A->ROW;
63 const INT nb = A->nb;
64 const INT nb2 = nb * nb;
65 const INT size = ROW * nb2;
66 const INT* IA = A->IA;
67 const INT* JA = A->JA;
68 const REAL* val = A->val;
69
70 // local variables
71 INT i, k;
72 SHORT nthreads = 1, use_openmp = FALSE;
73 REAL* diaginv = (REAL*)fasp_mem_calloc(size, sizeof(REAL));
74
75#ifdef _OPENMP
76 if (ROW > OPENMP_HOLDS) {
77 use_openmp = TRUE;
78 nthreads = fasp_get_num_threads();
79 }
80#endif
81
82 // get all the diagonal sub-blocks
83 if (use_openmp) {
84 INT mybegin, myend, myid;
85#ifdef _OPENMP
86#pragma omp parallel for private(myid, mybegin, myend, i, k)
87#endif
88 for (myid = 0; myid < nthreads; myid++) {
89 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
90 for (i = mybegin; i < myend; i++) {
91 for (k = IA[i]; k < IA[i + 1]; ++k)
92 if (JA[k] == i)
93 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
94 }
95 }
96 } else {
97 for (i = 0; i < ROW; ++i) {
98 for (k = IA[i]; k < IA[i + 1]; ++k) {
99 if (JA[k] == i)
100 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
101 }
102 }
103 }
104
105 // compute the inverses of all the diagonal sub-blocks
106 if (nb > 1) {
107 if (use_openmp) {
108 INT mybegin, myend, myid;
109#ifdef _OPENMP
110#pragma omp parallel for private(myid, mybegin, myend, i)
111#endif
112 for (myid = 0; myid < nthreads; myid++) {
113 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
114 for (i = mybegin; i < myend; i++) {
115 fasp_smat_inv(diaginv + i * nb2, nb);
116 }
117 }
118 } else {
119 for (i = 0; i < ROW; ++i) {
120 fasp_smat_inv(diaginv + i * nb2, nb);
121 }
122 }
123 } else {
124 if (use_openmp) {
125 INT mybegin, myend, myid;
126#ifdef _OPENMP
127#pragma omp parallel for private(myid, mybegin, myend, i)
128#endif
129 for (myid = 0; myid < nthreads; myid++) {
130 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
131 for (i = mybegin; i < myend; i++) {
132 diaginv[i] = 1.0 / diaginv[i];
133 }
134 }
135 } else {
136 for (i = 0; i < ROW; ++i) {
137 // zero-diagonal should be tested previously
138 diaginv[i] = 1.0 / diaginv[i];
139 }
140 }
141 }
142
143 fasp_smoother_dbsr_jacobi1(A, b, u, diaginv);
144
145 fasp_mem_free(diaginv);
146 diaginv = NULL;
147}
148
163void fasp_smoother_dbsr_jacobi_setup(dBSRmat* A, REAL* diaginv)
164{
165 // members of A
166 const INT ROW = A->ROW;
167 const INT nb = A->nb;
168 const INT nb2 = nb * nb;
169 const INT* IA = A->IA;
170 const INT* JA = A->JA;
171 const REAL* val = A->val;
172
173 // local variables
174 INT i, k;
175
176 SHORT nthreads = 1, use_openmp = FALSE;
177
178#ifdef _OPENMP
179 if (ROW > OPENMP_HOLDS) {
180 use_openmp = TRUE;
181 nthreads = fasp_get_num_threads();
182 }
183#endif
184
185 // get all the diagonal sub-blocks
186 if (use_openmp) {
187 INT mybegin, myend, myid;
188#ifdef _OPENMP
189#pragma omp parallel for private(myid, mybegin, myend, i, k)
190#endif
191 for (myid = 0; myid < nthreads; myid++) {
192 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
193 for (i = mybegin; i < myend; i++) {
194 for (k = IA[i]; k < IA[i + 1]; ++k)
195 if (JA[k] == i)
196 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
197 }
198 }
199 } else {
200 for (i = 0; i < ROW; ++i) {
201 for (k = IA[i]; k < IA[i + 1]; ++k) {
202 if (JA[k] == i)
203 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
204 }
205 }
206 }
207
208 // compute the inverses of all the diagonal sub-blocks
209 if (nb > 1) {
210 if (use_openmp) {
211 INT mybegin, myend, myid;
212#ifdef _OPENMP
213#pragma omp parallel for private(myid, mybegin, myend, i)
214#endif
215 for (myid = 0; myid < nthreads; myid++) {
216 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
217 for (i = mybegin; i < myend; i++) {
218 fasp_smat_inv(diaginv + i * nb2, nb);
219 }
220 }
221 } else {
222 for (i = 0; i < ROW; ++i) {
223 fasp_smat_inv(diaginv + i * nb2, nb);
224 }
225 }
226 } else {
227 if (use_openmp) {
228 INT mybegin, myend, myid;
229#ifdef _OPENMP
230#pragma omp parallel for private(myid, mybegin, myend, i)
231#endif
232 for (myid = 0; myid < nthreads; myid++) {
233 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
234 for (i = mybegin; i < myend; i++) {
235 diaginv[i] = 1.0 / diaginv[i];
236 }
237 }
238 } else {
239 for (i = 0; i < ROW; ++i) {
240 // zero-diagonal should be tested previously
241 diaginv[i] = 1.0 / diaginv[i];
242 }
243 }
244 }
245}
246
263void fasp_smoother_dbsr_jacobi1(dBSRmat* A, dvector* b, dvector* u, REAL* diaginv)
264{
265 // members of A
266 const INT ROW = A->ROW;
267 const INT nb = A->nb;
268 const INT nb2 = nb * nb;
269 const INT size = ROW * nb;
270 const INT* IA = A->IA;
271 const INT* JA = A->JA;
272 REAL* val = A->val;
273
274 SHORT nthreads = 1, use_openmp = FALSE;
275
276#ifdef _OPENMP
277 if (ROW > OPENMP_HOLDS) {
278 use_openmp = TRUE;
279 nthreads = fasp_get_num_threads();
280 }
281#endif
282
283 // values of dvector b and u
284 REAL* b_val = b->val;
285 REAL* u_val = u->val;
286
287 // auxiliary array
288 REAL* b_tmp = NULL;
289
290 // local variables
291 INT i, j, k;
292 INT pb;
293
294 // b_tmp = b_val
295 b_tmp = (REAL*)fasp_mem_calloc(size, sizeof(REAL));
296 memcpy(b_tmp, b_val, size * sizeof(REAL));
297
298 // No need to assign the smoothing order since the result doesn't depend on it
299 if (nb == 1) {
300 if (use_openmp) {
301 INT mybegin, myend, myid;
302#ifdef _OPENMP
303#pragma omp parallel for private(myid, mybegin, myend, i, j, k)
304#endif
305 for (myid = 0; myid < nthreads; myid++) {
306 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
307 for (i = mybegin; i < myend; i++) {
308 for (k = IA[i]; k < IA[i + 1]; ++k) {
309 j = JA[k];
310 if (j != i) b_tmp[i] -= val[k] * u_val[j];
311 }
312 }
313 }
314#ifdef _OPENMP
315#pragma omp parallel for private(myid, mybegin, myend, i)
316#endif
317 for (myid = 0; myid < nthreads; myid++) {
318 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
319 for (i = mybegin; i < myend; i++) {
320 u_val[i] = b_tmp[i] * diaginv[i];
321 }
322 }
323 } else {
324 for (i = 0; i < ROW; ++i) {
325 for (k = IA[i]; k < IA[i + 1]; ++k) {
326 j = JA[k];
327 if (j != i) b_tmp[i] -= val[k] * u_val[j];
328 }
329 }
330 for (i = 0; i < ROW; ++i) {
331 u_val[i] = b_tmp[i] * diaginv[i];
332 }
333 }
334
335 fasp_mem_free(b_tmp);
336 b_tmp = NULL;
337 } else if (nb > 1) {
338 if (use_openmp) {
339 INT mybegin, myend, myid;
340#ifdef _OPENMP
341#pragma omp parallel for private(myid, mybegin, myend, i, pb, k, j)
342#endif
343 for (myid = 0; myid < nthreads; myid++) {
344 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
345 for (i = mybegin; i < myend; i++) {
346 pb = i * nb;
347 for (k = IA[i]; k < IA[i + 1]; ++k) {
348 j = JA[k];
349 if (j != i)
350 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb,
351 b_tmp + pb, nb);
352 }
353 }
354 }
355#ifdef _OPENMP
356#pragma omp parallel for private(myid, mybegin, myend, i, pb)
357#endif
358 for (myid = 0; myid < nthreads; myid++) {
359 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
360 for (i = mybegin; i < myend; i++) {
361 pb = i * nb;
362 fasp_blas_smat_mxv(diaginv + nb2 * i, b_tmp + pb, u_val + pb, nb);
363 }
364 }
365 } else {
366 for (i = 0; i < ROW; ++i) {
367 pb = i * nb;
368 for (k = IA[i]; k < IA[i + 1]; ++k) {
369 j = JA[k];
370 if (j != i)
371 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp + pb,
372 nb);
373 }
374 }
375
376 for (i = 0; i < ROW; ++i) {
377 pb = i * nb;
378 fasp_blas_smat_mxv(diaginv + nb2 * i, b_tmp + pb, u_val + pb, nb);
379 }
380 }
381 fasp_mem_free(b_tmp);
382 b_tmp = NULL;
383 } else {
384 printf("### ERROR: nb is illegal! [%s:%d]\n", __FILE__, __LINE__);
385 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
386 }
387}
388
410void fasp_smoother_dbsr_gs(dBSRmat* A, dvector* b, dvector* u, INT order, INT* mark)
411{
412 // members of A
413 const INT ROW = A->ROW;
414 const INT nb = A->nb;
415 const INT nb2 = nb * nb;
416 const INT size = ROW * nb2;
417 const INT* IA = A->IA;
418 const INT* JA = A->JA;
419 const REAL* val = A->val;
420
421 // local variables
422 INT i, k;
423 SHORT nthreads = 1, use_openmp = FALSE;
424 REAL* diaginv = (REAL*)fasp_mem_calloc(size, sizeof(REAL));
425
426#ifdef _OPENMP
427 if (ROW > OPENMP_HOLDS) {
428 use_openmp = TRUE;
429 nthreads = fasp_get_num_threads();
430 }
431#endif
432
433 // get all the diagonal sub-blocks
434 if (use_openmp) {
435 INT mybegin, myend, myid;
436#ifdef _OPENMP
437#pragma omp parallel for private(myid, mybegin, myend, i, k)
438#endif
439 for (myid = 0; myid < nthreads; myid++) {
440 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
441 for (i = mybegin; i < myend; i++) {
442 for (k = IA[i]; k < IA[i + 1]; ++k)
443 if (JA[k] == i)
444 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
445 }
446 }
447 } else {
448 for (i = 0; i < ROW; ++i) {
449 for (k = IA[i]; k < IA[i + 1]; ++k) {
450 if (JA[k] == i)
451 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
452 }
453 }
454 }
455
456 // compute the inverses of all the diagonal sub-blocks
457 if (nb > 1) {
458 if (use_openmp) {
459 INT mybegin, myend, myid;
460#ifdef _OPENMP
461#pragma omp parallel for private(myid, mybegin, myend, i)
462#endif
463 for (myid = 0; myid < nthreads; myid++) {
464 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
465 for (i = mybegin; i < myend; i++) {
466 fasp_smat_inv(diaginv + i * nb2, nb);
467 }
468 }
469 } else {
470 for (i = 0; i < ROW; ++i) {
471 fasp_smat_inv(diaginv + i * nb2, nb);
472 }
473 }
474 } else {
475 if (use_openmp) {
476 INT mybegin, myend, myid;
477#ifdef _OPENMP
478#pragma omp parallel for private(myid, mybegin, myend, i)
479#endif
480 for (myid = 0; myid < nthreads; myid++) {
481 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
482 for (i = mybegin; i < myend; i++) {
483 diaginv[i] = 1.0 / diaginv[i];
484 }
485 }
486 } else {
487 for (i = 0; i < ROW; ++i) {
488 // zero-diagonal should be tested previously
489 diaginv[i] = 1.0 / diaginv[i];
490 }
491 }
492 }
493
494 fasp_smoother_dbsr_gs1(A, b, u, order, mark, diaginv);
495
496 fasp_mem_free(diaginv);
497 diaginv = NULL;
498}
499
520void fasp_smoother_dbsr_gs1(
521 dBSRmat* A, dvector* b, dvector* u, INT order, INT* mark, REAL* diaginv)
522{
523 if (!mark) {
524 if (order == ASCEND) // smooth ascendingly
525 {
526 fasp_smoother_dbsr_gs_ascend(A, b, u, diaginv);
527 } else if (order == DESCEND) // smooth descendingly
528 {
529 fasp_smoother_dbsr_gs_descend(A, b, u, diaginv);
530 }
531 }
532 // smooth according to the order 'mark' defined by user
533 else {
534 fasp_smoother_dbsr_gs_order1(A, b, u, diaginv, mark);
535 }
536}
537
552void fasp_smoother_dbsr_gs_ascend(dBSRmat* A, dvector* b, dvector* u, REAL* diaginv)
553{
554 // members of A
555 const INT ROW = A->ROW;
556 const INT nb = A->nb;
557 const INT nb2 = nb * nb;
558 const INT* IA = A->IA;
559 const INT* JA = A->JA;
560 REAL* val = A->val;
561
562 // values of dvector b and u
563 REAL* b_val = b->val;
564 REAL* u_val = u->val;
565
566 // local variables
567 INT i, j, k;
568 INT pb;
569 REAL rhs = 0.0;
570
571 if (nb == 1) {
572 for (i = 0; i < ROW; ++i) {
573 rhs = b_val[i];
574 for (k = IA[i]; k < IA[i + 1]; ++k) {
575 j = JA[k];
576 if (j != i) rhs -= val[k] * u_val[j];
577 }
578 u_val[i] = rhs * diaginv[i];
579 }
580 } else if (nb > 1) {
581 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
582
583 for (i = 0; i < ROW; ++i) {
584 pb = i * nb;
585 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
586 for (k = IA[i]; k < IA[i + 1]; ++k) {
587 j = JA[k];
588 if (j != i)
589 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
590 }
591 fasp_blas_smat_mxv(diaginv + nb2 * i, b_tmp, u_val + pb, nb);
592 }
593
594 fasp_mem_free(b_tmp);
595 b_tmp = NULL;
596 } else {
597 printf("### ERROR: nb is illegal! [%s:%d]\n", __FILE__, __LINE__);
598 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
599 }
600}
601
619void fasp_smoother_dbsr_gs_ascend1(dBSRmat* A, dvector* b, dvector* u)
620{
621 // members of A
622 const INT ROW = A->ROW;
623 const INT nb = A->nb;
624 const INT nb2 = nb * nb;
625 const INT* IA = A->IA;
626 const INT* JA = A->JA;
627 REAL* val = A->val;
628
629 // values of dvector b and u
630 REAL* b_val = b->val;
631 REAL* u_val = u->val;
632
633 // local variables
634 INT i, j, k;
635 INT pb;
636 REAL rhs = 0.0;
637
638 if (nb == 1) {
639 for (i = 0; i < ROW; ++i) {
640 rhs = b_val[i];
641 for (k = IA[i]; k < IA[i + 1]; ++k) {
642 j = JA[k];
643 if (j != i) rhs -= val[k] * u_val[j];
644 }
645 u_val[i] = rhs;
646 }
647 } else if (nb > 1) {
648 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
649
650 for (i = 0; i < ROW; ++i) {
651 pb = i * nb;
652 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
653 for (k = IA[i]; k < IA[i + 1]; ++k) {
654 j = JA[k];
655 if (j != i)
656 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
657 }
658 memcpy(u_val + pb, b_tmp, nb * sizeof(REAL));
659 }
660
661 fasp_mem_free(b_tmp);
662 b_tmp = NULL;
663 } else {
664 printf("### ERROR: nb is illegal! [%s:%d]\n", __FILE__, __LINE__);
665 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
666 }
667}
668
683void fasp_smoother_dbsr_gs_descend(dBSRmat* A, dvector* b, dvector* u, REAL* diaginv)
684{
685 // members of A
686 const INT ROW = A->ROW;
687 const INT nb = A->nb;
688 const INT nb2 = nb * nb;
689 const INT* IA = A->IA;
690 const INT* JA = A->JA;
691 REAL* val = A->val;
692
693 // values of dvector b and u
694 REAL* b_val = b->val;
695 REAL* u_val = u->val;
696
697 // local variables
698 INT i, j, k;
699 INT pb;
700 REAL rhs = 0.0;
701
702 if (nb == 1) {
703 for (i = ROW - 1; i >= 0; i--) {
704 rhs = b_val[i];
705 for (k = IA[i]; k < IA[i + 1]; ++k) {
706 j = JA[k];
707 if (j != i) rhs -= val[k] * u_val[j];
708 }
709 u_val[i] = rhs * diaginv[i];
710 }
711 } else if (nb > 1) {
712 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
713
714 for (i = ROW - 1; i >= 0; i--) {
715 pb = i * nb;
716 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
717 for (k = IA[i]; k < IA[i + 1]; ++k) {
718 j = JA[k];
719 if (j != i)
720 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
721 }
722 fasp_blas_smat_mxv(diaginv + nb2 * i, b_tmp, u_val + pb, nb);
723 }
724
725 fasp_mem_free(b_tmp);
726 b_tmp = NULL;
727 } else {
728 printf("### ERROR: nb is illegal! [%s:%d]\n", __FILE__, __LINE__);
729 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
730 }
731}
732
751void fasp_smoother_dbsr_gs_descend1(dBSRmat* A, dvector* b, dvector* u)
752{
753 // members of A
754 const INT ROW = A->ROW;
755 const INT nb = A->nb;
756 const INT nb2 = nb * nb;
757 const INT* IA = A->IA;
758 const INT* JA = A->JA;
759 REAL* val = A->val;
760
761 // values of dvector b and u
762 REAL* b_val = b->val;
763 REAL* u_val = u->val;
764
765 // local variables
766 INT i, j, k;
767 INT pb;
768 REAL rhs = 0.0;
769
770 if (nb == 1) {
771 for (i = ROW - 1; i >= 0; i--) {
772 rhs = b_val[i];
773 for (k = IA[i]; k < IA[i + 1]; ++k) {
774 j = JA[k];
775 if (j != i) rhs -= val[k] * u_val[j];
776 }
777 u_val[i] = rhs;
778 }
779 } else if (nb > 1) {
780 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
781
782 for (i = ROW - 1; i >= 0; i--) {
783 pb = i * nb;
784 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
785 for (k = IA[i]; k < IA[i + 1]; ++k) {
786 j = JA[k];
787 if (j != i)
788 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
789 }
790 memcpy(u_val + pb, b_tmp, nb * sizeof(REAL));
791 }
792
793 fasp_mem_free(b_tmp);
794 b_tmp = NULL;
795 } else {
796 printf("### ERROR: nb is illegal! [%s:%d]\n", __FILE__, __LINE__);
797 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
798 }
799}
800
816void fasp_smoother_dbsr_gs_order1(
817 dBSRmat* A, dvector* b, dvector* u, REAL* diaginv, INT* mark)
818{
819 // members of A
820 const INT ROW = A->ROW;
821 const INT nb = A->nb;
822 const INT nb2 = nb * nb;
823 const INT* IA = A->IA;
824 const INT* JA = A->JA;
825 REAL* val = A->val;
826
827 // values of dvector b and u
828 REAL* b_val = b->val;
829 REAL* u_val = u->val;
830
831 // local variables
832 INT i, j, k;
833 INT I, pb;
834 REAL rhs = 0.0;
835
836 if (nb == 1) {
837 for (I = 0; I < ROW; ++I) {
838 i = mark[I];
839 rhs = b_val[i];
840 for (k = IA[i]; k < IA[i + 1]; ++k) {
841 j = JA[k];
842 if (j != i) rhs -= val[k] * u_val[j];
843 }
844 u_val[i] = rhs * diaginv[i];
845 }
846 } else if (nb > 1) {
847 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
848
849 for (I = 0; I < ROW; ++I) {
850 i = mark[I];
851 pb = i * nb;
852 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
853 for (k = IA[i]; k < IA[i + 1]; ++k) {
854 j = JA[k];
855 if (j != i)
856 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
857 }
858 fasp_blas_smat_mxv(diaginv + nb2 * i, b_tmp, u_val + pb, nb);
859 }
860
861 fasp_mem_free(b_tmp);
862 b_tmp = NULL;
863 } else {
864 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
865 }
866}
867
888void fasp_smoother_dbsr_gs_order2(
889 dBSRmat* A, dvector* b, dvector* u, INT* mark, REAL* work)
890{
891 // members of A
892 const INT ROW = A->ROW;
893 const INT nb = A->nb;
894 const INT nb2 = nb * nb;
895 const INT* IA = A->IA;
896 const INT* JA = A->JA;
897 REAL* val = A->val;
898
899 // values of dvector b and u
900 REAL* b_val = b->val;
901 REAL* u_val = u->val;
902
903 // auxiliary array
904 REAL* b_tmp = work;
905
906 // local variables
907 INT i, j, k, I, pb;
908 REAL rhs = 0.0;
909
910 if (nb == 1) {
911 for (I = 0; I < ROW; ++I) {
912 i = mark[I];
913 rhs = b_val[i];
914 for (k = IA[i]; k < IA[i + 1]; ++k) {
915 j = JA[k];
916 if (j != i) rhs -= val[k] * u_val[j];
917 }
918 u_val[i] = rhs;
919 }
920 } else if (nb > 1) {
921 for (I = 0; I < ROW; ++I) {
922 i = mark[I];
923 pb = i * nb;
924 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
925 for (k = IA[i]; k < IA[i + 1]; ++k) {
926 j = JA[k];
927 if (j != i)
928 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
929 }
930 memcpy(u_val + pb, b_tmp, nb * sizeof(REAL));
931 }
932 } else {
933 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
934 }
935}
936
959void fasp_smoother_dbsr_sor(
960 dBSRmat* A, dvector* b, dvector* u, INT order, INT* mark, REAL weight)
961{
962 // members of A
963 const INT ROW = A->ROW;
964 const INT nb = A->nb;
965 const INT nb2 = nb * nb;
966 const INT size = ROW * nb2;
967 const INT* IA = A->IA;
968 const INT* JA = A->JA;
969 const REAL* val = A->val;
970
971 // local variables
972 INT i, k;
973 REAL* diaginv = NULL;
974
975 SHORT nthreads = 1, use_openmp = FALSE;
976
977#ifdef _OPENMP
978 if (ROW > OPENMP_HOLDS) {
979 use_openmp = TRUE;
980 nthreads = fasp_get_num_threads();
981 }
982#endif
983
984 // allocate memory
985 diaginv = (REAL*)fasp_mem_calloc(size, sizeof(REAL));
986
987 // get all the diagonal sub-blocks
988 if (use_openmp) {
989 INT mybegin, myend, myid;
990#ifdef _OPENMP
991#pragma omp parallel for private(myid, mybegin, myend, i, k)
992#endif
993 for (myid = 0; myid < nthreads; myid++) {
994 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
995 for (i = mybegin; i < myend; i++) {
996 for (k = IA[i]; k < IA[i + 1]; ++k)
997 if (JA[k] == i)
998 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
999 }
1000 }
1001 } else {
1002 for (i = 0; i < ROW; ++i) {
1003 for (k = IA[i]; k < IA[i + 1]; ++k) {
1004 if (JA[k] == i)
1005 memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
1006 }
1007 }
1008 }
1009
1010 // compute the inverses of all the diagonal sub-blocks
1011 if (nb > 1) {
1012 if (use_openmp) {
1013 INT mybegin, myend, myid;
1014#ifdef _OPENMP
1015#pragma omp parallel for private(myid, mybegin, myend, i)
1016#endif
1017 for (myid = 0; myid < nthreads; myid++) {
1018 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1019 for (i = mybegin; i < myend; i++) {
1020 fasp_smat_inv(diaginv + i * nb2, nb);
1021 }
1022 }
1023 } else {
1024 for (i = 0; i < ROW; ++i) {
1025 fasp_smat_inv(diaginv + i * nb2, nb);
1026 }
1027 }
1028 } else {
1029 if (use_openmp) {
1030 INT mybegin, myend, myid;
1031#ifdef _OPENMP
1032#pragma omp parallel for private(myid, mybegin, myend, i)
1033#endif
1034 for (myid = 0; myid < nthreads; myid++) {
1035 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1036 for (i = mybegin; i < myend; i++) {
1037 diaginv[i] = 1.0 / diaginv[i];
1038 }
1039 }
1040 } else {
1041 for (i = 0; i < ROW; ++i) {
1042 // zero-diagonal should be tested previously
1043 diaginv[i] = 1.0 / diaginv[i];
1044 }
1045 }
1046 }
1047
1048 fasp_smoother_dbsr_sor1(A, b, u, order, mark, diaginv, weight);
1049
1050 fasp_mem_free(diaginv);
1051 diaginv = NULL;
1052}
1053
1075void fasp_smoother_dbsr_sor1(dBSRmat* A,
1076 dvector* b,
1077 dvector* u,
1078 INT order,
1079 INT* mark,
1080 REAL* diaginv,
1081 REAL weight)
1082{
1083 if (!mark) {
1084 if (order == ASCEND) // smooth ascendingly
1085 {
1086 fasp_smoother_dbsr_sor_ascend(A, b, u, diaginv, weight);
1087 } else if (order == DESCEND) // smooth descendingly
1088 {
1089 fasp_smoother_dbsr_sor_descend(A, b, u, diaginv, weight);
1090 }
1091 }
1092 // smooth according to the order 'mark' defined by user
1093 else {
1094 fasp_smoother_dbsr_sor_order(A, b, u, diaginv, mark, weight);
1095 }
1096}
1097
1115void fasp_smoother_dbsr_sor_ascend(
1116 dBSRmat* A, dvector* b, dvector* u, REAL* diaginv, REAL weight)
1117{
1118 // members of A
1119 const INT ROW = A->ROW;
1120 const INT nb = A->nb;
1121 const INT* IA = A->IA;
1122 const INT* JA = A->JA;
1123 const REAL* val = A->val;
1124
1125 // values of dvector b and u
1126 const REAL* b_val = b->val;
1127 REAL* u_val = u->val;
1128
1129 // local variables
1130 const INT nb2 = nb * nb;
1131 INT i, j, k;
1132 INT pb;
1133 REAL rhs = 0.0;
1134 REAL one_minus_weight = 1.0 - weight;
1135
1136#ifdef _OPENMP
1137 // variables for OpenMP
1138 INT myid, mybegin, myend;
1139 INT nthreads = fasp_get_num_threads();
1140#endif
1141
1142 if (nb == 1) {
1143#ifdef _OPENMP
1144 if (ROW > OPENMP_HOLDS) {
1145#pragma omp parallel for private(myid, mybegin, myend, i, rhs, k, j)
1146 for (myid = 0; myid < nthreads; myid++) {
1147 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1148 for (i = mybegin; i < myend; i++) {
1149 rhs = b_val[i];
1150 for (k = IA[i]; k < IA[i + 1]; ++k) {
1151 j = JA[k];
1152 if (j != i) rhs -= val[k] * u_val[j];
1153 }
1154 u_val[i] =
1155 one_minus_weight * u_val[i] + weight * (rhs * diaginv[i]);
1156 }
1157 }
1158 } else {
1159#endif
1160 for (i = 0; i < ROW; ++i) {
1161 rhs = b_val[i];
1162 for (k = IA[i]; k < IA[i + 1]; ++k) {
1163 j = JA[k];
1164 if (j != i) rhs -= val[k] * u_val[j];
1165 }
1166 u_val[i] = one_minus_weight * u_val[i] + weight * (rhs * diaginv[i]);
1167 }
1168#ifdef _OPENMP
1169 }
1170#endif
1171 } else if (nb > 1) {
1172#ifdef _OPENMP
1173 if (ROW > OPENMP_HOLDS) {
1174 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb * nthreads, sizeof(REAL));
1175#pragma omp parallel for private(myid, mybegin, myend, i, pb, k, j)
1176 for (myid = 0; myid < nthreads; myid++) {
1177 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1178 for (i = mybegin; i < myend; i++) {
1179 pb = i * nb;
1180 memcpy(b_tmp + myid * nb, b_val + pb, nb * sizeof(REAL));
1181 for (k = IA[i]; k < IA[i + 1]; ++k) {
1182 j = JA[k];
1183 if (j != i)
1184 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp,
1185 nb);
1186 }
1187 fasp_blas_smat_aAxpby(weight, diaginv + nb2 * i, b_tmp + myid * nb,
1188 one_minus_weight, u_val + pb, nb);
1189 }
1190 }
1191 fasp_mem_free(b_tmp);
1192 b_tmp = NULL;
1193 } else {
1194#endif
1195 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
1196 for (i = 0; i < ROW; ++i) {
1197 pb = i * nb;
1198 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
1199 for (k = IA[i]; k < IA[i + 1]; ++k) {
1200 j = JA[k];
1201 if (j != i)
1202 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
1203 }
1204 fasp_blas_smat_aAxpby(weight, diaginv + nb2 * i, b_tmp,
1205 one_minus_weight, u_val + pb, nb);
1206 }
1207 fasp_mem_free(b_tmp);
1208 b_tmp = NULL;
1209#ifdef _OPENMP
1210 }
1211#endif
1212 } else {
1213 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
1214 }
1215}
1216
1234void fasp_smoother_dbsr_sor_descend(
1235 dBSRmat* A, dvector* b, dvector* u, REAL* diaginv, REAL weight)
1236{
1237 // members of A
1238 const INT ROW = A->ROW;
1239 const INT nb = A->nb;
1240 const INT nb2 = nb * nb;
1241 const INT* IA = A->IA;
1242 const INT* JA = A->JA;
1243 REAL* val = A->val;
1244 const REAL one_minus_weight = 1.0 - weight;
1245
1246 // values of dvector b and u
1247 REAL* b_val = b->val;
1248 REAL* u_val = u->val;
1249
1250 // local variables
1251 INT i, j, k;
1252 INT pb;
1253 REAL rhs = 0.0;
1254
1255#ifdef _OPENMP
1256 // variables for OpenMP
1257 INT myid, mybegin, myend;
1258 INT nthreads = fasp_get_num_threads();
1259#endif
1260
1261 if (nb == 1) {
1262#ifdef _OPENMP
1263 if (ROW > OPENMP_HOLDS) {
1264#pragma omp parallel for private(myid, mybegin, myend, i, rhs, k, j)
1265 for (myid = 0; myid < nthreads; myid++) {
1266 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1267 mybegin = ROW - 1 - mybegin;
1268 myend = ROW - 1 - myend;
1269 for (i = mybegin; i > myend; i--) {
1270 rhs = b_val[i];
1271 for (k = IA[i]; k < IA[i + 1]; ++k) {
1272 j = JA[k];
1273 if (j != i) rhs -= val[k] * u_val[j];
1274 }
1275 u_val[i] =
1276 one_minus_weight * u_val[i] + weight * (rhs * diaginv[i]);
1277 }
1278 }
1279 } else {
1280#endif
1281 for (i = ROW - 1; i >= 0; i--) {
1282 rhs = b_val[i];
1283 for (k = IA[i]; k < IA[i + 1]; ++k) {
1284 j = JA[k];
1285 if (j != i) rhs -= val[k] * u_val[j];
1286 }
1287 u_val[i] = one_minus_weight * u_val[i] + weight * (rhs * diaginv[i]);
1288 }
1289#ifdef _OPENMP
1290 }
1291#endif
1292 } else if (nb > 1) {
1293#ifdef _OPENMP
1294 if (ROW > OPENMP_HOLDS) {
1295 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb * nthreads, sizeof(REAL));
1296#pragma omp parallel for private(myid, mybegin, myend, i, pb, k, j)
1297 for (myid = 0; myid < nthreads; myid++) {
1298 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1299 mybegin = ROW - 1 - mybegin;
1300 myend = ROW - 1 - myend;
1301 for (i = mybegin; i > myend; i--) {
1302 pb = i * nb;
1303 memcpy(b_tmp + myid * nb, b_val + pb, nb * sizeof(REAL));
1304 for (k = IA[i]; k < IA[i + 1]; ++k) {
1305 j = JA[k];
1306 if (j != i)
1307 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb,
1308 b_tmp + myid * nb, nb);
1309 }
1310 fasp_blas_smat_aAxpby(weight, diaginv + nb2 * i, b_tmp + myid * nb,
1311 one_minus_weight, u_val + pb, nb);
1312 }
1313 }
1314 fasp_mem_free(b_tmp);
1315 b_tmp = NULL;
1316 } else {
1317#endif
1318 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
1319 for (i = ROW - 1; i >= 0; i--) {
1320 pb = i * nb;
1321 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
1322 for (k = IA[i]; k < IA[i + 1]; ++k) {
1323 j = JA[k];
1324 if (j != i)
1325 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
1326 }
1327 fasp_blas_smat_aAxpby(weight, diaginv + nb2 * i, b_tmp,
1328 one_minus_weight, u_val + pb, nb);
1329 }
1330 fasp_mem_free(b_tmp);
1331 b_tmp = NULL;
1332#ifdef _OPENMP
1333 }
1334#endif
1335 } else {
1336 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
1337 }
1338}
1339
1358void fasp_smoother_dbsr_sor_order(
1359 dBSRmat* A, dvector* b, dvector* u, REAL* diaginv, INT* mark, REAL weight)
1360{
1361 // members of A
1362 const INT ROW = A->ROW;
1363 const INT nb = A->nb;
1364 const INT nb2 = nb * nb;
1365 const INT* IA = A->IA;
1366 const INT* JA = A->JA;
1367 REAL* val = A->val;
1368 const REAL one_minus_weight = 1.0 - weight;
1369
1370 // values of dvector b and u
1371 REAL* b_val = b->val;
1372 REAL* u_val = u->val;
1373
1374 // local variables
1375 INT i, j, k;
1376 INT I, pb;
1377 REAL rhs = 0.0;
1378
1379#ifdef _OPENMP
1380 // variables for OpenMP
1381 INT myid, mybegin, myend;
1382 INT nthreads = fasp_get_num_threads();
1383#endif
1384
1385 if (nb == 1) {
1386#ifdef _OPENMP
1387 if (ROW > OPENMP_HOLDS) {
1388#pragma omp parallel for private(myid, mybegin, myend, I, i, rhs, k, j)
1389 for (myid = 0; myid < nthreads; myid++) {
1390 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1391 for (I = mybegin; I < myend; ++I) {
1392 i = mark[I];
1393 rhs = b_val[i];
1394 for (k = IA[i]; k < IA[i + 1]; ++k) {
1395 j = JA[k];
1396 if (j != i) rhs -= val[k] * u_val[j];
1397 }
1398 u_val[i] =
1399 one_minus_weight * u_val[i] + weight * (rhs * diaginv[i]);
1400 }
1401 }
1402 } else {
1403#endif
1404 for (I = 0; I < ROW; ++I) {
1405 i = mark[I];
1406 rhs = b_val[i];
1407 for (k = IA[i]; k < IA[i + 1]; ++k) {
1408 j = JA[k];
1409 if (j != i) rhs -= val[k] * u_val[j];
1410 }
1411 u_val[i] = one_minus_weight * u_val[i] + weight * (rhs * diaginv[i]);
1412 }
1413#ifdef _OPENMP
1414 }
1415#endif
1416 } else if (nb > 1) {
1417#ifdef _OPENMP
1418 if (ROW > OPENMP_HOLDS) {
1419 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb * nthreads, sizeof(REAL));
1420#pragma omp parallel for private(myid, mybegin, myend, I, i, pb, k, j)
1421 for (myid = 0; myid < nthreads; myid++) {
1422 fasp_get_start_end(myid, nthreads, ROW, &mybegin, &myend);
1423 for (I = mybegin; I < myend; ++I) {
1424 i = mark[I];
1425 pb = i * nb;
1426 memcpy(b_tmp + myid * nb, b_val + pb, nb * sizeof(REAL));
1427 for (k = IA[i]; k < IA[i + 1]; ++k) {
1428 j = JA[k];
1429 if (j != i)
1430 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb,
1431 b_tmp + myid * nb, nb);
1432 }
1433 fasp_blas_smat_aAxpby(weight, diaginv + nb2 * i, b_tmp + myid * nb,
1434 one_minus_weight, u_val + pb, nb);
1435 }
1436 }
1437 fasp_mem_free(b_tmp);
1438 b_tmp = NULL;
1439 } else {
1440#endif
1441 REAL* b_tmp = (REAL*)fasp_mem_calloc(nb, sizeof(REAL));
1442 for (I = 0; I < ROW; ++I) {
1443 i = mark[I];
1444 pb = i * nb;
1445 memcpy(b_tmp, b_val + pb, nb * sizeof(REAL));
1446 for (k = IA[i]; k < IA[i + 1]; ++k) {
1447 j = JA[k];
1448 if (j != i)
1449 fasp_blas_smat_ymAx(val + k * nb2, u_val + j * nb, b_tmp, nb);
1450 }
1451 fasp_blas_smat_aAxpby(weight, diaginv + nb2 * i, b_tmp,
1452 one_minus_weight, u_val + pb, nb);
1453 }
1454 fasp_mem_free(b_tmp);
1455 b_tmp = NULL;
1456#ifdef _OPENMP
1457 }
1458#endif
1459 } else {
1460 fasp_chkerr(ERROR_NUM_BLOCKS, __FUNCTION__);
1461 }
1462}
1463
1479void fasp_smoother_dbsr_ilu(dBSRmat* A, dvector* b, dvector* x, void* data)
1480{
1481 ILU_data* iludata = (ILU_data*)data;
1482 const INT nb = iludata->nb, m = A->ROW * nb, memneed = 5 * m;
1483
1484 REAL *xval = x->val, *bval = b->val;
1485 REAL* zr = iludata->work + 3 * m;
1486 REAL* z = zr + m;
1487
1488 double start, end;
1489
1490 if (iludata->nwork < memneed) goto MEMERR;
1491
1493 fasp_darray_cp(m, bval, zr);
1494 fasp_blas_dbsr_aAxpy(-1.0, A, xval, zr);
1495
1497#ifdef _OPENMP
1498
1499#if ILU_MC_OMP
1500 REAL* tz = (REAL*)fasp_mem_calloc(A->ROW * A->nb, sizeof(REAL));
1501 REAL* tzr = (REAL*)fasp_mem_calloc(A->ROW * A->nb, sizeof(REAL));
1502 perm(A->ROW, A->nb, zr, iludata->jlevL, tzr);
1503
1504 fasp_gettime(&start);
1505 fasp_precond_dbsr_ilu_mc_omp(tzr, tz, iludata);
1506 fasp_gettime(&end);
1507
1508 invperm(A->ROW, A->nb, tz, iludata->jlevL, z);
1509 fasp_mem_free(tzr);
1510 tzr = NULL;
1511 fasp_mem_free(tz);
1512 tz = NULL;
1513#else
1514 fasp_gettime(&start);
1515 fasp_precond_dbsr_ilu_ls_omp(zr, z, iludata);
1516 fasp_gettime(&end);
1517#endif
1518
1519 ilu_solve_time += end - start;
1520
1521#else
1522
1523 fasp_gettime(&start);
1524 fasp_precond_dbsr_ilu(zr, z, iludata);
1525 fasp_gettime(&end);
1526 ilu_solve_time += end - start;
1527
1528#endif
1529
1531 fasp_blas_darray_axpy(m, 1, z, xval);
1532
1533 return;
1534
1535MEMERR:
1536 printf("### ERROR: ILU needs %d memory, only %d available! [%s:%d]\n", memneed,
1537 iludata->nwork, __FILE__, __LINE__);
1538 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1539}
1540
1541/*---------------------------------*/
1542/*-- Private Functions --*/
1543/*---------------------------------*/
1544
1545#ifdef _OPENMP
1546
1547#if ILU_MC_OMP
1548
1564static inline void perm(const INT n, const INT nb, const REAL* x, const INT* p, REAL* y)
1565{
1566 INT i, j, indx, indy;
1567
1568#ifdef _OPENMP
1569#pragma omp parallel for private(i, j, indx, indy)
1570#endif
1571 for (i = 0; i < n; ++i) {
1572 indx = p[i] * nb;
1573 indy = i * nb;
1574 for (j = 0; j < nb; ++j) {
1575 y[indy + j] = x[indx + j];
1576 }
1577 }
1578}
1579
1595static inline void
1596invperm(const INT n, const INT nb, const REAL* x, const INT* p, REAL* y)
1597{
1598 INT i, j, indx, indy;
1599
1600#ifdef _OPENMP
1601#pragma omp parallel for private(i, j, indx, indy)
1602#endif
1603 for (i = 0; i < n; ++i) {
1604 indx = i * nb;
1605 indy = p[i] * nb;
1606 for (j = 0; j < nb; ++j) {
1607 y[indy + j] = x[indx + j];
1608 }
1609 }
1610}
1611
1612#endif // end of ILU_MC_OMP
1613
1614#endif // end of _OPENMP
1615
1616/*---------------------------------*/
1617/*-- End of File --*/
1618/*---------------------------------*/
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_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_gettime
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:37
fasp_blas_darray_axpy
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:90
fasp_smat_inv
SHORT fasp_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
Definition: BlaSmallMatInv.c:603
fasp_blas_smat_aAxpby
void fasp_blas_smat_aAxpby(const REAL alpha, const REAL *A, const REAL *x, const REAL beta, REAL *y, const INT n)
Compute y:=alpha*A*x + beta*y.
Definition: BlaSmallMat.c:1140
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_ymAx
void fasp_blas_smat_ymAx(const REAL *A, const REAL *x, REAL *y, const INT n)
Compute y := y - Ax, where 'A' is a n*n dense matrix.
Definition: BlaSmallMat.c:1028
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_smoother_dbsr_gs_descend
void fasp_smoother_dbsr_gs_descend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Gauss-Seidel relaxation in the descending order.
Definition: ItrSmootherBSR.c:683
fasp_smoother_dbsr_gs1
void fasp_smoother_dbsr_gs1(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL *diaginv)
Gauss-Seidel relaxation.
Definition: ItrSmootherBSR.c:520
fasp_smoother_dbsr_gs_order2
void fasp_smoother_dbsr_gs_order2(dBSRmat *A, dvector *b, dvector *u, INT *mark, REAL *work)
Gauss-Seidel relaxation in the user-defined order.
Definition: ItrSmootherBSR.c:888
fasp_smoother_dbsr_gs_descend1
void fasp_smoother_dbsr_gs_descend1(dBSRmat *A, dvector *b, dvector *u)
Gauss-Seidel relaxation in the descending order.
Definition: ItrSmootherBSR.c:751
fasp_smoother_dbsr_jacobi1
void fasp_smoother_dbsr_jacobi1(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Jacobi relaxation.
Definition: ItrSmootherBSR.c:263
fasp_smoother_dbsr_gs_ascend
void fasp_smoother_dbsr_gs_ascend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Gauss-Seidel relaxation in the ascending order.
Definition: ItrSmootherBSR.c:552
fasp_smoother_dbsr_gs_ascend1
void fasp_smoother_dbsr_gs_ascend1(dBSRmat *A, dvector *b, dvector *u)
Gauss-Seidel relaxation in the ascending order.
Definition: ItrSmootherBSR.c:619
fasp_smoother_dbsr_sor_ascend
void fasp_smoother_dbsr_sor_ascend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, REAL weight)
SOR relaxation in the ascending order.
Definition: ItrSmootherBSR.c:1115
fasp_smoother_dbsr_ilu
void fasp_smoother_dbsr_ilu(dBSRmat *A, dvector *b, dvector *x, void *data)
ILU method as the smoother in solving Au=b with multigrid method.
Definition: ItrSmootherBSR.c:1479
fasp_smoother_dbsr_jacobi
void fasp_smoother_dbsr_jacobi(dBSRmat *A, dvector *b, dvector *u)
Jacobi relaxation.
Definition: ItrSmootherBSR.c:59
fasp_smoother_dbsr_gs
void fasp_smoother_dbsr_gs(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark)
Gauss-Seidel relaxation.
Definition: ItrSmootherBSR.c:410
fasp_smoother_dbsr_sor1
void fasp_smoother_dbsr_sor1(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL *diaginv, REAL weight)
SOR relaxation.
Definition: ItrSmootherBSR.c:1075
ilu_solve_time
REAL ilu_solve_time
Definition: ItrSmootherBSR.c:39
fasp_smoother_dbsr_sor_order
void fasp_smoother_dbsr_sor_order(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark, REAL weight)
SOR relaxation in the user-defined order.
Definition: ItrSmootherBSR.c:1358
fasp_smoother_dbsr_sor_descend
void fasp_smoother_dbsr_sor_descend(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, REAL weight)
SOR relaxation in the descending order.
Definition: ItrSmootherBSR.c:1234
fasp_smoother_dbsr_jacobi_setup
void fasp_smoother_dbsr_jacobi_setup(dBSRmat *A, REAL *diaginv)
Setup for jacobi relaxation, fetch the diagonal sub-block matrixes and make them inverse first.
Definition: ItrSmootherBSR.c:163
fasp_smoother_dbsr_gs_order1
void fasp_smoother_dbsr_gs_order1(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv, INT *mark)
Gauss-Seidel relaxation in the user-defined order.
Definition: ItrSmootherBSR.c:816
fasp_smoother_dbsr_sor
void fasp_smoother_dbsr_sor(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL weight)
SOR relaxation.
Definition: ItrSmootherBSR.c:959
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_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_ilu
void fasp_precond_dbsr_ilu(REAL *r, REAL *z, void *data)
ILU preconditioner.
Definition: PreBSR.c:347
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
ERROR_NUM_BLOCKS
#define ERROR_NUM_BLOCKS
Definition: fasp_const.h:27
ASCEND
#define ASCEND
Definition: fasp_const.h:249
OPENMP_HOLDS
#define OPENMP_HOLDS
Definition: fasp_const.h:269
TRUE
#define TRUE
Definition of logic type.
Definition: fasp_const.h:61
FALSE
#define FALSE
Definition: fasp_const.h:62
DESCEND
#define DESCEND
Definition: fasp_const.h:250
ERROR_ALLOC_MEM
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
ILU_data
Data for ILU setup.
Definition: fasp.h:651
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::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
dBSRmat
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:34
dBSRmat::val
REAL * val
Definition: fasp_block.h:57
dBSRmat::nb
INT nb
dimension of each sub-block
Definition: fasp_block.h:46
dBSRmat::IA
INT * IA
integer array of row pointers, the size is ROW+1
Definition: fasp_block.h:60
dBSRmat::ROW
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:37
dBSRmat::JA
INT * JA
Definition: fasp_block.h:64
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360