Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
KryPgcr.c
Go to the documentation of this file.
1
17#include <math.h>
18
19#include "fasp.h"
20#include "fasp_functs.h"
21
22/*---------------------------------*/
23/*-- Declare Private Functions --*/
24/*---------------------------------*/
25
26#include "KryUtil.inl"
27
28static void dense_aAtxpby (INT, INT, REAL *, REAL, REAL *, REAL, REAL *);
29
56 dvector *b,
57 dvector *x,
58 precond *pc,
59 const REAL tol,
60 const INT MaxIt,
61 const SHORT restart,
62 const SHORT StopType,
63 const SHORT PrtLvl)
64{
65 const INT n = b->row;
66
67 // local variables
68 INT iter = 0;
69 int i, j, k, rst = -1; // must be signed! -zcs
70
71 REAL gamma, alpha, beta, checktol;
72 REAL absres0 = BIGREAL, absres = BIGREAL;
73 REAL relres = BIGREAL;
74
75 // allocate temp memory (need about (restart+4)*n REAL numbers)
76 REAL *c = NULL, *z = NULL, *alp = NULL, *tmpx = NULL;
77 REAL *norms = NULL, *r = NULL, *work = NULL;
78 REAL **h = NULL;
79
80 INT Restart = MIN(restart, MaxIt);
81 LONG worksize = n+2*Restart*n+Restart+Restart;
82
83 // Output some info for debugging
84 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GCR solver (CSR) ...\n");
85
86#if DEBUG_MODE > 0
87 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
88 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
89#endif
90
91 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
92
93 /* check whether memory is enough for GCR */
94 while ( (work == NULL) && (Restart > 5) ) {
95 Restart = Restart - 5;
96 worksize = n+2*Restart*n+Restart+Restart;
97 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
98 }
99
100 if ( work == NULL ) {
101 printf("### ERROR: No enough memory for GCR! [%s:%d]\n",
102 __FILE__, __LINE__ );
103 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
104 }
105
106 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
107 printf("### WARNING: GCR restart number set to %d!\n", Restart);
108 }
109
110 r = work; z = r+n; c = z + Restart*n; alp = c + Restart*n; tmpx = alp + Restart;
111
112 h = (REAL **)fasp_mem_calloc(Restart, sizeof(REAL *));
113 for (i = 0; i < Restart; i++) h[i] = (REAL*)fasp_mem_calloc(Restart, sizeof(REAL));
114
115 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
116
117 // r = b-A*x
118 fasp_darray_cp(n, b->val, r);
119 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
120
121 absres = fasp_blas_darray_dotprod(n, r, r);
122
123 absres0 = MAX(SMALLREAL,absres);
124
125 relres = absres/absres0;
126
127 // output iteration information if needed
128 fasp_itinfo(PrtLvl,StopType,0,relres,sqrt(absres0),0.0);
129
130 // store initial residual
131 norms[0] = relres;
132
133 checktol = MAX(tol*tol*absres0, absres*1.0e-4);
134
135 while ( iter < MaxIt && sqrt(relres) > tol ) {
136
137 i = -1; rst ++;
138
139 while ( i < Restart-1 && iter < MaxIt ) {
140
141 i++; iter++;
142
143 // z = B^-1r
144 if ( pc == NULL )
145 fasp_darray_cp(n, r, &z[i*n]);
146 else
147 pc->fct(r, &z[i*n], pc->data);
148
149 // c = Az
150 fasp_blas_dcsr_mxv(A, &z[i*n], &c[i*n]);
151
152 /* Modified Gram_Schmidt orthogonalization */
153 for ( j = 0; j < i; j++ ) {
154 gamma = fasp_blas_darray_dotprod(n, &c[j*n], &c[i*n]);
155 h[i][j] = gamma/h[j][j];
156 fasp_blas_darray_axpy(n, -h[i][j], &c[j*n], &c[i*n]);
157 }
158 // gamma = (c,c)
159 gamma = fasp_blas_darray_dotprod(n, &c[i*n], &c[i*n]);
160
161 h[i][i] = gamma;
162
163 // alpha = (c, r)
164 alpha = fasp_blas_darray_dotprod(n, &c[i*n], r);
165
166 beta = alpha/gamma;
167
168 alp[i] = beta;
169
170 // r = r - beta*c
171 fasp_blas_darray_axpy(n, -beta, &c[i*n], r);
172
173 // equivalent to ||r||_2
174 absres = absres - alpha*alpha/gamma;
175
176 if (absres < checktol) {
177 absres = fasp_blas_darray_dotprod(n, r, r);
178 checktol = MAX(tol*tol*absres0, absres*1.0e-4);
179 }
180
181 relres = absres / absres0;
182
183 norms[iter] = relres;
184
185 fasp_itinfo(PrtLvl, StopType, iter, sqrt(relres), sqrt(absres),
186 sqrt(norms[iter]/norms[iter-1]));
187
188 if (sqrt(relres) < tol) break;
189 }
190
191 for ( k = i; k >=0; k-- ) {
192 tmpx[k] = alp[k];
193 for (j=0; j<k; ++j) {
194 alp[j] -= h[k][j]*tmpx[k];
195 }
196 }
197
198 if (rst==0) dense_aAtxpby(n, i+1, z, 1.0, tmpx, 0.0, x->val);
199 else dense_aAtxpby(n, i+1, z, 1.0, tmpx, 1.0, x->val);
200
201 }
202
203 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,sqrt(relres));
204
205 // clean up memory
206 for (i = 0; i < Restart; i++) {
207 fasp_mem_free(h[i]); h[i] = NULL;
208 }
209 fasp_mem_free(h); h = NULL;
210
211 fasp_mem_free(work); work = NULL;
212 fasp_mem_free(norms); norms = NULL;
213
214#if DEBUG_MODE > 0
215 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
216#endif
217
218 if ( iter >= MaxIt )
219 return ERROR_SOLVER_MAXIT;
220 else
221 return iter;
222}
223
250 dvector *b,
251 dvector *x,
252 precond *pc,
253 const REAL tol,
254 const INT MaxIt,
255 const SHORT restart,
256 const SHORT StopType,
257 const SHORT PrtLvl)
258{
259 const INT n = b->row;
260
261 // local variables
262 INT iter = 0;
263 int i, j, k, rst = -1; // must be signed! -zcs
264
265 REAL gamma, alpha, beta, checktol;
266 REAL absres0 = BIGREAL, absres = BIGREAL;
267 REAL relres = BIGREAL;
268
269 // allocate temp memory (need about (restart+4)*n REAL numbers)
270 REAL *c = NULL, *z = NULL, *alp = NULL, *tmpx = NULL;
271 REAL *norms = NULL, *r = NULL, *work = NULL;
272 REAL **h = NULL;
273
274 INT Restart = MIN(restart, MaxIt);
275 LONG worksize = n+2*Restart*n+Restart+Restart;
276
277 // Output some info for debugging
278 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GCR solver (BLC) ...\n");
279
280#if DEBUG_MODE > 0
281 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
282 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
283#endif
284
285 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
286
287 /* check whether memory is enough for GCR */
288 while ( (work == NULL) && (Restart > 5) ) {
289 Restart = Restart - 5;
290 worksize = n+2*Restart*n+Restart+Restart;
291 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
292 }
293
294 if ( work == NULL ) {
295 printf("### ERROR: No enough memory for GCR! [%s:%d]\n",
296 __FILE__, __LINE__ );
297 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
298 }
299
300 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
301 printf("### WARNING: GCR restart number set to %d!\n", Restart);
302 }
303
304 r = work; z = r+n; c = z + Restart*n; alp = c + Restart*n; tmpx = alp + Restart;
305
306 h = (REAL **)fasp_mem_calloc(Restart, sizeof(REAL *));
307 for (i = 0; i < Restart; i++) h[i] = (REAL*)fasp_mem_calloc(Restart, sizeof(REAL));
308
309 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
310
311 // r = b-A*x
312 fasp_darray_cp(n, b->val, r);
313 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
314
315 absres = fasp_blas_darray_dotprod(n, r, r);
316
317 absres0 = MAX(SMALLREAL,absres);
318
319 relres = absres/absres0;
320
321 // output iteration information if needed
322 fasp_itinfo(PrtLvl,StopType,0,relres,sqrt(absres0),0.0);
323
324 // store initial residual
325 norms[0] = relres;
326
327 checktol = MAX(tol*tol*absres0, absres*1.0e-4);
328
329 while ( iter < MaxIt && sqrt(relres) > tol ) {
330
331 i = 0; rst ++;
332 while ( i < Restart && iter < MaxIt ) {
333
334 iter++;
335
336 // z = B^-1r
337 if ( pc == NULL )
338 fasp_darray_cp(n, r, &z[i*n]);
339 else
340 pc->fct(r, &z[i*n], pc->data);
341
342 // c = Az
343 fasp_blas_dblc_mxv(A, &z[i*n], &c[i*n]);
344
345 /* Modified Gram_Schmidt orthogonalization */
346 for ( j = 0; j < i; j++ ) {
347 gamma = fasp_blas_darray_dotprod(n, &c[j*n], &c[i*n]);
348 h[i][j] = gamma/h[j][j];
349 fasp_blas_darray_axpy(n, -h[i][j], &c[j*n], &c[i*n]);
350 }
351 // gamma = (c,c)
352 gamma = fasp_blas_darray_dotprod(n, &c[i*n], &c[i*n]);
353
354 h[i][i] = gamma;
355
356 // alpha = (c, r)
357 alpha = fasp_blas_darray_dotprod(n, &c[i*n], r);
358
359 beta = alpha/gamma;
360
361 alp[i] = beta;
362
363 // r = r - beta*c
364 fasp_blas_darray_axpy(n, -beta, &c[i*n], r);
365
366 // equivalent to ||r||_2
367 absres = absres - alpha*alpha/gamma;
368
369 if (absres < checktol) {
370 absres = fasp_blas_darray_dotprod(n, r, r);
371 checktol = MAX(tol*tol*absres0, absres*1.0e-4);
372 }
373
374 relres = absres / absres0;
375
376 norms[iter] = relres;
377
378 fasp_itinfo(PrtLvl, StopType, iter, sqrt(relres), sqrt(absres),
379 sqrt(norms[iter]/norms[iter-1]));
380
381 if (sqrt(relres) < tol) break;
382
383 i++;
384 }
385
386 for ( k = i; k >=0; k-- ) {
387 tmpx[k] = alp[k];
388 for (j=0; j<k; ++j) {
389 alp[j] -= h[k][j]*tmpx[k];
390 }
391 }
392
393 if (rst==0) dense_aAtxpby(n, i+1, z, 1.0, tmpx, 0.0, x->val);
394 else dense_aAtxpby(n, i+1, z, 1.0, tmpx, 1.0, x->val);
395 }
396
397 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,sqrt(relres));
398
399 // clean up memory
400 for (i = 0; i < Restart; i++) {
401 fasp_mem_free(h[i]); h[i] = NULL;
402 }
403 fasp_mem_free(h); h = NULL;
404
405 fasp_mem_free(work); work = NULL;
406 fasp_mem_free(norms); norms = NULL;
407
408#if DEBUG_MODE > 0
409 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
410#endif
411
412 if ( iter >= MaxIt )
413 return ERROR_SOLVER_MAXIT;
414 else
415 return iter;
416}
417
418/*---------------------------------*/
419/*-- Private Functions --*/
420/*---------------------------------*/
421
441static void dense_aAtxpby (INT n,
442 INT m,
443 REAL *A,
444 REAL alpha,
445 REAL *x,
446 REAL beta,
447 REAL *y)
448{
449 INT i, j;
450
451 for (i=0; i<m; i++) fasp_blas_darray_ax(n, x[i], &A[i*n]);
452
453 for (j=1; j<m; j++) {
454 for (i=0; i<n; i++) {
455 A[i] += A[i+j*n];
456 }
457 }
458
459 fasp_blas_darray_axpby(n, alpha, A, beta, y);
460}
461
462/*---------------------------------*/
463/*-- End of File --*/
464/*---------------------------------*/
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_itinfo(const INT ptrlvl, const INT stop_type, const INT iter, const REAL relres, const REAL absres, const REAL factor)
Print out iteration information for iterative solvers.
Definition: AuxMessage.c:41
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
REAL fasp_blas_darray_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies x and y.
Definition: BlaArray.c:741
void fasp_blas_darray_axpby(const INT n, const REAL a, const REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: BlaArray.c:580
void fasp_blas_darray_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: BlaArray.c:43
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:93
void fasp_blas_dblc_mxv(const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvBLC.c:164
void fasp_blas_dblc_aAxpy(const REAL alpha, const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvBLC.c:38
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:241
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
INT fasp_solver_dcsr_pgcr(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
A preconditioned GCR method for solving Au=b.
Definition: KryPgcr.c:55
INT fasp_solver_dblc_pgcr(dBLCmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
A preconditioned GCR method for solving Au=b.
Definition: KryPgcr.c:249
Main header file for the FASP project.
#define MIN(a, b)
Definition: fasp.h:75
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define LONG
Definition: fasp.h:65
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:74
#define INT
Definition: fasp.h:64
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:48
#define BIGREAL
Some global constants.
Definition: fasp_const.h:248
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
#define SMALLREAL
Definition: fasp_const.h:249
#define PRINT_MIN
Definition: fasp_const.h:74
Block REAL CSR matrix format.
Definition: fasp_block.h:74
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
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
Preconditioner data and action.
Definition: fasp.h:1081
void * data
data for preconditioner, void pointer
Definition: fasp.h:1084
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1087