Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
KryPgcg.c
Go to the documentation of this file.
1
22#include <math.h>
23
24#include "fasp.h"
25#include "fasp_functs.h"
26
27/*---------------------------------*/
28/*-- Declare Private Functions --*/
29/*---------------------------------*/
30
31#include "KryUtil.inl"
32
33/*---------------------------------*/
34/*-- Public Functions --*/
35/*---------------------------------*/
36
61 dvector *b,
62 dvector *u,
63 precond *pc,
64 const REAL tol,
65 const INT MaxIt,
66 const SHORT StopType,
67 const SHORT PrtLvl)
68{
69 INT iter=0, m=A->row, i;
70 REAL absres0 = BIGREAL, absres = BIGREAL;
71 REAL relres = BIGREAL, normb = BIGREAL;
72 REAL alpha, factor;
73
74 // allocate temp memory
75 REAL *work = (REAL *)fasp_mem_calloc(2*m+MaxIt+MaxIt*m,sizeof(REAL));
76
77 REAL *r, *Br, *beta, *p;
78 r = work; Br = r + m; beta = Br + m; p = beta + MaxIt;
79
80 // Output some info for debuging
81 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GCG solver (CSR) ...\n");
82
83#if DEBUG_MODE > 0
84 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
85 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
86#endif
87
88 normb=fasp_blas_darray_norm2(m,b->val);
89
90 // -------------------------------------
91 // 1st iteration (Steepest descent)
92 // -------------------------------------
93 // r = b-A*u
94 fasp_darray_cp(m,b->val,r);
95 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
96
97 // Br
98 if (pc != NULL)
99 pc->fct(r,p,pc->data); /* Preconditioning */
100 else
101 fasp_darray_cp(m,r,p); /* No preconditioner, B=I */
102
103 // alpha = (p'r)/(p'Ap)
104 alpha = fasp_blas_darray_dotprod (m,r,p) / fasp_blas_dcsr_vmv (A, p, p);
105
106 // u = u + alpha *p
107 fasp_blas_darray_axpy(m, alpha , p, u->val);
108
109 // r = r - alpha *Ap
110 fasp_blas_dcsr_aAxpy((-1.0*alpha),A,p,r);
111
112 // norm(r), factor
113 absres = fasp_blas_darray_norm2(m,r); factor = absres/absres0;
114
115 // compute relative residual
116 relres = absres/normb;
117
118 // output iteration information if needed
119 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
120
121 // update relative residual here
122 absres0 = absres;
123
124 for ( iter = 1; iter < MaxIt ; iter++) {
125
126 // Br
127 if (pc != NULL)
128 pc->fct(r, Br ,pc->data); // Preconditioning
129 else
130 fasp_darray_cp(m,r, Br); // No preconditioner, B=I
131
132 // form p
133 fasp_darray_cp(m, Br, p+iter*m);
134
135 for (i=0; i<iter; i++) {
136 beta[i] = (-1.0) * ( fasp_blas_dcsr_vmv (A, Br, p+i*m)
137 /fasp_blas_dcsr_vmv (A, p+i*m, p+i*m) );
138
139 fasp_blas_darray_axpy(m, beta[i], p+i*m, p+iter*m);
140 }
141
142 // -------------------------------------
143 // next iteration
144 // -------------------------------------
145
146 // alpha = (p'r)/(p'Ap)
147 alpha = fasp_blas_darray_dotprod(m,r,p+iter*m)
148 / fasp_blas_dcsr_vmv (A, p+iter*m, p+iter*m);
149
150 // u = u + alpha *p
151 fasp_blas_darray_axpy(m, alpha , p+iter*m, u->val);
152
153 // r = r - alpha *Ap
154 fasp_blas_dcsr_aAxpy((-1.0*alpha),A,p+iter*m,r);
155
156 // norm(r), factor
157 absres = fasp_blas_darray_norm2(m,r); factor = absres/absres0;
158
159 // compute relative residual
160 relres = absres/normb;
161
162 // output iteration information if needed
163 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
164
165 if (relres < tol) break;
166
167 // update relative residual here
168 absres0 = absres;
169
170 } // end of main GCG loop.
171
172 // finish iterative method
173 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
174
175 // clean up temp memory
176 fasp_mem_free(work); work = NULL;
177
178#if DEBUG_MODE > 0
179 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
180#endif
181
182 if (iter>MaxIt)
183 return ERROR_SOLVER_MAXIT;
184 else
185 return iter;
186}
187
214 dvector *b,
215 dvector *u,
216 precond *pc,
217 const REAL tol,
218 const INT MaxIt,
219 const SHORT StopType,
220 const SHORT PrtLvl)
221{
222 INT iter=0, m=b->row, i;
223 REAL absres0 = BIGREAL, absres = BIGREAL;
224 REAL relres = BIGREAL, normb = BIGREAL;
225 REAL alpha, factor, gama_1, gama_2;
226
227 // allocate temp memory
228 REAL *work = (REAL *)fasp_mem_calloc(3*m+MaxIt+MaxIt*m,sizeof(REAL));
229
230 REAL *r, *Br, *beta, *p, *q;
231 q = work; r = q + m; Br = r + m; beta = Br + m; p = beta + MaxIt;
232
233 // Output some info for debuging
234 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GCG solver (MatFree) ...\n");
235
236#if DEBUG_MODE > 0
237 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
238 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
239#endif
240
241 normb=fasp_blas_darray_norm2(m,b->val);
242
243 // -------------------------------------
244 // 1st iteration (Steepest descent)
245 // -------------------------------------
246 // r = b-A*u
247 mf->fct(mf->data, u->val, r);
248 fasp_blas_darray_axpby(m, 1.0, b->val, -1.0, r);
249
250 // Br
251 if (pc != NULL)
252 pc->fct(r,p,pc->data); /* Preconditioning */
253 else
254 fasp_darray_cp(m,r,p); /* No preconditioner, B=I */
255
256 // alpha = (p'r)/(p'Ap)
257 mf->fct(mf->data, p, q);
258 alpha = fasp_blas_darray_dotprod (m,r,p) / fasp_blas_darray_dotprod (m, p, q);
259
260 // u = u + alpha *p
261 fasp_blas_darray_axpy(m, alpha , p, u->val);
262
263 // r = r - alpha *Ap
264 mf->fct(mf->data, p, q);
265 fasp_blas_darray_axpby(m, (-1.0*alpha), q, 1.0, r);
266
267 // norm(r), factor
268 absres = fasp_blas_darray_norm2(m,r); factor = absres/absres0;
269
270 // compute relative residual
271 relres = absres/normb;
272
273 // output iteration information if needed
274 fasp_itinfo(PrtLvl,StopType,iter+1,relres,absres,factor);
275
276 // update relative residual here
277 absres0 = absres;
278
279 for ( iter = 1; iter < MaxIt ; iter++) {
280
281 // Br
282 if (pc != NULL)
283 pc->fct(r, Br ,pc->data); // Preconditioning
284 else
285 fasp_darray_cp(m,r, Br); // No preconditioner, B=I
286
287 // form p
288 fasp_darray_cp(m, Br, p+iter*m);
289
290 for (i=0; i<iter; i++) {
291 mf->fct(mf->data, Br, q);
292 gama_1 = fasp_blas_darray_dotprod(m, p+i*m, q);
293 mf->fct(mf->data, p+i*m, q);
294 gama_2 = fasp_blas_darray_dotprod(m, p+i*m, q);
295 beta[i] = (-1.0) * ( gama_1 / gama_2 );
296
297 fasp_blas_darray_axpy(m, beta[i], p+i*m, p+iter*m);
298 }
299
300 // -------------------------------------
301 // next iteration
302 // -------------------------------------
303
304 // alpha = (p'r)/(p'Ap)
305 mf->fct(mf->data, p+iter*m, q);
306 alpha = fasp_blas_darray_dotprod(m,r,p+iter*m)
307 / fasp_blas_darray_dotprod (m, q, p+iter*m);
308
309 // u = u + alpha *p
310 fasp_blas_darray_axpy(m, alpha , p+iter*m, u->val);
311
312 // r = r - alpha *Ap
313 mf->fct(mf->data, p+iter*m, q);
314 fasp_blas_darray_axpby(m, (-1.0*alpha), q, 1.0, r);
315
316 // norm(r), factor
317 absres = fasp_blas_darray_norm2(m,r); factor = absres/absres0;
318
319 // compute relative residual
320 relres = absres/normb;
321
322 // output iteration information if needed
323 fasp_itinfo(PrtLvl,StopType,iter+1,relres,absres,factor);
324
325 if (relres < tol) break;
326
327 // update relative residual here
328 absres0 = absres;
329
330 } // end of main GCG loop.
331
332 // finish iterative method
333 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
334
335 // clean up temp memory
336 fasp_mem_free(work); work = NULL;
337
338#if DEBUG_MODE > 0
339 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
340#endif
341
342 if (iter>MaxIt)
343 return ERROR_SOLVER_MAXIT;
344 else
345 return iter;
346}
347
348/*---------------------------------*/
349/*-- End of File --*/
350/*---------------------------------*/
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
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
REAL fasp_blas_darray_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: BlaArray.c:657
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:93
REAL fasp_blas_dcsr_vmv(const dCSRmat *A, const REAL *x, const REAL *y)
vector-Matrix-vector multiplication alpha = y'*A*x
Definition: BlaSpmvCSR.c:834
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_pgcg(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned generilzed conjugate gradient (GCG) method for solving Au=b.
Definition: KryPgcg.c:213
INT fasp_solver_dcsr_pgcg(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned generilzed conjugate gradient (GCG) method for solving Au=b.
Definition: KryPgcg.c:60
Main header file for the FASP project.
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#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 PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
INT row
row number of matrix A, m
Definition: fasp.h:146
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
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1095
void * data
data for MxV, can be a Matrix or something else
Definition: fasp.h:1098
void(* fct)(const void *, const REAL *, REAL *)
action for MxV, void function pointer
Definition: fasp.h:1101
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