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
24template <typename T>
25std::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
39enum AMG_type
40{
41 NONE_AMG,
42 HYPRE_AMG,
43 PETSC_AMG,
44 TRILINOS_ML
45};
46
47
48/*! \brief In case T does not match the PETSC precision compilation create a
49 * stub structure
50 *
51 * \param T precision
52 *
53 */
54template<typename T>
55class petsc_solver
56{
57public:
58
59 /*! \brief Solve the linear system.
60 *
61 * In this case just return an error
62 *
63 * \param A sparse matrix
64 * \param b right-hand-side
65 *
66 * \return the solution
67 *
68 */
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
79//! It contain statistic of the error of the calculated solution
80struct solError
81{
82 //! infinity norm of the error
83 PetscReal err_inf;
84
85 //! L1 norm of the error
86 PetscReal err_norm;
87
88 //! Number of iterations
89 PetscInt its;
90};
91
92
93/*! \brief This class is able to do Matrix inversion in parallel with PETSC solvers
94 *
95 * # Example of use of this class #
96 *
97 * \snippet eq_unit_test.hpp lid-driven cavity 2D
98 *
99 */
100template<>
101class petsc_solver<double>
102{
103 //! contain the infinity norm of the residual at each iteration
104 struct itError
105 {
106 //! Iteration
107 PetscInt it;
108
109 //! error norm at iteration it
110 PetscReal err_norm;
111 };
112
113 /*! \brief It contain the benchmark information for each solver
114 *
115 *
116 */
117 struct solv_bench_info
118 {
119 //! Solver Krylov name
120 std::string method;
121
122 //! Method name (short form)
123 std::string smethod;
124
125 //! time to converge in milliseconds
126 double time;
127
128 //! Solution error
129 solError err;
130
131 //! Convergence per iteration
132 openfpm::vector<itError> res;
133 };
134
135 //! indicate if the preconditioner is set
136 bool is_preconditioner_set = false;
137
138 //! KSP Maximum number of iterations
139 PetscInt maxits;
140
141 //! Main parallel solver
142 KSP ksp;
143
144 //! Temporal variable used for calculation of static members
145 size_t tmp;
146
147 //! The full set of solvers
148 openfpm::vector<std::string> solvs;
149
150
151 //! It contain the solver benchmark results
152 openfpm::vector<solv_bench_info> bench;
153
154 //! Type of the algebraic multi-grid preconditioner
155 AMG_type atype = NONE_AMG;
156
157 //! Block size
158 int block_sz = 0;
159
160 /*! \brief Calculate the residual error at time t for one method
161 *
162 * \param t time
163 * \param slv solver
164 *
165 * \return the residual
166 *
167 */
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
181 /*! \brief Here we write the benchmark report
182 *
183 *
184 */
185 void write_bench_report()
186 {
187 if (create_vcluster().getProcessUnitID() != 0)
188 return;
189
190 openfpm::vector<std::string> x;
191 openfpm::vector<openfpm::vector<double>> y;
192 openfpm::vector<std::string> yn;
193 openfpm::vector<double> xd;
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
293 /*! \brief It set up the solver based on the provided options
294 *
295 * \param A_ Matrix
296 * \param x_ solution
297 * \param b_ right-hand-side
298 *
299 */
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
316 /*! \brief Print a progress bar on standard out
317 *
318 *
319 */
320 void print_progress_bar()
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
329 /*! \brief procedure print the progress of the solver in benchmark mode
330 *
331 * \param ksp Solver
332 * \param it Iteration number
333 * \param res resudual
334 * \param data custom pointer to data
335 *
336 * \return always zero
337 *
338 */
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
375 /*! \brief This function print an "*" showing the progress of the solvers
376 *
377 * \param it iteration number
378 *
379 */
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
400 /*! \brief Print statistic about the solution error and method used
401 *
402 * \param err structure that contain the solution errors
403 *
404 */
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 /*! \brief It convert the KSP type into a human read-able string
420 *
421 * \param solv solver (short form)
422 *
423 * \return the name of the solver in long form
424 *
425 */
426 std::string to_string_method(const std::string & solv)
427 {
428 if (solv == std::string(KSPBCGS))
429 {
430 return std::string("BiCGStab (Stabilized version of BiConjugate Gradient Squared)");
431 }
432 else if (solv == std::string(KSPIBCGS))
433 {
434 return std::string("IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared)");
435 }
436 else if (solv == std::string(KSPFBCGS))
437 {
438 return std::string("FBiCGStab (Flexible Stabilized version of BiConjugate Gradient Squared)");
439 }
440 else if (solv == std::string(KSPFBCGSR))
441 {
442 return std::string("FBiCGStab-R (Modified Flexible Stabilized version of BiConjugate Gradient Squared)");
443 }
444 else if (solv == std::string(KSPBCGSL))
445 {
446 return std::string("BiCGStab(L) (Stabilized version of BiConjugate Gradient Squared with L search directions)");
447 }
448 else if (solv == std::string(KSPGMRES))
449 {
450 return std::string("GMRES(M) (Generalized Minimal Residual method with restart)");
451 }
452 else if (solv == std::string(KSPFGMRES))
453 {
454 return std::string("FGMRES(M) (Flexible Generalized Minimal Residual with restart)");
455 }
456 else if (solv == std::string(KSPLGMRES))
457 {
458 return std::string("LGMRES(M) (Augment Generalized Minimal Residual method with restart)");
459 }
460 else if (solv == std::string(KSPPGMRES))
461 {
462 return std::string("PGMRES(M) (Pipelined Generalized Minimal Residual method)");
463 }
464 else if (solv == std::string(KSPGCR))
465 {
466 return std::string("GCR (Generalized Conjugate Residual)");
467 }
468
469 return std::string("UNKNOWN");
470 }
471
472 /*! \brief Allocate a new benchmark slot for a method
473 *
474 * \param str Method name
475 *
476 */
477 void new_bench(const std::string & str)
478 {
479 // Add a new benchmark result
480 bench.add();
481 bench.last().method = to_string_method(str);
482 bench.last().smethod = std::string(str);
483 }
484
485 /*! \brief Copy the solution if better
486 *
487 * \param res Residual of the solution
488 * \param sol solution
489 * \param best_res residual of the best solution
490 * \param best_sol best solution
491 *
492 */
493 void copy_if_better(double res, Vec & sol, double & best_res, Vec & best_sol)
494 {
495 if (res < best_res)
496 {
497 VecCopy(sol,best_sol);
498 best_res = res;
499 }
500 }
501
502 /*! \brief Try to solve the system x=inv(A)*b using all the Krylov solvers with simple Jacobi pre-conditioner
503 *
504 * It try to solve the system using JACOBI pre-conditioner and all the Krylov solvers available at the end it write
505 * a report
506 *
507 * \param A_ Matrix
508 * \param b_ vector of coefficents
509 * \param x_ solution
510 *
511 */
512 void try_solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
513 {
514 Vec best_sol;
515 PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
516
517 // Best residual
518 double best_res = std::numeric_limits<double>::max();
519
520 // Create a new VCluster
521 auto & v_cl = create_vcluster();
522
523 destroyKSP();
524
525 // for each solver
526 for (size_t i = 0 ; i < solvs.size() ; i++)
527 {
528 initKSPForTest();
529
530 // Here we solve without preconditioner
531 PC pc;
532
533 // We set the pre-conditioner to none
534 PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
535 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCNONE));
536
537 setSolver(solvs.get(i).c_str());
538
539 // Setup for BCGSL, GMRES
540 if (solvs.get(i) == std::string(KSPBCGSL))
541 {
542 // we try from 2 to 6 as search direction
543 for (size_t j = 2 ; j < 6 ; j++)
544 {
545 new_bench(solvs.get(i));
546
547 if (v_cl.getProcessUnitID() == 0)
548 std::cout << "L = " << j << std::endl;
549 bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
550 searchDirections(j);
551 bench_solve_simple(A_,b_,x_,bench.last());
552
553 copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
554 }
555 }
556 else if (solvs.get(i) == std::string(KSPGMRES) ||
557 solvs.get(i) == std::string(std::string(KSPFGMRES)) ||
558 solvs.get(i) == std::string(std::string(KSPLGMRES)) )
559 {
560 // we try from 2 to 6 as search direction
561 for (size_t j = 50 ; j < 300 ; j += 50)
562 {
563 // Add a new benchmark result
564 new_bench(solvs.get(i));
565
566 if (v_cl.getProcessUnitID() == 0)
567 std::cout << "Restart = " << j << std::endl;
568 bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
569 setRestart(j);
570 bench_solve_simple(A_,b_,x_,bench.last());
571
572 copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
573 }
574 }
575 else
576 {
577 // Add a new benchmark result
578 new_bench(solvs.get(i));
579
580 bench_solve_simple(A_,b_,x_,bench.last());
581
582 copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
583 }
584
585 destroyKSP();
586 }
587
588 write_bench_report();
589
590 // Copy the best solution to x
591 PETSC_SAFE_CALL(VecCopy(best_sol,x_));
592 PETSC_SAFE_CALL(VecDestroy(&best_sol));
593 }
594
595
596 /*! \brief Benchmark solve simple solving x=inv(A)*b
597 *
598 * \param A_ Matrix A
599 * \param b_ vector b
600 * \param x_ solution x
601 * \param bench structure that store the benchmark information
602 *
603 */
604 void bench_solve_simple(const Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
605 {
606 // timer for benchmark
607 timer t;
608 t.start();
609 // Enable progress monitor
610 tmp = 0;
611 PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
612 PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress_residual,this,NULL));
613 print_progress_bar();
614 solve_simple(A_,b_,x_);
615
616 // New line
617 if (create_vcluster().getProcessUnitID() == 0)
618 std::cout << std::endl;
619
620 t.stop();
621 bench.time = t.getwct() * 1000.0;
622
623 // calculate error statistic about the solution and print
624 solError err = statSolutionError(A_,b_,x_,ksp);
625 print_stat(err);
626
627 bench.err = err;
628
629 // New line
630 if (create_vcluster().getProcessUnitID() == 0)
631 std::cout << std::endl;
632 }
633
634 /*! \brief solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
635 *
636 * \param A_ SparseMatrix
637 * \param b_ right-hand-side
638 * \param x_ solution
639 *
640 */
641 void solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
642 {
643// PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
644
645 // We set the size of x according to the Matrix A
646 PetscInt row;
647 PetscInt col;
648 PetscInt row_loc;
649 PetscInt col_loc;
650
651 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
652 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
653
654 // We set the Matrix operators
655 PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
656
657 // if we are on on best solve set-up a monitor function
658
659 PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
660 PETSC_SAFE_CALL(KSPSetUp(ksp));
661
662 // Solve the system
663 PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
664 }
665
666 /*! \brief solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
667 *
668 * \param b_ right hand side
669 * \param x_ solution
670 *
671 */
672 void solve_simple(const Vec & b_, Vec & x_)
673 {
674 // Solve the system
675 PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
676 }
677
678 /*! \brief Calculate statistic on the error solution
679 *
680 * \param A_ Matrix of the system
681 * \param b_ Right hand side of the matrix
682 * \param x_ Solution
683 * \param ksp Krylov solver
684 *
685 * \return the solution error
686 *
687 */
688 static solError statSolutionError(const Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
689 {
690 solError err;
691
692 err = getSolNormError(A_,b_,x_);
693
694 PetscInt its;
695 PETSC_SAFE_CALL(KSPGetIterationNumber(ksp,&its));
696 err.its = its;
697
698 return err;
699 }
700
701
702 /*! \brief initialize the KSP object
703 *
704 *
705 */
706 void initKSP()
707 {
708 PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
709 }
710
711 /*! \brief initialize the KSP object for solver testing
712 *
713 *
714 */
715 void initKSPForTest()
716 {
717 PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
718
719 setMaxIter(maxits);
720 }
721
722 /*! \brief Destroy the KSP object
723 *
724 *
725 */
726 void destroyKSP()
727 {
728 KSPDestroy(&ksp);
729 }
730
731 /*! \brief Return the norm error of the solution
732 *
733 * \param x_ the solution
734 * \param b_ the right-hand-side
735 * \param ksp Krylov solver
736 *
737 * \return the solution error
738 *
739 */
740 static solError getSolNormError(const Vec & b_, const Vec & x_,KSP ksp)
741 {
742 Mat A_;
743 Mat P_;
744 KSPGetOperators(ksp,&A_,&P_);
745
746 return getSolNormError(A_,b_,x_);
747 }
748
749 /*! \brief Return the norm error of the solution
750 *
751 * \param A_ the matrix that identity the linear system
752 * \param x_ the solution
753 * \param b_ the right-hand-side
754 *
755 * \return the solution error
756 *
757 */
758 static solError getSolNormError(const Mat & A_, const Vec & b_, const Vec & x_)
759 {
760 PetscScalar neg_one = -1.0;
761
762 // We set the size of x according to the Matrix A
763 PetscInt row;
764 PetscInt col;
765 PetscInt row_loc;
766 PetscInt col_loc;
767
768 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
769 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
770
771 /*
772 Here we calculate the residual error
773 */
774 PetscReal norm;
775 PetscReal norm_inf;
776
777 // Get a vector r for the residual
778 Vector<double,PETSC_BASE> r(row,row_loc);
779 Vec & r_ = r.getVec();
780
781 PETSC_SAFE_CALL(MatMult(A_,x_,r_));
782 PETSC_SAFE_CALL(VecAXPY(r_,neg_one,b_));
783
784 PETSC_SAFE_CALL(VecNorm(r_,NORM_1,&norm));
785 PETSC_SAFE_CALL(VecNorm(r_,NORM_INFINITY,&norm_inf));
786
787 solError err;
788 err.err_norm = norm / col;
789 err.err_inf = norm_inf;
790 err.its = 0;
791
792 return err;
793 }
794
795public:
796
797 //! Type of the solution object
798 typedef Vector<double,PETSC_BASE> return_type;
799
800 ~petsc_solver()
801 {
802 PETSC_SAFE_CALL(KSPDestroy(&ksp));
803 }
804
805 petsc_solver()
806 :maxits(300),tmp(0)
807 {
808 initKSP();
809
810 // Add the solvers
811
812 solvs.add(std::string(KSPBCGS));
813// solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
814// solvs.add(std::string(KSPFBCGS));
815// solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
816 solvs.add(std::string(KSPBCGSL));
817 solvs.add(std::string(KSPGMRES));
818 solvs.add(std::string(KSPFGMRES));
819 solvs.add(std::string(KSPLGMRES));
820// solvs.add(std::string(KSPPGMRES)); <--- Take forever
821// solvs.add(std::string(KSPGCR));
822
823 setSolver(KSPGMRES);
824 }
825
826 /*! \brief Add a test solver
827 *
828 * The try solve function use the most robust solvers in PETSC, if you want to add
829 * additionally solver like KSPIBCGS,KSPFBCGSR,KSPPGMRES, use addTestSolver(std::string(KSPIBCGS))
830 *
831 * \param solver additional solver solver to test
832 *
833 */
834 void addTestSolver(std::string & solver)
835 {
836 solvs.add(solver);
837 }
838
839 /*! \brief Remove a test solver
840 *
841 * The try solve function use the most robust solvers in PETSC, if you want to remove
842 * a solver like use removeTestSolver(std::string(KSPIBCGS))
843 *
844 * \param solver remove solver to test
845 *
846 */
847 void removeTestSolver(const std::string & solver)
848 {
849 // search for the test solver and remove it
850
851 for (size_t i = 0 ; i < solvs.size() ; i++)
852 {
853 if (solvs.get(i) == solver)
854 return;
855 }
856 }
857
858
859 /*! \brief Set the Petsc solver
860 *
861 * \see KSPType in PETSC manual for a list of all PETSC solvers
862 *
863 */
864 void log_monitor()
865 {
866 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-ksp_monitor",0));
867 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_print_statistics","2"));
868 }
869
870 /*! \brief Set the Petsc solver
871 *
872 * \see KSPType in PETSC manual for a list of all PETSC solvers
873 *
874 * \param type petsc solver type
875 *
876 */
877 void setSolver(KSPType type)
878 {
879 PetscOptionsSetValue(NULL,"-ksp_type",type);
880 }
881
882 /*! \brief Set the relative tolerance as stop criteria
883 *
884 * \see PETSC manual KSPSetTolerances for an explanation
885 *
886 * \param rtol_ Relative tolerance
887 *
888 */
889 void setRelTol(PetscReal rtol_)
890 {
891 PetscReal rtol;
892 PetscReal abstol;
893 PetscReal dtol;
894 PetscInt maxits;
895
896 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
897 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol_,abstol,dtol,maxits));
898 }
899
900 /*! \brief Set the absolute tolerance as stop criteria
901 *
902 * \see PETSC manual KSPSetTolerances for an explanation
903 *
904 * \param abstol_ Absolute tolerance
905 *
906 */
907 void setAbsTol(PetscReal abstol_)
908 {
909 PetscReal rtol;
910 PetscReal abstol;
911 PetscReal dtol;
912 PetscInt maxits;
913
914 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
915 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol_,dtol,maxits));
916 }
917
918 /*! \brief Set the divergence tolerance
919 *
920 * \see PETSC manual KSPSetTolerances for an explanation
921 *
922 * \param dtol_ tolerance
923 *
924 */
925 void setDivTol(PetscReal dtol_)
926 {
927 PetscReal rtol;
928 PetscReal abstol;
929 PetscReal dtol;
930 PetscInt maxits;
931
932 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
933 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol_,maxits));
934 }
935
936 /*! \brief Set the maximum number of iteration for Krylov solvers
937 *
938 * \param n maximum number of iterations
939 *
940 */
941 void setMaxIter(PetscInt n)
942 {
943 PetscReal rtol;
944 PetscReal abstol;
945 PetscReal dtol;
946 PetscInt maxits;
947
948 PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
949 PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,n));
950
951 this->maxits = n;
952 }
953
954 /*! For the BiCGStab(L) it define the number of search directions
955 *
956 * BiCG Methods can fail for base break-down (or near break-down). BiCGStab(L) try
957 * to avoid this problem. Such method has a parameter L (2 by default) Bigger is L
958 * more the system will try to avoid the breakdown
959 *
960 * \param l Increasing L should reduce the probability of failure of the solver because of break-down of the base
961 *
962 */
963 void searchDirections(PetscInt l)
964 {
965 PetscOptionsSetValue(NULL,"-ksp_bcgsl_ell",std::to_string(l).c_str());
966 }
967
968 /*! \brief For GMRES based method, the number of Krylov directions to orthogonalize against
969 *
970 * \param n number of directions
971 *
972 */
973 void setRestart(PetscInt n)
974 {
975 PetscOptionsSetValue(NULL,"-ksp_gmres_restart",std::to_string(n).c_str());
976 }
977
978 /*! \brief Set the preconditioner of the linear solver
979 *
980 * The preconditoner that can be set are the same as PETSC
981 *
982 * For a full list please visit
983 * Please visit: http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCType.html#PCType
984 *
985 * An exception is PCHYPRE with BOOMERAMG in this case use
986 * PCHYPRE_BOOMERAMG. Many preconditioners has default values, but many
987 * times the default values are not good. Here we list some interesting case
988 *
989 * ## Algebraic-multi-grid based preconditioners##
990 *
991 * Parameters for this type of preconditioner can be set using
992 * setPreconditionerAMG_* functions.
993 *
994 * * Number of levels set by setPreconditionerAMG_nl
995 * * Maximum number of cycles setPreconditionerAMG_maxit
996 * * Smooth or relax method setPreconditionerAMG_smooth,setPreconditionerAMG_relax
997 * * Cycle type and number of sweep (relaxation steps) when going up or down setPreconditionerAMG_cycleType
998 * * Coarsening method setPreconditionerAMG_coarsen
999 * * interpolation schemes setPreconditionerAMG_interp
1000 *
1001 *
1002 * \param type of the preconditioner
1003 *
1004 */
1005 void setPreconditioner(PCType type)
1006 {
1007 is_preconditioner_set = true;
1008
1009 if (std::string(type) == PCHYPRE_BOOMERAMG)
1010 {
1011 PC pc;
1012
1013 // We set the pre-conditioner
1014 PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
1015
1016 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCHYPRE));
1017
1018 PETSC_SAFE_CALL(PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO));
1019 PETSC_SAFE_CALL(PCFactorSetShiftAmount(pc, PETSC_DECIDE));
1020 PETSC_SAFE_CALL(PCHYPRESetType(pc, "boomeramg"));
1021 atype = HYPRE_AMG;
1022 }
1023 else
1024 {
1025 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCHYPRE));
1026 }
1027 }
1028
1029 /*! \brief Set the number of levels for the algebraic-multigrid preconditioner
1030 *
1031 * In case you select an algebraic preconditioner like PCHYPRE or PCGAMG you can
1032 * set the number of levels using this function
1033 *
1034 * \param nl number of levels
1035 *
1036 */
1037 void setPreconditionerAMG_nl(int nl)
1038 {
1039 if (atype == HYPRE_AMG)
1040 {
1041 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_max_levels",std::to_string(nl).c_str()));
1042 }
1043 else
1044 {
1045 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1046 }
1047 }
1048
1049 /*! \brief Set the maximum number of V or W cycle for algebraic-multi-grid
1050 *
1051 * \param nit number of levels
1052 *
1053 */
1054 void setPreconditionerAMG_maxit(int nit)
1055 {
1056 if (atype == HYPRE_AMG)
1057 {
1058 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_max_iter",std::to_string(nit).c_str()));
1059 }
1060 else
1061 {
1062 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1063 }
1064 }
1065
1066
1067 /*! \brief Set the relaxation method for the algebraic-multi-grid preconditioner
1068 *
1069 * Possible values for relazation can be
1070 * "Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel",
1071 * "SOR/Jacobi","backward-SOR/Jacobi",hybrid chaotic Gauss-Seidel (works only with OpenMP),
1072 * "symmetric-SOR/Jacobi","l1scaled-SOR/Jacobi","Gaussian-elimination","CG","Chebyshev",
1073 * "FCF-Jacobi","l1scaled-Jacobi"
1074 *
1075 *
1076 *
1077 *
1078 * Every smooth operator can have additional parameters to be set in order to
1079 * correctly work.
1080 *
1081 *
1082 * \param type of relax method
1083 * \param k where is applied REL_UP,REL_DOWN,REL_ALL (default is all)
1084 *
1085 */
1086 void setPreconditionerAMG_relax(const std::string & type, int k = REL_ALL)
1087 {
1088 if (atype == HYPRE_AMG)
1089 {
1090 if (k == REL_ALL)
1091 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_relax_type_all",type.c_str()));}
1092 else if (k == REL_UP)
1093 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_relax_type_up",type.c_str()));}
1094 else if (k == REL_DOWN)
1095 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_relax_type_down",type.c_str()));}
1096 }
1097 else
1098 {
1099 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1100 }
1101 }
1102
1103 /*! \brief It set the type of cycle and optionally the number of sweep
1104 *
1105 * This function set the cycle type for the multigrid methods.
1106 * Possible values are:
1107 * * V cycle
1108 * * W cycle
1109 *
1110 * Optionally you can set the number of sweep or relaxation steps
1111 * on each grid when going up and when going down
1112 *
1113 * \param cycle_type cycle type
1114 * \param sweep_up
1115 * \param sweep_dw
1116 * \param sweep_crs speep at the coarse level
1117 *
1118 */
1119 void setPreconditionerAMG_cycleType(const std::string & cycle_type, int sweep_up = -1, int sweep_dw = -1, int sweep_crs = -1)
1120 {
1121 if (atype == HYPRE_AMG)
1122 {
1123 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_cycle_type",cycle_type.c_str()));
1124
1125 if (sweep_dw != -1)
1126 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_down",std::to_string(sweep_up).c_str()));}
1127
1128 if (sweep_up != -1)
1129 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_up",std::to_string(sweep_dw).c_str()));}
1130
1131 if (sweep_crs != -1)
1132 {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_coarse",std::to_string(sweep_crs).c_str()));}
1133 }
1134 else
1135 {
1136 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1137 }
1138 }
1139
1140 /*! \brief Set the coarsening method for the algebraic-multi-grid preconditioner
1141 *
1142 * Possible values can be
1143 * "CLJP","Ruge-Stueben","modifiedRuge-Stueben","Falgout", "PMIS", "HMIS"
1144 *
1145 * \warning in case of big problem use PMIS or HMIS otherwise the AMG can hang (take a lot of time) in setup
1146 *
1147 * \param type of the preconditioner smoothing operator
1148 *
1149 */
1150 void setPreconditionerAMG_coarsen(const std::string & type)
1151 {
1152 if (atype == HYPRE_AMG)
1153 {
1154 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_coarsen_type",type.c_str()));
1155 }
1156 else
1157 {
1158 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1159 }
1160 }
1161
1162 /*! \brief Set the interpolation method for the algebraic-multi-grid preconditioner
1163 *
1164 * Possible values can be
1165 * "classical", "direct", "multipass", "multipass-wts", "ext+i",
1166 * "ext+i-cc", "standard", "standard-wts", "FF", "FF1"
1167 *
1168 *
1169 * \param type of the interpolation scheme
1170 *
1171 */
1172 void setPreconditionerAMG_interp(const std::string & type)
1173 {
1174 if (atype == HYPRE_AMG)
1175 {
1176 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_interp_type",type.c_str()));
1177 }
1178 else
1179 {
1180 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1181 }
1182 }
1183
1184 /*! \brief Set the block coarsening norm type.
1185 *
1186 * The use of this function make sanse if you specify
1187 * the degree of freedom for each node using setBlockSize
1188 *
1189 * * 0 (default each variable treat independently)
1190 * * 1 Frobenius norm
1191 * * 2 sum of absolute values of elements in each block
1192 * * 3 largest element in each block (not absolute value)
1193 * * 4 row-sum norm
1194 * * 6 sum of all values in each block
1195 *
1196 * In case the matrix represent a system of equations in general
1197 *
1198 * \param norm type
1199 *
1200 */
1201 void setPreconditionerAMG_coarsenNodalType(int norm)
1202 {
1203 if (atype == HYPRE_AMG)
1204 {
1205 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_nodal_coarsen",std::to_string(norm).c_str()));
1206 }
1207 else
1208 {
1209 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1210 }
1211 }
1212
1213 /*! \brief Indicate the number of levels in the ILU(k) for the Euclid smoother
1214 *
1215 * \param k number of levels for the Euclid smoother
1216 *
1217 */
1218 void setPreconditionerAMG_interp_eu_level(int k)
1219 {
1220 if (atype == HYPRE_AMG)
1221 {
1222 PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_eu_level",std::to_string(k).c_str()));
1223 }
1224 else
1225 {
1226 std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
1227 }
1228 }
1229
1230
1231
1232 /*! \brief Set how many degree of freedom each node has
1233 *
1234 * In case you are solving a system of equations this function,
1235 * help in setting the degree of freedom each grid point has.
1236 * Setting this parameter to the number of variables for each
1237 * grid point it should improve the convergenve of the solvers
1238 * in particular using algebraic-multi-grid
1239 *
1240 * \param block_sz number of degree of freedom
1241 *
1242 */
1243 void setBlockSize(int block_sz)
1244 {
1245 this->block_sz = block_sz;
1246 }
1247
1248 /*! \brief Here we invert the matrix and solve the system
1249 *
1250 * \warning umfpack is not a parallel solver, this function work only with one processor
1251 *
1252 * \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet
1253 *
1254 * \tparam impl Implementation of the SparseMatrix
1255 *
1256 * \param A sparse matrix
1257 * \param b vector
1258 * \param initial_guess true if x has the initial guess
1259 *
1260 * \return the solution
1261 *
1262 */
1263 Vector<double,PETSC_BASE> solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b, bool initial_guess = false)
1264 {
1265 Mat & A_ = A.getMat();
1266 const Vec & b_ = b.getVec();
1267
1268 // We set the size of x according to the Matrix A
1269 PetscInt row;
1270 PetscInt col;
1271 PetscInt row_loc;
1272 PetscInt col_loc;
1273
1274 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1275 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
1276 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
1277
1278 Vector<double,PETSC_BASE> x(row,row_loc);
1279 Vec & x_ = x.getVec();
1280
1281 pre_solve_impl(A_,b_,x_);
1282 solve_simple(A_,b_,x_);
1283
1284 x.update();
1285
1286 return x;
1287 }
1288
1289 /*! \brief Return the KSP solver
1290 *
1291 * In case you want to do fine tuning of the KSP solver before
1292 * solve your system whith this function you can retrieve the
1293 * KSP object
1294 *
1295 * \return the Krylov solver
1296 *
1297 */
1298 KSP getKSP()
1299 {
1300 return ksp;
1301 }
1302
1303 /*! \brief It return the resiual error
1304 *
1305 * \param A Sparse matrix
1306 * \param x solution
1307 * \param b right-hand-side
1308 *
1309 * \return the solution error norms
1310 *
1311 */
1312 solError get_residual_error(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
1313 {
1314 return getSolNormError(A.getMat(),b.getVec(),x.getVec());
1315 }
1316
1317 /*! \brief It return the resiual error
1318 *
1319 * \param x solution
1320 * \param b right-hand-side
1321 *
1322 * \return the solution error norms
1323 *
1324 */
1325 solError get_residual_error(const Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
1326 {
1327 return getSolNormError(b.getVec(),x.getVec(),ksp);
1328 }
1329
1330 /*! \brief Here we invert the matrix and solve the system
1331 *
1332 * \param A sparse matrix
1333 * \param b vector
1334 * \param x solution and initial guess
1335 *
1336 * \return true if succeed
1337 *
1338 */
1339 bool solve(SparseMatrix<double,int,PETSC_BASE> & A, Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
1340 {
1341 Mat & A_ = A.getMat();
1342 const Vec & b_ = b.getVec();
1343 Vec & x_ = x.getVec();
1344
1345 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_TRUE));
1346
1347 pre_solve_impl(A_,b_,x_);
1348 solve_simple(A_,b_,x_);
1349 x.update();
1350
1351 return true;
1352 }
1353
1354 /*! \brief Here we invert the matrix and solve the system
1355 *
1356 * \param b vector
1357 *
1358 * \return true if succeed
1359 *
1360 */
1361 Vector<double,PETSC_BASE> solve(const Vector<double,PETSC_BASE> & b)
1362 {
1363 const Vec & b_ = b.getVec();
1364
1365 // We set the size of x according to the Matrix A
1366 PetscInt row;
1367 PetscInt row_loc;
1368
1369 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1370 PETSC_SAFE_CALL(VecGetSize(b_,&row));
1371 PETSC_SAFE_CALL(VecGetLocalSize(b_,&row_loc));
1372
1373 Vector<double,PETSC_BASE> x(row,row_loc);
1374 Vec & x_ = x.getVec();
1375
1376 solve_simple(b_,x_);
1377 x.update();
1378
1379 return x;
1380 }
1381
1382 /*! \brief Here we invert the matrix and solve the system
1383 *
1384 * In this call we are interested in solving the system
1385 * with multiple right-hand-side and the same Matrix.
1386 * We do not set the Matrix again and this give us the
1387 * possibility to-skip the preconditioning setting that
1388 * in some case like Algebraic-multi-grid can be expensive
1389 *
1390 * \param b vector
1391 * \param x solution and initial guess
1392 *
1393 * \return true if succeed
1394 *
1395 */
1396 bool solve(Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
1397 {
1398 const Vec & b_ = b.getVec();
1399 Vec & x_ = x.getVec();
1400
1401 solve_simple(b_,x_);
1402 x.update();
1403
1404 return true;
1405 }
1406
1407 /*! \brief this function give you the possibility to set PETSC options
1408 *
1409 * this function call PetscOptionsSetValue
1410 *
1411 * \param name the name of the option
1412 * \param value the value of the option
1413 *
1414 */
1415 void setPetscOption(const char * name, const char * value)
1416 {
1417 PetscOptionsSetValue(NULL,name,value);
1418 }
1419
1420 /*! \brief Try to solve the system using all the solvers and generate a report
1421 *
1422 * In this mode the system will try different Solvers, Preconditioner and
1423 * combination of solvers in order to find the best solver in speed, and
1424 * precision. As output it will produce a performance report
1425 *
1426 * \param A Matrix to invert
1427 * \param b right hand side
1428 * \return the solution
1429 *
1430 */
1431 Vector<double,PETSC_BASE> try_solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b)
1432 {
1433 Mat & A_ = A.getMat();
1434 const Vec & b_ = b.getVec();
1435
1436 // We set the size of x according to the Matrix A
1437 PetscInt row;
1438 PetscInt col;
1439 PetscInt row_loc;
1440 PetscInt col_loc;
1441
1442 PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
1443 PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
1444 PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
1445
1446 Vector<double,PETSC_BASE> x(row,row_loc);
1447 Vec & x_ = x.getVec();
1448
1449 pre_solve_impl(A_,b_,x_);
1450 try_solve_simple(A_,b_,x_);
1451
1452 x.update();
1453
1454 return x;
1455 }
1456};
1457
1458#endif
1459
1460#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ */
1461