Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
SolSTR.c
Go to the documentation of this file.
1
16#include <math.h>
17#include <time.h>
18
19#include "fasp.h"
20#include "fasp_functs.h"
21
22/*---------------------------------*/
23/*-- Declare Private Functions --*/
24/*---------------------------------*/
25
26#include "KryUtil.inl"
27
28/*---------------------------------*/
29/*-- Public Functions --*/
30/*---------------------------------*/
31
52 dvector *b,
53 dvector *x,
54 precond *pc,
55 ITS_param *itparam)
56{
57 const SHORT prtlvl = itparam->print_level;
58 const SHORT itsolver_type = itparam->itsolver_type;
59 const SHORT stop_type = itparam->stop_type;
60 const SHORT restart = itparam->restart;
61 const INT MaxIt = itparam->maxit;
62 const REAL tol = itparam->tol;
63
64 // local variables
66 REAL solve_start, solve_end;
67
68#if DEBUG_MODE > 0
69 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
70 printf("### DEBUG: rhs/sol size: %d %d\n", b->row, x->row);
71#endif
72
73 fasp_gettime(&solve_start);
74
75 /* Safe-guard checks on parameters */
76 ITS_CHECK ( MaxIt, tol );
77
78 switch (itsolver_type) {
79
80 case SOLVER_CG:
81 iter=fasp_solver_dstr_pcg(A, b, x, pc, tol, MaxIt, stop_type, prtlvl);
82 break;
83
84 case SOLVER_BiCGstab:
85 iter=fasp_solver_dstr_pbcgs(A, b, x, pc, tol, MaxIt, stop_type, prtlvl);
86 break;
87
88 case SOLVER_GMRES:
89 iter=fasp_solver_dstr_pgmres(A, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
90 break;
91
92 case SOLVER_VGMRES:
93 iter=fasp_solver_dstr_pvgmres(A, b, x, pc, tol, MaxIt, restart, stop_type, prtlvl);
94 break;
95
96 default:
97 printf("### ERROR: Unknown iterative solver type %d! [%s]\n",
98 itsolver_type, __FUNCTION__);
99 return ERROR_SOLVER_TYPE;
100
101 }
102
103 if ( (prtlvl > PRINT_MIN) && (iter >= 0) ) {
104 fasp_gettime(&solve_end);
105 fasp_cputime("Iterative method", solve_end - solve_start);
106 }
107
108#if DEBUG_MODE > 0
109 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
110#endif
111
112 return iter;
113}
114
132 dvector *b,
133 dvector *x,
134 ITS_param *itparam)
135{
136 const SHORT prtlvl = itparam->print_level;
137 INT status = FASP_SUCCESS;
138 REAL solve_start, solve_end;
139
140#if DEBUG_MODE > 0
141 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
142#endif
143
144 // solver part
145 fasp_gettime(&solve_start);
146
147 status=fasp_solver_dstr_itsolver(A,b,x,NULL,itparam);
148
149 fasp_gettime(&solve_end);
150
151 if ( prtlvl >= PRINT_MIN )
152 fasp_cputime("Krylov method totally", solve_end - solve_start);
153
154#if DEBUG_MODE > 0
155 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
156#endif
157
158 return status;
159}
160
178 dvector *b,
179 dvector *x,
180 ITS_param *itparam)
181{
182 const SHORT prtlvl = itparam->print_level;
183 const INT ngrid = A->ngrid;
184
185 INT status = FASP_SUCCESS;
186 REAL solve_start, solve_end;
187 INT nc = A->nc, nc2 = nc*nc, i;
188
189 fasp_gettime(&solve_start);
190
191 // setup preconditioner
192 precond_diag_str diag;
193 fasp_dvec_alloc(ngrid*nc2, &(diag.diag));
194 fasp_darray_cp(ngrid*nc2, A->diag, diag.diag.val);
195
196 diag.nc = nc;
197
198 for (i=0;i<ngrid;++i) fasp_smat_inv(&(diag.diag.val[i*nc2]),nc);
199
200 precond *pc = (precond *)fasp_mem_calloc(1,sizeof(precond));
201
202 pc->data = &diag;
204
205#if DEBUG_MODE > 0
206 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
207#endif
208
209 // solver part
210 status=fasp_solver_dstr_itsolver(A,b,x,pc,itparam);
211
212 fasp_gettime(&solve_end);
213
214 if ( prtlvl >= PRINT_MIN )
215 fasp_cputime("Diag_Krylov method totally", solve_end - solve_start);
216
217#if DEBUG_MODE > 0
218 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
219#endif
220
221 return status;
222}
223
242 dvector *b,
243 dvector *x,
244 ITS_param *itparam,
245 ILU_param *iluparam)
246{
247 const SHORT prtlvl = itparam->print_level;
248 const INT ILU_lfil = iluparam->ILU_lfil;
249
250 INT status = FASP_SUCCESS;
251 REAL setup_start, setup_end, setup_time;
252 REAL solve_start, solve_end, solve_time;
253
254 //set up
255 dSTRmat LU;
256
257#if DEBUG_MODE > 0
258 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
259#endif
260
261 fasp_gettime(&setup_start);
262
263 if (ILU_lfil == 0) {
265 }
266 else if (ILU_lfil == 1) {
268 }
269 else {
270 printf("### ERROR: Illegal level of ILU fill-in (>1)! [%s]\n", __FUNCTION__);
271 return ERROR_MISC;
272 }
273
274 fasp_gettime(&setup_end);
275
276 setup_time = setup_end - setup_start;
277
278 if ( prtlvl > PRINT_NONE )
279 printf("Structrued ILU(%d) setup costs %f seconds.\n", ILU_lfil, setup_time);
280
281 precond pc; pc.data=&LU;
282 if (ILU_lfil == 0) {
284 }
285 else if (ILU_lfil == 1) {
287 }
288 else {
289 printf("### ERROR: Illegal level of ILU fill-in (>1)! [%s]\n", __FUNCTION__);
290 return ERROR_MISC;
291 }
292
293 // solver part
294 fasp_gettime(&solve_start);
295
296 status=fasp_solver_dstr_itsolver(A,b,x,&pc,itparam);
297
298 fasp_gettime(&solve_end);
299
300 if ( prtlvl >= PRINT_MIN ) {
301 solve_time = solve_end - solve_start;
302 printf("Iterative solver costs %f seconds.\n", solve_time);
303 fasp_cputime("ILU_Krylov method totally", setup_time+solve_time);
304 }
305
306 fasp_dstr_free(&LU);
307
308#if DEBUG_MODE > 0
309 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
310#endif
311
312 return status;
313}
314
335 dvector *b,
336 dvector *x,
337 ITS_param *itparam,
338 ivector *neigh,
339 ivector *order)
340{
341 // Parameter for iterative method
342 const SHORT prtlvl = itparam->print_level;
343
344 // Information about matrices
345 INT ngrid = A->ngrid;
346
347 // return parameter
348 INT status = FASP_SUCCESS;
349
350 // local parameter
351 REAL setup_start, setup_end, setup_time = 0;
352 REAL solve_start, solve_end, solve_time = 0;
353
354 dvector *diaginv;
355 ivector *pivot;
356
357#if DEBUG_MODE > 0
358 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
359#endif
360
361 // setup preconditioner
362 fasp_gettime(&setup_start);
363
364 diaginv = (dvector *)fasp_mem_calloc(ngrid, sizeof(dvector));
365 pivot = (ivector *)fasp_mem_calloc(ngrid, sizeof(ivector));
366 fasp_generate_diaginv_block(A, neigh, diaginv, pivot);
367
368 precond_data_str pcdata;
369 pcdata.A_str = A;
370 pcdata.diaginv = diaginv;
371 pcdata.pivot = pivot;
372 pcdata.order = order;
373 pcdata.neigh = neigh;
374
375 precond pc; pc.data = &pcdata; pc.fct = fasp_precond_dstr_blockgs;
376
377 fasp_gettime(&setup_end);
378
379 if ( prtlvl > PRINT_NONE ) {
380 setup_time = setup_end - setup_start;
381 printf("Preconditioner setup costs %f seconds.\n", setup_time);
382 }
383
384 // solver part
385 fasp_gettime(&solve_start);
386
387 status = fasp_solver_dstr_itsolver(A,b,x,&pc,itparam);
388
389 fasp_gettime(&solve_end);
390 solve_time = solve_end - solve_start;
391
392 if ( prtlvl >= PRINT_MIN ) {
393 fasp_cputime("Iterative solver", solve_time);
394 fasp_cputime("BlockGS_Krylov method totally", setup_time + solve_time);
395 }
396
397#if DEBUG_MODE > 0
398 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
399#endif
400
401 return status;
402}
403
404/*---------------------------------*/
405/*-- End of File --*/
406/*---------------------------------*/
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:164
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:36
void fasp_dvec_alloc(const INT m, dvector *u)
Create dvector data space of REAL type.
Definition: AuxVector.c:105
void fasp_ilu_dstr_setup0(dSTRmat *A, dSTRmat *LU)
Get ILU(0) decomposition of a structured matrix A.
void fasp_ilu_dstr_setup1(dSTRmat *A, dSTRmat *LU)
Get ILU(1) decoposition of a structured matrix A.
SHORT fasp_smat_inv(REAL *a, const INT n)
Compute the inverse matrix of a small full matrix a.
void fasp_dstr_free(dSTRmat *A)
Free STR sparse matrix data memeory space.
Definition: BlaSparseSTR.c:136
void fasp_generate_diaginv_block(dSTRmat *A, ivector *neigh, dvector *diaginv, ivector *pivot)
Generate inverse of diagonal block for block smoothers.
INT fasp_solver_dstr_pbcgs(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b for STR matrix.
Definition: KryPbcgs.c:1039
INT fasp_solver_dstr_pcg(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: KryPcg.c:978
INT fasp_solver_dstr_pgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:979
INT fasp_solver_dstr_pvgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during it...
Definition: KryPvgmres.c:1104
void fasp_precond_dstr_diag(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreSTR.c:44
void fasp_precond_dstr_blockgs(REAL *r, REAL *z, void *data)
CPR-type preconditioner (STR format)
Definition: PreSTR.c:1715
void fasp_precond_dstr_ilu1(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(1) decomposition.
Definition: PreSTR.c:349
void fasp_precond_dstr_ilu0(REAL *r, REAL *z, void *data)
Preconditioning using STR_ILU(0) decomposition.
Definition: PreSTR.c:71
INT fasp_solver_dstr_itsolver(dSTRmat *A, dvector *b, dvector *x, precond *pc, ITS_param *itparam)
Solve Ax=b by standard Krylov methods.
Definition: SolSTR.c:51
INT fasp_solver_dstr_krylov_diag(dSTRmat *A, dvector *b, dvector *x, ITS_param *itparam)
Solve Ax=b by diagonal preconditioned Krylov methods.
Definition: SolSTR.c:177
INT fasp_solver_dstr_krylov_blockgs(dSTRmat *A, dvector *b, dvector *x, ITS_param *itparam, ivector *neigh, ivector *order)
Solve Ax=b by diagonal preconditioned Krylov methods.
Definition: SolSTR.c:334
INT fasp_solver_dstr_krylov_ilu(dSTRmat *A, dvector *b, dvector *x, ITS_param *itparam, ILU_param *iluparam)
Solve Ax=b by structured ILU preconditioned Krylov methods.
Definition: SolSTR.c:241
INT fasp_solver_dstr_krylov(dSTRmat *A, dvector *b, dvector *x, ITS_param *itparam)
Solve Ax=b by standard Krylov methods.
Definition: SolSTR.c:131
Main header file for the FASP project.
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define INT
Definition: fasp.h:64
#define SOLVER_GMRES
Definition: fasp_const.h:106
#define ERROR_MISC
Definition: fasp_const.h:28
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
#define ERROR_SOLVER_TYPE
Definition: fasp_const.h:41
#define SOLVER_VGMRES
Definition: fasp_const.h:107
#define SOLVER_BiCGstab
Definition: fasp_const.h:104
#define SOLVER_CG
Definition: fasp_const.h:103
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
#define PRINT_MIN
Definition: fasp_const.h:74
Parameters for ILU.
Definition: fasp.h:396
INT ILU_lfil
level of fill-in for ILUk
Definition: fasp.h:405
Parameters for iterative solvers.
Definition: fasp.h:379
SHORT itsolver_type
Definition: fasp.h:382
SHORT print_level
Definition: fasp.h:381
REAL tol
Definition: fasp.h:388
INT restart
Definition: fasp.h:386
SHORT stop_type
Definition: fasp.h:385
INT maxit
Definition: fasp.h:387
Structure matrix of REAL type.
Definition: fasp.h:308
REAL * diag
diagonal entries (length is ngrid*(nc^2))
Definition: fasp.h:329
INT ngrid
number of grids
Definition: fasp.h:326
INT nc
size of each block (number of components)
Definition: fasp.h:323
Vector with n entries of REAL type.
Definition: fasp.h:346
REAL * val
actual vector entries
Definition: fasp.h:352
INT row
number of rows
Definition: fasp.h:349
Vector with n entries of INT type.
Definition: fasp.h:361
Data for preconditioners in dSTRmat format.
Definition: fasp.h:973
ivector * neigh
array to store neighbor information
Definition: fasp.h:1047
ivector * pivot
the pivot for the GS/block GS smoother (whole reservoir matrix)
Definition: fasp.h:1035
dSTRmat * A_str
store the whole reservoir block in STR format
Definition: fasp.h:1024
ivector * order
order for smoothing
Definition: fasp.h:1044
dvector * diaginv
the inverse of the diagonals for GS/block GS smoother (whole reservoir matrix)
Definition: fasp.h:1032
Data for diagonal preconditioners in dSTRmat format.
Definition: fasp.h:1065
INT nc
number of components
Definition: fasp.h:1068
dvector diag
diagonal elements
Definition: fasp.h:1071
Preconditioner data and action.
Definition: fasp.h:1081
void * data
data for preconditioner, void pointer
Definition: fasp.h:1084
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1087