Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pgcg.c
Go to the documentation of this file.
1 
11 #include <math.h>
12 
13 #include "fasp.h"
14 #include "fasp_functs.h"
15 #include "itsolver_util.inl"
16 
17 /*---------------------------------*/
18 /*-- Public Functions --*/
19 /*---------------------------------*/
20 
45  dvector *b,
46  dvector *u,
47  precond *pc,
48  const REAL tol,
49  const INT MaxIt,
50  const SHORT stop_type,
51  const SHORT prtlvl)
52 {
53  INT iter=0, m=A->row, i;
54  REAL absres0 = BIGREAL, absres = BIGREAL;
55  REAL relres = BIGREAL, normb = BIGREAL;
56  REAL alpha, factor;
57 
58 
59  // allocate temp memory
60  REAL *work = (REAL *)fasp_mem_calloc(2*m+MaxIt+MaxIt*m,sizeof(REAL));
61 
62  REAL *r, *Br, *beta, *p;
63  r = work; Br = r + m; beta = Br + m; p = beta + MaxIt;
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  normb=fasp_blas_array_norm2(m,b->val);
71 
72  // -------------------------------------
73  // 1st iteration (Steepest descent)
74  // -------------------------------------
75  // r = b-A*u
76  fasp_array_cp(m,b->val,r);
77  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
78 
79  // Br
80  if (pc != NULL)
81  pc->fct(r,p,pc->data); /* Preconditioning */
82  else
83  fasp_array_cp(m,r,p); /* No preconditioner, B=I */
84 
85  // alpha = (p'r)/(p'Ap)
86  alpha = fasp_blas_array_dotprod (m,r,p) / fasp_blas_dcsr_vmv (A, p, p);
87 
88  // u = u + alpha *p
89  fasp_blas_array_axpy(m, alpha , p, u->val);
90 
91  // r = r - alpha *Ap
92  fasp_blas_dcsr_aAxpy((-1.0*alpha),A,p,r);
93 
94  // norm(r), factor
95  absres = fasp_blas_array_norm2(m,r); factor = absres/absres0;
96 
97  // compute relative residual
98  relres = absres/normb;
99 
100  // output iteration information if needed
101  print_itinfo(prtlvl,stop_type,iter+1,relres,absres,factor);
102 
103  // update relative residual here
104  absres0 = absres;
105 
106  for ( iter = 1; iter < MaxIt ; iter++) {
107 
108  // Br
109  if (pc != NULL)
110  pc->fct(r, Br ,pc->data); // Preconditioning
111  else
112  fasp_array_cp(m,r, Br); // No preconditioner, B=I
113 
114  // form p
115  fasp_array_cp(m, Br, p+iter*m);
116 
117  for (i=0; i<iter; i++) {
118  beta[i] = (-1.0) * ( fasp_blas_dcsr_vmv (A, Br, p+i*m)
119  /fasp_blas_dcsr_vmv (A, p+i*m, p+i*m) );
120 
121  fasp_blas_array_axpy(m, beta[i], p+i*m, p+iter*m);
122  }
123 
124  // -------------------------------------
125  // next iteration
126  // -------------------------------------
127 
128  // alpha = (p'r)/(p'Ap)
129  alpha = fasp_blas_array_dotprod(m,r,p+iter*m)
130  / fasp_blas_dcsr_vmv (A, p+iter*m, p+iter*m);
131 
132  // u = u + alpha *p
133  fasp_blas_array_axpy(m, alpha , p+iter*m, u->val);
134 
135  // r = r - alpha *Ap
136  fasp_blas_dcsr_aAxpy((-1.0*alpha),A,p+iter*m,r);
137 
138  // norm(r), factor
139  absres = fasp_blas_array_norm2(m,r); factor = absres/absres0;
140 
141  // compute relative residual
142  relres = absres/normb;
143 
144  // output iteration information if needed
145  print_itinfo(prtlvl,stop_type,iter+1,relres,absres,factor);
146 
147  if (relres < tol) break;
148 
149  // update relative residual here
150  absres0 = absres;
151 
152  } // end of main GCG loop.
153 
154  // finish the iterative method
155  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
156 
157  // clean up temp memory
158  fasp_mem_free(work);
159 
160 #if DEBUG_MODE > 0
161  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
162 #endif
163 
164  if (iter>MaxIt)
165  return ERROR_SOLVER_MAXIT;
166  else
167  return iter;
168 }
169 /*---------------------------------*/
170 /*-- End of File --*/
171 /*---------------------------------*/
#define REAL
Definition: fasp.h:67
INT fasp_solver_dcsr_pgcg(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned generilzed conjugate gradient (GCG) method for solving Au=b.
Definition: pgcg.c:44
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_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.
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
INT row
row number of matrix A, m
Definition: fasp.h:151
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 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
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
REAL fasp_blas_dcsr_vmv(dCSRmat *A, REAL *x, REAL *y)
vector-Matrix-vector multiplication alpha = y'*A*x
Definition: blas_csr.c:704