Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
itsolver_str.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <time.h>
8 
9 #include "fasp.h"
10 #include "fasp_functs.h"
11 #include "itsolver_util.inl"
12 
13 /*---------------------------------*/
14 /*-- Public Functions --*/
15 /*---------------------------------*/
16 
35  dvector *b,
36  dvector *x,
37  precond *pc,
38  itsolver_param *itparam)
39 {
40  const SHORT prtlvl = itparam->print_level;
41  const SHORT itsolver_type = itparam->itsolver_type;
42  const SHORT stop_type = itparam->stop_type;
43  const SHORT restart = itparam->restart;
44  const INT MaxIt = itparam->maxit;
45  const REAL tol = itparam->tol;
46 
47  // local variables
48  INT iter = ERROR_SOLVER_TYPE;
49  REAL solver_start, solver_end, solver_duration;
50 
51 #if DEBUG_MODE > 0
52  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
53  printf("### DEBUG: rhs/sol size: %d %d\n", b->row, x->row);
54 #endif
55 
56  fasp_gettime(&solver_start);
57 
58  /* Safe-guard checks on parameters */
59  ITS_CHECK ( MaxIt, tol );
60 
61  switch (itsolver_type) {
62 
63  case SOLVER_CG:
64  if ( prtlvl > PRINT_NONE ) printf("\nCalling CG solver (STR) ...\n");
65  iter=fasp_solver_dstr_pcg(A, b, x, pc, tol, MaxIt, stop_type, prtlvl);
66  break;
67 
68  case SOLVER_BiCGstab:
69  if ( prtlvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (STR) ...\n");
70  iter=fasp_solver_dstr_pbcgs(A, b, x, pc, tol, MaxIt, stop_type, prtlvl);
71  break;
72 
73  case SOLVER_GMRES:
74  if ( prtlvl > PRINT_NONE ) printf("\nCalling GMRES solver (STR) ...\n");
75  iter=fasp_solver_dstr_pgmres(A, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
76  break;
77 
78  case SOLVER_VGMRES:
79  if ( prtlvl > PRINT_NONE ) printf("\nCalling vGMRES solver (STR) ...\n");
80  iter=fasp_solver_dstr_pvgmres(A, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
81  break;
82 
83  default:
84  printf("### ERROR: Unknown itertive solver type %d!\n", itsolver_type);
85 
86  }
87 
88  if ( (prtlvl > PRINT_MIN) && (iter >= 0) ) {
89  fasp_gettime(&solver_end);
90  solver_duration = solver_end - solver_start;
91  print_cputime("Iterative method", solver_duration);
92  }
93 
94 #if DEBUG_MODE > 0
95  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
96 #endif
97 
98  return iter;
99 }
100 
118  dvector *b,
119  dvector *x,
120  itsolver_param *itparam)
121 {
122  const SHORT prtlvl = itparam->print_level;
123  INT status = FASP_SUCCESS;
124  REAL solver_start, solver_end, solver_duration;
125 
126 #if DEBUG_MODE > 0
127  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
128 #endif
129 
130  // solver part
131  fasp_gettime(&solver_start);
132 
133  status=fasp_solver_dstr_itsolver(A,b,x,NULL,itparam);
134 
135  fasp_gettime(&solver_end);
136 
137  if ( prtlvl >= PRINT_MIN ) {
138  solver_duration = solver_end - solver_start;
139  print_cputime("Krylov method totally", solver_duration);
140  }
141 
142 #if DEBUG_MODE > 0
143  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
144 #endif
145 
146  return status;
147 }
148 
166  dvector *b,
167  dvector *x,
168  itsolver_param *itparam)
169 {
170  const SHORT prtlvl = itparam->print_level;
171  const INT ngrid = A->ngrid;
172 
173  INT status = FASP_SUCCESS;
174  REAL solver_start, solver_end, solver_duration;
175  INT nc = A->nc, nc2 = nc*nc, i;
176 
177  // setup preconditioner
178  precond_diagstr diag;
179  fasp_dvec_alloc(ngrid*nc2, &(diag.diag));
180  fasp_array_cp(ngrid*nc2, A->diag, diag.diag.val);
181 
182  diag.nc = nc;
183 
184  for (i=0;i<ngrid;++i) fasp_blas_smat_inv(&(diag.diag.val[i*nc2]),nc);
185 
186  precond *pc = (precond *)fasp_mem_calloc(1,sizeof(precond));
187 
188  pc->data = &diag;
190 
191 #if DEBUG_MODE > 0
192  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
193 #endif
194 
195  // solver part
196  fasp_gettime(&solver_start);
197 
198  status=fasp_solver_dstr_itsolver(A,b,x,pc,itparam);
199 
200  fasp_gettime(&solver_end);
201 
202  solver_duration = solver_end - solver_start;
203 
204  if ( prtlvl >= PRINT_MIN )
205  print_cputime("Diag_Krylov method totally", solver_duration);
206 
207 #if DEBUG_MODE > 0
208  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
209 #endif
210 
211  return status;
212 }
213 
232  dvector *b,
233  dvector *x,
234  itsolver_param *itparam,
235  ILU_param *iluparam)
236 {
237  const SHORT prtlvl = itparam->print_level;
238  const INT ILU_lfil = iluparam->ILU_lfil;
239 
240  INT status = FASP_SUCCESS;
241  REAL setup_start, setup_end, setup_duration;
242  REAL solver_start, solver_end, solver_duration;
243 
244  //set up
245  dSTRmat LU;
246 
247 #if DEBUG_MODE > 0
248  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
249 #endif
250 
251  fasp_gettime(&setup_start);
252 
253  if (ILU_lfil == 0) {
254  fasp_ilu_dstr_setup0(A,&LU);
255  }
256  else if (ILU_lfil == 1) {
257  fasp_ilu_dstr_setup1(A,&LU);
258  }
259  else {
260  printf("### ERROR: Illegal level of fill-in for ILU (lfil>=2)!\n");
261  return ERROR_MISC;
262  }
263 
264  fasp_gettime(&setup_end);
265 
266  setup_duration = setup_end - setup_start;
267 
268  if ( prtlvl > PRINT_NONE )
269  printf("structrued ILU(%d) setup costs %f seconds.\n", ILU_lfil, setup_duration);
270 
271  precond pc; pc.data=&LU;
272  if (ILU_lfil == 0) {
274  }
275  else if (ILU_lfil == 1) {
277  }
278  else {
279  printf("### ERROR: Illegal level of fill-in for ILU (lfil>=2)!\n");
280  return ERROR_MISC;
281  }
282 
283  // solver part
284  fasp_gettime(&solver_start);
285 
286  status=fasp_solver_dstr_itsolver(A,b,x,&pc,itparam);
287 
288  fasp_gettime(&solver_end);
289 
290  if ( prtlvl >= PRINT_MIN ) {
291  solver_duration = solver_end - solver_start;
292  printf("Iterative solver costs %f seconds.\n", solver_duration);
293  print_cputime("ILU_Krylov method totally", setup_duration+solver_duration);
294  }
295 
296  fasp_dstr_free(&LU);
297 
298 #if DEBUG_MODE > 0
299  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
300 #endif
301 
302  return status;
303 }
304 
325  dvector *b,
326  dvector *x,
327  itsolver_param *itparam,
328  ivector *neigh,
329  ivector *order)
330 {
331  // Parameter for iterative method
332  const SHORT prtlvl = itparam->print_level;
333 
334  // Information about matrices
335  INT ngrid = A->ngrid;
336 
337  // return parameter
338  INT status = FASP_SUCCESS;
339 
340  // local parameter
341  REAL setup_start, setup_end, setup_duration;
342  REAL solver_start, solver_end, solver_duration;
343 
344  dvector *diaginv;
345  ivector *pivot;
346 
347 #if DEBUG_MODE > 0
348  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
349 #endif
350 
351  // setup preconditioner
352  fasp_gettime(&setup_start);
353 
354  diaginv = (dvector *)fasp_mem_calloc(ngrid, sizeof(dvector));
355  pivot = (ivector *)fasp_mem_calloc(ngrid, sizeof(ivector));
356  fasp_generate_diaginv_block(A, neigh, diaginv, pivot);
357 
358  precond_data_str pcdata;
359  pcdata.A_str = A;
360  pcdata.diaginv = diaginv;
361  pcdata.pivot = pivot;
362  pcdata.order = order;
363  pcdata.neigh = neigh;
364 
365  precond pc; pc.data = &pcdata; pc.fct = fasp_precond_dstr_blockgs;
366 
367  fasp_gettime(&setup_end);
368 
369  if ( prtlvl > PRINT_NONE ) {
370  setup_duration = setup_end - setup_start;
371  printf("Preconditioner setup costs %f seconds.\n", setup_duration);
372  }
373 
374  // solver part
375  fasp_gettime(&solver_start);
376 
377  status=fasp_solver_dstr_itsolver(A,b,x,&pc,itparam);
378 
379  fasp_gettime(&solver_end);
380 
381  if ( prtlvl >= PRINT_MIN ) {
382  solver_duration = solver_end - solver_start;
383  printf("Iterative solver costs %f seconds.\n", solver_duration);
384  print_cputime("BlockGS_Krylov method totally", setup_duration+solver_duration);
385  }
386 
387 #if DEBUG_MODE > 0
388  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
389 #endif
390 
391  return status;
392 }
393 
394 /*---------------------------------*/
395 /*-- End of File --*/
396 /*---------------------------------*/
INT fasp_solver_dstr_pgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
Preconditioned GMRES method for solving Au=b.
Definition: pgmres.c:963
Data passed to the preconditioner for dSTRmat matrices.
Definition: fasp.h:891
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
void fasp_generate_diaginv_block(dSTRmat *A, ivector *neigh, dvector *diaginv, ivector *pivot)
Generate inverse of diagonal block for block smoothers.
INT ngrid
number of grids
Definition: fasp.h:322
INT fasp_blas_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
Definition: smat.c:554
INT nc
size of each block (number of components)
Definition: fasp.h:319
INT nc
number of components
Definition: fasp.h:986
Structure matrix of REAL type.
Definition: fasp.h:304
ivector * neigh
array to store neighbor information
Definition: fasp.h:965
void fasp_precond_dstr_diag(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: precond_str.c:27
REAL * val
actual vector entries
Definition: fasp.h:348
INT ILU_lfil
level of fill-in for ILUk
Definition: fasp.h:383
void fasp_ilu_dstr_setup1(dSTRmat *A, dSTRmat *LU)
Get ILU(1) decoposition of a structured matrix A.
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 fasp_precond_dstr_ilu1(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(1) decomposition.
Definition: precond_str.c:336
dvector * diaginv
the inverse of the diagonals for GS/block GS smoother (whole reservoir matrix)
Definition: fasp.h:950
INT fasp_solver_dstr_pbcgs(dSTRmat *A, 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.c:1117
#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 ERROR_SOLVER_TYPE
Definition: fasp_const.h:47
INT restart
Definition: fasp.h:1112
ivector * pivot
the pivot for the GS/block GS smoother (whole reservoir matrix)
Definition: fasp.h:953
SHORT stop_type
Definition: fasp.h:1109
void * data
data for preconditioner, void pointer
Definition: fasp.h:1002
dvector diag
diagonal elements
Definition: fasp.h:989
ivector * order
order for smoothing
Definition: fasp.h:962
#define SOLVER_BiCGstab
Definition: fasp_const.h:104
#define ERROR_MISC
Definition: fasp_const.h:35
Main header file for FASP.
#define SOLVER_CG
Definition: fasp_const.h:103
void fasp_ilu_dstr_setup0(dSTRmat *A, dSTRmat *LU)
Get ILU(0) decomposition of a structured matrix A.
Definition: ilu_setup_str.c:28
void fasp_precond_dstr_blockgs(REAL *r, REAL *z, void *data)
CPR-type preconditioner (STR format)
Definition: precond_str.c:1706
void fasp_gettime(REAL *time)
Get system time.
Definition: timing.c:28
#define SOLVER_GMRES
Definition: fasp_const.h:106
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
Parameters passed to iterative solvers.
Definition: fasp.h:1105
INT fasp_solver_dstr_pcg(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: pcg.c:957
INT fasp_solver_dstr_krylov_ilu(dSTRmat *A, dvector *b, dvector *x, itsolver_param *itparam, ILU_param *iluparam)
Solve Ax=b by structured ILU preconditioned Krylov methods.
Definition: itsolver_str.c:231
INT row
number of rows
Definition: fasp.h:345
REAL tol
Definition: fasp.h:1111
REAL * diag
diagonal entries (length is ngrid*(nc^2))
Definition: fasp.h:325
INT fasp_solver_dstr_krylov_diag(dSTRmat *A, dvector *b, dvector *x, itsolver_param *itparam)
Solve Ax=b by diagonal preconditioned Krylov methods.
Definition: itsolver_str.c:165
SHORT print_level
Definition: fasp.h:1113
Parameters for ILU.
Definition: fasp.h:374
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT fasp_solver_dstr_krylov_blockgs(dSTRmat *A, dvector *b, dvector *x, itsolver_param *itparam, ivector *neigh, ivector *order)
Solve Ax=b by diagonal preconditioned Krylov methods.
Definition: itsolver_str.c:324
INT fasp_solver_dstr_krylov(dSTRmat *A, dvector *b, dvector *x, itsolver_param *itparam)
Solve Ax=b by standard Krylov methods.
Definition: itsolver_str.c:117
void print_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: message.c:165
dSTRmat * A_str
store the whole reservoir block in STR format
Definition: fasp.h:942
void fasp_dvec_alloc(const INT m, dvector *u)
Create dvector data space of REAL type.
Definition: vec.c:99
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
void fasp_precond_dstr_ilu0(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(0) decomposition.
Definition: precond_str.c:54
void fasp_dstr_free(dSTRmat *A)
Free STR sparse matrix data memeory space.
Definition: sparse_str.c:152
#define SOLVER_VGMRES
Definition: fasp_const.h:107
INT fasp_solver_dstr_pvgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during th...
Definition: pvgmres.c:1083
INT fasp_solver_dstr_itsolver(dSTRmat *A, dvector *b, dvector *x, precond *pc, itsolver_param *itparam)
Solve Ax=b by standard Krylov methods.
Definition: itsolver_str.c:34
Data passed to diagonal preconditioner for dSTRmat matrices.
Definition: fasp.h:983
#define PRINT_MIN
Definition: fasp_const.h:80
SHORT itsolver_type
Definition: fasp.h:1107