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