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