Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
itsolver_mf.c
Go to the documentation of this file.
1 
6 #include <time.h>
7 
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11 
12 #include "fasp.h"
13 #include "fasp_functs.h"
14 #include "fasp_block.h"
15 
16 static void fasp_blas_mxv_csr (void *A, REAL *x, REAL *y);
17 static void fasp_blas_mxv_bsr (void *A, REAL *x, REAL *y);
18 static void fasp_blas_mxv_str (void *A, REAL *x, REAL *y);
19 static void fasp_blas_mxv_bcsr (void *A, REAL *x, REAL *y);
20 static void fasp_blas_mxv_bbsr (void *A, REAL *x, REAL *y);
21 static void fasp_blas_mxv_csrl (void *A, REAL *x, REAL *y);
22 
23 #include "itsolver_util.inl"
24 
25 /*---------------------------------*/
26 /*-- Public Functions --*/
27 /*---------------------------------*/
28 
51  dvector *b,
52  dvector *x,
53  precond *pc,
54  itsolver_param *itparam)
55 {
56  const SHORT prtlvl = itparam->print_level;
57  const SHORT itsolver_type = itparam->itsolver_type;
58  const SHORT stop_type = itparam->stop_type;
59  const SHORT restart = itparam->restart;
60  const INT MaxIt = itparam->maxit;
61  const REAL tol = itparam->tol;
62 
63  /* Local Variables */
64  REAL solver_start, solver_end, solver_duration;
65  INT iter = ERROR_SOLVER_TYPE;
66 
67  fasp_gettime(&solver_start);
68 
69 #if DEBUG_MODE > 0
70  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
71  printf("### DEBUG: rhs/sol size: %d %d\n", b->row, x->row);
72 #endif
73 
74  /* Safe-guard checks on parameters */
75  ITS_CHECK ( MaxIt, tol );
76 
77  /* Choose a desirable Krylov iterative solver */
78  switch ( itsolver_type ) {
79  case SOLVER_CG:
80  if ( prtlvl > PRINT_NONE ) printf("\nCalling PCG solver (matrix-free) ...\n");
81  iter = fasp_solver_pcg(mf, b, x, pc, tol, MaxIt, stop_type, prtlvl);
82  break;
83 
84  case SOLVER_BiCGstab:
85  if ( prtlvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (matrix-free) ...\n");
86  iter = fasp_solver_pbcgs(mf, b, x, pc, tol, MaxIt, stop_type, prtlvl);
87  break;
88 
89  case SOLVER_MinRes:
90  if ( prtlvl > PRINT_NONE ) printf("\nCalling MinRes solver (matrix-free) ...\n");
91  iter = fasp_solver_pminres(mf, b, x, pc, tol, MaxIt, stop_type, prtlvl);
92  break;
93 
94  case SOLVER_GMRES:
95  if ( prtlvl > PRINT_NONE ) printf("\nCalling GMRes solver (matrix-free) ...\n");
96  iter = fasp_solver_pgmres(mf, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
97  break;
98 
99  case SOLVER_VGMRES:
100  if ( prtlvl > PRINT_NONE ) printf("\nCalling vGMRes solver (matrix-free) ...\n");
101  iter = fasp_solver_pvgmres(mf, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
102  break;
103 
104  case SOLVER_VFGMRES:
105  if ( prtlvl > PRINT_NONE ) printf("\nCalling vFGMRes solver (matrix-free) ...\n");
106  iter = fasp_solver_pvfgmres(mf, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
107  break;
108 
109  case SOLVER_GCG:
110  if ( prtlvl > PRINT_NONE ) printf("\nCalling GCG solver (matrix-free) ...\n");
111  iter = fasp_solver_pgcg(mf, b, x, pc, tol, MaxIt, stop_type, prtlvl);
112  break;
113 
114  default:
115  printf("### ERROR: Unknown itertive solver type %d!\n", itsolver_type);
116 
117  }
118 
119  if ( (prtlvl >= PRINT_SOME) && (iter >= 0) ) {
120  fasp_gettime(&solver_end);
121  solver_duration = solver_end - solver_start;
122  print_cputime("Iterative method", solver_duration);
123  }
124 
125 #if DEBUG_MODE > 0
126  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
127 #endif
128 
129  return iter;
130 }
131 
151  dvector *b,
152  dvector *x,
153  itsolver_param *itparam)
154 {
155  const SHORT prtlvl = itparam->print_level;
156 
157  /* Local Variables */
158  INT status = FASP_SUCCESS;
159  REAL solver_start, solver_end, solver_duration;
160 
161 #if DEBUG_MODE > 0
162  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
163  printf("### DEBUG: rhs/sol size: %d %d\n", b->row, x->row);
164 #endif
165 
166  fasp_gettime(&solver_start);
167 
168  status = fasp_solver_itsolver(mf,b,x,NULL,itparam);
169 
170  if ( prtlvl >= PRINT_MIN ) {
171  fasp_gettime(&solver_end);
172  solver_duration = solver_end - solver_start;
173  print_cputime("Krylov method totally", solver_duration);
174  }
175 
176 #if DEBUG_MODE > 0
177  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
178 #endif
179 
180  return status;
181 }
182 
197 void fasp_solver_itsolver_init (INT matrix_format,
198  mxv_matfree *mf,
199  void *A)
200 {
201  switch ( matrix_format ) {
202 
203  case MAT_CSR:
204  mf->fct = fasp_blas_mxv_csr;
205  break;
206 
207  case MAT_BSR:
208  mf->fct = fasp_blas_mxv_bsr;
209  break;
210 
211  case MAT_STR:
212  mf->fct = fasp_blas_mxv_str;
213  break;
214 
215  case MAT_bCSR:
216  mf->fct = fasp_blas_mxv_bcsr;
217  break;
218 
219  case MAT_bBSR:
220  mf->fct = fasp_blas_mxv_bbsr;
221  break;
222 
223  case MAT_CSRL:
224  mf->fct = fasp_blas_mxv_csrl;
225  break;
226 
227  default:
228  printf("### ERROR: Wrong matrix format %d!\n", matrix_format);
229  exit(ERROR_DATA_STRUCTURE);
230 
231  }
232 
233  mf->data = A;
234 }
235 
236 /*---------------------------------*/
237 /*-- Private Functions --*/
238 /*---------------------------------*/
239 
252 static void fasp_blas_mxv_csr (void *A,
253  REAL *x,
254  REAL *y)
255 {
256  dCSRmat *a = (dCSRmat *)A;
257  fasp_blas_dcsr_mxv(a, x, y);
258 }
259 
272 static void fasp_blas_mxv_bsr (void *A,
273  REAL *x,
274  REAL *y)
275 {
276  dBSRmat *a = (dBSRmat *)A;
277  fasp_blas_dbsr_mxv(a, x, y);
278 }
279 
292 static void fasp_blas_mxv_str (void *A,
293  REAL *x,
294  REAL *y)
295 {
296  dSTRmat *a = (dSTRmat *)A;
297  fasp_blas_dstr_mxv(a, x, y);
298 }
299 
312 static void fasp_blas_mxv_bcsr (void *A,
313  REAL *x,
314  REAL *y)
315 {
316  block_dCSRmat *a = (block_dCSRmat *)A;
317  fasp_blas_bdcsr_mxv(a, x, y);
318 }
319 
332 static void fasp_blas_mxv_bbsr (void *A,
333  REAL *x,
334  REAL *y)
335 {
336  block_BSR *a = (block_BSR *)A;
337  fasp_blas_bdbsr_mxv(a, x, y);
338 }
339 
352 static void fasp_blas_mxv_csrl (void *A,
353  REAL *x,
354  REAL *y)
355 {
356  dCSRLmat *a = (dCSRLmat *)A;
357  fasp_blas_dcsrl_mxv(a, x, y);
358 }
359 
360 /*---------------------------------*/
361 /*-- End of File --*/
362 /*---------------------------------*/
INT fasp_solver_krylov(mxv_matfree *mf, dvector *b, dvector *x, itsolver_param *itparam)
Solve Ax=b by standard Krylov methods – without preconditioner.
Definition: itsolver_mf.c:150
INT fasp_solver_pcg(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient (CG) method for solving Au=b.
Definition: pcg_mf.c:86
#define REAL
Definition: fasp.h:67
void fasp_blas_bdcsr_mxv(block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_bcsr.c:155
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
Block REAL matrix format for reservoir simulation.
Definition: fasp_block.h:172
#define SOLVER_GCG
Definition: fasp_const.h:109
Structure matrix of REAL type.
Definition: fasp.h:304
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
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 MAT_CSR
Definition: fasp_const.h:90
#define PRINT_SOME
Definition: fasp_const.h:81
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
Preconditioner data and action.
Definition: fasp.h:999
#define MAT_BSR
Definition: fasp_const.h:91
#define ERROR_SOLVER_TYPE
Definition: fasp_const.h:47
INT restart
Definition: fasp.h:1112
#define MAT_bBSR
Definition: fasp_const.h:94
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
SHORT stop_type
Definition: fasp.h:1109
#define SOLVER_BiCGstab
Definition: fasp_const.h:104
INT fasp_solver_itsolver(mxv_matfree *mf, dvector *b, dvector *x, precond *pc, itsolver_param *itparam)
Solve Ax=b by preconditioned Krylov methods for CSR matrices.
Definition: itsolver_mf.c:50
#define SOLVER_VFGMRES
Definition: fasp_const.h:108
void fasp_blas_dstr_mxv(dSTRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_str.c:117
Main header file for FASP.
#define SOLVER_CG
Definition: fasp_const.h:103
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1013
void fasp_gettime(REAL *time)
Get system time.
Definition: timing.c:28
#define MAT_STR
Definition: fasp_const.h:92
#define SOLVER_GMRES
Definition: fasp_const.h:106
void fasp_blas_bdbsr_mxv(block_BSR *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_bcsr.c:326
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
INT fasp_solver_pbcgs(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned BiCGstab method for solving Au=b.
Definition: pbcgs_mf.c:91
#define MAT_bCSR
Definition: fasp_const.h:93
Parameters passed to iterative solvers.
Definition: fasp.h:1105
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT fasp_solver_pgcg(mxv_matfree *mf, 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_mf.c:47
INT row
number of rows
Definition: fasp.h:345
REAL tol
Definition: fasp.h:1111
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
void fasp_solver_itsolver_init(INT matrix_format, mxv_matfree *mf, void *A)
Initialize itsovlers.
Definition: itsolver_mf.c:197
SHORT print_level
Definition: fasp.h:1113
#define ERROR_DATA_STRUCTURE
Definition: fasp_const.h:38
Block REAL CSR matrix format.
Definition: fasp_block.h:84
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
void fasp_blas_dbsr_mxv(dBSRmat *A, REAL *x, REAL *y)
Compute y := A*x.
Definition: blas_bsr.c:895
void print_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: message.c:165
INT fasp_solver_pminres(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
A preconditioned minimal residual (Minres) method for solving Au=b.
Definition: pminres_mf.c:89
Header file for FASP block matrices.
Sparse matrix of REAL type in CSRL format.
Definition: fasp.h:265
void fasp_blas_dcsrl_mxv(dCSRLmat *A, REAL *x, REAL *y)
Compute y = A*x for a sparse matrix in CSRL format.
Definition: blas_csrl.c:28
#define SOLVER_VGMRES
Definition: fasp_const.h:107
void * data
data for MxV, can be a Matrix or something else
Definition: fasp.h:1016
#define SOLVER_MinRes
Definition: fasp_const.h:105
#define MAT_CSRL
Definition: fasp_const.h:95
#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
SHORT itsolver_type
Definition: fasp.h:1107