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