Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pcg_mf.c
Go to the documentation of this file.
1 
48 #include <math.h>
49 
50 #ifdef _OPENMP
51 #include <omp.h>
52 #endif
53 
54 #include "fasp.h"
55 #include "fasp_functs.h"
56 #include "itsolver_util.inl"
57 
58 /*---------------------------------*/
59 /*-- Public Functions --*/
60 /*---------------------------------*/
61 
87  dvector *b,
88  dvector *u,
89  precond *pc,
90  const REAL tol,
91  const INT MaxIt,
92  const SHORT stop_type,
93  const SHORT prtlvl)
94 {
95  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
96  const INT m=b->row;
97  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
98  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
99 
100  // local variables
101  INT iter = 0, stag, more_step, restart_step;
102  REAL absres0 = BIGREAL, absres = BIGREAL;
103  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
104  REAL reldiff, factor, infnormu;
105  REAL alpha, beta, temp1, temp2;
106 
107  // allocate temp memory (need 4*m REAL numbers)
108  REAL *work=(REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
109  REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
110 
111 #if DEBUG_MODE > 0
112  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
113  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
114 #endif
115 
116  // initialize counters
117  stag=1; more_step=1; restart_step=1;
118 
119  // r = b-A*u
120  mf->fct(mf->data, u->val, r);
121  fasp_blas_array_axpby(m, 1.0, b->val, -1.0, r);
122 
123  if (pc != NULL)
124  pc->fct(r,z,pc->data); /* Apply preconditioner */
125  else
126  fasp_array_cp(m,r,z); /* No preconditioner */
127 
128  // compute initial relative residual
129  switch (stop_type) {
130  case STOP_REL_PRECRES:
131  absres0=sqrt(fasp_blas_array_dotprod(m,r,z));
132  normr0=MAX(SMALLREAL,absres0);
133  relres=absres0/normr0;
134  break;
135  case STOP_MOD_REL_RES:
136  absres0=fasp_blas_array_norm2(m,r);
137  normu=MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
138  relres=absres0/normu;
139  break;
140  default:
141  absres0=fasp_blas_array_norm2(m,r);
142  normr0=MAX(SMALLREAL,absres0);
143  relres=absres0/normr0;
144  break;
145  }
146 
147  // if initial residual is small, no need to iterate!
148  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
149 
150  fasp_array_cp(m,z,p);
151  temp1=fasp_blas_array_dotprod(m,z,r);
152 
153  while ( iter++ < MaxIt ) {
154 
155  // t=A*p
156  mf->fct(mf->data, p, t);
157 
158  // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
159  temp2=fasp_blas_array_dotprod(m,t,p);
160  alpha=temp1/temp2;
161 
162  // u_k=u_{k-1} + alpha_k*p_{k-1}
163  fasp_blas_array_axpy(m,alpha,p,u->val);
164 
165  // r_k=r_{k-1} - alpha_k*A*p_{k-1}
166  fasp_blas_array_axpy(m,-alpha,t,r);
167  absres=fasp_blas_array_norm2(m,r);
168 
169  // compute reducation factor of residual ||r||
170  factor=absres/absres0;
171 
172  // compute relative residual
173  switch (stop_type) {
174  case STOP_REL_PRECRES:
175  // z = B(r)
176  if (pc != NULL)
177  pc->fct(r,z,pc->data); /* Apply preconditioner */
178  else
179  fasp_array_cp(m,r,z); /* No preconditioner */
180  temp2=fasp_blas_array_dotprod(m,z,r);
181  relres=sqrt(ABS(temp2))/normr0;
182  break;
183  case STOP_MOD_REL_RES:
184  relres=absres/normu;
185  break;
186  default:
187  relres=absres/normr0;
188  break;
189  }
190 
191  // output iteration information if needed
192  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
193 
194  // solution check, if soultion is too small, return ERROR_SOLVER_SOLSTAG.
195  infnormu = fasp_blas_array_norminf(m, u->val);
196  if ( infnormu <= sol_inf_tol ) {
197  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
198  iter = ERROR_SOLVER_SOLSTAG;
199  break;
200  }
201 
202  // compute relative difference
203  normu=fasp_blas_dvec_norm2(u);
204  reldiff=ABS(alpha)*fasp_blas_array_norm2(m,p)/normu;
205 
206  // stagnation check
207  if ( (stag<=MaxStag) & (reldiff<maxdiff) ) {
208 
209  if ( prtlvl >= PRINT_MORE ) {
210  ITS_DIFFRES(reldiff,relres);
211  ITS_RESTART;
212  }
213 
214  mf->fct(mf->data, u->val, r);
215  fasp_blas_array_axpby(m, 1.0, b->val, -1.0, r);
216  absres=fasp_blas_array_norm2(m,r);
217 
218  // relative residual
219  switch (stop_type) {
220  case STOP_REL_PRECRES:
221  // z = B(r)
222  if (pc != NULL)
223  pc->fct(r,z,pc->data); /* Apply preconditioner */
224  else
225  fasp_array_cp(m,r,z); /* No preconditioner */
226  temp2=fasp_blas_array_dotprod(m,z,r);
227  relres=sqrt(ABS(temp2))/normr0;
228  break;
229  case STOP_MOD_REL_RES:
230  relres=absres/normu;
231  break;
232  default:
233  relres=absres/normr0;
234  break;
235  }
236 
237  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
238 
239  if ( relres < tol )
240  break;
241  else {
242  if ( stag >= MaxStag ) {
243  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
244  iter = ERROR_SOLVER_STAG;
245  break;
246  }
247  fasp_array_set(m,p,0.0);
248  ++stag;
249  ++restart_step;
250  }
251  } // end of staggnation check!
252 
253  // safe-guard check
254  if ( relres < tol ) {
255  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
256 
257  mf->fct(mf->data, u->val, r);
258  fasp_blas_array_axpby(m, 1.0, b->val, -1.0, r);
259 
260  // relative residual
261  switch (stop_type) {
262  case STOP_REL_PRECRES:
263  // z = B(r)
264  if (pc != NULL)
265  pc->fct(r,z,pc->data); /* Apply preconditioner */
266  else
267  fasp_array_cp(m,r,z); /* No preconditioner */
268  temp2=fasp_blas_array_dotprod(m,z,r);
269  relres=sqrt(ABS(temp2))/normr0;
270  break;
271  case STOP_MOD_REL_RES:
272  absres=fasp_blas_array_norm2(m,r);
273  relres=absres/normu;
274  break;
275  default:
276  absres=fasp_blas_array_norm2(m,r);
277  relres=absres/normr0;
278  break;
279  }
280 
281  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
282 
283  // check convergence
284  if ( relres < tol ) break;
285 
286  if ( more_step >= MaxRestartStep ) {
287  if ( prtlvl > PRINT_MIN) ITS_ZEROTOL;
288  iter = ERROR_SOLVER_TOLSMALL;
289  break;
290  }
291 
292  // prepare for restarting the method
293  fasp_array_set(m,p,0.0);
294  ++more_step;
295  ++restart_step;
296 
297  } // end of safe-guard check!
298 
299  // update relative residual here
300  absres0 = absres;
301 
302  // compute z_k = B(r_k)
303  if ( stop_type != STOP_REL_PRECRES ) {
304  if ( pc != NULL )
305  pc->fct(r,z,pc->data); /* Apply preconditioner */
306  else
307  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
308  }
309 
310  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
311  temp2=fasp_blas_array_dotprod(m,z,r);
312  beta=temp2/temp1;
313  temp1=temp2;
314 
315  // compute p_k = z_k + beta_k*p_{k-1}
316  fasp_blas_array_axpby(m,1.0,z,beta,p);
317 
318  } // end of main PCG loop.
319 
320 FINISHED: // finish the iterative method
321  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
322 
323  // clean up temp memory
324  fasp_mem_free(work);
325 
326 #if DEBUG_MODE > 0
327  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
328 #endif
329 
330  if (iter>MaxIt)
331  return ERROR_SOLVER_MAXIT;
332  else
333  return iter;
334 }
335 
336 /*---------------------------------*/
337 /*-- End of File --*/
338 /*---------------------------------*/
INT fasp_solver_pcg(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient (CG) method for solving Au=b.
Definition: pcg_mf.c:86
#define REAL
Definition: fasp.h:67
#define STOP_MOD_REL_RES
Definition: fasp_const.h:133
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 fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
REAL fasp_blas_dvec_norm2(dvector *x)
L2 norm of dvector x.
Definition: blas_vec.c:265
void * data
data for preconditioner, void pointer
Definition: fasp.h:1002
#define STOP_REL_PRECRES
Definition: fasp_const.h:132
Main header file for FASP.
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1013
#define PRINT_MORE
Definition: fasp_const.h:82
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
#define ERROR_SOLVER_TOLSMALL
Definition: fasp_const.h:51
#define MAX_RESTART
Definition: fasp_const.h:245
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
#define ERROR_SOLVER_SOLSTAG
Definition: fasp_const.h:50
#define ERROR_SOLVER_STAG
Definition: fasp_const.h:49
REAL fasp_blas_array_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: blas_array.c:388
INT row
number of rows
Definition: fasp.h:345
#define ABS(a)
Definition: fasp.h:74
#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
#define MAX_STAG
Definition: fasp_const.h:246
#define STAG_RATIO
Definition: fasp_const.h:247
void * data
data for MxV, can be a Matrix or something else
Definition: fasp.h:1016
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
#define SMALLREAL
Definition: fasp_const.h:238
#define PRINT_MIN
Definition: fasp_const.h:80