8 #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ 9 #define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ 15 #include "Vector/Vector.hpp" 17 #include <petsctime.h> 18 #include "Plot/GoogleChart.hpp" 19 #include "Matrix/SparseMatrix.hpp" 20 #include "Vector/Vector.hpp" 27 std::ostringstream out;
28 out << std::setprecision(n) << a_value;
32 #define UMFPACK_NONE 0 37 #define PCHYPRE_BOOMERAMG "petsc_boomeramg" 71 std::cerr <<
"Error Petsc only suppor double precision" <<
"/n";
75 #define SOLVER_NOOPTION 0 76 #define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1 77 #define SOLVER_PRINT_DETERMINANT 2 117 struct solv_bench_info
136 bool is_preconditioner_set =
false;
155 AMG_type atype = NONE_AMG;
170 double s_int = slv.time / slv.res.size();
173 size_t pt = std::floor(t / s_int);
175 if (pt < slv.res.size())
176 return slv.res.get(pt).err_norm;
178 return slv.res.last().err_norm;
187 if (create_vcluster().getProcessUnitID() != 0)
195 for (
size_t i = 0 ; i < bench.
size() ; i++)
196 x.add(bench.get(i).smethod);
200 yn.add(
"Norm infinity");
201 yn.add(
"Norm average");
202 yn.add(
"Number of iterations");
206 for (
size_t i = 0 ; i < bench.
size() ; i++)
207 y.add({bench.get(i).err.err_inf,bench.get(i).err.err_norm,(double)bench.get(i).err.its});
212 options.
title = std::string(
"Residual error");
213 options.
yAxis = std::string(
"Error");
214 options.
xAxis = std::string(
"Method");
215 options.
stype = std::string(
"bars");
216 options.
more = std::string(
"vAxis: {scaleType: 'log'}");
220 std::stringstream ss;
222 cg.
addHTML(
"<h2>Robustness</h2>");
229 yn.add(
"Time in millseconds");
230 options.
title = std::string(
"Time");
232 for (
size_t i = 0 ; i < bench.
size() ; i++)
233 y.add({bench.get(i).time});
246 double max_time = 0.0;
248 for (
size_t i = 0 ; i < bench.
size() ; i++)
250 if (bench.get(i).time > max_time) {max_time = bench.get(i).time;}
251 yn.add(bench.get(i).smethod);
254 size_t n_int = maxits;
258 double dt = max_time / n_int;
264 for (
double t = dt ; t <= max_time + 0.05 * max_time ; t += dt)
268 for (
size_t j = 0 ; j < bench.
size() ; j++)
269 y.last().add(calculate_it(t,bench.get(j)));
272 std::stringstream ss2;
273 ss2 <<
"<h2>Convergence with time</h2><br>" << std::endl;
275 options.
title = std::string(
"Residual norm infinity");
276 options.
yAxis = std::string(
"residual");
277 options.
xAxis = std::string(
"Time in milliseconds");
279 options.
more = std::string(
"vAxis: {scaleType: 'log'},hAxis: {scaleType: 'log'}");
284 ss <<
"<h2>Solvers Tested</h2><br><br>" << std::endl;
285 for (
size_t i = 0 ; i < bench.
size() ; i++)
286 ss << bench.get(i).smethod <<
" = " << bench.get(i).method <<
"<br>" << std::endl;
290 cg.
write(
"gc_solver.html");
302 PETSC_SAFE_CALL(KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED));
304 if (atype == HYPRE_AMG)
309 PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
311 PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
312 PCFactorSetShiftAmount(pc, PETSC_DECIDE);
322 auto & v_cl = create_vcluster();
324 if (v_cl.getProcessUnitID() == 0)
325 std::cout <<
"-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
351 PETSC_SAFE_CALL(KSPGetOperators(ksp,&A,NULL));
352 PETSC_SAFE_CALL(MatCreateVecs(A,&w,&v));
353 PETSC_SAFE_CALL(KSPBuildResidual(ksp,v,w,&Br));
354 PETSC_SAFE_CALL(KSPBuildResidual(ksp,v,w,&Br));
355 PETSC_SAFE_CALL(VecNorm(Br,NORM_INFINITY,&erri.err_norm));
358 size_t old_size = pts->
bench.last().res.
size();
359 pts->
bench.last().res.resize(it+1);
362 err_fill = pts->
bench.last().res.get(old_size-1);
366 for (
long int i = old_size ; i < (
long int)it ; i++)
367 pts->
bench.last().res.get(i) = err_fill;
370 pts->
bench.last().res.get(it) = erri;
384 PETSC_SAFE_CALL(KSPGetTolerances(ksp,NULL,NULL,NULL,&n_max_it));
386 auto & v_cl = create_vcluster();
388 if (std::floor(it * 100.0 / n_max_it) > tmp)
390 size_t diff = it * 100.0 / n_max_it - tmp;
391 tmp = it * 100.0 / n_max_it;
392 if (v_cl.getProcessUnitID() == 0)
394 for (
size_t k = 0 ; k < diff ; k++) {std::cout <<
"*";}
395 std::cout << std::flush;
407 auto & v_cl = create_vcluster();
409 if (v_cl.getProcessUnitID() == 0)
412 KSPGetType(ksp,&typ);
414 std::cout <<
"Method: " << typ <<
" " <<
" pre-conditoner: " << PCJACOBI <<
" iterations: " << err.
its << std::endl;
415 std::cout <<
"Norm of error: " << err.
err_norm <<
" Norm infinity: " << err.
err_inf << std::endl;
429 if (solv == std::string(KSPBCGS))
431 return std::string(
"BiCGStab (Stabilized version of BiConjugate Gradient Squared)");
433 else if (solv == std::string(KSPIBCGS))
435 return std::string(
"IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared)");
437 else if (solv == std::string(KSPFBCGS))
439 return std::string(
"FBiCGStab (Flexible Stabilized version of BiConjugate Gradient Squared)");
441 else if (solv == std::string(KSPFBCGSR))
443 return std::string(
"FBiCGStab-R (Modified Flexible Stabilized version of BiConjugate Gradient Squared)");
445 else if (solv == std::string(KSPBCGSL))
447 return std::string(
"BiCGStab(L) (Stabilized version of BiConjugate Gradient Squared with L search directions)");
449 else if (solv == std::string(KSPGMRES))
451 return std::string(
"GMRES(M) (Generalized Minimal Residual method with restart)");
453 else if (solv == std::string(KSPFGMRES))
455 return std::string(
"FGMRES(M) (Flexible Generalized Minimal Residual with restart)");
457 else if (solv == std::string(KSPLGMRES))
459 return std::string(
"LGMRES(M) (Augment Generalized Minimal Residual method with restart)");
461 else if (solv == std::string(KSPPGMRES))
463 return std::string(
"PGMRES(M) (Pipelined Generalized Minimal Residual method)");
465 else if (solv == std::string(KSPGCR))
467 return std::string(
"GCR (Generalized Conjugate Residual)");
470 return std::string(
"UNKNOWN");
482 bench.last().method = to_string_method(str);
483 bench.last().smethod = std::string(str);
498 VecCopy(sol,best_sol);
516 PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
519 double best_res = std::numeric_limits<double>::max();
522 auto & v_cl = create_vcluster();
527 for (
size_t i = 0 ; i < solvs.
size() ; i++)
535 PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
536 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_type",PCNONE));
538 setSolver(solvs.get(i).c_str());
541 if (solvs.get(i) == std::string(KSPBCGSL))
544 for (
size_t j = 2 ; j < 6 ; j++)
546 new_bench(solvs.get(i));
548 if (v_cl.getProcessUnitID() == 0)
549 std::cout <<
"L = " << j << std::endl;
550 bench.last().smethod += std::string(
"(") + std::to_string(j) + std::string(
")");
552 bench_solve_simple(A_,b_,x_,bench.last());
554 copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
557 else if (solvs.get(i) == std::string(KSPGMRES) ||
558 solvs.get(i) == std::string(std::string(KSPFGMRES)) ||
559 solvs.get(i) == std::string(std::string(KSPLGMRES)) )
562 for (
size_t j = 50 ; j < 300 ; j += 50)
565 new_bench(solvs.get(i));
567 if (v_cl.getProcessUnitID() == 0)
568 std::cout <<
"Restart = " << j << std::endl;
569 bench.last().smethod += std::string(
"(") + std::to_string(j) + std::string(
")");
571 bench_solve_simple(A_,b_,x_,bench.last());
573 copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
579 new_bench(solvs.get(i));
581 bench_solve_simple(A_,b_,x_,bench.last());
583 copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
589 write_bench_report();
592 PETSC_SAFE_CALL(VecCopy(best_sol,x_));
593 PETSC_SAFE_CALL(VecDestroy(&best_sol));
612 PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
613 PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress_residual,
this,NULL));
614 print_progress_bar();
615 solve_simple(A_,b_,x_);
618 if (create_vcluster().getProcessUnitID() == 0)
619 std::cout << std::endl;
622 bench.time = t.
getwct() * 1000.0;
625 solError err = statSolutionError(A_,b_,x_,ksp);
631 if (create_vcluster().getProcessUnitID() == 0)
632 std::cout << std::endl;
652 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
653 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
656 PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
660 PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
661 PETSC_SAFE_CALL(KSPSetUp(ksp));
664 PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
676 PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
693 err = getSolNormError(A_,b_,x_);
696 PETSC_SAFE_CALL(KSPGetIterationNumber(ksp,&its));
709 PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
718 PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
745 KSPGetOperators(ksp,&A_,&P_);
747 return getSolNormError(A_,b_,x_);
761 PetscScalar neg_one = -1.0;
769 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
770 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
780 Vec & r_ = r.getVec();
782 PETSC_SAFE_CALL(MatMult(A_,x_,r_));
783 PETSC_SAFE_CALL(VecAXPY(r_,neg_one,b_));
785 PETSC_SAFE_CALL(VecNorm(r_,NORM_1,&norm));
786 PETSC_SAFE_CALL(VecNorm(r_,NORM_INFINITY,&norm_inf));
803 PETSC_SAFE_CALL(KSPDestroy(&ksp));
813 solvs.add(std::string(KSPBCGS));
817 solvs.add(std::string(KSPBCGSL));
818 solvs.add(std::string(KSPGMRES));
819 solvs.add(std::string(KSPFGMRES));
820 solvs.add(std::string(KSPLGMRES));
833 auto & v_cl = create_vcluster();
835 if (v_cl.getProcessUnitID() == 0)
841 PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
842 PETSC_SAFE_CALL(PCGetType(pc,&type_pc));
844 std::cout <<
"The precoditioner is: " << type_pc << std::endl;
854 auto & v_cl = create_vcluster();
856 if (v_cl.getProcessUnitID() == 0)
861 PETSC_SAFE_CALL(KSPGetType(ksp,&type_ksp));
863 std::cout <<
"The ksp_type is: " << type_ksp << std::endl;
892 for (
size_t i = 0 ; i < solvs.
size() ; i++)
894 if (solvs.get(i) == solver)
907 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-ksp_monitor",0));
908 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_print_statistics",
"2"));
920 PetscOptionsSetValue(NULL,
"-ksp_type",type);
937 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
938 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol_,abstol,dtol,maxits));
955 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
956 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol_,dtol,maxits));
973 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
974 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol_,maxits));
989 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
990 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,n));
1006 PetscOptionsSetValue(NULL,
"-ksp_bcgsl_ell",std::to_string(l).c_str());
1016 PetscOptionsSetValue(NULL,
"-ksp_gmres_restart",std::to_string(n).c_str());
1048 is_preconditioner_set =
true;
1050 if (std::string(type) == PCHYPRE_BOOMERAMG)
1055 PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
1057 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_type",PCHYPRE));
1059 PETSC_SAFE_CALL(PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO));
1060 PETSC_SAFE_CALL(PCFactorSetShiftAmount(pc, PETSC_DECIDE));
1061 #ifdef PETSC_HAVE_HYPRE 1062 PETSC_SAFE_CALL(PCHYPRESetType(pc,
"boomeramg"));
1068 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_type",type));
1082 if (atype == HYPRE_AMG)
1084 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_max_levels",std::to_string(nl).c_str()));
1088 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1099 if (atype == HYPRE_AMG)
1101 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_max_iter",std::to_string(nit).c_str()));
1105 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1131 if (atype == HYPRE_AMG)
1134 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_relax_type_all",type.c_str()));}
1135 else if (k == REL_UP)
1136 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_relax_type_up",type.c_str()));}
1137 else if (k == REL_DOWN)
1138 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_relax_type_down",type.c_str()));}
1142 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1164 if (atype == HYPRE_AMG)
1166 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_cycle_type",cycle_type.c_str()));
1169 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_grid_sweeps_down",std::to_string(sweep_up).c_str()));}
1172 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_grid_sweeps_up",std::to_string(sweep_dw).c_str()));}
1174 if (sweep_crs != -1)
1175 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_grid_sweeps_coarse",std::to_string(sweep_crs).c_str()));}
1179 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1195 if (atype == HYPRE_AMG)
1197 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_coarsen_type",type.c_str()));
1201 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1217 if (atype == HYPRE_AMG)
1219 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_interp_type",type.c_str()));
1223 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1246 if (atype == HYPRE_AMG)
1248 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_nodal_coarsen",std::to_string(norm).c_str()));
1252 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1263 if (atype == HYPRE_AMG)
1265 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,
"-pc_hypre_boomeramg_eu_level",std::to_string(k).c_str()));
1269 std::cout << __FILE__ <<
":" << __LINE__ <<
"Warning: the selected preconditioner does not support this option" << std::endl;
1288 this->block_sz = block_sz;
1308 Mat & A_ = A.getMat();
1309 const Vec & b_ = b.getVec();
1317 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1318 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
1319 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
1322 Vec & x_ = x.getVec();
1324 pre_solve_impl(A_,b_,x_);
1325 solve_simple(A_,b_,x_);
1349 Mat & A_ = A.getMat();
1350 const Vec & b_ = b.getVec();
1357 MatNullSpace nullspace;
1360 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1361 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
1362 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
1366 Vec & x_ = x.getVec();
1369 PETSC_SAFE_CALL(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace));
1370 PETSC_SAFE_CALL(MatNullSpaceRemove(nullspace,b_));
1371 PETSC_SAFE_CALL(MatNullSpaceDestroy(&nullspace));
1373 pre_solve_impl(A_,b_,x_);
1374 solve_simple(A_,b_,x_);
1406 return getSolNormError(A.getMat(),b.getVec(),x.getVec());
1419 return getSolNormError(b.getVec(),x.getVec(),ksp);
1433 Mat & A_ = A.getMat();
1434 const Vec & b_ = b.getVec();
1435 Vec & x_ = x.getVec();
1437 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_TRUE));
1439 pre_solve_impl(A_,b_,x_);
1440 solve_simple(A_,b_,x_);
1455 const Vec & b_ = b.getVec();
1461 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1462 PETSC_SAFE_CALL(VecGetSize(b_,&row));
1463 PETSC_SAFE_CALL(VecGetLocalSize(b_,&row_loc));
1466 Vec & x_ = x.getVec();
1468 solve_simple(b_,x_);
1490 const Vec & b_ = b.getVec();
1491 Vec & x_ = x.getVec();
1493 solve_simple(b_,x_);
1509 PetscOptionsSetValue(NULL,name,value);
1525 Mat & A_ = A.getMat();
1526 const Vec & b_ = b.getVec();
1534 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1535 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
1536 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
1539 Vec & x_ = x.getVec();
1541 pre_solve_impl(A_,b_,x_);
1542 try_solve_simple(A_,b_,x_);
void setRelTol(PetscReal rtol_)
Set the relative tolerance as stop criteria.
solError err
Solution error.
Vector< double, PETSC_BASE > return_type
Type of the solution object.
void setPreconditioner(PCType type)
Set the preconditioner of the linear solver.
openfpm::vector< itError > res
Convergence per iteration.
void log_monitor()
Set the Petsc solver.
bool solve(Vector< double, PETSC_BASE > &x, const Vector< double, PETSC_BASE > &b)
Here we invert the matrix and solve the system.
void addHTML(const std::string &html)
Add HTML text.
static Vector< T > solve(const SparseMatrix< T, impl > &A, const Vector< T > &b)
Solve the linear system.
void print_stat(solError &err)
Print statistic about the solution error and method used.
double time
time to converge in milliseconds
void AddLinesGraph(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add a simple lines graph.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
std::string smethod
Method name (short form)
std::string title
Title of the chart.
void setMaxIter(PetscInt n)
Set the maximum number of iteration for Krylov solvers.
KSP getKSP()
Return the KSP solver.
void destroyKSP()
Destroy the KSP object.
In case T does not match the PETSC precision compilation create a stub structure.
Vector< double, PETSC_BASE > solve(const Vector< double, PETSC_BASE > &b)
Here we invert the matrix and solve the system.
It contain statistic of the error of the calculated solution.
std::string method
Solver Krylov name.
void bench_solve_simple(const Mat &A_, const Vec &b_, Vec &x_, solv_bench_info &bench)
Benchmark solve simple solving x=inv(A)*b.
void setPreconditionerAMG_cycleType(const std::string &cycle_type, int sweep_up=-1, int sweep_dw=-1, int sweep_crs=-1)
It set the type of cycle and optionally the number of sweep.
void initKSPForTest()
initialize the KSP object for solver testing
void setPreconditionerAMG_interp_eu_level(int k)
Indicate the number of levels in the ILU(k) for the Euclid smoother.
KSP ksp
Main parallel solver.
void setAbsTol(PetscReal abstol_)
Set the absolute tolerance as stop criteria.
Vector< double, PETSC_BASE > solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, bool initial_guess=false)
Here we invert the matrix and solve the system.
void setPreconditionerAMG_interp(const std::string &type)
Set the interpolation method for the algebraic-multi-grid preconditioner.
void solve_simple(const Vec &b_, Vec &x_)
solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
double getwct()
Return the elapsed real time.
Vector< double, PETSC_BASE > try_solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b)
Try to solve the system using all the solvers and generate a report.
std::string to_string_with_precision(const T myValue, const size_t n=6)
Converts value into string maintaining a desired precision.
void setPreconditionerAMG_maxit(int nit)
Set the maximum number of V or W cycle for algebraic-multi-grid.
Vector< double, PETSC_BASE > with_constant_nullspace_solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, bool initial_guess=false)
Here we invert the matrix and solve the system using a Nullspace for Neumann BC.
void setRestart(PetscInt n)
For GMRES based method, the number of Krylov directions to orthogonalize against.
std::string yAxis
Y axis name.
solError get_residual_error(const Vector< double, PETSC_BASE > &x, const Vector< double, PETSC_BASE > &b)
It return the resiual error.
void print_preconditioner()
Print the preconditioner used by the solver.
void progress(PetscInt it)
This function print an "*" showing the progress of the solvers.
openfpm::vector< solv_bench_info > bench
It contain the solver benchmark results.
void print_ksptype()
Print the ksp_type used by the solver.
void solve_simple(const Mat &A_, const Vec &b_, Vec &x_)
solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
Small class to produce graph with Google chart in HTML.
void copy_if_better(double res, Vec &sol, double &best_res, Vec &best_sol)
Copy the solution if better.
void try_solve_simple(const Mat &A_, const Vec &b_, Vec &x_)
Try to solve the system x=inv(A)*b using all the Krylov solvers with simple Jacobi pre-conditioner.
Sparse Matrix implementation.
void addTestSolver(std::string &solver)
Add a test solver.
void new_bench(const std::string &str)
Allocate a new benchmark slot for a method.
void setPreconditionerAMG_nl(int nl)
Set the number of levels for the algebraic-multigrid preconditioner.
This class is able to do Matrix inversion in parallel with PETSC solvers.
void setSolver(KSPType type)
Set the Petsc solver.
PetscReal err_norm
error norm at iteration it
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
void start()
Start the timer.
void pre_solve_impl(const Mat &A_, const Vec &b_, Vec &x_)
It set up the solver based on the provided options.
openfpm::vector< std::string > solvs
The full set of solvers.
size_t lineWidth
Width of the line.
static solError getSolNormError(const Mat &A_, const Vec &b_, const Vec &x_)
Return the norm error of the solution.
void setDivTol(PetscReal dtol_)
Set the divergence tolerance.
void print_progress_bar()
Print a progress bar on standard out.
static double calculate_it(double t, solv_bench_info &slv)
Calculate the residual error at time t for one method.
PetscInt its
Number of iterations.
void removeTestSolver(const std::string &solver)
Remove a test solver.
PetscInt maxits
KSP Maximum number of iterations.
void write(std::string file)
It write the graphs on file in html format using Google charts.
void setBlockSize(int block_sz)
Set how many degree of freedom each node has.
PetscReal err_norm
L1 norm of the error.
void write_bench_report()
Here we write the benchmark report.
void setPetscOption(const char *name, const char *value)
this function give you the possibility to set PETSC options
void setPreconditionerAMG_coarsen(const std::string &type)
Set the coarsening method for the algebraic-multi-grid preconditioner.
static solError getSolNormError(const Vec &b_, const Vec &x_, KSP ksp)
Return the norm error of the solution.
static PetscErrorCode monitor_progress_residual(KSP ksp, PetscInt it, PetscReal res, void *data)
procedure print the progress of the solver in benchmark mode
void initKSP()
initialize the KSP object
void setPreconditionerAMG_relax(const std::string &type, int k=REL_ALL)
Set the relaxation method for the algebraic-multi-grid preconditioner.
void AddHistGraph(openfpm::vector< Y > &y)
Add an histogram graph.
void searchDirections(PetscInt l)
solError get_residual_error(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &x, const Vector< double, PETSC_BASE > &b)
It return the resiual error.
bool solve(SparseMatrix< double, int, PETSC_BASE > &A, Vector< double, PETSC_BASE > &x, const Vector< double, PETSC_BASE > &b)
Here we invert the matrix and solve the system.
static solError statSolutionError(const Mat &A_, const Vec &b_, Vec &x_, KSP ksp)
Calculate statistic on the error solution.
std::string xAxis
X axis name.
size_t tmp
Temporal variable used for calculation of static members.
Class for cpu time benchmarking.
void setPreconditionerAMG_coarsenNodalType(int norm)
Set the block coarsening norm type.
void stop()
Stop the timer.
PetscReal err_inf
infinity norm of the error
std::string to_string_method(const std::string &solv)
It convert the KSP type into a human read-able string.