Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pbcgs_mf.c
Go to the documentation of this file.
1 
56 #include <math.h>
57 
58 #include "fasp.h"
59 #include "fasp_functs.h"
60 #include "itsolver_util.inl"
61 
62 /*---------------------------------*/
63 /*-- Public Functions --*/
64 /*---------------------------------*/
65 
92  dvector *b,
93  dvector *u,
94  precond *pc,
95  const REAL tol,
96  const INT MaxIt,
97  const SHORT stop_type,
98  const SHORT prtlvl)
99 {
100  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
101  const INT m=b->row;
102  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
103  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
104  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
105 
106  // local variables
107  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
108  REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
109  REAL normd = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
110  REAL reldiff, factor, infnormu;
111  REAL alpha, beta, omega, temp1, temp2, tempr;
112 
113  REAL *uval=u->val, *bval=b->val;
114 
115  // allocate temp memory (need 8*m REAL)
116  REAL *work=(REAL *)fasp_mem_calloc(8*m,sizeof(REAL));
117  REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
118  REAL *rho=t+m, *pp=rho+m, *s=pp+m, *sp=s+m;
119 
120 #if DEBUG_MODE > 0
121  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
122  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
123 #endif
124 
125  // initialize counters
126  stag=more_step=restart_step=1;
127 
128  // r = b-A*u
129  mf->fct(mf->data, uval, r);
130  fasp_blas_array_axpby(m, 1.0, bval, -1.0, r);
131 
132  // pp=precond(p)
133  fasp_array_cp(m,r,p);
134  if (pc != NULL)
135  pc->fct(p,pp,pc->data); /* Apply preconditioner */
136  else
137  fasp_array_cp(m,p,pp); /* No preconditioner */
138 
139  // compute initial relative residual
140  switch (stop_type) {
141  case STOP_REL_PRECRES:
142  absres0=sqrt(ABS(fasp_blas_array_dotprod(m,r,pp)));
143  normr0=MAX(SMALLREAL,absres0);
144  relres=absres0/normr0;
145  break;
146  case STOP_MOD_REL_RES:
147  absres0=fasp_blas_array_norm2(m,r);
148  normu=MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
149  relres=absres0/normu;
150  break;
151  default:
152  absres0=fasp_blas_array_norm2(m,r);
153  normr0=MAX(SMALLREAL,absres0);
154  relres=absres0/normr0;
155  break;
156  }
157 
158  // if initial residual is small, no need to iterate!
159  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
160 
161  // rho = r* := r
162  fasp_array_cp(m,r,rho);
163  temp1=fasp_blas_array_dotprod(m,r,rho);
164 
165  while ( iter++ < MaxIt ) {
166 
167  // z = A*pp
168  mf->fct(mf->data, pp, z);
169 
170  // alpha = (r,rho)/(A*p,rho)
171  temp2=fasp_blas_array_dotprod(m,z,rho);
172  if (ABS(temp2)>SMALLREAL) {
173  alpha=temp1/temp2;
174  }
175  else {
176  ITS_DIVZERO;
177  return ERROR_SOLVER_MISC;
178  }
179 
180  // s = r - alpha z
181  fasp_array_cp(m,r,s);
182  fasp_blas_array_axpy(m,-alpha,z,s);
183 
184  // sp = precond(s)
185  if (pc != NULL)
186  pc->fct(s,sp,pc->data); /* Apply preconditioner */
187  else
188  fasp_array_cp(m,s,sp); /* No preconditioner */
189 
190  // t = A*sp;
191  mf->fct(mf->data, sp, t);
192 
193  // omega = (t,s)/(t,t)
194  tempr=fasp_blas_array_dotprod(m,t,t);
195 
196  if (ABS(tempr)>SMALLREAL) {
197  omega=fasp_blas_array_dotprod(m,s,t)/tempr;
198  }
199  else {
200  if ( prtlvl >= PRINT_SOME ) ITS_DIVZERO;
201  omega=0.0;
202  }
203 
204  // delu = alpha pp + omega sp
205  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
206 
207  // u = u + delu
208  fasp_blas_array_axpy(m,1.0,sp,uval);
209 
210  // r = s - omega t
211  fasp_blas_array_axpy(m,-omega,t,s);
212  fasp_array_cp(m,s,r);
213 
214  // beta = (r,rho)/(rp,rho)
215  temp2=temp1;
216  temp1=fasp_blas_array_dotprod(m,r,rho);
217 
218  if (ABS(temp2)>SMALLREAL) {
219  beta=(temp1*alpha)/(temp2*omega);
220  }
221  else {
222  ITS_DIVZERO;
223  return ERROR_SOLVER_MISC;
224  }
225 
226  // p = p - omega z
227  fasp_blas_array_axpy(m,-omega,z,p);
228 
229  // p = r + beta p
230  fasp_blas_array_axpby(m,1.0,r,beta,p);
231 
232  // pp = precond(p)
233  if (pc != NULL)
234  pc->fct(p,pp,pc->data); /* Apply preconditioner */
235  else
236  fasp_array_cp(m,p,pp); /* No preconditioner */
237 
238  // compute reducation factor of residual ||r||
239  absres=fasp_blas_array_norm2(m,r);
240  factor=absres/absres0;
241 
242  // relative difference
243  normd = fasp_blas_array_norm2(m,sp);
244  normu = fasp_blas_array_norm2(m,uval);
245  reldiff = normd/normu;
246 
247  if ( normd<TOL_s ) {
248  ITS_SMALLSP; goto FINISHED;
249  }
250 
251  // relative residual
252  switch (stop_type) {
253  case STOP_REL_PRECRES:
254  if (pc == NULL) fasp_array_cp(m,r,z);
255  else pc->fct(r,z,pc->data);
256  tempr=sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
257  relres=tempr/normr0;
258  break;
259  case STOP_MOD_REL_RES:
260  relres=absres/normu;
261  break;
262  default:
263  relres=absres/normr0;
264  break;
265  }
266 
267  // output iteration information if needed
268  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
269 
270  // solution check, if soultion is too small, return ERROR_SOLVER_SOLSTAG.
271  infnormu = fasp_blas_array_norminf(m, uval);
272  if ( infnormu <= sol_inf_tol ) {
273  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
274  iter = ERROR_SOLVER_SOLSTAG;
275  goto FINISHED;
276  }
277 
278  // stagnation check
279  if ( (stag<=MaxStag) && (reldiff<maxdiff) ) {
280 
281  if ( prtlvl >= PRINT_MORE ) {
282  ITS_DIFFRES(reldiff,relres);
283  ITS_RESTART;
284  }
285 
286  // re-init iteration param
287  mf->fct(mf->data, uval, r);
288  fasp_blas_array_axpby(m, 1.0, bval, -1.0, r);
289 
290  // pp=precond(p)
291  fasp_array_cp(m,r,p);
292  if (pc != NULL)
293  pc->fct(p,pp,pc->data); /* Apply preconditioner */
294  else
295  fasp_array_cp(m,p,pp); /* No preconditioner */
296  // rho = r* := r
297  fasp_array_cp(m,r,rho);
298  temp1=fasp_blas_array_dotprod(m,r,rho);
299  absres=fasp_blas_array_norm2(m,r);
300 
301  // relative residual
302  switch (stop_type) {
303  case STOP_REL_PRECRES:
304  if (pc != NULL)
305  pc->fct(r,z,pc->data);
306  else
307  fasp_array_cp(m,r,z);
308  tempr=sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
309  relres=tempr/normr0;
310  break;
311  case STOP_MOD_REL_RES:
312  relres=absres/normu;
313  break;
314  default:
315  relres=absres/normr0;
316  break;
317  }
318 
319  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
320 
321  if ( relres < tol )
322  break;
323  else {
324  if ( stag >= MaxStag ) {
325  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
326  iter = ERROR_SOLVER_STAG;
327  goto FINISHED;
328  }
329  ++stag;
330  ++restart_step;
331  }
332 
333  } // end of stagnation check!
334 
335  // safe-guard check
336  if ( relres < tol ) {
337  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
338 
339  // re-init iteration param
340  mf->fct(mf->data, uval, r);
341  fasp_blas_array_axpby(m, 1.0, bval, -1.0, r);
342 
343  // pp=precond(p)
344  fasp_array_cp(m,r,p);
345  if (pc != NULL)
346  pc->fct(p,pp,pc->data); /* Apply preconditioner */
347  else
348  fasp_array_cp(m,p,pp); /* No preconditioner */
349  // rho = r* := r
350  fasp_array_cp(m,r,rho);
351  temp1=fasp_blas_array_dotprod(m,r,rho);
352  absres=fasp_blas_array_norm2(m,r);
353 
354  // relative residual
355  switch (stop_type) {
356  case STOP_REL_PRECRES:
357  if (pc != NULL)
358  pc->fct(r,z,pc->data);
359  else
360  fasp_array_cp(m,r,z);
361  tempr=sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
362  relres=tempr/normr0;
363  break;
364  case STOP_MOD_REL_RES:
365  relres=absres/normu;
366  break;
367  default:
368  relres=absres/normr0;
369  break;
370  }
371 
372  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
373 
374  // check convergence
375  if ( relres < tol ) break;
376 
377  if ( more_step >= MaxRestartStep ) {
378  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
379  iter = ERROR_SOLVER_TOLSMALL;
380  goto FINISHED;
381  }
382  else {
383  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
384  }
385 
386  ++more_step;
387  ++restart_step;
388  } // end if safe guard
389 
390  absres0=absres;
391  }
392 
393 FINISHED: // finish the iterative method
394  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
395 
396  // clean up temp memory
397  fasp_mem_free(work);
398 
399 #if DEBUG_MODE > 0
400  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
401 #endif
402 
403  if (iter>MaxIt)
404  return ERROR_SOLVER_MAXIT;
405  else
406  return iter;
407 }
408 
409 /*---------------------------------*/
410 /*-- End of File --*/
411 /*---------------------------------*/
#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_MISC
Definition: fasp_const.h:53
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:54
#define PRINT_SOME
Definition: fasp_const.h:81
#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
#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
INT fasp_solver_pbcgs(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned BiCGstab method for solving Au=b.
Definition: pbcgs_mf.c:91
#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