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