Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pgcg_mf.c
Go to the documentation of this file.
1 
11 #include <math.h>
12 
13 #include "fasp.h"
14 #include "fasp_functs.h"
15 #include "itsolver_util.inl"
16 
17 /*---------------------------------*/
18 /*-- Public Functions --*/
19 /*---------------------------------*/
20 
48  dvector *b,
49  dvector *u,
50  precond *pc,
51  const REAL tol,
52  const INT MaxIt,
53  const SHORT stop_type,
54  const SHORT prtlvl)
55 {
56  INT iter=0, m=b->row, i;
57  REAL absres0 = BIGREAL, absres = BIGREAL;
58  REAL relres = BIGREAL, normb = BIGREAL;
59  REAL alpha, factor, gama_1, gama_2;
60 
61  // allocate temp memory
62  REAL *work = (REAL *)fasp_mem_calloc(3*m+MaxIt+MaxIt*m,sizeof(REAL));
63 
64  REAL *r, *Br, *beta, *p, *q;
65  q = work; r = q + m; Br = r + m; beta = Br + m; p = beta + MaxIt;
66 
67 #if DEBUG_MODE > 0
68  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
69  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
70 #endif
71 
72  normb=fasp_blas_array_norm2(m,b->val);
73 
74  // -------------------------------------
75  // 1st iteration (Steepest descent)
76  // -------------------------------------
77  // r = b-A*u
78  mf->fct(mf->data, u->val, r);
79  fasp_blas_array_axpby(m, 1.0, b->val, -1.0, r);
80 
81  // Br
82  if (pc != NULL)
83  pc->fct(r,p,pc->data); /* Preconditioning */
84  else
85  fasp_array_cp(m,r,p); /* No preconditioner, B=I */
86 
87  // alpha = (p'r)/(p'Ap)
88  mf->fct(mf->data, p, q);
89  alpha = fasp_blas_array_dotprod (m,r,p) / fasp_blas_array_dotprod (m, p, q);
90 
91  // u = u + alpha *p
92  fasp_blas_array_axpy(m, alpha , p, u->val);
93 
94  // r = r - alpha *Ap
95  mf->fct(mf->data, p, q);
96  fasp_blas_array_axpby(m, (-1.0*alpha), q, 1.0, r);
97 
98  // norm(r), factor
99  absres = fasp_blas_array_norm2(m,r); factor = absres/absres0;
100 
101  // compute relative residual
102  relres = absres/normb;
103 
104  // output iteration information if needed
105  print_itinfo(prtlvl,stop_type,iter+1,relres,absres,factor);
106 
107  // update relative residual here
108  absres0 = absres;
109 
110  for ( iter = 1; iter < MaxIt ; iter++) {
111 
112  // Br
113  if (pc != NULL)
114  pc->fct(r, Br ,pc->data); // Preconditioning
115  else
116  fasp_array_cp(m,r, Br); // No preconditioner, B=I
117 
118  // form p
119  fasp_array_cp(m, Br, p+iter*m);
120 
121  for (i=0; i<iter; i++) {
122  mf->fct(mf->data, Br, q);
123  gama_1 = fasp_blas_array_dotprod(m, p+i*m, q);
124  mf->fct(mf->data, p+i*m, q);
125  gama_2 = fasp_blas_array_dotprod(m, p+i*m, q);
126  beta[i] = (-1.0) * ( gama_1 / gama_2 );
127 
128  fasp_blas_array_axpy(m, beta[i], p+i*m, p+iter*m);
129  }
130 
131  // -------------------------------------
132  // next iteration
133  // -------------------------------------
134 
135  // alpha = (p'r)/(p'Ap)
136  mf->fct(mf->data, p+iter*m, q);
137  alpha = fasp_blas_array_dotprod(m,r,p+iter*m)
138  / fasp_blas_array_dotprod (m, q, p+iter*m);
139 
140  // u = u + alpha *p
141  fasp_blas_array_axpy(m, alpha , p+iter*m, u->val);
142 
143  // r = r - alpha *Ap
144  mf->fct(mf->data, p+iter*m, q);
145  fasp_blas_array_axpby(m, (-1.0*alpha), q, 1.0, r);
146 
147  // norm(r), factor
148  absres = fasp_blas_array_norm2(m,r); factor = absres/absres0;
149 
150  // compute relative residual
151  relres = absres/normb;
152 
153  // output iteration information if needed
154  print_itinfo(prtlvl,stop_type,iter+1,relres,absres,factor);
155 
156  if (relres < tol) break;
157 
158  // update relative residual here
159  absres0 = absres;
160 
161  } // end of main GCG loop.
162 
163  // finish the iterative method
164  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
165 
166  // clean up temp memory
167  fasp_mem_free(work);
168 
169 #if DEBUG_MODE > 0
170  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
171 #endif
172 
173  if (iter>MaxIt)
174  return ERROR_SOLVER_MAXIT;
175  else
176  return iter;
177 }
178 
179 /*---------------------------------*/
180 /*-- End of File --*/
181 /*---------------------------------*/
#define REAL
Definition: fasp.h:67
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1005
Vector with n entries of REAL type.
Definition: fasp.h:342
void(* fct)(void *, REAL *, REAL *)
action for MxV, void function pointer
Definition: fasp.h:1019
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:54
#define INT
Definition: fasp.h:64
void print_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: message.c:36
Preconditioner data and action.
Definition: fasp.h:999
void fasp_blas_array_axpby(const INT n, const REAL a, REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: blas_array.c:218
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void * data
data for preconditioner, void pointer
Definition: fasp.h:1002
Main header file for FASP.
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1013
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
INT fasp_solver_pgcg(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned generilzed conjugate gradient (GCG) method for solving Au=b.
Definition: pgcg_mf.c:47
INT row
number of rows
Definition: fasp.h:345
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
REAL fasp_blas_array_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies (x,y)
Definition: blas_array.c:267
REAL fasp_blas_array_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: blas_array.c:347
#define BIGREAL
Some global constants.
Definition: fasp_const.h:237
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
void * data
data for MxV, can be a Matrix or something else
Definition: fasp.h:1016