Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pvfgmres_mf.c
Go to the documentation of this file.
1 
17 #include <math.h>
18 
19 #include "fasp.h"
20 #include "fasp_functs.h"
21 #include "itsolver_util.inl"
22 
23 /*---------------------------------*/
24 /*-- Public Functions --*/
25 /*---------------------------------*/
26 
56  dvector *b,
57  dvector *x,
58  precond *pc,
59  const REAL tol,
60  const INT MaxIt,
61  const SHORT restart,
62  const SHORT stop_type,
63  const SHORT prtlvl)
64 {
65  const INT n = b->row;
66  const INT min_iter = 0;
67 
68  //--------------------------------------------//
69  // Newly added parameters to monitor when //
70  // to change the restart parameter //
71  //--------------------------------------------//
72  const REAL cr_max = 0.99; // = cos(8^o) (experimental)
73  const REAL cr_min = 0.174; // = cos(80^o) (experimental)
74 
75  // local variables
76  INT iter = 0;
77  INT i,j,k;
78 
79  REAL epsmac = SMALLREAL;
80  REAL r_norm, b_norm, den_norm;
81  REAL epsilon, gamma, t;
82 
83  REAL *c = NULL, *s = NULL, *rs = NULL;
84  REAL *norms = NULL, *r = NULL;
85  REAL **p = NULL, **hh = NULL, **z=NULL;
86  REAL *work = NULL;
87 
88  REAL cr = 1.0; // convergence rate
89  REAL r_norm_old = 0.0; // save the residual norm of the previous restart cycle
90  INT d = 3; // reduction for the restart parameter
91  INT restart_max = restart; // upper bound for restart in each restart cycle
92  INT restart_min = 3; // lower bound for restart in each restart cycle
93 
94  unsigned INT Restart = restart; // the real restart in some fixed restarted cycle
95  unsigned INT Restart1 = Restart + 1;
96  unsigned LONG worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
97 
98 #if DEBUG_MODE > 0
99  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
100  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
101 #endif
102 
103  /* allocate memory and setup temp work space */
104  work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
105 
106  /* check whether memory is enough for GMRES */
107  while ( (work == NULL) && (Restart > 5) ) {
108  Restart = Restart - 5;
109  worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
110  work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
111  Restart1 = Restart + 1;
112  }
113 
114  if ( work == NULL ) {
115  printf("### ERROR: No enough memory for vFGMRES %s : %s : %d!\n",
116  __FILE__, __FUNCTION__, __LINE__ );
117  exit(ERROR_ALLOC_MEM);
118  }
119 
120  if ( prtlvl > PRINT_MIN && Restart < restart ) {
121  printf("### WARNING: vFGMRES restart number set to %d!\n", Restart);
122  }
123 
124  p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
125  hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
126  z = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
127  norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
128 
129  r = work; rs = r + n; c = rs + Restart1; s = c + Restart;
130  for (i = 0; i < Restart1; i ++) p[i] = s + Restart + i*n;
131  for (i = 0; i < Restart1; i ++) hh[i] = p[Restart] + n + i*Restart;
132  for (i = 0; i < Restart1; i ++) z[i] = hh[Restart] + Restart + i*n;
133 
134  /* initialization */
135  mf->fct(mf->data, x->val, p[0]);
136  fasp_blas_array_axpby(n, 1.0, b->val, -1.0, p[0]);
137 
138  b_norm = fasp_blas_array_norm2(n, b->val);
139  r_norm = fasp_blas_array_norm2(n, p[0]);
140  norms[0] = r_norm;
141 
142  if ( prtlvl >= PRINT_SOME) {
143  ITS_PUTNORM("right-hand side", b_norm);
144  ITS_PUTNORM("residual", r_norm);
145  }
146 
147  if (b_norm > 0.0) den_norm = b_norm;
148  else den_norm = r_norm;
149 
150  epsilon = tol*den_norm;
151 
152  /* outer iteration cycle */
153  while (iter < MaxIt) {
154  rs[0] = r_norm;
155  r_norm_old = r_norm;
156  if (r_norm == 0.0) {
157  fasp_mem_free(work);
158  fasp_mem_free(p);
159  fasp_mem_free(hh);
160  fasp_mem_free(norms);
161  fasp_mem_free(z);
162  return iter;
163  }
164 
165  //-----------------------------------//
166  // adjust the restart parameter //
167  //-----------------------------------//
168 
169  if (cr > cr_max || iter == 0) {
170  Restart = restart_max;
171  }
172  else if (cr < cr_min) {
173  // Restart = Restart;
174  }
175  else {
176  if (Restart - d > restart_min) {
177  Restart -= d;
178  }
179  else {
180  Restart = restart_max;
181  }
182  }
183 
184  if (r_norm <= epsilon && iter >= min_iter) {
185  mf->fct(mf->data, x->val, r);
186  fasp_blas_array_axpby(n, 1.0, b->val, -1.0, r);
187  r_norm = fasp_blas_array_norm2(n, r);
188 
189  if (r_norm <= epsilon) {
190  break;
191  }
192  else {
193  if (prtlvl >= PRINT_SOME) ITS_FACONV;
194  }
195  }
196 
197  t = 1.0 / r_norm;
198  //for (j = 0; j < n; j ++) p[0][j] *= t;
199  fasp_blas_array_ax(n, t, p[0]);
200 
201  /* RESTART CYCLE (right-preconditioning) */
202  i = 0;
203  while (i < Restart && iter < MaxIt) {
204 
205  i ++; iter ++;
206 
207  /* apply the preconditioner */
208  if (pc == NULL)
209  fasp_array_cp(n, p[i-1], z[i-1]);
210  else
211  pc->fct(p[i-1], z[i-1], pc->data);
212 
213  mf->fct(mf->data, z[i-1], p[i]);
214 
215  /* modified Gram_Schmidt */
216  for (j = 0; j < i; j ++) {
217  hh[j][i-1] = fasp_blas_array_dotprod(n, p[j], p[i]);
218  fasp_blas_array_axpy(n, -hh[j][i-1], p[j], p[i]);
219  }
220  t = fasp_blas_array_norm2(n, p[i]);
221  hh[i][i-1] = t;
222  if (t != 0.0) {
223  t = 1.0/t;
224  //for (j = 0; j < n; j ++) p[i][j] *= t;
225  fasp_blas_array_ax(n, t, p[i]);
226  }
227 
228  for (j = 1; j < i; ++j) {
229  t = hh[j-1][i-1];
230  hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
231  hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
232  }
233  t= hh[i][i-1]*hh[i][i-1];
234  t+= hh[i-1][i-1]*hh[i-1][i-1];
235  gamma = sqrt(t);
236  if (gamma == 0.0) gamma = epsmac;
237  c[i-1] = hh[i-1][i-1] / gamma;
238  s[i-1] = hh[i][i-1] / gamma;
239  rs[i] = -s[i-1]*rs[i-1];
240  rs[i-1] = c[i-1]*rs[i-1];
241  hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
242 
243  r_norm = fabs(rs[i]);
244  norms[iter] = r_norm;
245 
246  if (b_norm > 0 ) {
247  print_itinfo(prtlvl,stop_type,iter,norms[iter]/b_norm,
248  norms[iter],norms[iter]/norms[iter-1]);
249  }
250  else {
251  print_itinfo(prtlvl,stop_type,iter,norms[iter],norms[iter],
252  norms[iter]/norms[iter-1]);
253  }
254 
255  /* should we exit the restart cycle? */
256  if (r_norm <= epsilon && iter >= min_iter) break;
257 
258  } /* end of restart cycle */
259 
260  /* now compute solution, first solve upper triangular system */
261 
262  rs[i-1] = rs[i-1] / hh[i-1][i-1];
263  for (k = i-2; k >= 0; k --) {
264  t = 0.0;
265  for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
266 
267  t += rs[k];
268  rs[k] = t / hh[k][k];
269  }
270 
271  fasp_array_cp(n, z[i-1], r);
272  //for (j = 0; j < n; j ++) r[j] *= rs[i-1];
273  fasp_blas_array_ax(n, rs[i-1], r);
274  for (j = i-2; j >= 0; j --) fasp_blas_array_axpy(n, rs[j], z[j], r);
275 
276  fasp_blas_array_axpy(n, 1.0, r, x->val);
277 
278  if (r_norm <= epsilon && iter >= min_iter) {
279  mf->fct(mf->data, x->val, r);
280  fasp_blas_array_axpby(n, 1.0, b->val, -1.0, r);
281  r_norm = fasp_blas_array_norm2(n, r);
282 
283  if (r_norm <= epsilon) {
284  break;
285  }
286  else {
287  if (prtlvl >= PRINT_SOME) ITS_FACONV;
288  fasp_array_cp(n, r, p[0]); i = 0;
289  }
290  } /* end of convergence check */
291 
292 
293  /* compute residual vector and continue loop */
294  for (j = i; j > 0; j--) {
295  rs[j-1] = -s[j-1]*rs[j];
296  rs[j] = c[j-1]*rs[j];
297  }
298 
299  if (i) fasp_blas_array_axpy(n, rs[i]-1.0, p[i], p[i]);
300 
301  for (j = i-1 ; j > 0; j --) fasp_blas_array_axpy(n, rs[j], p[j], p[i]);
302 
303  if (i) {
304  fasp_blas_array_axpy(n, rs[0]-1.0, p[0], p[0]);
305  fasp_blas_array_axpy(n, 1.0, p[i], p[0]);
306  }
307 
308  //-----------------------------------//
309  // compute the convergence rate //
310  //-----------------------------------//
311  cr = r_norm / r_norm_old;
312 
313  } /* end of iteration while loop */
314 
315  if (prtlvl > PRINT_NONE) ITS_FINAL(iter,MaxIt,r_norm);
316 
317  /*-------------------------------------------
318  * Free some stuff
319  *------------------------------------------*/
320  fasp_mem_free(work);
321  fasp_mem_free(p);
322  fasp_mem_free(hh);
323  fasp_mem_free(norms);
324  fasp_mem_free(z);
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 /*---------------------------------*/
#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 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
Main header file for FASP.
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1013
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:37
#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 row
number of rows
Definition: fasp.h:345
#define LONG
Definition: fasp.h:65
INT fasp_solver_pvfgmres(mxv_matfree *mf, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
Solve "Ax=b" using PFGMRES(right preconditioned) iterative method in which the restart parameter can ...
Definition: pvfgmres_mf.c:55
#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
void fasp_blas_array_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: blas_array.c:35
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
#define SMALLREAL
Definition: fasp_const.h:238
#define PRINT_MIN
Definition: fasp_const.h:80