OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
petsc_solver.hpp
1 /*
2  * petsc_solver.hpp
3  *
4  * Created on: Apr 26, 2016
5  * Author: i-bird
6  */
7 
8 #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_
9 #define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_
10 
11 #include "config.h"
12 
13 #ifdef HAVE_PETSC
14 
15 #include "Vector/Vector.hpp"
16 #include <petscksp.h>
17 #include <petsctime.h>
18 #include "Plot/GoogleChart.hpp"
19 #include "Matrix/SparseMatrix.hpp"
20 #include "Vector/Vector.hpp"
21 #include <sstream>
22 #include <iomanip>
23 
24 template <typename T>
25 std::string to_string_with_precision(const T a_value, const int n = 6)
26 {
27  std::ostringstream out;
28  out << std::setprecision(n) << a_value;
29  return out.str();
30 }
31 
32 #define UMFPACK_NONE 0
33 #define REL_DOWN 1
34 #define REL_UP 2
35 #define REL_ALL 3
36 
37 #define PCHYPRE_BOOMERAMG "petsc_boomeramg"
38 
39 enum AMG_type
40 {
41  NONE_AMG,
42  HYPRE_AMG,
43  PETSC_AMG,
44  TRILINOS_ML
45 };
46 
47 
54 template<typename T>
56 {
57 public:
58 
69  template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b)
70  {
71  std::cerr << "Error Petsc only suppor double precision" << "/n";
72  }
73 };
74 
75 #define SOLVER_NOOPTION 0
76 #define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
77 #define SOLVER_PRINT_DETERMINANT 2
78 
80 struct solError
81 {
83  PetscReal err_inf;
84 
86  PetscReal err_norm;
87 
89  PetscInt its;
90 };
91 
92 
100 template<>
101 class petsc_solver<double>
102 {
104  struct itError
105  {
107  PetscInt it;
108 
110  PetscReal err_norm;
111  };
112 
117  struct solv_bench_info
118  {
120  std::string method;
121 
123  std::string smethod;
124 
126  double time;
127 
130 
133  };
134 
136  bool is_preconditioner_set = false;
137 
139  PetscInt maxits;
140 
142  KSP ksp;
143 
145  size_t tmp;
146 
149 
150 
153 
155  AMG_type atype = NONE_AMG;
156 
158  int block_sz = 0;
159 
168  static double calculate_it(double t, solv_bench_info & slv)
169  {
170  double s_int = slv.time / slv.res.size();
171 
172  // Calculate the discrete point in time
173  size_t pt = std::floor(t / s_int);
174 
175  if (pt < slv.res.size())
176  return slv.res.get(pt).err_norm;
177 
178  return slv.res.last().err_norm;
179  }
180 
186  {
187  if (create_vcluster().getProcessUnitID() != 0)
188  return;
189 
194 
195  for (size_t i = 0 ; i < bench.size() ; i++)
196  x.add(bench.get(i).smethod);
197 
198  // Each colum can have multiple data set (in this case 4 dataset)
199  // Each dataset can have a name
200  yn.add("Norm infinity");
201  yn.add("Norm average");
202  yn.add("Number of iterations");
203 
204  // Each colums can have multiple data-set
205 
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});
208 
209  // Google charts options
210  GCoptions options;
211 
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'}");
217 
218  GoogleChart cg;
219 
220  std::stringstream ss;
221 
222  cg.addHTML("<h2>Robustness</h2>");
223  cg.AddHistGraph(x,y,yn,options);
224  cg.addHTML("<h2>Speed</h2>");
225 
226  y.clear();
227  yn.clear();
228 
229  yn.add("Time in millseconds");
230  options.title = std::string("Time");
231 
232  for (size_t i = 0 ; i < bench.size() ; i++)
233  y.add({bench.get(i).time});
234 
235  cg.AddHistGraph(x,y,yn,options);
236 
237  // Convergence in time
238 
239  x.clear();
240  y.clear();
241  yn.clear();
242  xd.clear();
243 
244  // Get the maximum in time across all the solvers
245 
246  double max_time = 0.0;
247 
248  for (size_t i = 0 ; i < bench.size() ; i++)
249  {
250  if (bench.get(i).time > max_time) {max_time = bench.get(i).time;}
251  yn.add(bench.get(i).smethod);
252  }
253 
254  size_t n_int = maxits;
255 
256  // calculate dt
257 
258  double dt = max_time / n_int;
259 
260  // Add
261 
262  // For each solver we have a convergence plot
263 
264  for (double t = dt ; t <= max_time + 0.05 * max_time ; t += dt)
265  {
266  y.add();
267  xd.add(t);
268  for (size_t j = 0 ; j < bench.size() ; j++)
269  y.last().add(calculate_it(t,bench.get(j)));
270  }
271 
272  std::stringstream ss2;
273  ss2 << "<h2>Convergence with time</h2><br>" << std::endl;
274 
275  options.title = std::string("Residual norm infinity");
276  options.yAxis = std::string("residual");
277  options.xAxis = std::string("Time in milliseconds");
278  options.lineWidth = 1.0;
279  options.more = std::string("vAxis: {scaleType: 'log'},hAxis: {scaleType: 'log'}");
280 
281  cg.addHTML(ss2.str());
282  cg.AddLinesGraph(xd,y,yn,options);
283 
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;
287 
288  cg.addHTML(ss.str());
289 
290  cg.write("gc_solver.html");
291  }
292 
300  void pre_solve_impl(const Mat & A_, const Vec & b_, Vec & x_)
301  {
302  PETSC_SAFE_CALL(KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED));
303 
304  if (atype == HYPRE_AMG)
305  {
306  PC pc;
307 
308  // We set the pre-conditioner
309  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
310 
311  PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
312  PCFactorSetShiftAmount(pc, PETSC_DECIDE);
313  }
314  }
315 
321  {
322  auto & v_cl = create_vcluster();
323 
324  if (v_cl.getProcessUnitID() == 0)
325  std::cout << "-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
326  }
327 
328 
339  static PetscErrorCode monitor_progress_residual(KSP ksp,PetscInt it,PetscReal res,void* data)
340  {
341  petsc_solver<double> * pts = (petsc_solver *)data;
342 
343  pts->progress(it);
344 
345  itError erri;
346  erri.it = it;
347 
348  Mat A;
349  Vec Br,v,w;
350 
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));
356  itError err_fill;
357 
358  size_t old_size = pts->bench.last().res.size();
359  pts->bench.last().res.resize(it+1);
360 
361  if (old_size > 0)
362  err_fill = pts->bench.last().res.get(old_size-1);
363  else
364  err_fill = erri;
365 
366  for (long int i = old_size ; i < (long int)it ; i++)
367  pts->bench.last().res.get(i) = err_fill;
368 
369  // Add the error per iteration
370  pts->bench.last().res.get(it) = erri;
371 
372  return 0;
373  }
374 
380  void progress(PetscInt it)
381  {
382  PetscInt n_max_it;
383 
384  PETSC_SAFE_CALL(KSPGetTolerances(ksp,NULL,NULL,NULL,&n_max_it));
385 
386  auto & v_cl = create_vcluster();
387 
388  if (std::floor(it * 100.0 / n_max_it) > tmp)
389  {
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)
393  {
394  for (size_t k = 0 ; k < diff ; k++) {std::cout << "*";}
395  std::cout << std::flush;
396  }
397  }
398  }
399 
405  void print_stat(solError & err)
406  {
407  auto & v_cl = create_vcluster();
408 
409  if (v_cl.getProcessUnitID() == 0)
410  {
411  KSPType typ;
412  KSPGetType(ksp,&typ);
413 
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;
416  }
417  }
418 
419 
427  std::string to_string_method(const std::string & solv)
428  {
429  if (solv == std::string(KSPBCGS))
430  {
431  return std::string("BiCGStab (Stabilized version of BiConjugate Gradient Squared)");
432  }
433  else if (solv == std::string(KSPIBCGS))
434  {
435  return std::string("IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared)");
436  }
437  else if (solv == std::string(KSPFBCGS))
438  {
439  return std::string("FBiCGStab (Flexible Stabilized version of BiConjugate Gradient Squared)");
440  }
441  else if (solv == std::string(KSPFBCGSR))
442  {
443  return std::string("FBiCGStab-R (Modified Flexible Stabilized version of BiConjugate Gradient Squared)");
444  }
445  else if (solv == std::string(KSPBCGSL))
446  {
447  return std::string("BiCGStab(L) (Stabilized version of BiConjugate Gradient Squared with L search directions)");
448  }
449  else if (solv == std::string(KSPGMRES))
450  {
451  return std::string("GMRES(M) (Generalized Minimal Residual method with restart)");
452  }
453  else if (solv == std::string(KSPFGMRES))
454  {
455  return std::string("FGMRES(M) (Flexible Generalized Minimal Residual with restart)");
456  }
457  else if (solv == std::string(KSPLGMRES))
458  {
459  return std::string("LGMRES(M) (Augment Generalized Minimal Residual method with restart)");
460  }
461  else if (solv == std::string(KSPPGMRES))
462  {
463  return std::string("PGMRES(M) (Pipelined Generalized Minimal Residual method)");
464  }
465  else if (solv == std::string(KSPGCR))
466  {
467  return std::string("GCR (Generalized Conjugate Residual)");
468  }
469 
470  return std::string("UNKNOWN");
471  }
472 
478  void new_bench(const std::string & str)
479  {
480  // Add a new benchmark result
481  bench.add();
482  bench.last().method = to_string_method(str);
483  bench.last().smethod = std::string(str);
484  }
485 
494  void copy_if_better(double res, Vec & sol, double & best_res, Vec & best_sol)
495  {
496  if (res < best_res)
497  {
498  VecCopy(sol,best_sol);
499  best_res = res;
500  }
501  }
502 
513  void try_solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
514  {
515  Vec best_sol;
516  PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
517 
518  // Best residual
519  double best_res = std::numeric_limits<double>::max();
520 
521  // Create a new VCluster
522  auto & v_cl = create_vcluster();
523 
524  destroyKSP();
525 
526  // for each solver
527  for (size_t i = 0 ; i < solvs.size() ; i++)
528  {
529  initKSPForTest();
530 
531  // Here we solve without preconditioner
532  PC pc;
533 
534  // We set the pre-conditioner to none
535  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
536  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCNONE));
537 
538  setSolver(solvs.get(i).c_str());
539 
540  // Setup for BCGSL, GMRES
541  if (solvs.get(i) == std::string(KSPBCGSL))
542  {
543  // we try from 2 to 6 as search direction
544  for (size_t j = 2 ; j < 6 ; j++)
545  {
546  new_bench(solvs.get(i));
547 
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(")");
551  searchDirections(j);
552  bench_solve_simple(A_,b_,x_,bench.last());
553 
554  copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
555  }
556  }
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)) )
560  {
561  // we try from 2 to 6 as search direction
562  for (size_t j = 50 ; j < 300 ; j += 50)
563  {
564  // Add a new benchmark result
565  new_bench(solvs.get(i));
566 
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(")");
570  setRestart(j);
571  bench_solve_simple(A_,b_,x_,bench.last());
572 
573  copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
574  }
575  }
576  else
577  {
578  // Add a new benchmark result
579  new_bench(solvs.get(i));
580 
581  bench_solve_simple(A_,b_,x_,bench.last());
582 
583  copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
584  }
585 
586  destroyKSP();
587  }
588 
589  write_bench_report();
590 
591  // Copy the best solution to x
592  PETSC_SAFE_CALL(VecCopy(best_sol,x_));
593  PETSC_SAFE_CALL(VecDestroy(&best_sol));
594  }
595 
596 
605  void bench_solve_simple(const Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
606  {
607  // timer for benchmark
608  timer t;
609  t.start();
610  // Enable progress monitor
611  tmp = 0;
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_);
616 
617  // New line
618  if (create_vcluster().getProcessUnitID() == 0)
619  std::cout << std::endl;
620 
621  t.stop();
622  bench.time = t.getwct() * 1000.0;
623 
624  // calculate error statistic about the solution and print
625  solError err = statSolutionError(A_,b_,x_,ksp);
626  print_stat(err);
627 
628  bench.err = err;
629 
630  // New line
631  if (create_vcluster().getProcessUnitID() == 0)
632  std::cout << std::endl;
633  }
634 
642  void solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
643  {
644 // PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
645 
646  // We set the size of x according to the Matrix A
647  PetscInt row;
648  PetscInt col;
649  PetscInt row_loc;
650  PetscInt col_loc;
651 
652  PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
653  PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
654 
655  // We set the Matrix operators
656  PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
657 
658  // if we are on on best solve set-up a monitor function
659 
660  PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
661  PETSC_SAFE_CALL(KSPSetUp(ksp));
662 
663  // Solve the system
664  PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
665  }
666 
673  void solve_simple(const Vec & b_, Vec & x_)
674  {
675  // Solve the system
676  PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
677  }
678 
689  static solError statSolutionError(const Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
690  {
691  solError err;
692 
693  err = getSolNormError(A_,b_,x_);
694 
695  PetscInt its;
696  PETSC_SAFE_CALL(KSPGetIterationNumber(ksp,&its));
697  err.its = its;
698 
699  return err;
700  }
701 
702 
707  void initKSP()
708  {
709  PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
710  }
711 
717  {
718  PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
719 
720  setMaxIter(maxits);
721  }
722 
727  void destroyKSP()
728  {
729  KSPDestroy(&ksp);
730  }
731 
741  static solError getSolNormError(const Vec & b_, const Vec & x_,KSP ksp)
742  {
743  Mat A_;
744  Mat P_;
745  KSPGetOperators(ksp,&A_,&P_);
746 
747  return getSolNormError(A_,b_,x_);
748  }
749 
759  static solError getSolNormError(const Mat & A_, const Vec & b_, const Vec & x_)
760  {
761  PetscScalar neg_one = -1.0;
762 
763  // We set the size of x according to the Matrix A
764  PetscInt row;
765  PetscInt col;
766  PetscInt row_loc;
767  PetscInt col_loc;
768 
769  PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
770  PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
771 
772  /*
773  Here we calculate the residual error
774  */
775  PetscReal norm;
776  PetscReal norm_inf;
777 
778  // Get a vector r for the residual
779  Vector<double,PETSC_BASE> r(row,row_loc);
780  Vec & r_ = r.getVec();
781 
782  PETSC_SAFE_CALL(MatMult(A_,x_,r_));
783  PETSC_SAFE_CALL(VecAXPY(r_,neg_one,b_));
784 
785  PETSC_SAFE_CALL(VecNorm(r_,NORM_1,&norm));
786  PETSC_SAFE_CALL(VecNorm(r_,NORM_INFINITY,&norm_inf));
787 
788  solError err;
789  err.err_norm = norm / col;
790  err.err_inf = norm_inf;
791  err.its = 0;
792 
793  return err;
794  }
795 
796 public:
797 
800 
801  ~petsc_solver()
802  {
803  PETSC_SAFE_CALL(KSPDestroy(&ksp));
804  }
805 
806  petsc_solver()
807  :maxits(300),tmp(0)
808  {
809  initKSP();
810 
811  // Add the solvers
812 
813  solvs.add(std::string(KSPBCGS));
814 // solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
815 // solvs.add(std::string(KSPFBCGS));
816 // solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
817  solvs.add(std::string(KSPBCGSL));
818  solvs.add(std::string(KSPGMRES));
819  solvs.add(std::string(KSPFGMRES));
820  solvs.add(std::string(KSPLGMRES));
821 // solvs.add(std::string(KSPPGMRES)); <--- Take forever
822 // solvs.add(std::string(KSPGCR));
823 
824  //setSolver(KSPGMRES);
825  }
826 
832  {
833  auto & v_cl = create_vcluster();
834 
835  if (v_cl.getProcessUnitID() == 0)
836  {
837  PC pc;
838  PCType type_pc;
839 
840  // We set the pre-conditioner
841  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
842  PETSC_SAFE_CALL(PCGetType(pc,&type_pc));
843 
844  std::cout << "The precoditioner is: " << type_pc << std::endl;
845  }
846  }
847 
853  {
854  auto & v_cl = create_vcluster();
855 
856  if (v_cl.getProcessUnitID() == 0)
857  {
858  KSPType type_ksp;
859 
860  // We get the KSP_type
861  PETSC_SAFE_CALL(KSPGetType(ksp,&type_ksp));
862 
863  std::cout << "The ksp_type is: " << type_ksp << std::endl;
864  }
865  }
866 
875  void addTestSolver(std::string & solver)
876  {
877  solvs.add(solver);
878  }
879 
888  void removeTestSolver(const std::string & solver)
889  {
890  // search for the test solver and remove it
891 
892  for (size_t i = 0 ; i < solvs.size() ; i++)
893  {
894  if (solvs.get(i) == solver)
895  return;
896  }
897  }
898 
899 
905  void log_monitor()
906  {
907  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-ksp_monitor",0));
908  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_print_statistics","2"));
909  }
910 
918  void setSolver(KSPType type)
919  {
920  PetscOptionsSetValue(NULL,"-ksp_type",type);
921  }
922 
930  void setRelTol(PetscReal rtol_)
931  {
932  PetscReal rtol;
933  PetscReal abstol;
934  PetscReal dtol;
935  PetscInt maxits;
936 
937  PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
938  PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol_,abstol,dtol,maxits));
939  }
940 
948  void setAbsTol(PetscReal abstol_)
949  {
950  PetscReal rtol;
951  PetscReal abstol;
952  PetscReal dtol;
953  PetscInt maxits;
954 
955  PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
956  PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol_,dtol,maxits));
957  }
958 
966  void setDivTol(PetscReal dtol_)
967  {
968  PetscReal rtol;
969  PetscReal abstol;
970  PetscReal dtol;
971  PetscInt maxits;
972 
973  PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
974  PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol_,maxits));
975  }
976 
982  void setMaxIter(PetscInt n)
983  {
984  PetscReal rtol;
985  PetscReal abstol;
986  PetscReal dtol;
987  PetscInt maxits;
988 
989  PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
990  PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,n));
991 
992  this->maxits = n;
993  }
994 
1004  void searchDirections(PetscInt l)
1005  {
1006  PetscOptionsSetValue(NULL,"-ksp_bcgsl_ell",std::to_string(l).c_str());
1007  }
1008 
1014  void setRestart(PetscInt n)
1015  {
1016  PetscOptionsSetValue(NULL,"-ksp_gmres_restart",std::to_string(n).c_str());
1017  }
1018 
1046  void setPreconditioner(PCType type)
1047  {
1048  is_preconditioner_set = true;
1049 
1050  if (std::string(type) == PCHYPRE_BOOMERAMG)
1051  {
1052  PC pc;
1053 
1054  // We set the pre-conditioner
1055  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
1056 
1057  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCHYPRE));
1058 
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"));
1063 #endif
1064  atype = HYPRE_AMG;
1065  }
1066  else
1067  {
1068  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",type));
1069  }
1070  }
1071 
1081  {
1082  if (atype == HYPRE_AMG)
1083  {
1084  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_max_levels",std::to_string(nl).c_str()));
1085  }
1086  else
1087  {
1088  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1089  }
1090  }
1091 
1098  {
1099  if (atype == HYPRE_AMG)
1100  {
1101  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_max_iter",std::to_string(nit).c_str()));
1102  }
1103  else
1104  {
1105  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1106  }
1107  }
1108 
1109 
1129  void setPreconditionerAMG_relax(const std::string & type, int k = REL_ALL)
1130  {
1131  if (atype == HYPRE_AMG)
1132  {
1133  if (k == REL_ALL)
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()));}
1139  }
1140  else
1141  {
1142  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1143  }
1144  }
1145 
1162  void setPreconditionerAMG_cycleType(const std::string & cycle_type, int sweep_up = -1, int sweep_dw = -1, int sweep_crs = -1)
1163  {
1164  if (atype == HYPRE_AMG)
1165  {
1166  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_cycle_type",cycle_type.c_str()));
1167 
1168  if (sweep_dw != -1)
1169  {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_down",std::to_string(sweep_up).c_str()));}
1170 
1171  if (sweep_up != -1)
1172  {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_up",std::to_string(sweep_dw).c_str()));}
1173 
1174  if (sweep_crs != -1)
1175  {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_coarse",std::to_string(sweep_crs).c_str()));}
1176  }
1177  else
1178  {
1179  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1180  }
1181  }
1182 
1193  void setPreconditionerAMG_coarsen(const std::string & type)
1194  {
1195  if (atype == HYPRE_AMG)
1196  {
1197  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_coarsen_type",type.c_str()));
1198  }
1199  else
1200  {
1201  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1202  }
1203  }
1204 
1215  void setPreconditionerAMG_interp(const std::string & type)
1216  {
1217  if (atype == HYPRE_AMG)
1218  {
1219  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_interp_type",type.c_str()));
1220  }
1221  else
1222  {
1223  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1224  }
1225  }
1226 
1245  {
1246  if (atype == HYPRE_AMG)
1247  {
1248  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_nodal_coarsen",std::to_string(norm).c_str()));
1249  }
1250  else
1251  {
1252  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1253  }
1254  }
1255 
1262  {
1263  if (atype == HYPRE_AMG)
1264  {
1265  PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_eu_level",std::to_string(k).c_str()));
1266  }
1267  else
1268  {
1269  std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1270  }
1271  }
1272 
1273 
1274 
1286  void setBlockSize(int block_sz)
1287  {
1288  this->block_sz = block_sz;
1289  }
1290 
1307  {
1308  Mat & A_ = A.getMat();
1309  const Vec & b_ = b.getVec();
1310 
1311  // We set the size of x according to the Matrix A
1312  PetscInt row;
1313  PetscInt col;
1314  PetscInt row_loc;
1315  PetscInt col_loc;
1316 
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));
1320 
1321  Vector<double,PETSC_BASE> x(row,row_loc);
1322  Vec & x_ = x.getVec();
1323 
1324  pre_solve_impl(A_,b_,x_);
1325  solve_simple(A_,b_,x_);
1326 
1327  x.update();
1328 
1329  return x;
1330  }
1331 
1348  {
1349  Mat & A_ = A.getMat();
1350  const Vec & b_ = b.getVec();
1351 
1352  // We set the size of x according to the Matrix A
1353  PetscInt row;
1354  PetscInt col;
1355  PetscInt row_loc;
1356  PetscInt col_loc;
1357  MatNullSpace nullspace;
1358 
1359 
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));
1363 
1364 
1365  Vector<double,PETSC_BASE> x(row,row_loc);
1366  Vec & x_ = x.getVec();
1367 
1368  //Removing Null Space from RHS
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));
1372 
1373  pre_solve_impl(A_,b_,x_);
1374  solve_simple(A_,b_,x_);
1375 
1376  x.update();
1377 
1378  return x;
1379  }
1380 
1390  KSP getKSP()
1391  {
1392  return ksp;
1393  }
1394 
1405  {
1406  return getSolNormError(A.getMat(),b.getVec(),x.getVec());
1407  }
1408 
1418  {
1419  return getSolNormError(b.getVec(),x.getVec(),ksp);
1420  }
1421 
1432  {
1433  Mat & A_ = A.getMat();
1434  const Vec & b_ = b.getVec();
1435  Vec & x_ = x.getVec();
1436 
1437  PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_TRUE));
1438 
1439  pre_solve_impl(A_,b_,x_);
1440  solve_simple(A_,b_,x_);
1441  x.update();
1442 
1443  return true;
1444  }
1445 
1454  {
1455  const Vec & b_ = b.getVec();
1456 
1457  // We set the size of x according to the Matrix A
1458  PetscInt row;
1459  PetscInt row_loc;
1460 
1461  PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1462  PETSC_SAFE_CALL(VecGetSize(b_,&row));
1463  PETSC_SAFE_CALL(VecGetLocalSize(b_,&row_loc));
1464 
1465  Vector<double,PETSC_BASE> x(row,row_loc);
1466  Vec & x_ = x.getVec();
1467 
1468  solve_simple(b_,x_);
1469  x.update();
1470 
1471  return x;
1472  }
1473 
1489  {
1490  const Vec & b_ = b.getVec();
1491  Vec & x_ = x.getVec();
1492 
1493  solve_simple(b_,x_);
1494  x.update();
1495 
1496  return true;
1497  }
1498 
1507  void setPetscOption(const char * name, const char * value)
1508  {
1509  PetscOptionsSetValue(NULL,name,value);
1510  }
1511 
1524  {
1525  Mat & A_ = A.getMat();
1526  const Vec & b_ = b.getVec();
1527 
1528  // We set the size of x according to the Matrix A
1529  PetscInt row;
1530  PetscInt col;
1531  PetscInt row_loc;
1532  PetscInt col_loc;
1533 
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));
1537 
1538  Vector<double,PETSC_BASE> x(row,row_loc);
1539  Vec & x_ = x.getVec();
1540 
1541  pre_solve_impl(A_,b_,x_);
1542  try_solve_simple(A_,b_,x_);
1543 
1544  x.update();
1545 
1546  return x;
1547  }
1548 };
1549 
1550 #endif
1551 
1552 #endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ */
void setRelTol(PetscReal rtol_)
Set the relative tolerance as stop criteria.
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.
Definition: Vector.hpp:39
std::string smethod
Method name (short form)
std::string title
Title of the chart.
Definition: GoogleChart.hpp:28
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.
std::string more
more
Definition: GoogleChart.hpp:67
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.
size_t size()
Stub size.
Definition: map_vector.hpp:211
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.
Definition: timer.hpp:130
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.
Definition: GoogleChart.hpp:30
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.
Definition: timer.hpp:90
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.
Definition: GoogleChart.hpp:56
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.
std::string stype
Definition: GoogleChart.hpp:38
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.
Definition: GoogleChart.hpp:32
Google chart options.
Definition: GoogleChart.hpp:25
size_t tmp
Temporal variable used for calculation of static members.
Class for cpu time benchmarking.
Definition: timer.hpp:27
void setPreconditionerAMG_coarsenNodalType(int norm)
Set the block coarsening norm type.
void stop()
Stop the timer.
Definition: timer.hpp:119
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.