petsc_solver.hpp

/*
 * petsc_solver.hpp
 *
 *  Created on: Apr 26, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_
#define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_

#include "config.h"

#ifdef HAVE_PETSC

#include "Vector/Vector.hpp"
#include <petscksp.h>
#include <petsctime.h>
#include "Plot/GoogleChart.hpp"
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include <sstream>
#include <iomanip>

template <typename T>
std::string to_string_with_precision(const T a_value, const int n = 6)
{
    std::ostringstream out;
    out << std::setprecision(n) << a_value;
    return out.str();
}

#define UMFPACK_NONE 0
#define REL_DOWN 1
#define REL_UP 2
#define REL_ALL 3

#define PCHYPRE_BOOMERAMG "petsc_boomeramg"

enum AMG_type
{
    NONE_AMG,
    HYPRE_AMG,
    PETSC_AMG,
    TRILINOS_ML
};


template<typename T>
class petsc_solver
{
public:

    template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b)
    {
        std::cerr << "Error Petsc only suppor double precision" << "/n";
    }
};

#define SOLVER_NOOPTION 0
#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
#define SOLVER_PRINT_DETERMINANT 2

struct solError
{
    PetscReal err_inf;

    PetscReal err_norm;

    PetscInt its;
};


template<>
class petsc_solver<double>
{
    struct itError
    {
        PetscInt it;

        PetscReal err_norm;
    };

    struct solv_bench_info
    {
        std::string method;

        std::string smethod;

        double time;

        solError err;

        openfpm::vector<itError> res;
    };

    bool is_preconditioner_set = false;

    PetscInt maxits;

    KSP ksp;

    size_t tmp;

    openfpm::vector<std::string> solvs;


    openfpm::vector<solv_bench_info> bench;

    AMG_type atype = NONE_AMG;

    int block_sz = 0;

    static double calculate_it(double t, solv_bench_info & slv)
    {
        double s_int = slv.time / slv.res.size();

        // Calculate the discrete point in time
        size_t pt = std::floor(t / s_int);

        if (pt < slv.res.size())
            return slv.res.get(pt).err_norm;

        return slv.res.last().err_norm;
    }

    void write_bench_report()
    {
        if (create_vcluster().getProcessUnitID() != 0)
            return;

        openfpm::vector<std::string> x;
        openfpm::vector<openfpm::vector<double>> y;
        openfpm::vector<std::string> yn;
        openfpm::vector<double> xd;

        for (size_t i = 0 ; i < bench.size() ; i++)
            x.add(bench.get(i).smethod);

        // Each colum can have multiple data set (in this case 4 dataset)
        // Each dataset can have a name
        yn.add("Norm infinity");
        yn.add("Norm average");
        yn.add("Number of iterations");

        // Each colums can have multiple data-set

        for (size_t i = 0 ; i < bench.size() ; i++)
            y.add({bench.get(i).err.err_inf,bench.get(i).err.err_norm,(double)bench.get(i).err.its});

        // Google charts options
        GCoptions options;

        options.title = std::string("Residual error");
        options.yAxis = std::string("Error");
        options.xAxis = std::string("Method");
        options.stype = std::string("bars");
        options.more = std::string("vAxis: {scaleType: 'log'}");

        GoogleChart cg;

        std::stringstream ss;

        cg.addHTML("<h2>Robustness</h2>");
        cg.AddHistGraph(x,y,yn,options);
        cg.addHTML("<h2>Speed</h2>");

        y.clear();
        yn.clear();

        yn.add("Time in millseconds");
        options.title = std::string("Time");

        for (size_t i = 0 ; i < bench.size() ; i++)
            y.add({bench.get(i).time});

        cg.AddHistGraph(x,y,yn,options);

        // Convergence in time

        x.clear();
        y.clear();
        yn.clear();
        xd.clear();

        // Get the maximum in time across all the solvers

        double max_time = 0.0;

        for (size_t i = 0 ; i < bench.size() ; i++)
        {
            if (bench.get(i).time > max_time)   {max_time = bench.get(i).time;}
            yn.add(bench.get(i).smethod);
        }

        size_t n_int = maxits;

        // calculate dt

        double dt = max_time / n_int;

        // Add

        // For each solver we have a convergence plot

        for (double t = dt ; t <= max_time + 0.05 * max_time ; t += dt)
        {
            y.add();
            xd.add(t);
            for (size_t j = 0 ; j < bench.size() ; j++)
                y.last().add(calculate_it(t,bench.get(j)));
        }

        std::stringstream ss2;
        ss2 << "<h2>Convergence with time</h2><br>" << std::endl;

        options.title = std::string("Residual norm infinity");
        options.yAxis = std::string("residual");
        options.xAxis = std::string("Time in milliseconds");
        options.lineWidth = 1.0;
        options.more = std::string("vAxis: {scaleType: 'log'},hAxis: {scaleType: 'log'}");

        cg.addHTML(ss2.str());
        cg.AddLinesGraph(xd,y,yn,options);

        ss << "<h2>Solvers Tested</h2><br><br>" << std::endl;
        for (size_t i = 0 ; i < bench.size() ; i++)
            ss << bench.get(i).smethod << " = " << bench.get(i).method << "<br>" << std::endl;

        cg.addHTML(ss.str());

        cg.write("gc_solver.html");
    }

    void pre_solve_impl(const Mat & A_, const Vec & b_, Vec & x_)
    {
        PETSC_SAFE_CALL(KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED));

        if (atype == HYPRE_AMG)
        {
            PC pc;

            // We set the pre-conditioner
            PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));

            PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
            PCFactorSetShiftAmount(pc, PETSC_DECIDE);
        }
    }

    void print_progress_bar()
    {
        auto & v_cl = create_vcluster();

        if (v_cl.getProcessUnitID() == 0)
            std::cout << "-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
    }


    static PetscErrorCode monitor_progress_residual(KSP ksp,PetscInt it,PetscReal res,void* data)
    {
        petsc_solver<double> * pts = (petsc_solver *)data;

        pts->progress(it);

        itError erri;
        erri.it = it;

        Mat A;
        Vec Br,v,w;

        PETSC_SAFE_CALL(KSPGetOperators(ksp,&A,NULL));
        PETSC_SAFE_CALL(MatCreateVecs(A,&w,&v));
        PETSC_SAFE_CALL(KSPBuildResidual(ksp,v,w,&Br));
        PETSC_SAFE_CALL(KSPBuildResidual(ksp,v,w,&Br));
        PETSC_SAFE_CALL(VecNorm(Br,NORM_INFINITY,&erri.err_norm));
        itError err_fill;

        size_t old_size = pts->bench.last().res.size();
        pts->bench.last().res.resize(it+1);

        if (old_size > 0)
            err_fill = pts->bench.last().res.get(old_size-1);
        else
            err_fill = erri;

        for (long int i = old_size ; i < (long int)it ; i++)
            pts->bench.last().res.get(i) = err_fill;

        // Add the error per iteration
        pts->bench.last().res.get(it) = erri;

        return 0;
    }

    void progress(PetscInt it)
    {
        PetscInt n_max_it;

        PETSC_SAFE_CALL(KSPGetTolerances(ksp,NULL,NULL,NULL,&n_max_it));

        auto & v_cl = create_vcluster();

        if (std::floor(it * 100.0 / n_max_it) > tmp)
        {
            size_t diff = it * 100.0 / n_max_it - tmp;
            tmp = it * 100.0 / n_max_it;
            if (v_cl.getProcessUnitID() == 0)
            {
                for (size_t k = 0 ; k <  diff ; k++)    {std::cout << "*";}
                std::cout << std::flush;
            }
        }
    }

    void print_stat(solError & err)
    {
        auto & v_cl = create_vcluster();

        if (v_cl.getProcessUnitID() == 0)
        {
            KSPType typ;
            KSPGetType(ksp,&typ);

            std::cout << "Method: " << typ << " " << " pre-conditoner: " << PCJACOBI << "  iterations: " << err.its << std::endl;
            std::cout << "Norm of error: " << err.err_norm << "   Norm infinity: " << err.err_inf << std::endl;
        }
    }

    std::string to_string_method(const std::string & solv)
    {
        if (solv == std::string(KSPBCGS))
        {
            return std::string("BiCGStab (Stabilized version of BiConjugate Gradient Squared)");
        }
        else if (solv == std::string(KSPIBCGS))
        {
            return std::string("IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared)");
        }
        else if (solv == std::string(KSPFBCGS))
        {
            return std::string("FBiCGStab (Flexible Stabilized version of BiConjugate Gradient Squared)");
        }
        else if (solv == std::string(KSPFBCGSR))
        {
            return std::string("FBiCGStab-R (Modified Flexible Stabilized version of BiConjugate Gradient Squared)");
        }
        else if (solv == std::string(KSPBCGSL))
        {
            return std::string("BiCGStab(L) (Stabilized version of BiConjugate Gradient Squared with L search directions)");
        }
        else if (solv == std::string(KSPGMRES))
        {
            return std::string("GMRES(M) (Generalized Minimal Residual method with restart)");
        }
        else if (solv == std::string(KSPFGMRES))
        {
            return std::string("FGMRES(M) (Flexible Generalized Minimal Residual with restart)");
        }
        else if (solv == std::string(KSPLGMRES))
        {
            return std::string("LGMRES(M) (Augment Generalized Minimal Residual method with restart)");
        }
        else if (solv == std::string(KSPPGMRES))
        {
            return std::string("PGMRES(M) (Pipelined Generalized Minimal Residual method)");
        }
        else if (solv == std::string(KSPGCR))
        {
            return std::string("GCR (Generalized Conjugate Residual)");
        }

        return std::string("UNKNOWN");
    }

    void new_bench(const std::string & str)
    {
        // Add a new benchmark result
        bench.add();
        bench.last().method = to_string_method(str);
        bench.last().smethod = std::string(str);
    }

    void copy_if_better(double res, Vec & sol, double & best_res, Vec & best_sol)
    {
        if (res < best_res)
        {
            VecCopy(sol,best_sol);
            best_res = res;
        }
    }

    void try_solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
    {
        Vec best_sol;
        PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));

        // Best residual
        double best_res = std::numeric_limits<double>::max();

        // Create a new VCluster
        auto & v_cl = create_vcluster();

        destroyKSP();

        // for each solver
        for (size_t i = 0 ; i < solvs.size() ; i++)
        {
            initKSPForTest();

            // Here we solve without preconditioner
            PC pc;

            // We set the pre-conditioner to none
            PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCNONE));

            setSolver(solvs.get(i).c_str());

            // Setup for BCGSL, GMRES
            if (solvs.get(i) == std::string(KSPBCGSL))
            {
                // we try from 2 to 6 as search direction
                for (size_t j = 2 ; j < 6 ; j++)
                {
                    new_bench(solvs.get(i));

                    if (v_cl.getProcessUnitID() == 0)
                        std::cout << "L = " << j << std::endl;
                    bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
                    searchDirections(j);
                    bench_solve_simple(A_,b_,x_,bench.last());

                    copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
                }
            }
            else if (solvs.get(i) == std::string(KSPGMRES) ||
                     solvs.get(i) == std::string(std::string(KSPFGMRES)) ||
                     solvs.get(i) == std::string(std::string(KSPLGMRES)) )
            {
                // we try from 2 to 6 as search direction
                for (size_t j = 50 ; j < 300 ; j += 50)
                {
                    // Add a new benchmark result
                    new_bench(solvs.get(i));

                    if (v_cl.getProcessUnitID() == 0)
                        std::cout << "Restart = " << j << std::endl;
                    bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
                    setRestart(j);
                    bench_solve_simple(A_,b_,x_,bench.last());

                    copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
                }
            }
            else
            {
                // Add a new benchmark result
                new_bench(solvs.get(i));

                bench_solve_simple(A_,b_,x_,bench.last());

                copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
            }

            destroyKSP();
        }

        write_bench_report();

        // Copy the best solution to x
        PETSC_SAFE_CALL(VecCopy(best_sol,x_));
        PETSC_SAFE_CALL(VecDestroy(&best_sol));
    }


    void bench_solve_simple(const Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
    {
        // timer for benchmark
        timer t;
        t.start();
        // Enable progress monitor
        tmp = 0;
        PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
        PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress_residual,this,NULL));
        print_progress_bar();
        solve_simple(A_,b_,x_);

        // New line
        if (create_vcluster().getProcessUnitID() == 0)
            std::cout << std::endl;

        t.stop();
        bench.time = t.getwct() * 1000.0;

        // calculate error statistic about the solution and print
        solError err = statSolutionError(A_,b_,x_,ksp);
        print_stat(err);

        bench.err = err;

        // New line
        if (create_vcluster().getProcessUnitID() == 0)
            std::cout << std::endl;
    }

    void solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
    {
//      PETSC_SAFE_CALL(KSPSetType(ksp,s_type));

        // We set the size of x according to the Matrix A
        PetscInt row;
        PetscInt col;
        PetscInt row_loc;
        PetscInt col_loc;

        PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
        PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

        // We set the Matrix operators
        PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));

        // if we are on on best solve set-up a monitor function

        PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
        PETSC_SAFE_CALL(KSPSetUp(ksp));

        // Solve the system
        PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
    }

    void solve_simple(const Vec & b_, Vec & x_)
    {
        // Solve the system
        PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
    }

    static solError statSolutionError(const Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
    {
        solError err;

        err = getSolNormError(A_,b_,x_);

        PetscInt its;
        PETSC_SAFE_CALL(KSPGetIterationNumber(ksp,&its));
        err.its = its;

        return err;
    }


    void initKSP()
    {
        PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
    }

    void initKSPForTest()
    {
        PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));

        setMaxIter(maxits);
    }

    void destroyKSP()
    {
        KSPDestroy(&ksp);
    }

    static solError getSolNormError(const Vec & b_, const Vec & x_,KSP ksp)
    {
        Mat A_;
        Mat P_;
        KSPGetOperators(ksp,&A_,&P_);

        return getSolNormError(A_,b_,x_);
    }

    static solError getSolNormError(const Mat & A_, const Vec & b_, const Vec & x_)
    {
        PetscScalar neg_one = -1.0;

        // We set the size of x according to the Matrix A
        PetscInt row;
        PetscInt col;
        PetscInt row_loc;
        PetscInt col_loc;

        PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
        PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

        /*
              Here we calculate the residual error
        */
        PetscReal norm;
        PetscReal norm_inf;

        // Get a vector r for the residual
        Vector<double,PETSC_BASE> r(row,row_loc);
        Vec & r_ = r.getVec();

        PETSC_SAFE_CALL(MatMult(A_,x_,r_));
        PETSC_SAFE_CALL(VecAXPY(r_,neg_one,b_));

        PETSC_SAFE_CALL(VecNorm(r_,NORM_1,&norm));
        PETSC_SAFE_CALL(VecNorm(r_,NORM_INFINITY,&norm_inf));

        solError err;
        err.err_norm = norm / col;
        err.err_inf = norm_inf;
        err.its = 0;

        return err;
    }

