Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
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
51 const INT i_1,
52 const INT i_n,
53 const INT s,
54 dCSRmat *A,
55 dvector *b,
56 INT L,
57 const REAL w)
58{
59 const INT N = ABS(i_n - i_1) + 1;
60 const INT *ia = A->IA, *ja = A->JA;
61 const REAL *aval = A->val, *bval = b->val;
62 REAL *uval = u->val;
63
64 // local variables
65 INT i,j,k,begin_row,end_row;
66
67 // OpenMP variables
68#ifdef _OPENMP
69 INT myid, mybegin, myend;
70 INT nthreads = fasp_get_num_threads();
71#endif
72
73 REAL *t = (REAL *)fasp_mem_calloc(N,sizeof(REAL));
74 REAL *d = (REAL *)fasp_mem_calloc(N,sizeof(REAL));
75
76 while (L--) {
77
78 if (s>0) {
79#ifdef _OPENMP
80 if (N > OPENMP_HOLDS) {
81#pragma omp parallel for private(myid, mybegin, myend, begin_row, end_row, i, k, j)
82 for (myid=0; myid<nthreads; ++myid) {
83 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
84 mybegin += i_1; myend += i_1;
85 for (i=mybegin; i<myend; i+=s) {
86 t[i]=bval[i];
87 begin_row=ia[i],end_row=ia[i+1];
88 for (k=begin_row; k<end_row; ++k) {
89 j=ja[k];
90 if (i!=j) t[i]-=aval[k]*uval[j];
91 else d[i]=aval[k];
92 }
93 }
94 }
95 }
96 else {
97#endif
98 for (i=i_1;i<=i_n;i+=s) {
99 t[i]=bval[i];
100 begin_row=ia[i]; end_row=ia[i+1];
101 for (k=begin_row;k<end_row;++k) {
102 j=ja[k];
103 if (i!=j) t[i]-=aval[k]*uval[j];
104 else d[i]=aval[k];
105 }
106 }
107#ifdef _OPENMP
108 }
109#endif
110
111#ifdef _OPENMP
112#pragma omp parallel for private (i)
113#endif
114 for (i=i_1;i<=i_n;i+=s) {
115 if (ABS(d[i])>SMALLREAL) uval[i]=(1-w)*uval[i]+ w*t[i]/d[i];
116 }
117
118 }
119
120 else {
121
122#ifdef _OPENMP
123 if (N > OPENMP_HOLDS) {
124#pragma omp parallel for private(myid, mybegin, myend, i, begin_row, end_row, k, j)
125 for (myid=0; myid<nthreads; myid++) {
126 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
127 mybegin = i_1-mybegin; myend = i_1-myend;
128 for (i=mybegin; i>myend; i+=s) {
129 t[i]=bval[i];
130 begin_row=ia[i],end_row=ia[i+1];
131 for (k=begin_row; k<end_row; ++k) {
132 j=ja[k];
133 if (i!=j) t[i]-=aval[k]*uval[j];
134 else d[i]=aval[k];
135 }
136 }
137 }
138 }
139 else {
140#endif
141 for (i=i_1;i>=i_n;i+=s) {
142 t[i]=bval[i];
143 begin_row=ia[i]; end_row=ia[i+1];
144 for (k=begin_row;k<end_row;++k) {
145 j=ja[k];
146 if (i!=j) t[i]-=aval[k]*uval[j];
147 else d[i]=aval[k];
148 }
149 }
150#ifdef _OPENMP
151 }
152#endif
153
154#ifdef _OPENMP
155#pragma omp parallel for private(i)
156#endif
157 for (i=i_1;i>=i_n;i+=s) {
158 if (ABS(d[i])>SMALLREAL) uval[i]=(1-w)*uval[i]+ w*t[i]/d[i];
159 }
160
161 }
162
163 } // end while
164
165 fasp_mem_free(t); t = NULL;
166 fasp_mem_free(d); d = NULL;
167
168 return;
169}
170
191 const INT i_1,
192 const INT i_n,
193 const INT s,
194 dCSRmat *A,
195 dvector *b,
196 INT L)
197{
198 const INT *ia = A->IA, *ja = A->JA;
199 const REAL *aval = A->val, *bval = b->val;
200 REAL *uval = u->val;
201
202 // local variables
203 INT i,j,k,begin_row,end_row;
204 REAL t,d=0.0;
205
206#ifdef _OPENMP
207 const INT N = ABS(i_n - i_1)+1;
208 INT myid, mybegin, myend;
209 INT nthreads = fasp_get_num_threads();
210#endif
211
212 if (s > 0) {
213
214 while (L--) {
215#ifdef _OPENMP
216 if (N >OPENMP_HOLDS) {
217#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, d, k, j)
218 for (myid=0; myid<nthreads; myid++) {
219 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
220 mybegin += i_1, myend += i_1;
221 for (i=mybegin; i<myend; i+=s) {
222 t = bval[i];
223 begin_row=ia[i],end_row=ia[i+1];
224#if DIAGONAL_PREF // diagonal first
225 d=aval[begin_row];
226 if (ABS(d)>SMALLREAL) {
227 for (k=begin_row+1;k<end_row;++k) {
228 j=ja[k];
229 t-=aval[k]*uval[j];
230 }
231 uval[i]=t/d;
232 }
233#else // general order
234 for (k=begin_row;k<end_row;++k) {
235 j=ja[k];
236 if (i!=j)
237 t-=aval[k]*uval[j];
238 else if (ABS(aval[k])>SMALLREAL) d=1.e+0/aval[k];
239 }
240 uval[i]=t*d;
241#endif // end DIAGONAL_PREF
242 } // end for i
243 }
244
245 }
246
247 else {
248#endif
249 for (i=i_1;i<=i_n;i+=s) {
250 t = bval[i];
251 begin_row=ia[i]; end_row=ia[i+1];
252
253#if DIAGONAL_PREF // diagonal first
254 d=aval[begin_row];
255 if (ABS(d)>SMALLREAL) {
256 for (k=begin_row+1;k<end_row;++k) {
257 j=ja[k];
258 t-=aval[k]*uval[j];
259 }
260 uval[i]=t/d;
261 }
262#else // general order
263 for (k=begin_row;k<end_row;++k) {
264 j=ja[k];
265 if (i!=j)
266 t-=aval[k]*uval[j];
267 else if (ABS(aval[k])>SMALLREAL) d=1.e+0/aval[k];
268 }
269 uval[i]=t*d;
270#endif
271 } // end for i
272#ifdef _OPENMP
273 }
274#endif
275 } // end while
276
277 } // if s
278 else {
279
280 while (L--) {
281#ifdef _OPENMP
282 if (N > OPENMP_HOLDS) {
283#pragma omp parallel for private(myid, mybegin, myend, i, begin_row, end_row, d, k, j, t)
284 for (myid=0; myid<nthreads; myid++) {
285 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
286 mybegin = i_1 - mybegin; myend = i_1 - myend;
287 for (i=mybegin; i>myend; i+=s) {
288 t=bval[i];
289 begin_row=ia[i],end_row=ia[i+1];
290#if DIAGONAL_PREF // diagonal first
291 d=aval[begin_row];
292 if (ABS(d)>SMALLREAL) {
293 for (k=begin_row+1;k<end_row;++k) {
294 j=ja[k];
295 t-=aval[k]*uval[j];
296 }
297 uval[i]=t/d;
298 }
299#else // general order
300 for (k=begin_row;k<end_row;++k) {
301 j=ja[k];
302 if (i!=j)
303 t-=aval[k]*uval[j];
304 else if (ABS(aval[k])>SMALLREAL) d=1.0/aval[k];
305 }
306 uval[i]=t*d;
307#endif
308 } // end for i
309 }
310 }
311 else {
312#endif
313 for (i=i_1;i>=i_n;i+=s) {
314 t=bval[i];
315 begin_row=ia[i]; end_row=ia[i+1];
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) d=1.0/aval[k];
331 }
332 uval[i]=t*d;
333#endif
334 } // end for i
335#ifdef _OPENMP
336 }
337#endif
338 } // end while
339
340 } // end if
341
342 return;
343}
344
364 dCSRmat *A,
365 dvector *b,
366 INT L,
367 INT *mark,
368 const INT order)
369{
370 const INT nrow = b->row; // number of rows
371 const INT *ia = A->IA, *ja = A->JA;
372 const REAL *aval = A->val, *bval = b->val;
373 REAL *uval = u->val;
374
375 INT i,j,k,begin_row,end_row;
376 REAL t,d=0.0;
377
378#ifdef _OPENMP
379 INT myid,mybegin,myend;
380 INT nthreads = fasp_get_num_threads();
381#endif
382
383 // F-point first, C-point second
384 if (order == FPFIRST) {
385
386 while (L--) {
387
388#ifdef _OPENMP
389 if (nrow > OPENMP_HOLDS) {
390#pragma omp parallel for private(myid, mybegin, myend, i,t,begin_row,end_row,k,j,d)
391 for (myid = 0; myid < nthreads; myid ++) {
392 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
393 for (i=mybegin; i<myend; i++) {
394 if (mark[i] != 1) {
395 t = bval[i];
396 begin_row = ia[i], end_row = ia[i+1];
397#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
398 d = aval[begin_row];
399 for (k = begin_row+1; k < end_row; k ++) {
400 j = ja[k];
401 t -= aval[k]*uval[j];
402 } // end for k
403#else
404 for (k = begin_row; k < end_row; k ++) {
405 j = ja[k];
406 if (i!=j) t -= aval[k]*uval[j];
407 else d = aval[k];
408 } // end for k
409#endif // end if DIAG_PREF
410 if (ABS(d) > SMALLREAL) uval[i] = t/d;
411 }
412 } // end for i
413 }
414 }
415 else {
416#endif
417 for (i = 0; i < nrow; i ++) {
418 if (mark[i] != 1) {
419 t = bval[i];
420 begin_row = ia[i]; end_row = ia[i+1];
421#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
422 d = aval[begin_row];
423 for (k = begin_row+1; k < end_row; k ++) {
424 j = ja[k];
425 t -= aval[k]*uval[j];
426 } // end for k
427#else
428 for (k = begin_row; k < end_row; k ++) {
429 j = ja[k];
430 if (i!=j) t -= aval[k]*uval[j];
431 else d = aval[k];
432 } // end for k
433#endif // end if DIAG_PREF
434 if (ABS(d) > SMALLREAL) uval[i] = t/d;
435 }
436 } // end for i
437#ifdef _OPENMP
438 }
439#endif
440
441#ifdef _OPENMP
442 if (nrow > OPENMP_HOLDS) {
443#pragma omp parallel for private(myid,mybegin,myend,i,t,begin_row,end_row,k,j,d)
444 for (myid = 0; myid < nthreads; myid ++) {
445 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
446 for (i=mybegin; i<myend; i++) {
447 if (mark[i] == 1) {
448 t = bval[i];
449 begin_row = ia[i], end_row = ia[i+1];
450#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
451 d = aval[begin_row];
452 for (k = begin_row+1; k < end_row; k ++) {
453 j = ja[k];
454 t -= aval[k]*uval[j];
455 } // end for k
456#else
457 for (k = begin_row; k < end_row; k ++) {
458 j = ja[k];
459 if (i!=j) t -= aval[k]*uval[j];
460 else d = aval[k];
461 } // end for k
462#endif // end if DIAG_PREF
463 if (ABS(d) > SMALLREAL) uval[i] = t/d;
464 }
465 } // end for i
466 }
467 }
468 else {
469#endif
470 for (i = 0; i < nrow; i ++) {
471 if (mark[i] == 1) {
472 t = bval[i];
473 begin_row = ia[i]; end_row = ia[i+1];
474#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
475 d = aval[begin_row];
476 for (k = begin_row+1; k < end_row; k ++) {
477 j = ja[k];
478 t -= aval[k]*uval[j];
479 } // end for k
480#else
481 for (k = begin_row; k < end_row; k ++) {
482 j = ja[k];
483 if (i!=j) t -= aval[k]*uval[j];
484 else d = aval[k];
485 } // end for k
486#endif // end if DIAG_PREF
487 if (ABS(d) > SMALLREAL) uval[i] = t/d;
488 }
489 } // end for i
490#ifdef _OPENMP
491 }
492#endif
493 } // end while
494
495 }
496
497 // C-point first, F-point second
498 else {
499
500 while (L--) {
501#ifdef _OPENMP
502 if (nrow > OPENMP_HOLDS) {
503#pragma omp parallel for private(myid,mybegin,myend,t,i,begin_row,end_row,k,j,d)
504 for (myid = 0; myid < nthreads; myid ++) {
505 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
506 for (i=mybegin; i<myend; i++) {
507 if (mark[i] == 1) {
508 t = bval[i];
509 begin_row = ia[i],end_row = ia[i+1];
510#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
511 d = aval[begin_row];
512 for (k = begin_row+1; k < end_row; k ++) {
513 j = ja[k];
514 t -= aval[k]*uval[j];
515 } // end for k
516#else
517 for (k = begin_row; k < end_row; k ++) {
518 j = ja[k];
519 if (i!=j) t -= aval[k]*uval[j];
520 else d = aval[k];
521 } // end for k
522#endif // end if DIAG_PREF
523 if (ABS(d) > SMALLREAL) uval[i] = t/d;
524 }
525 } // end for i
526 }
527 }
528 else {
529#endif
530 for (i = 0; i < nrow; i ++) {
531 if (mark[i] == 1) {
532 t = bval[i];
533 begin_row = ia[i]; end_row = ia[i+1];
534#if DIAGONAL_PREF // Added by Chensong on 09/22/2012
535 d = aval[begin_row];
536 for (k = begin_row+1; k < end_row; k ++) {
537 j = ja[k];
538 t -= aval[k]*uval[j];
539 } // end for k
540#else
541 for (k = begin_row; k < end_row; k ++) {
542 j = ja[k];
543 if (i!=j) t -= aval[k]*uval[j];
544 else d = aval[k];
545 } // end for k
546#endif // end if DIAG_PREF
547 if (ABS(d) > SMALLREAL) uval[i] = t/d;
548 }
549 } // end for i
550#ifdef _OPENMP
551 }
552#endif
553
554#ifdef _OPENMP
555 if (nrow > OPENMP_HOLDS) {
556#pragma omp parallel for private(myid, mybegin, myend, i,t,begin_row,end_row,k,j,d)
557 for (myid = 0; myid < nthreads; myid ++) {
558 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
559 for (i=mybegin; i<myend; i++) {
560 if (mark[i] != 1) {
561 t = bval[i];
562 begin_row = ia[i],end_row = ia[i+1];
563#if DIAGONAL_PREF // Added by Chensong on 01/17/2013
564 d = aval[begin_row];
565 for (k = begin_row+1; k < end_row; k ++) {
566 j = ja[k];
567 t -= aval[k]*uval[j];
568 } // end for k
569#else
570 for (k = begin_row; k < end_row; k ++) {
571 j = ja[k];
572 if (i!=j) t -= aval[k]*uval[j];
573 else d = aval[k];
574 } // end for k
575#endif // end if DIAG_PREF
576 if (ABS(d) > SMALLREAL) uval[i] = t/d;
577 }
578 } // end for i
579 }
580 }
581 else {
582#endif
583 for (i = 0; i < nrow; i ++) {
584 if (mark[i] != 1) {
585 t = bval[i];
586 begin_row = ia[i]; end_row = ia[i+1];
587#if DIAGONAL_PREF // Added by Chensong on 09/22/2012
588 d = aval[begin_row];
589 for (k = begin_row+1; k < end_row; k ++) {
590 j = ja[k];
591 t -= aval[k]*uval[j];
592 } // end for k
593#else
594 for (k = begin_row; k < end_row; k ++) {
595 j = ja[k];
596 if (i!=j) t -= aval[k]*uval[j];
597 else 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#ifdef _OPENMP
604 }
605#endif
606 } // end while
607
608 } // end if order
609
610 return;
611}
612
629 dCSRmat *A,
630 dvector *b,
631 INT L)
632{
633 const INT nm1=b->row-1;
634 const INT *ia=A->IA,*ja=A->JA;
635 const REAL *aval=A->val,*bval=b->val;
636 REAL *uval=u->val;
637
638 // local variables
639 INT i,j,k,begin_row,end_row;
640 REAL t,d=0;
641
642#ifdef _OPENMP
643 INT myid, mybegin, myend, up;
644 INT nthreads = fasp_get_num_threads();
645#endif
646
647 while (L--) {
648 // forward sweep
649#ifdef _OPENMP
650 up = nm1 + 1;
651 if (up > OPENMP_HOLDS) {
652#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, j, k, d)
653 for (myid=0; myid<nthreads; myid++) {
654 fasp_get_start_end(myid, nthreads, up, &mybegin, &myend);
655 for (i=mybegin; i<myend; i++) {
656 t=bval[i];
657 begin_row=ia[i], end_row=ia[i+1];
658 for (k=begin_row;k<end_row;++k) {
659 j=ja[k];
660 if (i!=j) t-=aval[k]*uval[j];
661 else d=aval[k];
662 } // end for k
663 if (ABS(d)>SMALLREAL) uval[i]=t/d;
664 } // end for i
665 }
666 }
667 else {
668#endif
669 for (i=0;i<=nm1;++i) {
670 t=bval[i];
671 begin_row=ia[i]; end_row=ia[i+1];
672 for (k=begin_row;k<end_row;++k) {
673 j=ja[k];
674 if (i!=j) t-=aval[k]*uval[j];
675 else d=aval[k];
676 } // end for k
677 if (ABS(d)>SMALLREAL) uval[i]=t/d;
678 } // end for i
679#ifdef _OPENMP
680 }
681#endif
682
683 // backward sweep
684#ifdef _OPENMP
685 up = nm1;
686 if (up > OPENMP_HOLDS) {
687#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, d)
688 for (myid=0; myid<nthreads; myid++) {
689 fasp_get_start_end(myid, nthreads, up, &mybegin, &myend);
690 mybegin = nm1-1-mybegin; myend = nm1-1-myend;
691 for (i=mybegin; i>myend; i--) {
692 t=bval[i];
693 begin_row=ia[i], end_row=ia[i+1];
694 for (k=begin_row; k<end_row; k++) {
695 j=ja[k];
696 if (i!=j) t-=aval[k]*uval[j];
697 else d=aval[k];
698 } // end for k
699 if (ABS(d)>SMALLREAL) uval[i]=t/d;
700 } // end for i
701 }
702 }
703 else {
704#endif
705 for (i=nm1-1;i>=0;--i) {
706 t=bval[i];
707 begin_row=ia[i]; end_row=ia[i+1];
708 for (k=begin_row;k<end_row;++k) {
709 j=ja[k];
710 if (i!=j) t-=aval[k]*uval[j];
711 else d=aval[k];
712 } // end for k
713 if (ABS(d)>SMALLREAL) uval[i]=t/d;
714 } // end for i
715#ifdef _OPENMP
716 }
717#endif
718 } // end while
719
720 return;
721}
722
745 const INT i_1,
746 const INT i_n,
747 const INT s,
748 dCSRmat *A,
749 dvector *b,
750 INT L,
751 const REAL w)
752{
753 const INT *ia=A->IA,*ja=A->JA;
754 const REAL *aval=A->val,*bval=b->val;
755 REAL *uval=u->val;
756
757 // local variables
758 INT i,j,k,begin_row,end_row;
759 REAL t, d=0;
760
761#ifdef _OPENMP
762 const INT N = ABS(i_n - i_1)+1;
763 INT myid, mybegin, myend;
764 INT nthreads = fasp_get_num_threads();
765#endif
766
767 while (L--) {
768 if (s>0) {
769#ifdef _OPENMP
770 if (N > OPENMP_HOLDS) {
771#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, d)
772 for (myid=0; myid<nthreads; myid++) {
773 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
774 mybegin += i_1, myend += i_1;
775 for (i=mybegin; i<myend; i+=s) {
776 t=bval[i];
777 begin_row=ia[i], end_row=ia[i+1];
778 for (k=begin_row; k<end_row; k++) {
779 j=ja[k];
780 if (i!=j)
781 t-=aval[k]*uval[j];
782 else
783 d=aval[k];
784 }
785 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
786 }
787 }
788
789 }
790 else {
791#endif
792 for (i=i_1; i<=i_n; i+=s) {
793 t=bval[i];
794 begin_row=ia[i]; end_row=ia[i+1];
795 for (k=begin_row; k<end_row; ++k) {
796 j=ja[k];
797 if (i!=j)
798 t-=aval[k]*uval[j];
799 else
800 d=aval[k];
801 }
802 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
803 }
804#ifdef _OPENMP
805 }
806#endif
807 }
808 else {
809#ifdef _OPENMP
810 if (N > OPENMP_HOLDS) {
811#pragma omp parallel for private(myid, mybegin, myend, i, t, begin_row, end_row, k, j, d)
812 for (myid=0; myid<nthreads; myid++) {
813 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
814 mybegin = i_1 - mybegin, myend = i_1 - myend;
815 for (i=mybegin; i>myend; i+=s) {
816 t=bval[i];
817 begin_row=ia[i],end_row=ia[i+1];
818 for (k=begin_row;k<end_row;++k) {
819 j=ja[k];
820 if (i!=j)
821 t-=aval[k]*uval[j];
822 else
823 d=aval[k];
824 }
825 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
826 }
827 }
828 }
829 else {
830#endif
831 for (i=i_1;i>=i_n;i+=s) {
832 t=bval[i];
833 begin_row=ia[i]; end_row=ia[i+1];
834 for (k=begin_row;k<end_row;++k) {
835 j=ja[k];
836 if (i!=j)
837 t-=aval[k]*uval[j];
838 else
839 d=aval[k];
840 }
841 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
842 }
843#ifdef _OPENMP
844 }
845#endif
846 }
847 } // end while
848
849 return;
850}
851
872 dCSRmat *A,
873 dvector *b,
874 INT L,
875 const REAL w,
876 INT *mark,
877 const INT order )
878{
879 const INT nrow = b->row; // number of rows
880 const INT *ia = A->IA, *ja=A->JA;
881 const REAL *aval = A->val,*bval=b->val;
882 REAL *uval = u->val;
883
884 // local variables
885 INT i,j,k,begin_row,end_row;
886 REAL t,d=0.0;
887
888#ifdef _OPENMP
889 INT myid, mybegin, myend;
890 INT nthreads = fasp_get_num_threads();
891#endif
892
893 // F-point first
894 if (order == -1) {
895 while (L--) {
896#ifdef _OPENMP
897 if (nrow > OPENMP_HOLDS) {
898#pragma omp parallel for private (myid, mybegin, myend, i, t, begin_row, end_row, k, j, d)
899 for (myid = 0; myid < nthreads; myid++) {
900 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
901 for (i = mybegin; i < myend; i ++) {
902 if (mark[i] == 0 || mark[i] == 2) {
903 t = bval[i];
904 begin_row = ia[i], end_row = ia[i+1];
905 for (k = begin_row; k < end_row; k ++) {
906 j = ja[k];
907 if (i!=j) t -= aval[k]*uval[j];
908 else d = aval[k];
909 } // end for k
910 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
911 }
912 }
913 }
914 } // end for i
915 else {
916#endif
917 for (i = 0; i < nrow; i ++) {
918 if (mark[i] == 0 || mark[i] == 2) {
919 t = bval[i];
920 begin_row = ia[i]; end_row = ia[i+1];
921 for (k = begin_row; k < end_row; k ++) {
922 j = ja[k];
923 if (i!=j) t -= aval[k]*uval[j];
924 else d = aval[k];
925 } // end for k
926 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
927 }
928 } // end for i
929#ifdef _OPENMP
930 }
931#endif
932
933#ifdef _OPENMP
934 if (nrow > OPENMP_HOLDS) {
935#pragma omp parallel for private(myid, i, mybegin, myend, t, begin_row, end_row, k, j, d)
936 for (myid = 0; myid < nthreads; myid++) {
937 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
938 for (i = mybegin; i < myend; i++) {
939 if (mark[i] == 1) {
940 t = bval[i];
941 begin_row = ia[i], end_row = ia[i+1];
942 for (k = begin_row; k < end_row; k ++) {
943 j = ja[k];
944 if (i!=j) t -= aval[k]*uval[j];
945 else d = aval[k];
946 } // end for k
947 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
948 }
949 } // end for i
950 }
951 }
952 else {
953#endif
954 for (i = 0; i < nrow; i ++) {
955 if (mark[i] == 1) {
956 t = bval[i];
957 begin_row = ia[i]; end_row = ia[i+1];
958 for (k = begin_row; k < end_row; k ++) {
959 j = ja[k];
960 if (i!=j) t -= aval[k]*uval[j];
961 else d = aval[k];
962 } // end for k
963 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
964 }
965 } // end for i
966#ifdef _OPENMP
967 }
968#endif
969 } // end while
970 }
971 else {
972 while (L--) {
973#ifdef _OPENMP
974 if (nrow > OPENMP_HOLDS) {
975#pragma omp parallel for private(myid, mybegin, myend, i, t, k, j, d, begin_row, end_row)
976 for (myid = 0; myid < nthreads; myid++) {
977 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
978 for (i = mybegin; i < myend; i++) {
979 if (mark[i] == 1) {
980 t = bval[i];
981 begin_row = ia[i], end_row = ia[i+1];
982 for (k = begin_row; k < end_row; k ++) {
983 j = ja[k];
984 if (i!=j) t -= aval[k]*uval[j];
985 else d = aval[k];
986 } // end for k
987 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
988 }
989 } // end for i
990 }
991 }
992 else {
993#endif
994 for (i = 0; i < nrow; i ++) {
995 if (mark[i] == 1) {
996 t = bval[i];
997 begin_row = ia[i]; end_row = ia[i+1];
998 for (k = begin_row; k < end_row; k ++) {
999 j = ja[k];
1000 if (i!=j) t -= aval[k]*uval[j];
1001 else d = aval[k];
1002 } // end for k
1003 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
1004 }
1005 } // end for i
1006#ifdef _OPENMP
1007 }
1008#endif
1009
1010#ifdef _OPENMP
1011 if (nrow > OPENMP_HOLDS) {
1012#pragma omp parallel for private (myid, mybegin, myend, i, t, begin_row, end_row, k, j, d)
1013 for (myid = 0; myid < nthreads; myid++) {
1014 fasp_get_start_end(myid, nthreads, nrow, &mybegin, &myend);
1015 for (i = mybegin; i < myend; i++) {
1016 if (mark[i] != 1) {
1017 t = bval[i];
1018 begin_row = ia[i], end_row = ia[i+1];
1019 for (k = begin_row; k < end_row; k ++) {
1020 j = ja[k];
1021 if (i!=j) t -= aval[k]*uval[j];
1022 else d = aval[k];
1023 } // end for k
1024 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
1025 }
1026 }
1027 }
1028 } // end for i
1029 else {
1030#endif
1031 for (i = 0; i < nrow; i ++) {
1032 if (mark[i] != 1) {
1033 t = bval[i];
1034 begin_row = ia[i]; end_row = ia[i+1];
1035 for (k = begin_row; k < end_row; k ++) {
1036 j = ja[k];
1037 if (i!=j) t -= aval[k]*uval[j];
1038 else d = aval[k];
1039 } // end for k
1040 if (ABS(d)>SMALLREAL) uval[i]=w*(t/d)+(1-w)*uval[i];
1041 }
1042 } // end for i
1043#ifdef _OPENMP
1044 }
1045#endif
1046 } // end while
1047 }
1048
1049 return;
1050}
1051
1066 dvector *b,
1067 dvector *x,
1068 void *data)
1069{
1070 const INT m=A->row, m2=2*m, memneed=3*m;
1071 const ILU_data *iludata=(ILU_data *)data;
1072
1073 REAL *zz = iludata->work;
1074 REAL *zr = iludata->work+m;
1075 REAL *z = iludata->work+m2;
1076
1077 if (iludata->nwork<memneed) goto MEMERR;
1078
1079 {
1080 INT i, j, jj, begin_row, end_row;
1081 REAL *lu = iludata->luval;
1082 INT *ijlu = iludata->ijlu;
1083 REAL *xval = x->val, *bval = b->val;
1084
1086 fasp_darray_cp(m,bval,zr); fasp_blas_dcsr_aAxpy(-1.0,A,xval,zr);
1087
1088 // forward sweep: solve unit lower matrix equation L*zz=zr
1089 zz[0]=zr[0];
1090 for (i=1;i<m;++i) {
1091 begin_row=ijlu[i]; end_row=ijlu[i+1];
1092 for (j=begin_row;j<end_row;++j) {
1093 jj=ijlu[j];
1094 if (jj<i) zr[i]-=lu[j]*zz[jj];
1095 else break;
1096 }
1097 zz[i]=zr[i];
1098 }
1099
1100 // backward sweep: solve upper matrix equation U*z=zz
1101 z[m-1]=zz[m-1]*lu[m-1];
1102 for (i=m-2;i>=0;--i) {
1103 begin_row=ijlu[i]; end_row=ijlu[i+1]-1;
1104 for (j=end_row;j>=begin_row;--j) {
1105 jj=ijlu[j];
1106 if (jj>i) zz[i]-=lu[j]*z[jj];
1107 else break;
1108 }
1109 z[i]=zz[i]*lu[i];
1110 }
1111
1112 fasp_blas_darray_axpy(m,1,z,xval);
1113 }
1114
1115 return;
1116
1117MEMERR:
1118 printf("### ERROR: ILU needs %d memory, only %d available! [%s:%d]\n",
1119 memneed, iludata->nwork, __FILE__, __LINE__);
1120 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1121}
1122
1145 const INT i_1,
1146 const INT i_n,
1147 const INT s,
1148 dCSRmat *A,
1149 dvector *b,
1150 INT L,
1151 const REAL w)
1152{
1153 const INT *ia=A->IA,*ja=A->JA;
1154 const REAL *aval=A->val,*bval=b->val;
1155 REAL *uval=u->val;
1156
1157 // local variables
1158 INT i,j,k,begin_row,end_row;
1159 REAL temp1,temp2,alpha;
1160
1161#ifdef _OPENMP
1162 const INT N = ABS(i_n - i_1)+1;
1163 INT myid, mybegin, myend;
1164 INT nthreads = fasp_get_num_threads();
1165#endif
1166
1167 if (s > 0) {
1168
1169 while (L--) {
1170#ifdef _OPENMP
1171 if (N > OPENMP_HOLDS) {
1172#pragma omp parallel for private(myid, mybegin, myend, i, temp1, temp2, begin_row, end_row, k, alpha, j)
1173 for (myid=0; myid<nthreads; myid++) {
1174 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1175 mybegin += i_1, myend += i_1;
1176 for (i=mybegin; i<myend; i+=s) {
1177 temp1 = 0; temp2 = 0;
1178 begin_row=ia[i], end_row=ia[i+1];
1179 for (k=begin_row; k<end_row; k++) {
1180 j=ja[k];
1181 temp1 += aval[k]*aval[k];
1182 temp2 += aval[k]*uval[j];
1183 } // end for k
1184 }
1185 alpha = (bval[i] - temp2)/temp1;
1186 for (k=begin_row; k<end_row; ++k){
1187 j = ja[k];
1188 uval[j] += w*alpha*aval[k];
1189 }// end for k
1190 } // end for i
1191 }
1192 else {
1193#endif
1194 for (i=i_1;i<=i_n;i+=s) {
1195 temp1 = 0; temp2 = 0;
1196 begin_row=ia[i]; end_row=ia[i+1];
1197 for (k=begin_row;k<end_row;++k) {
1198 j=ja[k];
1199 temp1 += aval[k]*aval[k];
1200 temp2 += aval[k]*uval[j];
1201 } // end for k
1202 alpha = (bval[i] - temp2)/temp1;
1203 for (k=begin_row;k<end_row;++k){
1204 j = ja[k];
1205 uval[j] += w*alpha*aval[k];
1206 }// end for k
1207 } // end for i
1208#ifdef _OPENMP
1209 }
1210#endif
1211 } // end while
1212
1213 } // if s
1214
1215 else {
1216 while (L--) {
1217#ifdef _OPENMP
1218 if (N > OPENMP_HOLDS) {
1219#pragma omp parallel for private(myid, mybegin, myend, i, temp1, temp2, begin_row, end_row, k, alpha, j)
1220 for (myid=0; myid<nthreads; myid++) {
1221 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1222 mybegin = i_1 - mybegin, myend = i_1 - myend;
1223 for (i=mybegin; i>myend; i+=s) {
1224 temp1 = 0; temp2 = 0;
1225 begin_row=ia[i], end_row=ia[i+1];
1226 for (k=begin_row;k<end_row;++k) {
1227 j=ja[k];
1228 temp1 += aval[k]*aval[k];
1229 temp2 += aval[k]*uval[j];
1230 } // end for k
1231 alpha = (bval[i] - temp2)/temp1;
1232 for (k=begin_row;k<end_row;++k){
1233 j = ja[k];
1234 uval[j] += w*alpha*aval[k];
1235 }// end for k
1236 } // end for i
1237 }
1238 }
1239 else {
1240#endif
1241 for (i=i_1;i>=i_n;i+=s) {
1242 temp1 = 0; temp2 = 0;
1243 begin_row=ia[i]; end_row=ia[i+1];
1244 for (k=begin_row;k<end_row;++k) {
1245 j=ja[k];
1246 temp1 += aval[k]*aval[k];
1247 temp2 += aval[k]*uval[j];
1248 } // end for k
1249 alpha = (bval[i] - temp2)/temp1;
1250 for (k=begin_row;k<end_row;++k){
1251 j = ja[k];
1252 uval[j] += w*alpha*aval[k];
1253 }// end for k
1254 } // end for i
1255#ifdef _OPENMP
1256 }
1257#endif
1258 } // end while
1259
1260 } // end if
1261
1262 return;
1263}
1264
1285 const INT i_1,
1286 const INT i_n,
1287 const INT s,
1288 dCSRmat *A,
1289 dvector *b,
1290 INT L)
1291{
1292 const INT N = ABS(i_n - i_1)+1;
1293 const INT *ia=A->IA, *ja=A->JA;
1294 const REAL *aval=A->val,*bval=b->val;
1295 REAL *uval=u->val;
1296
1297 // local variables
1298 INT i,j,k,begin_row,end_row;
1299
1300#ifdef _OPENMP
1301 INT myid, mybegin, myend;
1302 INT nthreads = fasp_get_num_threads();
1303#endif
1304
1305 // Checks should be outside of for; t,d can be allocated before calling!!! --Chensong
1306 REAL *t = (REAL *)fasp_mem_calloc(N,sizeof(REAL));
1307 REAL *d = (REAL *)fasp_mem_calloc(N,sizeof(REAL));
1308
1309 while (L--) {
1310 if (s>0) {
1311#ifdef _OPENMP
1312 if (N > OPENMP_HOLDS) {
1313#pragma omp parallel for private(myid, mybegin, myend, i, begin_row, end_row, k, j)
1314 for (myid=0; myid<nthreads; myid++) {
1315 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1316 mybegin += i_1, myend += i_1;
1317 for (i=mybegin; i<myend; i+=s) {
1318 t[i]=bval[i]; d[i]=0.0;
1319 begin_row=ia[i], end_row=ia[i+1];
1320 for (k=begin_row; k<end_row; k++) {
1321 j=ja[k];
1322 t[i]-=aval[k]*uval[j];
1323 d[i]+=ABS(aval[k]);
1324 }
1325 }
1326 }
1327#pragma omp parallel for private(i)
1328 for (i=i_1;i<=i_n;i+=s) {
1329 if (ABS(d[i])>SMALLREAL) u->val[i]+=t[i]/d[i];
1330 }
1331 }
1332 else {
1333#endif
1334 for (i=i_1;i<=i_n;i+=s) {
1335 t[i]=bval[i]; d[i]=0.0;
1336 begin_row=ia[i]; end_row=ia[i+1];
1337 for (k=begin_row;k<end_row;++k) {
1338 j=ja[k];
1339 t[i]-=aval[k]*uval[j];
1340 d[i]+=ABS(aval[k]);
1341 }
1342 }
1343
1344 for (i=i_1;i<=i_n;i+=s) {
1345 if (ABS(d[i])>SMALLREAL) u->val[i]+=t[i]/d[i];
1346 }
1347#ifdef _OPENMP
1348 }
1349#endif
1350 }
1351 else {
1352#ifdef _OPENMP
1353 if (N > OPENMP_HOLDS) {
1354#pragma omp parallel for private(myid, mybegin, myend, i, k, j, begin_row, end_row)
1355 for (myid=0; myid<nthreads; myid++) {
1356 fasp_get_start_end(myid, nthreads, N, &mybegin, &myend);
1357 mybegin = i_1 - mybegin, myend = i_1 - myend;
1358 for (i=mybegin; i>myend; i+=s) {
1359 t[i]=bval[i]; d[i]=0.0;
1360 begin_row=ia[i], end_row=ia[i+1];
1361 for (k=begin_row; k<end_row; k++) {
1362 j=ja[k];
1363 t[i]-=aval[k]*uval[j];
1364 d[i]+=ABS(aval[k]);
1365 }
1366 }
1367 }
1368#pragma omp parallel for private(i)
1369 for (i=i_1;i>=i_n;i+=s) {
1370 if (ABS(d[i])>SMALLREAL) u->val[i]+=t[i]/d[i];
1371 }
1372 }
1373 else {
1374#endif
1375 for (i=i_1;i>=i_n;i+=s) {
1376 t[i]=bval[i];d[i]=0.0;
1377 begin_row=ia[i]; end_row=ia[i+1];
1378 for (k=begin_row;k<end_row;++k) {
1379 j=ja[k];
1380 t[i]-=aval[k]*uval[j];
1381 d[i]+=ABS(aval[k]);
1382 }
1383 }
1384
1385 for (i=i_1;i>=i_n;i+=s) {
1386 if (ABS(d[i])>SMALLREAL) u->val[i]+=t[i]/d[i];
1387 }
1388#ifdef _OPENMP
1389 }
1390#endif
1391 }
1392
1393 } // end while
1394
1395 fasp_mem_free(t); t = NULL;
1396 fasp_mem_free(d); d = NULL;
1397
1398 return;
1399}
1400
1401#if 0
1422static dCSRmat form_contractor (dCSRmat *A,
1423 const INT smoother,
1424 const INT steps,
1425 const INT ndeg,
1426 const REAL relax,
1427 const REAL dtol)
1428{
1429 const INT n=A->row;
1430 INT i;
1431
1432 REAL *work = (REAL *)fasp_mem_calloc(2*n,sizeof(REAL));
1433
1434 dvector b, x;
1435 b.row=x.row=n;
1436 b.val=work; x.val=work+n;
1437
1438 INT *index = (INT *)fasp_mem_calloc(n,sizeof(INT));
1439
1440 for (i=0; i<n; ++i) index[i]=i;
1441
1442 dCSRmat B = fasp_dcsr_create(n, n, n*n); // too much memory required, need to change!!
1443
1444 dCSRmat C, D;
1445
1446 for (i=0; i<n; ++i){
1447
1448 // get i-th column
1449 fasp_dcsr_getcol(i, A, b.val);
1450
1451 // set x =0.0
1452 fasp_dvec_set(n, &x, 0.0);
1453
1454 // smooth
1455 switch (smoother) {
1456 case GS:
1457 fasp_smoother_dcsr_gs(&x, 0, n-1, 1, A, &b, steps);
1458 break;
1459 case POLY:
1460 fasp_smoother_dcsr_poly(A, &b, &x, n, ndeg, steps);
1461 break;
1462 case JACOBI:
1463 fasp_smoother_dcsr_jacobi(&x, 0, n-1, 1, A, &b, steps);
1464 break;
1465 case SGS:
1466 fasp_smoother_dcsr_sgs(&x, A, &b, steps);
1467 break;
1468 case SOR:
1469 fasp_smoother_dcsr_sor(&x, 0, n-1, 1, A, &b, steps, relax);
1470 break;
1471 case SSOR:
1472 fasp_smoother_dcsr_sor(&x, 0, n-1, 1, A, &b, steps, relax);
1473 fasp_smoother_dcsr_sor(&x, n-1, 0,-1, A, &b, steps, relax);
1474 break;
1475 case GSOR:
1476 fasp_smoother_dcsr_gs(&x, 0, n-1, 1, A, &b, steps);
1477 fasp_smoother_dcsr_sor(&x, n-1, 0, -1, A, &b, steps, relax);
1478 break;
1479 case SGSOR:
1480 fasp_smoother_dcsr_gs(&x, 0, n-1, 1, A, &b, steps);
1481 fasp_smoother_dcsr_gs(&x, n-1, 0,-1, A, &b, steps);
1482 fasp_smoother_dcsr_sor(&x, 0, n-1, 1, A, &b, steps, relax);
1483 fasp_smoother_dcsr_sor(&x, n-1, 0,-1, A, &b, steps, relax);
1484 break;
1485 default:
1486 printf("### ERROR: Unknown smoother type! [%s:%d]\n",
1487 __FILE__, __LINE__);
1488 fasp_chkerr(ERROR_INPUT_PAR, __FUNCTION__);
1489 }
1490
1491 // store to B
1492 B.IA[i] = i*n;
1493 memcpy(&(B.JA[i*n]), index, n*sizeof(INT));
1494 memcpy(&(B.val[i*n]), x.val, x.row*sizeof(REAL));
1495
1496 }
1497
1498 B.IA[n] = n*n;
1499
1500 // drop small entries
1501 compress_dCSRmat(&B, &D, dtol);
1502
1503 // get contractor
1504 fasp_dcsr_trans(&D, &C);
1505
1506 // clean up
1507 fasp_mem_free(work); work = NULL;
1508 fasp_dcsr_free(&B);
1509 fasp_dcsr_free(&D);
1510
1511 return C;
1512}
1513#endif
1514
1515/*---------------------------------*/
1516/*-- End of File --*/
1517/*---------------------------------*/
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:164
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:155
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
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:92
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
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:93
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:47
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:184
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
INT fasp_dcsr_trans(const dCSRmat *A, dCSRmat *AT)
Find transpose of dCSRmat matrix A.
Definition: BlaSparseCSR.c:952
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:493
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.
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.
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.
void fasp_smoother_dcsr_sgs(dvector *u, dCSRmat *A, dvector *b, INT L)
Symmetric Gauss-Seidel method as a smoother.
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.
void fasp_smoother_dcsr_ilu(dCSRmat *A, dvector *b, dvector *x, void *data)
ILU method as a smoother.
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.
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.
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.
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
Main header file for the FASP project.
#define REAL
Definition: fasp.h:67
#define ABS(a)
Definition: fasp.h:76
#define INT
Definition: fasp.h:64
#define FPFIRST
Definition: fasp_const.h:241
#define OPENMP_HOLDS
Definition: fasp_const.h:260
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
#define ERROR_INPUT_PAR
Definition: fasp_const.h:24
#define SMALLREAL
Definition: fasp_const.h:249
Data for ILU setup.
Definition: fasp.h:637
INT * ijlu
integer array of row pointers and column indexes, the size is nzlu
Definition: fasp.h:655
REAL * luval
nonzero entries of LU
Definition: fasp.h:658
INT nwork
work space size
Definition: fasp.h:664
REAL * work
work space
Definition: fasp.h:667
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
REAL * val
nonzero entries of A
Definition: fasp.h:161
INT row
row number of matrix A, m
Definition: fasp.h:146
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:155
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:158
Vector with n entries of REAL type.
Definition: fasp.h:346
REAL * val
actual vector entries
Definition: fasp.h:352
INT row
number of rows
Definition: fasp.h:349