Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pgcr.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #include "fasp.h"
9 #include "fasp_functs.h"
10 #include "itsolver_util.inl"
11 
12 static void dense_aAtxpby (INT, INT, REAL *, REAL, REAL *, REAL, REAL *);
13 
38  dvector *b,
39  dvector *x,
40  precond *pc,
41  const REAL tol,
42  const INT MaxIt,
43  const SHORT restart,
44  const SHORT stop_type,
45  const SHORT prtlvl)
46 {
47  const INT n = b->row;
48 
49  // local variables
50  INT iter = 0, rst = -1;
51  INT i, j, k;
52 
53  REAL gamma, alpha, beta, checktol;
54  REAL absres0 = BIGREAL, absres = BIGREAL;
55  REAL relres = BIGREAL;
56 
57  // allocate temp memory (need about (restart+4)*n REAL numbers)
58  REAL *c = NULL, *z = NULL, *alp = NULL, *tmpx = NULL;
59  REAL *norms = NULL, *r = NULL, *work = NULL;
60  REAL **h = NULL;
61 
62  INT Restart = MIN(restart, MaxIt);
63  LONG worksize = n+2*Restart*n+Restart+Restart;
64 
65 #if DEBUG_MODE > 0
66  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
67  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
68 #endif
69 
70  work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
71 
72  /* check whether memory is enough for GMRES */
73  while ( (work == NULL) && (Restart > 5) ) {
74  Restart = Restart - 5;
75  worksize = n+2*Restart*n+Restart+Restart;
76  work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
77  }
78 
79  if ( work == NULL ) {
80  printf("### ERROR: No enough memory for GMRES %s : %s : %d!\n",
81  __FILE__, __FUNCTION__, __LINE__ );
82  exit(ERROR_ALLOC_MEM);
83  }
84 
85  if ( prtlvl > PRINT_MIN && Restart < restart ) {
86  printf("### WARNING: GMRES restart number set to %d!\n", Restart);
87  }
88 
89  r = work, z = r+n, c = z + Restart*n, alp = c + Restart*n, tmpx = alp + Restart;
90 
91  h = (REAL **)fasp_mem_calloc(Restart, sizeof(REAL *));
92  for (i = 0; i < Restart; i++) h[i] = (REAL*)fasp_mem_calloc(Restart, sizeof(REAL));
93 
94  norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
95 
96  // r = b-A*x
97  fasp_array_cp(n, b->val, r);
98  fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
99 
100  absres = fasp_blas_array_dotprod(n, r, r);
101 
102  absres0 = MAX(SMALLREAL,absres);
103 
104  relres = absres/absres0;
105 
106  // output iteration information if needed
107  print_itinfo(prtlvl,stop_type,0,relres,sqrt(absres0),0.0);
108 
109  // store initial residual
110  norms[0] = relres;
111 
112  checktol = MAX(tol*tol*absres0, absres*1.0e-4);
113 
114  while ( iter < MaxIt && sqrt(relres) > tol ) {
115 
116  i = -1; rst ++;
117 
118  while ( i < Restart-1 && iter < MaxIt ) {
119 
120  i++; iter++;
121 
122  // z = B^-1r
123  if ( pc == NULL )
124  fasp_array_cp(n, r, &z[i*n]);
125  else
126  pc->fct(r, &z[i*n], pc->data);
127 
128  // c = Az
129  fasp_blas_dcsr_mxv(A, &z[i*n], &c[i*n]);
130 
131  /* Modified Gram_Schmidt orthogonalization */
132  for ( j = 0; j < i; j++ ) {
133  gamma = fasp_blas_array_dotprod(n, &c[j*n], &c[i*n]);
134  h[i][j] = gamma/h[j][j];
135  fasp_blas_array_axpy(n, -h[i][j], &c[j*n], &c[i*n]);
136  }
137  // gamma = (c,c)
138  gamma = fasp_blas_array_dotprod(n, &c[i*n], &c[i*n]);
139 
140  h[i][i] = gamma;
141 
142  // alpha = (c, r)
143  alpha = fasp_blas_array_dotprod(n, &c[i*n], r);
144 
145  beta = alpha/gamma;
146 
147  alp[i] = beta;
148 
149  // r = r - beta*c
150  fasp_blas_array_axpy(n, -beta, &c[i*n], r);
151 
152  // equivalent to ||r||_2
153  absres = absres - alpha*alpha/gamma;
154 
155  if (absres < checktol) {
156  absres = fasp_blas_array_dotprod(n, r, r);
157  checktol = MAX(tol*tol*absres0, absres*1.0e-4);
158  }
159 
160  relres = absres / absres0;
161 
162  norms[iter] = relres;
163 
164  print_itinfo(prtlvl, stop_type, iter, sqrt(relres), sqrt(absres),
165  sqrt(norms[iter]/norms[iter-1]));
166 
167  if (sqrt(relres) < tol) break;
168  }
169 
170  for ( k = i; k >=0; k-- ) {
171  tmpx[k] = alp[k];
172  for (j=0; j<k; ++j) {
173  alp[j] -= h[k][j]*tmpx[k];
174  }
175  }
176 
177  if (rst==0) dense_aAtxpby(n, i+1, z, 1.0, tmpx, 0.0, x->val);
178  else dense_aAtxpby(n, i+1, z, 1.0, tmpx, 1.0, x->val);
179 
180  }
181 
182  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,sqrt(relres));
183 
184  // clean up memory
185  for (i = 0; i < Restart; i++) fasp_mem_free(h[i]);
186  fasp_mem_free(h);
187 
188  fasp_mem_free(work);
189  fasp_mem_free(norms);
190 
191 #if DEBUG_MODE > 0
192  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
193 #endif
194 
195  if ( iter >= MaxIt )
196  return ERROR_SOLVER_MAXIT;
197  else
198  return iter;
199 }
200 
227  dvector *b,
228  dvector *x,
229  precond *pc,
230  const REAL tol,
231  const INT MaxIt,
232  const SHORT restart,
233  const SHORT stop_type,
234  const SHORT prtlvl)
235 {
236  INT i, j;
237  INT iter = 0;
238  INT m = restart;
239  REAL alpha, beta, gamma, tempe, tempb=.0, tempu, temp2;
240  REAL absres0 = BIGREAL, absres, relres1 = BIGREAL, infnormu, factor;
241 
242  const INT nrow = b->row;
243  const REAL sol_inf_tol = 1e-16;
244 
245  // default restart number
246  if ((m < 1)||(m > nrow)||(m >150)) m=10;
247 
248  dvector *v = (dvector *)fasp_mem_calloc(m,sizeof(dvector));
249 
250  for (i=0; i<=m-1; ++i) {
251  v[i].row = nrow;
252  v[i].val = (REAL*)fasp_mem_calloc(nrow,sizeof(REAL));
253  }
254 
255  dvector *s = (dvector *)fasp_mem_calloc(m,sizeof(dvector));
256 
257  for (i=0; i<=m-1;++i) {
258  s[i].row = nrow;
259  s[i].val = (REAL*)fasp_mem_calloc(nrow,sizeof(REAL));
260  }
261 
262  REAL *r = (REAL*)fasp_mem_calloc(nrow,sizeof(REAL));
263 
264  // compute norm for right hand side
265  switch (stop_type) {
266  case STOP_REL_RES:
267  tempb = fasp_blas_array_norm2(nrow,b->val);
268  break;
269 
270  case STOP_REL_PRECRES:
271  if (pc != NULL)
272  pc->fct(b->val,s[0].val,pc->data);
273  else
274  fasp_array_cp(nrow,b->val,s[0].val);
275  tempb = sqrt(ABS(fasp_blas_array_dotprod(nrow,b->val,s[0].val)));
276  break;
277 
278  case STOP_MOD_REL_RES:
279  break;
280 
281  default:
282  printf("### ERROR: Unknown stopping criteria!\n");
283  iter = ERROR_INPUT_PAR;
284  goto FINISHED;
285  }
286 
287  tempb = MAX(SMALLREAL,tempb);
288  tempu = MAX(SMALLREAL,fasp_blas_array_norm2(nrow,x->val));
289 
290  // r = b-A*u
291  fasp_array_cp(nrow, b->val, r);
292  fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
293  tempe = fasp_blas_array_norm2(nrow, r);
294  tempb = MAX(SMALLREAL, tempe);
295 
296  switch (stop_type) {
297  case STOP_REL_RES:
298  relres1 = tempe/tempb;
299  break;
300 
301  case STOP_REL_PRECRES:
302  if (pc == NULL)
303  fasp_array_cp(nrow, r, s[0].val);
304  else
305  pc->fct(r, s[0].val, pc->data);
306  temp2 = sqrt(ABS(fasp_blas_array_dotprod(nrow, r, s[0].val)));
307  relres1 = temp2/tempb;
308  break;
309 
310  default: // case STOP_MOD_REL_RES
311  relres1 = tempe/tempu;
312  break;
313  }
314 
315  if (relres1<tol) { fasp_mem_free(r); goto FINISHED; }
316 
317  if (iter < 0) goto FINISHED;
318 
319  while (iter++ < MaxIt)
320  {
321  for (j=0; j<m; ++j)
322  {
323  if (pc == NULL) {
324  fasp_array_cp(nrow, r, s[j].val);
325  }
326  else {
327  pc->fct(r, s[j].val, pc->data);
328  }
329 
330  fasp_blas_dcsr_aAxpy(1.0, A, s[j].val, v[j].val);
331 
332  for (i=0; i<j; ++i)
333  {
334  alpha = fasp_blas_array_dotprod(nrow, v[j].val, v[i].val);
335  fasp_blas_array_axpy(nrow, -alpha, v[i].val, v[j].val);
336  fasp_blas_array_axpy(nrow, -alpha, s[i].val, s[j].val);
337  }
338 
339  beta = fasp_blas_array_norm2(nrow, v[j].val);
340  fasp_blas_array_ax(nrow, 1.0/beta, v[j].val);
341  fasp_blas_array_ax(nrow, 1.0/beta, s[j].val);
342 
343  gamma = fasp_blas_array_dotprod(nrow, v[j].val, r);
344  fasp_blas_array_axpy(nrow, gamma, s[j].val, x->val);
345 
346  // r = b-A*u
347  fasp_array_cp(nrow, b->val, r);
348  fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
349 
350  // absolute and relative residuals
351  absres = sqrt(fasp_blas_array_dotprod(nrow, r, r));
352  tempu = sqrt(fasp_blas_dvec_dotprod(x, x));
353 
354  switch (stop_type) {
355  case STOP_REL_RES:
356  relres1 = absres/tempb;
357  break;
358 
359  case STOP_REL_PRECRES:
360  if (pc == NULL)
361  fasp_array_cp(nrow, r, s[j].val);
362  else
363  pc->fct(r, s[j].val, pc->data);
364  temp2 = sqrt(ABS(fasp_blas_array_dotprod(nrow, r, s[j].val)));
365  relres1 = temp2/tempb;
366  break;
367 
368  default: // case STOP_MOD_REL_RES
369  relres1 = absres/tempu;
370  break;
371  }
372 
373  // contraction factor
374  factor = absres/absres0;
375 
376  // output iteration information if needed
377  print_itinfo(prtlvl, stop_type, (iter-1)*m+j+1, relres1, absres, factor);
378 
379  absres0 = absres;
380 
381  // solution check, if soultion is too small, return ERROR_SOLVER_SOLSTAG.
382  infnormu = fasp_blas_array_norminf(nrow, x->val);
383 
384  if (infnormu <= sol_inf_tol)
385  {
386  print_message(prtlvl, "Inf_norm of the solution is too small!\n");
387  iter = ERROR_SOLVER_SOLSTAG;
388  goto FINISHED;
389  }
390 
391  if (relres1<tol) goto FINISHED;
392  }
393  }
394 
395 FINISHED:
396 
397  if (prtlvl > 0) {
398  if (iter > MaxIt){
399  printf("Maximal iteration %d exceeded with relative residual %e.\n", MaxIt, relres1);
400  iter = ERROR_SOLVER_MAXIT;
401  }
402  else
403  printf("Number of iterations = %d with relative residual %e.\n", iter, relres1);
404  }
405 
406  fasp_mem_free(r);
407  for (i=0; i<=m-1; ++i) fasp_mem_free(v[i].val);
408  fasp_mem_free(v);
409  for (i=0; i<=m-1; ++i) fasp_mem_free(s[i].val);
410  fasp_mem_free(s);
411 
412  return iter;
413 }
414 
415 /*---------------------------------*/
416 /*-- Private Functions --*/
417 /*---------------------------------*/
418 
438 static void dense_aAtxpby (INT n,
439  INT m,
440  REAL *A,
441  REAL alpha,
442  REAL *x,
443  REAL beta,
444  REAL *y)
445 {
446  INT i, j;
447 
448  for (i=0; i<m; i++) fasp_blas_array_ax(n, x[i], &A[i*n]);
449 
450  for (j=1; j<m; j++) {
451  for (i=0; i<n; i++) {
452  A[i] += A[i+j*n];
453  }
454  }
455 
456  fasp_blas_array_axpby(n, alpha, A, beta, y);
457 }
458 
459 /*---------------------------------*/
460 /*-- End of File --*/
461 /*---------------------------------*/
#define ERROR_INPUT_PAR
Definition: fasp_const.h:31
#define REAL
Definition: fasp.h:67
#define STOP_MOD_REL_RES
Definition: fasp_const.h:133
INT fasp_solver_dcsr_pgcr1(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
A preconditioned GCR method for solving Au=b.
Definition: pgcr.c:226
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
#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 * data
data for preconditioner, void pointer
Definition: fasp.h:1002
#define STOP_REL_PRECRES
Definition: fasp_const.h:132
REAL fasp_blas_dvec_dotprod(dvector *x, dvector *y)
Inner product of two vectors (x,y)
Definition: blas_vec.c:121
Main header file for FASP.
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:37
INT fasp_solver_dcsr_pgcr(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
A preconditioned GCR method for solving Au=b.
Definition: pgcr.c:37
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
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 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 LONG
Definition: fasp.h:65
#define ABS(a)
Definition: fasp.h:74
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define MIN(a, b)
Definition: fasp.h:73
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 print_message(const INT ptrlvl, const char *message)
Print output information if necessary.
Definition: message.c:182
#define BIGREAL
Some global constants.
Definition: fasp_const.h:237
void fasp_blas_array_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: blas_array.c:35
void fasp_blas_dcsr_aAxpy(const REAL alpha, dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_csr.c:479
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(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
void fasp_blas_dcsr_mxv(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_csr.c:225