public:

    typedef Vector<double,PETSC_BASE> return_type;

    ~petsc_solver()
    {
        PETSC_SAFE_CALL(KSPDestroy(&ksp));
    }

    petsc_solver()
    :maxits(300),tmp(0)
    {
        initKSP();

        // Add the solvers

        solvs.add(std::string(KSPBCGS));
//      solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
//      solvs.add(std::string(KSPFBCGS));
//      solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
        solvs.add(std::string(KSPBCGSL));
        solvs.add(std::string(KSPGMRES));
        solvs.add(std::string(KSPFGMRES));
        solvs.add(std::string(KSPLGMRES));
//      solvs.add(std::string(KSPPGMRES)); <--- Take forever
//      solvs.add(std::string(KSPGCR));

        setSolver(KSPGMRES);
    }

    void addTestSolver(std::string & solver)
    {
        solvs.add(solver);
    }

    void removeTestSolver(const std::string & solver)
    {
        // search for the test solver and remove it

        for (size_t i = 0 ; i < solvs.size() ; i++)
        {
            if (solvs.get(i) == solver)
                return;
        }
    }


    void log_monitor()
    {
        PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-ksp_monitor",0));
        PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_print_statistics","2"));
    }

    void setSolver(KSPType type)
    {
        PetscOptionsSetValue(NULL,"-ksp_type",type);
    }

    void setRelTol(PetscReal rtol_)
    {
        PetscReal rtol;
        PetscReal abstol;
        PetscReal dtol;
        PetscInt maxits;

        PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
        PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol_,abstol,dtol,maxits));
    }

    void setAbsTol(PetscReal abstol_)
    {
        PetscReal rtol;
        PetscReal abstol;
        PetscReal dtol;
        PetscInt maxits;

        PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
        PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol_,dtol,maxits));
    }

    void setDivTol(PetscReal dtol_)
    {
        PetscReal rtol;
        PetscReal abstol;
        PetscReal dtol;
        PetscInt maxits;

        PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
        PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol_,maxits));
    }

    void setMaxIter(PetscInt n)
    {
        PetscReal rtol;
        PetscReal abstol;
        PetscReal dtol;
        PetscInt maxits;

        PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
        PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,n));

        this->maxits = n;
    }

    void searchDirections(PetscInt l)
    {
        PetscOptionsSetValue(NULL,"-ksp_bcgsl_ell",std::to_string(l).c_str());
    }

    void setRestart(PetscInt n)
    {
        PetscOptionsSetValue(NULL,"-ksp_gmres_restart",std::to_string(n).c_str());
    }

    void setPreconditioner(PCType type)
    {
        is_preconditioner_set = true;

        if (std::string(type) == PCHYPRE_BOOMERAMG)
        {
            PC pc;

            // We set the pre-conditioner
            PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));

            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCHYPRE));

            PETSC_SAFE_CALL(PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO));
            PETSC_SAFE_CALL(PCFactorSetShiftAmount(pc, PETSC_DECIDE));
            PETSC_SAFE_CALL(PCHYPRESetType(pc, "boomeramg"));
            atype = HYPRE_AMG;
        }
        else
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",PCHYPRE));
        }
    }

    void setPreconditionerAMG_nl(int nl)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_max_levels",std::to_string(nl).c_str()));
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }

    void setPreconditionerAMG_maxit(int nit)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_max_iter",std::to_string(nit).c_str()));
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }


    void setPreconditionerAMG_relax(const std::string & type, int k = REL_ALL)
    {
        if (atype == HYPRE_AMG)
        {
            if (k == REL_ALL)
            {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_relax_type_all",type.c_str()));}
            else if (k == REL_UP)
            {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_relax_type_up",type.c_str()));}
            else if (k == REL_DOWN)
            {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_relax_type_down",type.c_str()));}
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }

    void setPreconditionerAMG_cycleType(const std::string & cycle_type, int sweep_up = -1, int sweep_dw = -1, int sweep_crs = -1)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_cycle_type",cycle_type.c_str()));

            if (sweep_dw != -1)
            {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_down",std::to_string(sweep_up).c_str()));}

            if (sweep_up != -1)
            {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_up",std::to_string(sweep_dw).c_str()));}

            if (sweep_crs != -1)
            {PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_grid_sweeps_coarse",std::to_string(sweep_crs).c_str()));}
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }

    void setPreconditionerAMG_coarsen(const std::string & type)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_coarsen_type",type.c_str()));
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }

    void setPreconditionerAMG_interp(const std::string & type)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_interp_type",type.c_str()));
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }

    void setPreconditionerAMG_coarsenNodalType(int norm)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_nodal_coarsen",std::to_string(norm).c_str()));
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }

    void setPreconditionerAMG_interp_eu_level(int k)
    {
        if (atype == HYPRE_AMG)
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_hypre_boomeramg_eu_level",std::to_string(k).c_str()));
        }
        else
        {
            std::cout << __FILE__ << ":" << __LINE__ << "Warning: the selected preconditioner does not support this option" << std::endl;
        }
    }



    void setBlockSize(int block_sz)
    {
        this->block_sz = block_sz;
    }

    Vector<double,PETSC_BASE> solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b, bool initial_guess = false)
    {
        Mat & A_ = A.getMat();
        const Vec & b_ = b.getVec();

        // We set the size of x according to the Matrix A
        PetscInt row;
        PetscInt col;
        PetscInt row_loc;
        PetscInt col_loc;

        PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
        PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
        PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

        Vector<double,PETSC_BASE> x(row,row_loc);
        Vec & x_ = x.getVec();

        pre_solve_impl(A_,b_,x_);
        solve_simple(A_,b_,x_);

        x.update();

        return x;
    }

    KSP getKSP()
    {
        return ksp;
    }

    solError get_residual_error(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
    {
        return getSolNormError(A.getMat(),b.getVec(),x.getVec());
    }

    solError get_residual_error(const Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
    {
        return getSolNormError(b.getVec(),x.getVec(),ksp);
    }

    bool solve(SparseMatrix<double,int,PETSC_BASE> & A, Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
    {
        Mat & A_ = A.getMat();
        const Vec & b_ = b.getVec();
        Vec & x_ = x.getVec();

        PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_TRUE));

        pre_solve_impl(A_,b_,x_);
        solve_simple(A_,b_,x_);
        x.update();

        return true;
    }

    Vector<double,PETSC_BASE> solve(const Vector<double,PETSC_BASE> & b)
    {
        const Vec & b_ = b.getVec();

        // We set the size of x according to the Matrix A
        PetscInt row;
        PetscInt row_loc;

        PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
        PETSC_SAFE_CALL(VecGetSize(b_,&row));
        PETSC_SAFE_CALL(VecGetLocalSize(b_,&row_loc));

        Vector<double,PETSC_BASE> x(row,row_loc);
        Vec & x_ = x.getVec();

        solve_simple(b_,x_);
        x.update();

        return x;
    }

    bool solve(Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
    {
        const Vec & b_ = b.getVec();
        Vec & x_ = x.getVec();

        solve_simple(b_,x_);
        x.update();

        return true;
    }

    void setPetscOption(const char * name, const char * value)
    {
        PetscOptionsSetValue(NULL,name,value);
    }

    Vector<double,PETSC_BASE> try_solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b)
    {
        Mat & A_ = A.getMat();
        const Vec & b_ = b.getVec();

        // We set the size of x according to the Matrix A
        PetscInt row;
        PetscInt col;
        PetscInt row_loc;
        PetscInt col_loc;

        PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
        PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
        PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

        Vector<double,PETSC_BASE> x(row,row_loc);
        Vec & x_ = x.getVec();

        pre_solve_impl(A_,b_,x_);
        try_solve_simple(A_,b_,x_);

        x.update();

        return x;
    }
};

#endif

#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ */