Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
ItrSmootherCSRpoly.c
Go to the documentation of this file.
1
23#include <math.h>
24#include <time.h>
25#include <float.h>
26#include <limits.h>
27
28#ifdef _OPENMP
29#include <omp.h>
30#endif
31
32#include "fasp.h"
33#include "fasp_functs.h"
34
35/*---------------------------------*/
36/*-- Declare Private Functions --*/
37/*---------------------------------*/
38
39static void bminax (REAL *,INT *,INT *, REAL *, REAL *,INT *, REAL *);
40static void Diaginv (dCSRmat *, REAL *);
41static REAL DinvAnorminf (dCSRmat *, REAL *);
42static void Diagx (REAL *, INT, REAL *, REAL *);
43static void Rr (dCSRmat *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, INT);
44static void fasp_aux_uuplv0_ (REAL *, REAL *, INT *);
45static void fasp_aux_norm1_ (INT *, INT *, REAL *, INT *, REAL *);
46
47/*---------------------------------*/
48/*-- Public Function --*/
49/*---------------------------------*/
50
67void fasp_smoother_dcsr_poly (dCSRmat *Amat,
68 dvector *brhs,
69 dvector *usol,
70 INT n,
71 INT ndeg,
72 INT L)
73{
74 // local variables
75 INT i;
76 REAL *b = brhs->val, *u = usol->val;
77 REAL *Dinv = NULL, *r = NULL, *rbar = NULL, *v0 = NULL, *v1 = NULL;
78 REAL *error = NULL, *k = NULL;
79 REAL mu0, mu1, smu0, smu1;
80
81 /* allocate memory */
82 Dinv = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
83 r = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
84 rbar = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
85 v0 = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
86 v1 = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
87 error = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
88 k = (REAL *) fasp_mem_calloc(6,sizeof(REAL)); // coefficients for calculation
89
90 // get the inverse of the diagonal of A
91 Diaginv(Amat, Dinv);
92
93 // set up parameter
94 mu0 = DinvAnorminf(Amat, Dinv); // get the inf norm of Dinv*A;
95
96 mu0 = 1.0/mu0; mu1 = 4.0*mu0; // default set 8;
97 smu0 = sqrt(mu0); smu1 = sqrt(mu1);
98
99 k[1] = (mu0+mu1)/2.0;
100 k[2] = (smu0 + smu1)*(smu0 + smu1)/2.0;
101 k[3] = mu0 * mu1;
102
103 // 4.0*mu0*mu1/(sqrt(mu0)+sqrt(mu1))/(sqrt(mu0)+sqrt(mu1));
104 k[4] = 2.0*k[3]/k[2];
105
106 // square of (sqrt(kappa)-1)/(sqrt(kappa)+1);
107 k[5] = (mu1-2.0*smu0*smu1+mu0)/(mu1+2.0*smu0*smu1+mu0);
108
109#if DEBUG_MODE > 0
110 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
111#endif
112
113 // Update
114 for ( i=0; i<L; i++ ) {
115 // get residual
116 fasp_blas_dcsr_mxv(Amat, u, r);// r= Amat*u;
117 fasp_blas_darray_axpyz(n, -1, r, b, r);// r= -r+b;
118
119 // Get correction error = R*r
120 Rr(Amat, Dinv, r, rbar, v0, v1, error, k, ndeg);
121
122 // update solution
123 fasp_blas_darray_axpy(n, 1, error, u);
124
125 }
126
127#if DEBUG_MODE > 1
128 printf("### DEBUG: Degree of polysmoothing is: %d\n",ndeg);
129#endif
130
131 // free memory
132 fasp_mem_free(Dinv); Dinv = NULL;
133 fasp_mem_free(r); r = NULL;
134 fasp_mem_free(rbar); rbar = NULL;
135 fasp_mem_free(v0); v0 = NULL;
136 fasp_mem_free(v1); v1 = NULL;
137 fasp_mem_free(error); error = NULL;
138 fasp_mem_free(k); k = NULL;
139
140#if DEBUG_MODE > 0
141 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
142#endif
143
144 return;
145}
146
165void fasp_smoother_dcsr_poly_old (dCSRmat *Amat,
166 dvector *brhs,
167 dvector *usol,
168 INT n,
169 INT ndeg,
170 INT L)
171{
172 INT *ia=Amat->IA,*ja=Amat->JA;
173 INT i,j,k,it,jk,iaa,iab,ndeg0; // id and ij for scaling of A
174
175 REAL *a=Amat->val, *b=brhs->val, *u=usol->val;
176 REAL *v,*v0,*r,*vsave; // one can get away without r as well;
177 REAL smaxa,smina,delinv,s,smu0,smu1,skappa,th,th1,sq;
178 REAL ri,ari,vj,ravj,snj,sm,sm01,smsqrt,delta,delta2,chi;
179
180#ifdef _OPENMP
181 // variables for OpenMP
182 INT myid, mybegin, myend;
183 INT nthreads = fasp_get_num_threads();
184#endif
185
186#if DEBUG_MODE > 0
187 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
188#endif
189
190 /* WORKING MEM */
191 v = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
192 v0 = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
193 vsave = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
194 r = (REAL *) fasp_mem_calloc(n,sizeof(REAL));
195
196 /* COMPUTE PARAMS*/
197 // min INT for approx -- could be done upfront
198 // i.e., only once per level... only norm1 ...
199 fasp_aux_norm1_(ia,ja,a,&n,&smaxa);
200 smina=smaxa/8;
201 delinv=(smaxa+smina)/(smaxa-smina);
202 th=delinv+sqrt(delinv*delinv-1e+00);
203 th1=1e+00/th;
204 sq=(th-th1)*(th-th1);
205 //
206 ndeg0=(int)floor(log(2*(2e0+th+th1)/sq)/log(th)+1e0);
207 if (ndeg0 < ndeg) ndeg0=ndeg;
208 //
209 smu0=1e+00/smaxa;
210 smu1=1e+00/smina;
211 skappa=sqrt(smaxa/smina);
212 delta=(skappa-1e+00)/(skappa+1);
213 delta2=delta*delta;
214 s=sqrt(smu0)+sqrt(smu1);
215 s=s*s;
216 smsqrt=0.5e+00*s;
217 chi=4e+00*smu0*smu1/s;
218 sm=0.5e+00*(smu0+smu1);
219 sm01=smu0*smu1;
220
221#if DEBUG_MODE > 1
222 printf("### DEBUG: Degree of polysmoothing is: %d\n",ndeg);
223#endif
224
225 /* BEGIN POLY ITS */
226
227 /* auv_(ia,ja,a,u,u,&n,&err0); NA: u = 0 */
228 //bminax(b,ia,ja,a,u,&n,r);
229 //for (i=0; i < n; ++i) {res0 += r[i]*r[i];}
230 //res0=sqrt(res0);
231
232 for (it = 0 ; it < L; it++) {
233 bminax(b,ia,ja,a,u,&n,r);
234#ifdef _OPENMP
235#pragma omp parallel for private(myid,mybegin,myend,i,iaa,iab,ari,jk,j,ri) if(n>OPENMP_HOLDS)
236 for (myid=0; myid<nthreads; ++myid) {
237 fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
238 for (i=mybegin; i<myend; ++i) {
239#else
240 for (i=0; i < n ; ++i) {
241#endif
242 iaa = ia[i];
243 iab = ia[i+1];
244 ari=0e+00; /* ari is (A*r)[i] */
245 if(iab > iaa) {
246 for (jk = iaa; jk < iab; jk++) {
247 j=ja[jk];
248 ari += a[jk] * r[j];
249 }
250 }
251 ri=r[i];
252 v0[i]=sm*ri;
253 v[i]=smsqrt*ri-sm01*ari;
254 }
255#ifdef _OPENMP
256 }
257#endif
258 for (i=1; i < ndeg0; ++i) {
259 //for (j=0; j < n ; ++j) vsave[j]=v[j];
260 fasp_darray_cp(n, v, vsave);
261
262#ifdef _OPENMP
263#pragma omp parallel for private(myid,mybegin,myend,j,ravj,iaa,iab,jk,k,vj,snj) if(n>OPENMP_HOLDS)
264 for (myid=0; myid<nthreads; ++myid) {
265 fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
266 for (j=mybegin; j<myend; ++j) {
267#else
268 for (j=0; j < n ; ++j) {
269#endif
270 /* ravj = (r- A*v)[j] */
271 ravj= r[j];
272 iaa = ia[j];
273 iab = ia[j+1];
274 if(iab > iaa) {
275 for (jk = iaa; jk < iab; jk++) {
276 k=ja[jk];
277 ravj -= a[jk] * vsave[k];
278 }
279 }
280 vj=v[j];
281 snj = chi*ravj+delta2*(vj-v0[j]);
282 v0[j]=vj;
283 v[j]=vj+snj;
284 }
285 }
286#ifdef _OPENMP
287 }
288#endif
289 fasp_aux_uuplv0_(u,v,&n);
290 //bminax(b,ia,ja,a,u,&n,r);
291 //for (i=0; i < n ; ++i)
292 //resk += r[i]*r[i];
293 //resk=sqrt(resk);
294 //fprintf("\nres0=%12.5g\n",res0);
295 //fprintf("\nresk=%12.5g\n",resk);
296 //res0=resk;
297 //resk=0.0e0;
298 }
299
300 fasp_mem_free(v); v = NULL;
301 fasp_mem_free(v0); v0 = NULL;
302 fasp_mem_free(r); r = NULL;
303 fasp_mem_free(vsave); vsave = NULL;
304
305#if DEBUG_MODE > 0
306 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
307#endif
308
309 return;
310}
311
312/*---------------------------------*/
313/*-- Private Functions --*/
314/*---------------------------------*/
315
334static void bminax (REAL *b,
335 INT *ia,
336 INT *ja,
337 REAL *a,
338 REAL *x,
339 INT *nn,
340 REAL *res)
341{
342 /* Computes b-A*x */
343
344 INT i,j,jk,iaa,iab;
345 INT n;
346 REAL u;
347 n=*nn;
348
349#ifdef _OPENMP
350 // variables for OpenMP
351 INT myid, mybegin, myend;
352 INT nthreads = fasp_get_num_threads();
353#endif
354
355#ifdef _OPENMP
356#pragma omp parallel for private(myid,mybegin,myend,i,iaa,iab,u,jk,j) if(n>OPENMP_HOLDS)
357 for (myid=0; myid<nthreads; ++myid) {
358 fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
359 for (i=mybegin; i<myend; ++i) {
360#else
361 for (i=0; i < n ; ++i) {
362#endif
363 iaa = ia[i];
364 iab = ia[i+1];
365 u = b[i];
366 if(iab > iaa)
367 for (jk = iaa; jk < iab; jk++) {
368 j=ja[jk];
369 u -= a[jk] * x[j];
370 }
371 res[i] = u;
372 }
373#ifdef _OPENMP
374 }
375#endif
376 return;
377}
378
392static void Diaginv (dCSRmat *Amat,
393 REAL *Dinv)
394{
395 const INT n = Amat->row;
396 const INT *ia = Amat->IA, *ja = Amat->JA;
397 const REAL *a = Amat->val;
398 INT i,j;
399
400#ifdef _OPENMP
401#pragma omp parallel for private(j) if(n>OPENMP_HOLDS)
402#endif
403 for (i=0; i<n; i++) {
404 for(j=ia[i]; j<ia[i+1]; j++) {
405 if(i==ja[j]) // find the diagonal
406 break;
407 }
408 Dinv[i] = 1.0/a[j];
409 }
410 return;
411}
412
428static REAL DinvAnorminf (dCSRmat *Amat,
429 REAL *Dinv)
430{
431 //local variable
432 const INT n = Amat->row;
433 const INT *ia = Amat->IA;
434 const REAL *a = Amat->val;
435
436 INT i,j;
437 REAL norm, temp;
438
439#ifdef _OPENMP
440 // variables for OpenMP
441 INT myid, mybegin, myend;
442 REAL sub_norm = 0.0;
443 INT nthreads = fasp_get_num_threads();
444#endif
445
446 norm = 0.0;
447
448 // get the infinity norm of Dinv*A
449#ifdef _OPENMP
450#pragma omp parallel for private(myid,mybegin,myend,i,temp,sub_norm) if(n>OPENMP_HOLDS)
451 for (myid=0; myid<nthreads; ++myid) {
452 fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
453 sub_norm = 0.0;
454 for (i=mybegin; i<myend; ++i) {
455#else
456 for (i=0; i<n; i++) {
457#endif
458 temp = 0.0;
459 for (j=ia[i]; j<ia[i+1]; j++) {
460 temp += ABS(a[j]);
461 }
462 temp *= Dinv[i]; // temp is the L1 norm of the ith row of Dinv*A;
463#ifdef _OPENMP
464 sub_norm = MAX(sub_norm, temp);
465#else
466 norm = MAX(norm, temp);
467#endif
468 }
469#ifdef _OPENMP
470#pragma omp critical(norm)
471 norm = MAX(norm, sub_norm);
472 }
473#endif
474
475 return norm;
476}
477
493static void Diagx (REAL *Dinv,
494 INT n,
495 REAL *x,
496 REAL *b)
497{
498 INT i;
499
500 // Variables for OpenMP
501 SHORT nthreads = 1, use_openmp = FALSE;
502 INT myid, mybegin, myend;
503
504#ifdef _OPENMP
505 if (n > OPENMP_HOLDS) {
506 use_openmp = TRUE;
507 nthreads = fasp_get_num_threads();
508 }
509#endif
510
511 if (use_openmp) {
512#ifdef _OPENMP
513#pragma omp parallel for private(myid, mybegin, myend, i)
514#endif
515 for (myid = 0; myid < nthreads; myid++) {
516 fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
517 for (i = mybegin; i < myend; i++) {
518 b[i] = Dinv[i] * x[i];
519 }
520 }
521 }
522 else {
523 for (i=0; i<n; i++) {
524 b[i] = Dinv[i] * x[i];
525 }
526 }
527 return;
528}
529
551static void Rr (dCSRmat *Amat,
552 REAL *Dinv,
553 REAL *r,
554 REAL *rbar,
555 REAL *v0,
556 REAL *v1,
557 REAL *vnew,
558 REAL *k,
559 INT m)
560{
561 // local variables
562 const INT n = Amat->row;
563 INT i,j;
564
565#ifdef _OPENMP
566 // variables for OpenMP
567 INT myid, mybegin, myend;
568 INT nthreads = fasp_get_num_threads();
569#endif
570
571 //1 set up rbar
572 Diagx(Dinv, n, r, rbar);// rbar = Dinv *r;
573
574 //2 set up v0, v1;
575 fasp_blas_dcsr_mxv(Amat, rbar, v1);//v1= A*rbar;
576 Diagx(Dinv, n, v1, v1); // v1=Dinv *v1;
577
578#ifdef _OPENMP
579#pragma omp parallel for if(n>OPENMP_HOLDS)
580#endif
581 for(i=0;i<n;i++) {
582 v0[i] = k[1] * rbar[i];
583 v1[i] = k[2] * rbar[i] - k[3] * v1[i];
584 }
585
586 //3 iterate to get v_(j+1)
587
588 for (j=1;j<m;j++) {
589 fasp_blas_dcsr_mxv(Amat, v1, rbar);//rbar= A*v_(j);
590
591#ifdef _OPENMP
592#pragma omp parallel for private(myid,mybegin,myend,i) if(n>OPENMP_HOLDS)
593 for (myid=0; myid<nthreads; ++myid) {
594 fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
595 for (i=mybegin; i<myend; ++i) {
596#else
597 for(i=0;i<n;i++) {
598#endif
599 rbar[i] = (r[i] - rbar[i])*Dinv[i];// indeed rbar=Dinv*(r-A*v_(j));
600 vnew[i] = v1[i] + k[5] *(v1[i] - v0[i]) + k[4] * rbar[i];// compute v_(j+1)
601 // prepare for next cycle
602 v0[i]=v1[i];
603 v1[i]=vnew[i];
604 }
605#ifdef _OPENMP
606 }
607#endif
608 }
609}
610
611static void fasp_aux_uuplv0_ (REAL *u,
612 REAL *v,
613 INT *n)
614{
615 /*
616 This computes y = y + x.
617 */
618 INT i;
619 for ( i=0; i < *n; i++ ) u[i] += v[i];
620 return;
621}
622
623static void fasp_aux_norm1_ (INT *ia,
624 INT *ja,
625 REAL *a,
626 INT *nn,
627 REAL *a1norm)
628{
629 INT n,i,jk,iaa,iab;
630 REAL sum,s;
631 /* computes one norm of a matrix a and stores it in the variable
632 pointed to by *a1norm*/
633 n = *nn;
634 s = 0.0;
635 for ( i=0; i < n ; i++ ) {
636 iaa = ia[i];
637 iab = ia[i+1];
638 sum = 0e+00;
639 for ( jk = iaa; jk < iab; jk++ ) sum += fabs(a[jk]);
640 if ( sum > s ) s = sum;
641 }
642 *a1norm=s;
643}
644
645/*---------------------------------*/
646/*-- End of File --*/
647/*---------------------------------*/
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_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_blas_darray_axpyz
void fasp_blas_darray_axpyz(const INT n, const REAL a, const REAL *x, const REAL *y, REAL *z)
z = a*x + y
Definition: BlaArray.c:403
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_blas_dcsr_mxv
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:242
fasp_smoother_dcsr_poly_old
void fasp_smoother_dcsr_poly_old(dCSRmat *Amat, dvector *brhs, dvector *usol, INT n, INT ndeg, INT L)
poly approx to A^{-1} as MG smoother: JK&LTZ2010
Definition: ItrSmootherCSRpoly.c:165
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
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
ABS
#define ABS(a)
Definition: fasp.h:84
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
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
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
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