25 cout <<
"The input file " << myfile <<
" is missing!" << endl;
26 myfile =
solveDir +
"../conf/csr.fasp";
29 cout <<
"The input file " << myfile <<
" is missing!" << endl;
30 cout <<
"Using the default parameters of FASP" << endl;
33 cout <<
"Using the input file " << myfile << endl;
34 fasp_param_input(myfile.data(), &
inParam);
38 fasp_param_input(myfile.data(), &
inParam);
43 void ScalarFaspSolver::Allocate(
const vector<USI>& rowCapacity,
48 for (
OCP_USI n = 0; n < maxDim; n++) {
49 nnz += rowCapacity[n];
51 A = fasp_dcsr_create(maxDim, maxDim, nnz);
54 void ScalarFaspSolver::InitParam()
57 inParam.print_level = PRINT_MIN;
61 inParam.solver_type = SOLVER_VFGMRES;
62 inParam.precond_type = PREC_AMG;
63 inParam.stop_type = STOP_REL_RES;
81 inParam.SWZ_blksolver = SOLVER_DEFAULT;
86 inParam.AMG_cycle_type = V_CYCLE;
87 inParam.AMG_smoother = SMOOTHER_GS;
88 inParam.AMG_smooth_order = CF_ORDER;
90 inParam.AMG_postsmooth_iter = 1;
98 inParam.AMG_coarse_scaling = OFF;
100 inParam.AMG_nl_amli_krylov_type = 2;
103 inParam.AMG_coarsening_type = 1;
104 inParam.AMG_interpolation_type = 1;
106 inParam.AMG_strong_threshold = 0.3;
107 inParam.AMG_truncation_threshold = 0.2;
108 inParam.AMG_aggressive_level = 0;
109 inParam.AMG_aggressive_path = 1;
112 inParam.AMG_aggregation_type = PAIRWISE;
113 inParam.AMG_quality_bound = 8.0;
115 inParam.AMG_strong_coupled = 0.25;
116 inParam.AMG_max_aggregation = 9;
117 inParam.AMG_tentative_smooth = 0.67;
118 inParam.AMG_smooth_filter = ON;
119 inParam.AMG_smooth_restriction = ON;
122 void ScalarFaspSolver::AssembleMat(
const vector<vector<USI>>& colId,
123 const vector<vector<OCP_DBL>>& val,
126 vector<OCP_DBL>& rhs,
136 for (
OCP_USI i = 0; i < dim; i++) {
137 nnz += colId[i].size();
147 for (
OCP_USI i = 1; i < dim + 1; i++) {
148 USI nnz_Row = colId[i - 1].size();
149 A.IA[i] = A.IA[i - 1] + nnz_Row;
151 for (
USI j = 0; j < nnz_Row; j++) {
152 A.JA[count] = colId[i - 1][j];
153 A.val[count] = val[i - 1][j];
160 for (
int i = 0; i < dim; i++) {
161 if (!isfinite(b.val[i]))
OCP_ABORT(
"sFasp b is infinite!");
162 if (!isfinite(x.val[i]))
OCP_ABORT(
"sFasp x is infinite!");
165 for (
int i = 0; i < nnz; i++) {
166 if (!isfinite(A.val[i]))
OCP_ABORT(
"sFasp A is infinite!");
171 OCP_INT ScalarFaspSolver::Solve()
181 const char* outputfile =
"out/Solver.out";
182 printf(
"Redirecting outputs to file: %s ...\n", outputfile);
183 freopen(outputfile,
"w", stdout);
187 if (solver_type >= 1 && solver_type <= 20) {
190 if (precond_type == PREC_NULL) {
191 status = fasp_solver_dcsr_krylov(&A, &b, &x, &
itParam);
195 else if (precond_type == PREC_DIAG) {
196 status = fasp_solver_dcsr_krylov_diag(&A, &b, &x, &
itParam);
200 else if (precond_type == PREC_AMG || precond_type == PREC_FMG) {
201 if (print_level > PRINT_NONE) fasp_param_amg_print(&
amgParam);
206 else if (precond_type == PREC_ILU) {
207 if (print_level > PRINT_NONE) fasp_param_ilu_print(&
iluParam);
213 printf(
"### ERROR: Wrong preconditioner type %d!!!\n", precond_type);
214 status = ERROR_SOLVER_PRECTYPE;
220 else if (solver_type == SOLVER_AMG) {
221 if (print_level > PRINT_NONE) fasp_param_amg_print(&
amgParam);
222 fasp_solver_amg(&A, &b, &x, &
amgParam);
226 else if (solver_type == SOLVER_FMG) {
227 if (print_level > PRINT_NONE) fasp_param_amg_print(&
amgParam);
228 fasp_solver_famg(&A, &b, &x, &
amgParam);
232 else if (solver_type == SOLVER_MUMPS) {
233 status = fasp_solver_mumps(&A, &b, &x, print_level);
234 if (status >= 0) status = 1;
239 else if (solver_type == SOLVER_SUPERLU) {
240 status = fasp_solver_superlu(&A, &b, &x, print_level);
241 if (status >= 0) status = 1;
246 else if (solver_type == SOLVER_UMFPACK) {
248 dCSRmat A_trans = fasp_dcsr_create(A.row, A.col, A.nnz);
249 fasp_dcsr_transz(&A, NULL, &A_trans);
250 fasp_dcsr_sort(&A_trans);
251 status = fasp_solver_umfpack(&A_trans, &b, &x, print_level);
252 fasp_dcsr_free(&A_trans);
253 if (status >= 0) status = 1;
258 else if (solver_type == SOLVER_PARDISO) {
260 status = fasp_solver_pardiso(&A, &b, &x, print_level);
261 if (status >= 0) status = 1;
266 printf(
"### ERROR: Wrong solver type %d!!!\n", solver_type);
267 status = ERROR_SOLVER_TYPE;
271 printf(
"### ERROR: Solver failed! Exit status = %d.\n\n", status);
274 if (output_type) fclose(stdout);
278 void VectorFaspSolver::Allocate(
const vector<USI>& rowCapacity,
283 for (
OCP_USI n = 0; n < maxDim; n++) {
284 nnz += rowCapacity[n];
286 A = fasp_dbsr_create(maxDim, maxDim, nnz, blockDim, 0);
287 Asc = fasp_dbsr_create(maxDim, maxDim, nnz, blockDim, 0);
288 fsc = fasp_dvec_create(maxDim * blockDim);
289 order = fasp_ivec_create(maxDim);
290 Dmat.resize(maxDim * blockDim * blockDim);
293 void VectorFaspSolver::InitParam()
296 inParam.print_level = PRINT_MIN;
300 inParam.solver_type = SOLVER_VFGMRES;
303 inParam.stop_type = STOP_REL_RES;
321 inParam.SWZ_blksolver = SOLVER_DEFAULT;
324 inParam.AMG_type = CLASSIC_AMG;
326 inParam.AMG_cycle_type = V_CYCLE;
327 inParam.AMG_smoother = SMOOTHER_GS;
328 inParam.AMG_smooth_order = CF_ORDER;
329 inParam.AMG_presmooth_iter = 1;
330 inParam.AMG_postsmooth_iter = 1;
338 inParam.AMG_coarse_scaling = OFF;
340 inParam.AMG_nl_amli_krylov_type = 2;
343 inParam.AMG_coarsening_type = 1;
344 inParam.AMG_interpolation_type = 1;
346 inParam.AMG_strong_threshold = 0.3;
347 inParam.AMG_truncation_threshold = 0.2;
348 inParam.AMG_aggressive_level = 0;
349 inParam.AMG_aggressive_path = 1;
352 inParam.AMG_aggregation_type = PAIRWISE;
353 inParam.AMG_quality_bound = 8.0;
355 inParam.AMG_strong_coupled = 0.25;
356 inParam.AMG_max_aggregation = 9;
357 inParam.AMG_tentative_smooth = 0.67;
358 inParam.AMG_smooth_filter = ON;
359 inParam.AMG_smooth_restriction = ON;
362 void VectorFaspSolver::AssembleMat(
const vector<vector<USI>>& colId,
363 const vector<vector<OCP_DBL>>& val,
366 vector<OCP_DBL>& rhs,
369 const OCP_USI nrow = dim * blockDim;
382 for (
OCP_USI i = 0; i < dim; i++) {
383 nnz += colId[i].size();
402 for (
OCP_USI i = 1; i < dim + 1; i++) {
403 USI nnb_Row = colId[i - 1].size();
404 A.IA[i] = A.IA[i - 1] + nnb_Row;
406 for (
USI j = 0; j < nnb_Row; j++) {
407 A.JA[count1] = colId[i - 1][j];
410 size_row = nnb_Row * blockDim * blockDim;
411 for (
USI k = 0; k < size_row; k++) {
412 A.val[count2] = val[i - 1][k];
419 for (
int i = 0; i < nrow; i++) {
420 if (!isfinite(b.val[i]))
OCP_ABORT(
"vFasp b is infinite!");
421 if (!isfinite(x.val[i]))
OCP_ABORT(
"vFasp x is infinite!");
424 const OCP_USI len = A.NNZ * blockDim * blockDim;
425 for (
int i = 0; i < len; i++) {
426 if (!isfinite(A.val[i]))
OCP_ABORT(
"vFasp A is infinite!");
431 OCP_INT VectorFaspSolver::Solve()
441 #if WITH_FASP4BLKOIL || WITH_FASPCPR
446 const char* outputfile =
"../output/test.out";
447 printf(
"Redirecting outputs to file: %s ...\n", outputfile);
448 freopen(outputfile,
"w", stdout);
451 fasp_dvec_set(x.row, &x, 0);
454 if (solver_type >= 1 && solver_type <= 10) {
457 switch (precond_type) {
459 status = fasp_solver_dbsr_krylov(&A, &b, &x, &
itParam);
462 status = fasp_solver_dbsr_krylov_diag(&A, &b, &x, &
itParam);
471 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
477 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
485 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
487 status = fasp_solver_dbsr_krylov_FASP1_cuda_interface(
490 status = fasp_solver_dbsr_krylov_FASP1a(
495 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
497 status = fasp_solver_dbsr_krylov_FASP1_cuda_share_interface(
501 status = fasp_solver_dbsr_krylov_FASP1a_share_interface(
507 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
508 status = fasp_solver_dbsr_krylov_FASP2(
512 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
513 status = fasp_solver_dbsr_krylov_FASP3(
517 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
519 status = fasp_solver_dbsr_krylov_FASP4_cuda(
522 status = fasp_solver_dbsr_krylov_FASP4(
527 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
529 status = fasp_solver_dbsr_krylov_FASP4_cuda_share_interface(
533 status = fasp_solver_dbsr_krylov_FASP4_share_interface(
539 Decoupling(&A, &b, &Asc, &fsc, &order, Dmat.data(), decoup_type);
540 status = fasp_solver_dbsr_krylov_FASP5(
545 OCP_WARNING(
"fasp4blkoil was not linked correctly!");
546 OCP_ABORT(
"Preconditioner type " + to_string(precond_type) +
549 fill(Dmat.begin(), Dmat.end(), 0.0);
553 else if (solver_type == SOLVER_MUMPS) {
554 dCSRmat Acsr = fasp_format_dbsr_dcsr(&A);
555 status = fasp_solver_mumps(&Acsr, &b, &x, print_level);
556 fasp_dcsr_free(&Acsr);
557 if (status >= 0) status = 1;
562 else if (solver_type == SOLVER_SUPERLU) {
563 dCSRmat Acsr = fasp_format_dbsr_dcsr(&A);
564 status = fasp_solver_superlu(&Acsr, &b, &x, print_level);
565 fasp_dcsr_free(&Acsr);
566 if (status >= 0) status = 1;
571 else if (solver_type == SOLVER_UMFPACK) {
573 dCSRmat Acsr = fasp_format_dbsr_dcsr(&A);
574 dCSRmat A_trans = fasp_dcsr_create(Acsr.row, Acsr.col, Acsr.nnz);
575 fasp_dcsr_transz(&Acsr, NULL, &A_trans);
576 fasp_dcsr_sort(&A_trans);
577 status = fasp_solver_umfpack(&A_trans, &b, &x, print_level);
578 fasp_dcsr_free(&A_trans);
579 fasp_dcsr_free(&Acsr);
580 if (status >= 0) status = 1;
585 else if (solver_type == SOLVER_PARDISO) {
586 dCSRmat Acsr = fasp_format_dbsr_dcsr(&A);
587 fasp_dcsr_sort(&Acsr);
588 status = fasp_solver_pardiso(&Acsr, &b, &x, print_level);
589 fasp_dcsr_free(&Acsr);
590 if (status >= 0) status = 1;
595 printf(
"### ERROR: Wrong solver type %d!!!\n", solver_type);
596 status = ERROR_SOLVER_TYPE;
599 if (print_level > PRINT_MIN) {
601 cout <<
"\n### WARNING: Solver does not converge!\n" << endl;
603 cout <<
"\nSolver converges successfully!\n" << endl;
607 if (output_type) fclose(stdout);
Declaration of classes interfacing to the FASP solvers.
#define PC_BILU
BILU: block ILU preconditioner.
#define PC_FASP2
FASP2: MSP, experimental only.
#define PC_FASP3
FASP3: MSP, monolithic preconditioner.
#define PC_FASP5
FASP5: MSP, experimental only.
#define RESET_CONST
Sharing threshold for PC_FASP1_SHARE, PC_FASP4_SHARE.
#define PC_NULL
None: no preconditioner.
#define PC_FASP4
FASP4: MSP, default for FIM from 2015.
#define PC_DIAG
DIAG: diagonal preconditioner.
#define PC_FASP1_SHARE
Sharing setup stage for PC_FASP1, use with caution.
#define PC_FASP4_SHARE
Sharing setup stage for PC_FASP4, use with caution.
#define PC_FASP1
FASP1: MSP, default for FIM from 2020.
unsigned int USI
Generic unsigned integer.
unsigned int OCP_USI
Long unsigned integer.
#define OCP_WARNING(msg)
Log warning messages.
#define OCP_ABORT(msg)
Abort if critical error happens.
AMG_param amgParam
Parameters for AMG method.
input_param inParam
Parameters from input files.
ILU_param iluParam
Parameters for ILU method.
string solveFile
Relative path of fasp file.
SWZ_param swzParam
Parameters for Schwarz method.
string solveDir
Current work dir.
ITS_param itParam
Parameters for iterative method.
void SetupParam(const string &dir, const string &file) override
Set FASP parameters.
virtual void InitParam()=0
Initialize the params for linear solvers.