Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
ItrSmootherCSR.c
Go to the documentation of this file.
1
15#include <math.h>
16
17#ifdef _OPENMP
18#include <omp.h>
19#endif
20
21#include "fasp.h"
22#include "fasp_functs.h"
23
24/*---------------------------------*/
25/*-- Public Functions --*/
26/*---------------------------------*/
27
28// fasp_smoother_dcsr_jacobi_ff(x, A, b, nsweeps, ordering, relax);
34void fasp_smoother_dcsr_jacobi_ff(dvector* x,
35 dCSRmat* A,
36 dvector* b,
37 const INT nsweeps,
38 INT* ordering,
39 const REAL relax)
40{
41 REAL* r = (REAL*)fasp_mem_calloc(A->col, sizeof(REAL));
42 REAL* xval = (REAL*)fasp_mem_calloc(A->col, sizeof(REAL));
43 INT i, j, k, diagptr;
44 for (i = 0; i < A->col; ++i) {
45 xval[i] = x->val[i];
46 }
47 INT sweep_now;
48 for (sweep_now = 0; sweep_now < nsweeps; ++sweep_now) {
49 for (i = 0; i < A->row; ++i) {
50 if (ordering[i] == FGPT) {
51 diagptr = -1;
52 r[i] = b->val[i];
53 INT row_start = A->IA[i];
54 INT row_end = A->IA[i + 1];
55 for (j = row_start; j < row_end; ++j) {
56 k = A->JA[j];
57 // find diagonal entry and compute residual
58 if (k == i) diagptr = j;
59 r[i] -= A->val[j] * x->val[k]; // use old x entry
60 }
61 // update solution
62 xval[i] += relax * r[i] / A->val[diagptr]; // compute new x entry
63 }
64 }
65 for (i = 0; i < A->col; ++i) {
66 if (ordering[i] == FGPT) {
67 x->val[i] = xval[i];
68 }
69 }
70 }
71 fasp_mem_free(r);
72 fasp_mem_free(xval);
73 // printf("###DEBUG: fasp_smoother_dcsr_jacobi_ff() -- Done!\n");
74}
97void fasp_smoother_dcsr_jacobi(dvector* u,
98 const INT i_1,
99 const INT i_n,
100 const INT s,
101 dCSRmat* A,
102 dvector* b,
103 INT L,
104 const REAL w)
105{
106 const INT N = ABS(i_n - i_1) + 1;
107 const INT * ia = A->IA, *ja = A->JA;
108 const REAL *aval = A->val, *bval = b->val;
109 REAL* uval = u->val;
110
111 // local variables
112 INT i, j, k, begin_row, end_row;
113
114 // OpenMP variables
115#ifdef _OPENMP
116 INT myid, mybegin, myend;
117 INT nthreads = fasp_get_num_threads();
118#endif
119
120 REAL* t = (REAL*)fasp_mem_calloc(N, sizeof(REAL));
121 REAL* d = (REAL*)fasp_mem_calloc(N, sizeof(REAL));
122
123 while (L--) {
124
125 if (s > 0) {
126#ifdef _OPENMP
127 if (N > OPENMP_HOLDS) {
128#pragma omp parallel for private(myid, mybegin, myend, begin_row, end_row, i, k, j)
129 for (myid = 0; myid < nthreads; ++myid) {
130 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
131 mybegin += i_1;
132 myend += i_1;
133 for (i = mybegin; i < myend; i += s) {
134 t[i] = bval[i];
135 begin_row = ia[i], end_row = ia[i + 1];
136 for (k = begin_row; k < end_row; ++k) {
137 j = ja[k];
138 if (i != j)
139 t[i] -= aval[k] * uval[j];
140 else
141 d[i] = aval[k];
142 }
143 }
144 }
145 } else {
146#endif
147 for (i = i_1; i <= i_n; i += s) {
148 t[i] = bval[i];
149 begin_row = ia[i];
150 end_row = ia[i + 1];
151 for (k = begin_row; k < end_row; ++k) {
152 j = ja[k];
153 if (i != j)
154 t[i] -= aval[k] * uval[j];
155 else
156 d[i] = aval[k];
157 }
158 }
159#ifdef _OPENMP
160 }
161#endif
162
163#ifdef _OPENMP
164#pragma omp parallel for private(i)
165#endif
166 for (i = i_1; i <= i_n; i += s) {
167 if (ABS(d[i]) > SMALLREAL)
168 uval[i] = (1 - w) * uval[i] + w * t[i] / d[i];
169 }
170
171 }
172
173 else {
174
175#ifdef _OPENMP
176 if (N > OPENMP_HOLDS) {
177#pragma omp parallel for private(myid, mybegin, myend, i, begin_row, end_row, k, j)
178 for (myid = 0; myid < nthreads; myid++) {
179 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
180 mybegin = i_1 - mybegin;
181 myend = i_1 - myend;
182 for (i = mybegin; i > myend; i += s) {
183 t[i] = bval[i];
184 begin_row = ia[i], end_row = ia[i + 1];
185 for (k = begin_row; k < end_row; ++k) {
186 j = ja[k];
187 if (i != j)
188 t[i] -= aval[k] * uval[j];
189 else
190 d[i] = aval[k];
191 }
192 }
193 }
194 } else {
195#endif
196 for (i = i_1; i >= i_n; i += s) {
197 t[i] = bval[i];
198 begin_row = ia[i];
199 end_row = ia[i + 1];
200 for (k = begin_row; k < end_row; ++k) {
201 j = ja[k];
202 if (i != j)
203 t[i] -= aval[k] * uval[j];
204 else
205 d[i] = aval[k];
206 }
207 }
208#ifdef _OPENMP
209 }
210#endif
211
212#ifdef _OPENMP
213#pragma omp parallel for private(i)
214#endif
215 for (i = i_1; i >= i_n; i += s) {
216 if (ABS(d[i]) > SMALLREAL)
217 uval[i] = (1 - w) * uval[i] + w * t[i] / d[i];
218 }
219 }
220
221 } // end while
222
223 fasp_mem_free(t);
224 t = NULL;
225 fasp_mem_free(d);
226 d = NULL;
227
228 return;
229}
230
250void fasp_smoother_dcsr_gs(dvector* u,
251 const INT i_1,
252 const INT i_n,
253 const INT s,
254 dCSRmat* A,
255 dvector* b,
256 INT L)
257{
258 const INT * ia = A->IA, *ja = A->JA;
259 const REAL *aval = A->val, *bval = b->val;
260 REAL* uval = u->val;
261
262 // local variables
263 INT i, j, k, begin_row, end_row;
264 REAL t, d = 0.0;
265
266#ifdef _OPENMP
267 const INT N = ABS(i_n - i_1) + 1;
268 INT myid, mybegin, myend;
269 INT nthreads = fasp_get_num_threads();
270#endif
271
272 if (s > 0) {
273
274 while (L--) {
275#ifdef _OPENMP
276 if (N > OPENMP_HOLDS) {
277#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, d, k, \
278 j)
279 for (myid = 0; myid < nthreads; myid++) {
280 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
281 mybegin += i_1, myend += i_1;
282 for (i = mybegin; i < myend; i += s) {
283 t = bval[i];
284 begin_row = ia[i], end_row = ia[i + 1];
285#if DIAGONAL_PREF // diagonal first
286 d = aval[begin_row];
287 if (ABS(d) > SMALLREAL) {
288 for (k = begin_row + 1; k < end_row; ++k) {
289 j = ja[k];
290 t -= aval[k] * uval[j];
291 }
292 uval[i] = t / d;
293 }
294#else // general order
295 for (k = begin_row; k < end_row; ++k) {
296 j = ja[k];
297 if (i != j)
298 t -= aval[k] * uval[j];
299 else if (ABS(aval[k]) > SMALLREAL)
300 d = 1.e+0 / aval[k];
301 }
302 uval[i] = t * d;
303#endif // end DIAGONAL_PREF
304 } // end for i
305 }
306
307 }
308
309 else {
310#endif
311 for (i = i_1; i <= i_n; i += s) {
312 t = bval[i];
313 begin_row = ia[i];
314 end_row = ia[i + 1];
315
316#if DIAGONAL_PREF // diagonal first
317 d = aval[begin_row];
318 if (ABS(d) > SMALLREAL) {
319 for (k = begin_row + 1; k < end_row; ++k) {
320 j = ja[k];
321 t -= aval[k] * uval[j];
322 }
323 uval[i] = t / d;
324 }
325#else // general order
326 for (k = begin_row; k < end_row; ++k) {
327 j = ja[k];
328 if (i != j)
329 t -= aval[k] * uval[j];
330 else if (ABS(aval[k]) > SMALLREAL)
331 d = 1.e+0 / aval[k];
332 }
333 uval[i] = t * d;
334#endif
335 } // end for i
336#ifdef _OPENMP
337 }
338#endif
339 } // end while
340
341 } // if s
342 else {
343
344 while (L--) {
345#ifdef _OPENMP
346 if (N > OPENMP_HOLDS) {
347#pragma omp parallel for private(myid, mybegin, myend, i, begin_row, end_row, d, k, j, \
348 t)
349 for (myid = 0; myid < nthreads; myid++) {
350 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
351 mybegin = i_1 - mybegin;
352 myend = i_1 - myend;
353 for (i = mybegin; i > myend; i += s) {
354 t = bval[i];
355 begin_row = ia[i], end_row = ia[i + 1];
356#if DIAGONAL_PREF // diagonal first
357 d = aval[begin_row];
358 if (ABS(d) > SMALLREAL) {
359 for (k = begin_row + 1; k < end_row; ++k) {
360 j = ja[k];
361 t -= aval[k] * uval[j];
362 }
363 uval[i] = t / d;
364 }
365#else // general order
366 for (k = begin_row; k < end_row; ++k) {
367 j = ja[k];
368 if (i != j)
369 t -= aval[k] * uval[j];
370 else if (ABS(aval[k]) > SMALLREAL)
371 d = 1.0 / aval[k];
372 }
373 uval[i] = t * d;
374#endif
375 } // end for i
376 }
377 } else {
378#endif
379 for (i = i_1; i >= i_n; i += s) {
380 t = bval[i];
381 begin_row = ia[i];
382 end_row = ia[i + 1];
383#if DIAGONAL_PREF // diagonal first
384 d = aval[begin_row];
385 if (ABS(d) > SMALLREAL) {
386 for (k = begin_row + 1; k < end_row; ++k) {
387 j = ja[k];
388 t -= aval[k] * uval[j];
389 }
390 uval[i] = t / d;
391 }
392#else // general order
393 for (k = begin_row; k < end_row; ++k) {
394 j = ja[k];
395 if (i != j)
396 t -= aval[k] * uval[j];
397 else if (ABS(aval[k]) > SMALLREAL)
398 d = 1.0 / aval[k];
399 }
400 uval[i] = t * d;
401#endif
402 } // end for i
403#ifdef _OPENMP
404 }
405#endif
406 } // end while
407
408 } // end if
409
410 return;
411}
412
431void fasp_smoother_dcsr_gs_cf(
432 dvector* u, dCSRmat* A, dvector* b, INT L, INT* mark, const INT order)
433{
434 const INT nrow = b->row; // number of rows
435 const INT * ia = A->IA, *ja = A->JA;
436 const REAL *aval = A->val, *bval = b->val;
437 REAL* uval = u->val;
438
439 INT i, j, k, begin_row, end_row;
440 REAL t, d = 0.0;
441
442#ifdef _OPENMP
443 INT myid, mybegin, myend;
444 INT nthreads = fasp_get_num_threads();
445#endif
446
447 // F-point first, C-point second
448 if (order == FPFIRST) {
449
450 while (L--) {
451
452#ifdef _OPENMP
453 if (nrow > OPENMP_HOLDS) {
454#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
455 d)
456 for (myid = 0; myid < nthreads; myid++) {
457 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
458 for (i = mybegin; i < myend; i++) {
459 if (mark[i] != 1) {
460 t = bval[i];
461 begin_row = ia[i], end_row = ia[i + 1];
462#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
463 d = aval[begin_row];
464 for (k = begin_row + 1; k < end_row; k++) {
465 j = ja[k];
466 t -= aval[k] * uval[j];
467 } // end for k
468#else
469 for (k = begin_row; k < end_row; k++) {
470 j = ja[k];
471 if (i != j)
472 t -= aval[k] * uval[j];
473 else
474 d = aval[k];
475 } // end for k
476#endif // end if DIAG_PREF
477 if (ABS(d) > SMALLREAL) uval[i] = t / d;
478 }
479 } // end for i
480 }
481 } else {
482#endif
483 for (i = 0; i < nrow; i++) {
484 if (mark[i] != 1) {
485 t = bval[i];
486 begin_row = ia[i];
487 end_row = ia[i + 1];
488#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
489 d = aval[begin_row];
490 for (k = begin_row + 1; k < end_row; k++) {
491 j = ja[k];
492 t -= aval[k] * uval[j];
493 } // end for k
494#else
495 for (k = begin_row; k < end_row; k++) {
496 j = ja[k];
497 if (i != j)
498 t -= aval[k] * uval[j];
499 else
500 d = aval[k];
501 } // end for k
502#endif // end if DIAG_PREF
503 if (ABS(d) > SMALLREAL) uval[i] = t / d;
504 }
505 } // end for i
506#ifdef _OPENMP
507 }
508#endif
509
510#ifdef _OPENMP
511 if (nrow > OPENMP_HOLDS) {
512#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
513 d)
514 for (myid = 0; myid < nthreads; myid++) {
515 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
516 for (i = mybegin; i < myend; i++) {
517 if (mark[i] == 1) {
518 t = bval[i];
519 begin_row = ia[i], end_row = ia[i + 1];
520#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
521 d = aval[begin_row];
522 for (k = begin_row + 1; k < end_row; k++) {
523 j = ja[k];
524 t -= aval[k] * uval[j];
525 } // end for k
526#else
527 for (k = begin_row; k < end_row; k++) {
528 j = ja[k];
529 if (i != j)
530 t -= aval[k] * uval[j];
531 else
532 d = aval[k];
533 } // end for k
534#endif // end if DIAG_PREF
535 if (ABS(d) > SMALLREAL) uval[i] = t / d;
536 }
537 } // end for i
538 }
539 } else {
540#endif
541 for (i = 0; i < nrow; i++) {
542 if (mark[i] == 1) {
543 t = bval[i];
544 begin_row = ia[i];
545 end_row = ia[i + 1];
546#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
547 d = aval[begin_row];
548 for (k = begin_row + 1; k < end_row; k++) {
549 j = ja[k];
550 t -= aval[k] * uval[j];
551 } // end for k
552#else
553 for (k = begin_row; k < end_row; k++) {
554 j = ja[k];
555 if (i != j)
556 t -= aval[k] * uval[j];
557 else
558 d = aval[k];
559 } // end for k
560#endif // end if DIAG_PREF
561 if (ABS(d) > SMALLREAL) uval[i] = t / d;
562 }
563 } // end for i
564#ifdef _OPENMP
565 }
566#endif
567 } // end while
568
569 }
570
571 // C-point first, F-point second
572 else {
573
574 while (L--) {
575#ifdef _OPENMP
576 if (nrow > OPENMP_HOLDS) {
577#pragma omp parallel for private(myid, mybegin, myend, t, i, begin_row, end_row, k, j, \
578 d)
579 for (myid = 0; myid < nthreads; myid++) {
580 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
581 for (i = mybegin; i < myend; i++) {
582 if (mark[i] == 1) {
583 t = bval[i];
584 begin_row = ia[i], end_row = ia[i + 1];
585#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
586 d = aval[begin_row];
587 for (k = begin_row + 1; k < end_row; k++) {
588 j = ja[k];
589 t -= aval[k] * uval[j];
590 } // end for k
591#else
592 for (k = begin_row; k < end_row; k++) {
593 j = ja[k];
594 if (i != j)
595 t -= aval[k] * uval[j];
596 else
597 d = aval[k];
598 } // end for k
599#endif // end if DIAG_PREF
600 if (ABS(d) > SMALLREAL) uval[i] = t / d;
601 }
602 } // end for i
603 }
604 } else {
605#endif
606 for (i = 0; i < nrow; i++) {
607 if (mark[i] == 1) {
608 t = bval[i];
609 begin_row = ia[i];
610 end_row = ia[i + 1];
611#if DIAGONAL_PREF // Added by Chensong on 09/22/2012
612 d = aval[begin_row];
613 for (k = begin_row + 1; k < end_row; k++) {
614 j = ja[k];
615 t -= aval[k] * uval[j];
616 } // end for k
617#else
618 for (k = begin_row; k < end_row; k++) {
619 j = ja[k];
620 if (i != j)
621 t -= aval[k] * uval[j];
622 else
623 d = aval[k];
624 } // end for k
625#endif // end if DIAG_PREF
626 if (ABS(d) > SMALLREAL) uval[i] = t / d;
627 }
628 } // end for i
629#ifdef _OPENMP
630 }
631#endif
632
633#ifdef _OPENMP
634 if (nrow > OPENMP_HOLDS) {
635#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
636 d)
637 for (myid = 0; myid < nthreads; myid++) {
638 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
639 for (i = mybegin; i < myend; i++) {
640 if (mark[i] != 1) {
641 t = bval[i];
642 begin_row = ia[i], end_row = ia[i + 1];
643#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
644 d = aval[begin_row];
645 for (k = begin_row + 1; k < end_row; k++) {
646 j = ja[k];
647 t -= aval[k] * uval[j];
648 } // end for k
649#else
650 for (k = begin_row; k < end_row; k++) {
651 j = ja[k];
652 if (i != j)
653 t -= aval[k] * uval[j];
654 else
655 d = aval[k];
656 } // end for k
657#endif // end if DIAG_PREF
658 if (ABS(d) > SMALLREAL) uval[i] = t / d;
659 }
660 } // end for i
661 }
662 } else {
663#endif
664 for (i = 0; i < nrow; i++) {
665 if (mark[i] != 1) {
666 t = bval[i];
667 begin_row = ia[i];
668 end_row = ia[i + 1];
669#if DIAGONAL_PREF // Added by Chensong on 09/22/2012
670 d = aval[begin_row];
671 for (k = begin_row + 1; k < end_row; k++) {
672 j = ja[k];
673 t -= aval[k] * uval[j];
674 } // end for k
675#else
676 for (k = begin_row; k < end_row; k++) {
677 j = ja[k];
678 if (i != j)
679 t -= aval[k] * uval[j];
680 else
681 d = aval[k];
682 } // end for k
683#endif // end if DIAG_PREF
684 if (ABS(d) > SMALLREAL) uval[i] = t / d;
685 }
686 } // end for i
687#ifdef _OPENMP
688 }
689#endif
690 } // end while
691
692 } // end if order
693
694 return;
695}
696
712void fasp_smoother_dcsr_gs_ff(dvector* u, dCSRmat* A, dvector* b, INT L, INT* mark)
713{
714 const INT nrow = b->row; // number of rows
715 const INT * ia = A->IA, *ja = A->JA;
716 const REAL *aval = A->val, *bval = b->val;
717 REAL* uval = u->val;
718
719 INT i, j, k, begin_row, end_row;
720 REAL t, d = 0.0;
721
722#ifdef _OPENMP
723 INT myid, mybegin, myend;
724 INT nthreads = fasp_get_num_threads();
725#endif
726
727 while (L--) {
728
729#ifdef _OPENMP
730 if (nrow > OPENMP_HOLDS) {
731#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
732 d)
733 for (myid = 0; myid < nthreads; myid++) {
734 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
735 for (i = mybegin; i < myend; i++) {
736 if (mark[i] != 1) {
737 t = bval[i];
738 begin_row = ia[i], end_row = ia[i + 1];
739#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
740 d = aval[begin_row];
741 for (k = begin_row + 1; k < end_row; k++) {
742 j = ja[k];
743 t -= aval[k] * uval[j];
744 } // end for k
745#else
746 for (k = begin_row; k < end_row; k++) {
747 j = ja[k];
748 if (i != j)
749 t -= aval[k] * uval[j];
750 else
751 d = aval[k];
752 } // end for k
753#endif // end if DIAG_PREF
754 if (ABS(d) > SMALLREAL) uval[i] = t / d;
755 }
756 } // end for i
757 }
758 } else {
759#endif
760 for (i = 0; i < nrow; i++) {
761 if (mark[i] != 1) {
762 t = bval[i];
763 begin_row = ia[i];
764 end_row = ia[i + 1];
765#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
766 d = aval[begin_row];
767 for (k = begin_row + 1; k < end_row; k++) {
768 j = ja[k];
769 t -= aval[k] * uval[j];
770 } // end for k
771#else
772 for (k = begin_row; k < end_row; k++) {
773 j = ja[k];
774 if (i != j)
775 t -= aval[k] * uval[j];
776 else
777 d = aval[k];
778 } // end for k
779#endif // end if DIAG_PREF
780 if (ABS(d) > SMALLREAL) uval[i] = t / d;
781 }
782 } // end for i
783#ifdef _OPENMP
784 }
785#endif
786
787 } // end while
788
789 return;
790}
791
807void fasp_smoother_dcsr_sgs(dvector* u, dCSRmat* A, dvector* b, INT L)
808{
809 const INT nm1 = b->row - 1;
810 const INT * ia = A->IA, *ja = A->JA;
811 const REAL *aval = A->val, *bval = b->val;
812 REAL* uval = u->val;
813
814 // local variables
815 INT i, j, k, begin_row, end_row;
816 REAL t, d = 0;
817
818#ifdef _OPENMP
819 INT myid, mybegin, myend, up;
820 INT nthreads = fasp_get_num_threads();
821#endif
822
823 while (L--) {
824 // forward sweep
825#ifdef _OPENMP
826 up = nm1 + 1;
827 if (up > OPENMP_HOLDS) {
828#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, j, k, \
829 d)
830 for (myid = 0; myid < nthreads; myid++) {
831 fasp_get_start_end(myid, nthreads, up, &mybegin, &myend);
832 for (i = mybegin; i < myend; i++) {
833 t = bval[i];
834 begin_row = ia[i], end_row = ia[i + 1];
835 for (k = begin_row; k < end_row; ++k) {
836 j = ja[k];
837 if (i != j)
838 t -= aval[k] * uval[j];
839 else
840 d = aval[k];
841 } // end for k
842 if (ABS(d) > SMALLREAL) uval[i] = t / d;
843 } // end for i
844 }
845 } else {
846#endif
847 for (i = 0; i <= nm1; ++i) {
848 t = bval[i];
849 begin_row = ia[i];
850 end_row = ia[i + 1];
851 for (k = begin_row; k < end_row; ++k) {
852 j = ja[k];
853 if (i != j)
854 t -= aval[k] * uval[j];
855 else
856 d = aval[k];
857 } // end for k
858 if (ABS(d) > SMALLREAL) uval[i] = t / d;
859 } // end for i
860#ifdef _OPENMP
861 }
862#endif
863
864 // backward sweep
865#ifdef _OPENMP
866 up = nm1;
867 if (up > OPENMP_HOLDS) {
868#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
869 d)
870 for (myid = 0; myid < nthreads; myid++) {
871 fasp_get_start_end(myid, nthreads, up, &mybegin, &myend);
872 mybegin = nm1 - 1 - mybegin;
873 myend = nm1 - 1 - myend;
874 for (i = mybegin; i > myend; i--) {
875 t = bval[i];
876 begin_row = ia[i], end_row = ia[i + 1];
877 for (k = begin_row; k < end_row; k++) {
878 j = ja[k];
879 if (i != j)
880 t -= aval[k] * uval[j];
881 else
882 d = aval[k];
883 } // end for k
884 if (ABS(d) > SMALLREAL) uval[i] = t / d;
885 } // end for i
886 }
887 } else {
888#endif
889 for (i = nm1 - 1; i >= 0; --i) {
890 t = bval[i];
891 begin_row = ia[i];
892 end_row = ia[i + 1];
893 for (k = begin_row; k < end_row; ++k) {
894 j = ja[k];
895 if (i != j)
896 t -= aval[k] * uval[j];
897 else
898 d = aval[k];
899 } // end for k
900 if (ABS(d) > SMALLREAL) uval[i] = t / d;
901 } // end for i
902#ifdef _OPENMP
903 }
904#endif
905 } // end while
906
907 return;
908}
909
931void fasp_smoother_dcsr_sor(dvector* u,
932 const INT i_1,
933 const INT i_n,
934 const INT s,
935 dCSRmat* A,
936 dvector* b,
937 INT L,
938 const REAL w)
939{
940 const INT * ia = A->IA, *ja = A->JA;
941 const REAL *aval = A->val, *bval = b->val;
942 REAL* uval = u->val;
943
944 // local variables
945 INT i, j, k, begin_row, end_row;
946 REAL t, d = 0;
947
948#ifdef _OPENMP
949 const INT N = ABS(i_n - i_1) + 1;
950 INT myid, mybegin, myend;
951 INT nthreads = fasp_get_num_threads();
952#endif
953
954 while (L--) {
955 if (s > 0) {
956#ifdef _OPENMP
957 if (N > OPENMP_HOLDS) {
958#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
959 d)
960 for (myid = 0; myid < nthreads; myid++) {
961 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
962 mybegin += i_1, myend += i_1;
963 for (i = mybegin; i < myend; i += s) {
964 t = bval[i];
965 begin_row = ia[i], end_row = ia[i + 1];
966 for (k = begin_row; k < end_row; k++) {
967 j = ja[k];
968 if (i != j)
969 t -= aval[k] * uval[j];
970 else
971 d = aval[k];
972 }
973 if (ABS(d) > SMALLREAL)
974 uval[i] = w * (t / d) + (1 - w) * uval[i];
975 }
976 }
977
978 } else {
979#endif
980 for (i = i_1; i <= i_n; i += s) {
981 t = bval[i];
982 begin_row = ia[i];
983 end_row = ia[i + 1];
984 for (k = begin_row; k < end_row; ++k) {
985 j = ja[k];
986 if (i != j)
987 t -= aval[k] * uval[j];
988 else
989 d = aval[k];
990 }
991 if (ABS(d) > SMALLREAL) uval[i] = w * (t / d) + (1 - w) * uval[i];
992 }
993#ifdef _OPENMP
994 }
995#endif
996 } else {
997#ifdef _OPENMP
998 if (N > OPENMP_HOLDS) {
999#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
1000 d)
1001 for (myid = 0; myid < nthreads; myid++) {
1002 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1003 mybegin = i_1 - mybegin, myend = i_1 - myend;
1004 for (i = mybegin; i > myend; i += s) {
1005 t = bval[i];
1006 begin_row = ia[i], end_row = ia[i + 1];
1007 for (k = begin_row; k < end_row; ++k) {
1008 j = ja[k];
1009 if (i != j)
1010 t -= aval[k] * uval[j];
1011 else
1012 d = aval[k];
1013 }
1014 if (ABS(d) > SMALLREAL)
1015 uval[i] = w * (t / d) + (1 - w) * uval[i];
1016 }
1017 }
1018 } else {
1019#endif
1020 for (i = i_1; i >= i_n; i += s) {
1021 t = bval[i];
1022 begin_row = ia[i];
1023 end_row = ia[i + 1];
1024 for (k = begin_row; k < end_row; ++k) {
1025 j = ja[k];
1026 if (i != j)
1027 t -= aval[k] * uval[j];
1028 else
1029 d = aval[k];
1030 }
1031 if (ABS(d) > SMALLREAL) uval[i] = w * (t / d) + (1 - w) * uval[i];
1032 }
1033#ifdef _OPENMP
1034 }
1035#endif
1036 }
1037 } // end while
1038
1039 return;
1040}
1041
1061void fasp_smoother_dcsr_sor_cf(
1062 dvector* u, dCSRmat* A, dvector* b, INT L, const REAL w, INT* mark, const INT order)
1063{
1064 const INT nrow = b->row; // number of rows
1065 const INT * ia = A->IA, *ja = A->JA;
1066 const REAL *aval = A->val, *bval = b->val;
1067 REAL* uval = u->val;
1068
1069 // local variables
1070 INT i, j, k, begin_row, end_row;
1071 REAL t, d = 0.0;
1072
1073#ifdef _OPENMP
1074 INT myid, mybegin, myend;
1075 INT nthreads = fasp_get_num_threads();
1076#endif
1077
1078 // F-point first
1079 if (order == -1) {
1080 while (L--) {
1081#ifdef _OPENMP
1082 if (nrow > OPENMP_HOLDS) {
1083#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
1084 d)
1085 for (myid = 0; myid < nthreads; myid++) {
1086 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
1087 for (i = mybegin; i < myend; i++) {
1088 if (mark[i] == 0 || mark[i] == 2) {
1089 t = bval[i];
1090 begin_row = ia[i], end_row = ia[i + 1];
1091 for (k = begin_row; k < end_row; k++) {
1092 j = ja[k];
1093 if (i != j)
1094 t -= aval[k] * uval[j];
1095 else
1096 d = aval[k];
1097 } // end for k
1098 if (ABS(d) > SMALLREAL)
1099 uval[i] = w * (t / d) + (1 - w) * uval[i];
1100 }
1101 }
1102 }
1103 } // end for i
1104 else {
1105#endif
1106 for (i = 0; i < nrow; i++) {
1107 if (mark[i] == 0 || mark[i] == 2) {
1108 t = bval[i];
1109 begin_row = ia[i];
1110 end_row = ia[i + 1];
1111 for (k = begin_row; k < end_row; k++) {
1112 j = ja[k];
1113 if (i != j)
1114 t -= aval[k] * uval[j];
1115 else
1116 d = aval[k];
1117 } // end for k
1118 if (ABS(d) > SMALLREAL)
1119 uval[i] = w * (t / d) + (1 - w) * uval[i];
1120 }
1121 } // end for i
1122#ifdef _OPENMP
1123 }
1124#endif
1125
1126#ifdef _OPENMP
1127 if (nrow > OPENMP_HOLDS) {
1128#pragma omp parallel for private(myid, i, mybegin, myend, t, begin_row, end_row, k, j, \
1129 d)
1130 for (myid = 0; myid < nthreads; myid++) {
1131 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
1132 for (i = mybegin; i < myend; i++) {
1133 if (mark[i] == 1) {
1134 t = bval[i];
1135 begin_row = ia[i], end_row = ia[i + 1];
1136 for (k = begin_row; k < end_row; k++) {
1137 j = ja[k];
1138 if (i != j)
1139 t -= aval[k] * uval[j];
1140 else
1141 d = aval[k];
1142 } // end for k
1143 if (ABS(d) > SMALLREAL)
1144 uval[i] = w * (t / d) + (1 - w) * uval[i];
1145 }
1146 } // end for i
1147 }
1148 } else {
1149#endif
1150 for (i = 0; i < nrow; i++) {
1151 if (mark[i] == 1) {
1152 t = bval[i];
1153 begin_row = ia[i];
1154 end_row = ia[i + 1];
1155 for (k = begin_row; k < end_row; k++) {
1156 j = ja[k];
1157 if (i != j)
1158 t -= aval[k] * uval[j];
1159 else
1160 d = aval[k];
1161 } // end for k
1162 if (ABS(d) > SMALLREAL)
1163 uval[i] = w * (t / d) + (1 - w) * uval[i];
1164 }
1165 } // end for i
1166#ifdef _OPENMP
1167 }
1168#endif
1169 } // end while
1170 } else {
1171 while (L--) {
1172#ifdef _OPENMP
1173 if (nrow > OPENMP_HOLDS) {
1174#pragma omp parallel for private(myid, mybegin, myend, i, t, k, j, d, begin_row, \
1175 end_row)
1176 for (myid = 0; myid < nthreads; myid++) {
1177 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
1178 for (i = mybegin; i < myend; i++) {
1179 if (mark[i] == 1) {
1180 t = bval[i];
1181 begin_row = ia[i], end_row = ia[i + 1];
1182 for (k = begin_row; k < end_row; k++) {
1183 j = ja[k];
1184 if (i != j)
1185 t -= aval[k] * uval[j];
1186 else
1187 d = aval[k];
1188 } // end for k
1189 if (ABS(d) > SMALLREAL)
1190 uval[i] = w * (t / d) + (1 - w) * uval[i];
1191 }
1192 } // end for i
1193 }
1194 } else {
1195#endif
1196 for (i = 0; i < nrow; i++) {
1197 if (mark[i] == 1) {
1198 t = bval[i];
1199 begin_row = ia[i];
1200 end_row = ia[i + 1];
1201 for (k = begin_row; k < end_row; k++) {
1202 j = ja[k];
1203 if (i != j)
1204 t -= aval[k] * uval[j];
1205 else
1206 d = aval[k];
1207 } // end for k
1208 if (ABS(d) > SMALLREAL)
1209 uval[i] = w * (t / d) + (1 - w) * uval[i];
1210 }
1211 } // end for i
1212#ifdef _OPENMP
1213 }
1214#endif
1215
1216#ifdef _OPENMP
1217 if (nrow > OPENMP_HOLDS) {
1218#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, \
1219 d)
1220 for (myid = 0; myid < nthreads; myid++) {
1221 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
1222 for (i = mybegin; i < myend; i++) {
1223 if (mark[i] != 1) {
1224 t = bval[i];
1225 begin_row = ia[i], end_row = ia[i + 1];
1226 for (k = begin_row; k < end_row; k++) {
1227 j = ja[k];
1228 if (i != j)
1229 t -= aval[k] * uval[j];
1230 else
1231 d = aval[k];
1232 } // end for k
1233 if (ABS(d) > SMALLREAL)
1234 uval[i] = w * (t / d) + (1 - w) * uval[i];
1235 }
1236 }
1237 }
1238 } // end for i
1239 else {
1240#endif
1241 for (i = 0; i < nrow; i++) {
1242 if (mark[i] != 1) {
1243 t = bval[i];
1244 begin_row = ia[i];
1245 end_row = ia[i + 1];
1246 for (k = begin_row; k < end_row; k++) {
1247 j = ja[k];
1248 if (i != j)
1249 t -= aval[k] * uval[j];
1250 else
1251 d = aval[k];
1252 } // end for k
1253 if (ABS(d) > SMALLREAL)
1254 uval[i] = w * (t / d) + (1 - w) * uval[i];
1255 }
1256 } // end for i
1257#ifdef _OPENMP
1258 }
1259#endif
1260 } // end while
1261 }
1262
1263 return;
1264}
1265
1279void fasp_smoother_dcsr_ilu(dCSRmat* A, dvector* b, dvector* x, void* data)
1280{
1281 const INT m = A->row, m2 = 2 * m, memneed = 3 * m;
1282 const ILU_data* iludata = (ILU_data*)data;
1283
1284 REAL* zz = iludata->work;
1285 REAL* zr = iludata->work + m;
1286 REAL* z = iludata->work + m2;
1287
1288 if (iludata->nwork < memneed) goto MEMERR;
1289
1290 {
1291 INT i, j, jj, begin_row, end_row;
1292 REAL* lu = iludata->luval;
1293 INT* ijlu = iludata->ijlu;
1294 REAL *xval = x->val, *bval = b->val;
1295
1297 fasp_darray_cp(m, bval, zr);
1298 fasp_blas_dcsr_aAxpy(-1.0, A, xval, zr);
1299
1300 // forward sweep: solve unit lower matrix equation L*zz=zr
1301 zz[0] = zr[0];
1302 for (i = 1; i < m; ++i) {
1303 begin_row = ijlu[i];
1304 end_row = ijlu[i + 1];
1305 for (j = begin_row; j < end_row; ++j) {
1306 jj = ijlu[j];
1307 if (jj < i)
1308 zr[i] -= lu[j] * zz[jj];
1309 else
1310 break;
1311 }
1312 zz[i] = zr[i];
1313 }
1314
1315 // backward sweep: solve upper matrix equation U*z=zz
1316 z[m - 1] = zz[m - 1] * lu[m - 1];
1317 for (i = m - 2; i >= 0; --i) {
1318 begin_row = ijlu[i];
1319 end_row = ijlu[i + 1] - 1;
1320 for (j = end_row; j >= begin_row; --j) {
1321 jj = ijlu[j];
1322 if (jj > i)
1323 zz[i] -= lu[j] * z[jj];
1324 else
1325 break;
1326 }
1327 z[i] = zz[i] * lu[i];
1328 }
1329
1330 fasp_blas_darray_axpy(m, 1, z, xval);
1331 }
1332
1333 return;
1334
1335MEMERR:
1336 printf("### ERROR: ILU needs %d memory, only %d available! [%s:%d]\n", memneed,
1337 iludata->nwork, __FILE__, __LINE__);
1338 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1339}
1340
1362void fasp_smoother_dcsr_kaczmarz(dvector* u,
1363 const INT i_1,
1364 const INT i_n,
1365 const INT s,
1366 dCSRmat* A,
1367 dvector* b,
1368 INT L,
1369 const REAL w)
1370{
1371 const INT * ia = A->IA, *ja = A->JA;
1372 const REAL *aval = A->val, *bval = b->val;
1373 REAL* uval = u->val;
1374
1375 // local variables
1376 INT i, j, k, begin_row, end_row;
1377 REAL temp1, temp2, alpha;
1378
1379#ifdef _OPENMP
1380 const INT N = ABS(i_n - i_1) + 1;
1381 INT myid, mybegin, myend;
1382 INT nthreads = fasp_get_num_threads();
1383#endif
1384
1385 if (s > 0) {
1386
1387 while (L--) {
1388#ifdef _OPENMP
1389 if (N > OPENMP_HOLDS) {
1390#pragma omp parallel for private(myid, mybegin, myend, i, temp1, temp2, begin_row, \
1391 end_row, k, alpha, j)
1392 for (myid = 0; myid < nthreads; myid++) {
1393 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1394 mybegin += i_1, myend += i_1;
1395 for (i = mybegin; i < myend; i += s) {
1396 temp1 = 0;
1397 temp2 = 0;
1398 begin_row = ia[i], end_row = ia[i + 1];
1399 for (k = begin_row; k < end_row; k++) {
1400 j = ja[k];
1401 temp1 += aval[k] * aval[k];
1402 temp2 += aval[k] * uval[j];
1403 } // end for k
1404 }
1405 alpha = (bval[i] - temp2) / temp1;
1406 for (k = begin_row; k < end_row; ++k) {
1407 j = ja[k];
1408 uval[j] += w * alpha * aval[k];
1409 } // end for k
1410 } // end for i
1411 } else {
1412#endif
1413 for (i = i_1; i <= i_n; i += s) {
1414 temp1 = 0;
1415 temp2 = 0;
1416 begin_row = ia[i];
1417 end_row = ia[i + 1];
1418 for (k = begin_row; k < end_row; ++k) {
1419 j = ja[k];
1420 temp1 += aval[k] * aval[k];
1421 temp2 += aval[k] * uval[j];
1422 } // end for k
1423 alpha = (bval[i] - temp2) / temp1;
1424 for (k = begin_row; k < end_row; ++k) {
1425 j = ja[k];
1426 uval[j] += w * alpha * aval[k];
1427 } // end for k
1428 } // end for i
1429#ifdef _OPENMP
1430 }
1431#endif
1432 } // end while
1433
1434 } // if s
1435
1436 else {
1437 while (L--) {
1438#ifdef _OPENMP
1439 if (N > OPENMP_HOLDS) {
1440#pragma omp parallel for private(myid, mybegin, myend, i, temp1, temp2, begin_row, \
1441 end_row, k, alpha, j)
1442 for (myid = 0; myid < nthreads; myid++) {
1443 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1444 mybegin = i_1 - mybegin, myend = i_1 - myend;
1445 for (i = mybegin; i > myend; i += s) {
1446 temp1 = 0;
1447 temp2 = 0;
1448 begin_row = ia[i], end_row = ia[i + 1];
1449 for (k = begin_row; k < end_row; ++k) {
1450 j = ja[k];
1451 temp1 += aval[k] * aval[k];
1452 temp2 += aval[k] * uval[j];
1453 } // end for k
1454 alpha = (bval[i] - temp2) / temp1;
1455 for (k = begin_row; k < end_row; ++k) {
1456 j = ja[k];
1457 uval[j] += w * alpha * aval[k];
1458 } // end for k
1459 } // end for i
1460 }
1461 } else {
1462#endif
1463 for (i = i_1; i >= i_n; i += s) {
1464 temp1 = 0;
1465 temp2 = 0;
1466 begin_row = ia[i];
1467 end_row = ia[i + 1];
1468 for (k = begin_row; k < end_row; ++k) {
1469 j = ja[k];
1470 temp1 += aval[k] * aval[k];
1471 temp2 += aval[k] * uval[j];
1472 } // end for k
1473 alpha = (bval[i] - temp2) / temp1;
1474 for (k = begin_row; k < end_row; ++k) {
1475 j = ja[k];
1476 uval[j] += w * alpha * aval[k];
1477 } // end for k
1478 } // end for i
1479#ifdef _OPENMP
1480 }
1481#endif
1482 } // end while
1483
1484 } // end if
1485
1486 return;
1487}
1488
1508void fasp_smoother_dcsr_L1diag(dvector* u,
1509 const INT i_1,
1510 const INT i_n,
1511 const INT s,
1512 dCSRmat* A,
1513 dvector* b,
1514 INT L)
1515{
1516 const INT N = ABS(i_n - i_1) + 1;
1517 const INT * ia = A->IA, *ja = A->JA;
1518 const REAL *aval = A->val, *bval = b->val;
1519 REAL* uval = u->val;
1520
1521 // local variables
1522 INT i, j, k, begin_row, end_row;
1523
1524#ifdef _OPENMP
1525 INT myid, mybegin, myend;
1526 INT nthreads = fasp_get_num_threads();
1527#endif
1528
1529 // Checks should be outside of for; t,d can be allocated before calling!!!
1530 // --Chensong
1531 REAL* t = (REAL*)fasp_mem_calloc(N, sizeof(REAL));
1532 REAL* d = (REAL*)fasp_mem_calloc(N, sizeof(REAL));
1533
1534 while (L--) {
1535 if (s > 0) {
1536#ifdef _OPENMP
1537 if (N > OPENMP_HOLDS) {
1538#pragma omp parallel for private(myid, mybegin, myend, i, begin_row, end_row, k, j)
1539 for (myid = 0; myid < nthreads; myid++) {
1540 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1541 mybegin += i_1, myend += i_1;
1542 for (i = mybegin; i < myend; i += s) {
1543 t[i] = bval[i];
1544 d[i] = 0.0;
1545 begin_row = ia[i], end_row = ia[i + 1];
1546 for (k = begin_row; k < end_row; k++) {
1547 j = ja[k];
1548 t[i] -= aval[k] * uval[j];
1549 d[i] += ABS(aval[k]);
1550 }
1551 }
1552 }
1553#pragma omp parallel for private(i)
1554 for (i = i_1; i <= i_n; i += s) {
1555 if (ABS(d[i]) > SMALLREAL) u->val[i] += t[i] / d[i];
1556 }
1557 } else {
1558#endif
1559 for (i = i_1; i <= i_n; i += s) {
1560 t[i] = bval[i];
1561 d[i] = 0.0;
1562 begin_row = ia[i];
1563 end_row = ia[i + 1];
1564 for (k = begin_row; k < end_row; ++k) {
1565 j = ja[k];
1566 t[i] -= aval[k] * uval[j];
1567 d[i] += ABS(aval[k]);
1568 }
1569 }
1570
1571 for (i = i_1; i <= i_n; i += s) {
1572 if (ABS(d[i]) > SMALLREAL) u->val[i] += t[i] / d[i];
1573 }
1574#ifdef _OPENMP
1575 }
1576#endif
1577 } else {
1578#ifdef _OPENMP
1579 if (N > OPENMP_HOLDS) {
1580#pragma omp parallel for private(myid, mybegin, myend, i, k, j, begin_row, end_row)
1581 for (myid = 0; myid < nthreads; myid++) {
1582 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1583 mybegin = i_1 - mybegin, myend = i_1 - myend;
1584 for (i = mybegin; i > myend; i += s) {
1585 t[i] = bval[i];
1586 d[i] = 0.0;
1587 begin_row = ia[i], end_row = ia[i + 1];
1588 for (k = begin_row; k < end_row; k++) {
1589 j = ja[k];
1590 t[i] -= aval[k] * uval[j];
1591 d[i] += ABS(aval[k]);
1592 }
1593 }
1594 }
1595#pragma omp parallel for private(i)
1596 for (i = i_1; i >= i_n; i += s) {
1597 if (ABS(d[i]) > SMALLREAL) u->val[i] += t[i] / d[i];
1598 }
1599 } else {
1600#endif
1601 for (i = i_1; i >= i_n; i += s) {
1602 t[i] = bval[i];
1603 d[i] = 0.0;
1604 begin_row = ia[i];
1605 end_row = ia[i + 1];
1606 for (k = begin_row; k < end_row; ++k) {
1607 j = ja[k];
1608 t[i] -= aval[k] * uval[j];
1609 d[i] += ABS(aval[k]);
1610 }
1611 }
1612
1613 for (i = i_1; i >= i_n; i += s) {
1614 if (ABS(d[i]) > SMALLREAL) u->val[i] += t[i] / d[i];
1615 }
1616#ifdef _OPENMP
1617 }
1618#endif
1619 }
1620
1621 } // end while
1622
1623 fasp_mem_free(t);
1624 t = NULL;
1625 fasp_mem_free(d);
1626 d = NULL;
1627
1628 return;
1629}
1630
1631#if 0
1652static dCSRmat form_contractor (dCSRmat *A,
1653 const INT smoother,
1654 const INT steps,
1655 const INT ndeg,
1656 const REAL relax,
1657 const REAL dtol)
1658{
1659 const INT n=A->row;
1660 INT i;
1661
1662 REAL *work = (REAL *)fasp_mem_calloc(2*n,sizeof(REAL));
1663
1664 dvector b, x;
1665 b.row=x.row=n;
1666 b.val=work; x.val=work+n;
1667
1668 INT *index = (INT *)fasp_mem_calloc(n,sizeof(INT));
1669
1670 for (i=0; i<n; ++i) index[i]=i;
1671
1672 dCSRmat B = fasp_dcsr_create(n, n, n*n); // too much memory required, need to change!!
1673
1674 dCSRmat C, D;
1675
1676 for (i=0; i<n; ++i){
1677
1678 // get i-th column
1679 fasp_dcsr_getcol(i, A, b.val);
1680
1681 // set x =0.0
1682 fasp_dvec_set(n, &x, 0.0);
1683
1684 // smooth
1685 switch (smoother) {
1686 case GS:
1687 fasp_smoother_dcsr_gs(&x, 0, n-1, 1, A, &b, steps);
1688 break;
1689 case POLY:
1690 fasp_smoother_dcsr_poly(A, &b, &x, n, ndeg, steps);
1691 break;
1692 case JACOBI:
1693 fasp_smoother_dcsr_jacobi(&x, 0, n-1, 1, A, &b, steps);
1694 break;
1695 case SGS:
1696 fasp_smoother_dcsr_sgs(&x, A, &b, steps);
1697 break;
1698 case SOR:
1699 fasp_smoother_dcsr_sor(&x, 0, n-1, 1, A, &b, steps, relax);
1700 break;
1701 case SSOR:
1702 fasp_smoother_dcsr_sor(&x, 0, n-1, 1, A, &b, steps, relax);
1703 fasp_smoother_dcsr_sor(&x, n-1, 0,-1, A, &b, steps, relax);
1704 break;
1705 case GSOR:
1706 fasp_smoother_dcsr_gs(&x, 0, n-1, 1, A, &b, steps);
1707 fasp_smoother_dcsr_sor(&x, n-1, 0, -1, A, &b, steps, relax);
1708 break;
1709 case SGSOR:
1710 fasp_smoother_dcsr_gs(&x, 0, n-1, 1, A, &b, steps);
1711 fasp_smoother_dcsr_gs(&x, n-1, 0,-1, A, &b, steps);
1712 fasp_smoother_dcsr_sor(&x, 0, n-1, 1, A, &b, steps, relax);
1713 fasp_smoother_dcsr_sor(&x, n-1, 0,-1, A, &b, steps, relax);
1714 break;
1715 default:
1716 printf("### ERROR: Unknown smoother type! [%s:%d]\n",
1717 __FILE__, __LINE__);
1718 fasp_chkerr(ERROR_INPUT_PAR, __FUNCTION__);
1719 }
1720
1721 // store to B
1722 B.IA[i] = i*n;
1723 memcpy(&(B.JA[i*n]), index, n*sizeof(INT));
1724 memcpy(&(B.val[i*n]), x.val, x.row*sizeof(REAL));
1725
1726 }
1727
1728 B.IA[n] = n*n;
1729
1730 // drop small entries
1731 compress_dCSRmat(&B, &D, dtol);
1732
1733 // get contractor
1734 fasp_dcsr_trans(&D, &C);
1735
1736 // clean up
1737 fasp_mem_free(work); work = NULL;
1738 fasp_dcsr_free(&B);
1739 fasp_dcsr_free(&D);
1740
1741 return C;
1742}
1743#endif
1744
1745/*---------------------------------*/
1746/*-- End of File --*/
1747/*---------------------------------*/
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_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_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_dcsr_create
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:47
fasp_dcsr_free
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:184
fasp_dcsr_getcol
void fasp_dcsr_getcol(const INT n, const dCSRmat *A, REAL *col)
Get the n-th column of a CSR matrix A.
Definition: BlaSparseCSR.c:602
fasp_dcsr_trans
INT fasp_dcsr_trans(const dCSRmat *A, dCSRmat *AT)
Find transpose of dCSRmat matrix A.
Definition: BlaSparseCSR.c:952
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_smoother_dcsr_sor_cf
void fasp_smoother_dcsr_sor_cf(dvector *u, dCSRmat *A, dvector *b, INT L, const REAL w, INT *mark, const INT order)
SOR smoother with C/F ordering for Au=b.
Definition: ItrSmootherCSR.c:1061
fasp_smoother_dcsr_sor
void fasp_smoother_dcsr_sor(dvector *u, const INT i_1, const INT i_n, const INT s, dCSRmat *A, dvector *b, INT L, const REAL w)
SOR method as a smoother.
Definition: ItrSmootherCSR.c:931
fasp_smoother_dcsr_kaczmarz
void fasp_smoother_dcsr_kaczmarz(dvector *u, const INT i_1, const INT i_n, const INT s, dCSRmat *A, dvector *b, INT L, const REAL w)
Kaczmarz method as a smoother.
Definition: ItrSmootherCSR.c:1362
fasp_smoother_dcsr_jacobi_ff
void fasp_smoother_dcsr_jacobi_ff(dvector *x, dCSRmat *A, dvector *b, const INT nsweeps, INT *ordering, const REAL relax)
Weighted Jacobi method as a smoother only for the fine points.
Definition: ItrSmootherCSR.c:34
fasp_smoother_dcsr_sgs
void fasp_smoother_dcsr_sgs(dvector *u, dCSRmat *A, dvector *b, INT L)
Symmetric Gauss-Seidel method as a smoother.
Definition: ItrSmootherCSR.c:807
fasp_smoother_dcsr_gs
void fasp_smoother_dcsr_gs(dvector *u, const INT i_1, const INT i_n, const INT s, dCSRmat *A, dvector *b, INT L)
Gauss-Seidel method as a smoother.
Definition: ItrSmootherCSR.c:250
fasp_smoother_dcsr_gs_ff
void fasp_smoother_dcsr_gs_ff(dvector *u, dCSRmat *A, dvector *b, INT L, INT *mark)
Gauss-Seidel smoother with on F-points only for Au=b.
Definition: ItrSmootherCSR.c:712
fasp_smoother_dcsr_ilu
void fasp_smoother_dcsr_ilu(dCSRmat *A, dvector *b, dvector *x, void *data)
ILU method as a smoother.
Definition: ItrSmootherCSR.c:1279
fasp_smoother_dcsr_L1diag
void fasp_smoother_dcsr_L1diag(dvector *u, const INT i_1, const INT i_n, const INT s, dCSRmat *A, dvector *b, INT L)
Diagonal scaling (using L1 norm) as a smoother.
Definition: ItrSmootherCSR.c:1508
fasp_smoother_dcsr_jacobi
void fasp_smoother_dcsr_jacobi(dvector *u, const INT i_1, const INT i_n, const INT s, dCSRmat *A, dvector *b, INT L, const REAL w)
Weighted Jacobi method as a smoother.
Definition: ItrSmootherCSR.c:97
fasp_smoother_dcsr_gs_cf
void fasp_smoother_dcsr_gs_cf(dvector *u, dCSRmat *A, dvector *b, INT L, INT *mark, const INT order)
Gauss-Seidel smoother with C/F ordering for Au=b.
Definition: ItrSmootherCSR.c:431
fasp_smoother_dcsr_poly
void fasp_smoother_dcsr_poly(dCSRmat *Amat, dvector *brhs, dvector *usol, INT n, INT ndeg, INT L)
poly approx to A^{-1} as MG smoother
Definition: ItrSmootherCSRpoly.c:67
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
ABS
#define ABS(a)
Definition: fasp.h:84
INT
#define INT
Definition: fasp.h:72
FPFIRST
#define FPFIRST
Definition: fasp_const.h:248
FGPT
#define FGPT
Definition: fasp_const.h:233
OPENMP_HOLDS
#define OPENMP_HOLDS
Definition: fasp_const.h:269
ERROR_ALLOC_MEM
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
ERROR_INPUT_PAR
#define ERROR_INPUT_PAR
Definition: fasp_const.h:24
SMALLREAL
#define SMALLREAL
Definition: fasp_const.h:256
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::nwork
INT nwork
work space size
Definition: fasp.h:678
ILU_data::work
REAL * work
work space
Definition: fasp.h:681
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::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
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