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