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 print_preconditioner()
    {
        auto & v_cl = create_vcluster();
 
        if (v_cl.getProcessUnitID() == 0)
        {
            PC pc;
            PCType type_pc;
 
            // We set the pre-conditioner
            PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
            PETSC_SAFE_CALL(PCGetType(pc,&type_pc));
 
            std::cout << "The precoditioner is: " << type_pc << std::endl;
        }
    }
 
    void print_ksptype()
    {
        auto & v_cl = create_vcluster();
 
        if (v_cl.getProcessUnitID() == 0)
        {
            KSPType type_ksp;
 
            // We get the KSP_type
            PETSC_SAFE_CALL(KSPGetType(ksp,&type_ksp));
 
            std::cout << "The ksp_type is: " << type_ksp << std::endl;
        }
    }
 
    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));
#ifdef PETSC_HAVE_HYPRE
            PETSC_SAFE_CALL(PCHYPRESetType(pc, "boomeramg"));
#endif
            atype = HYPRE_AMG;
        }
        else
        {
            PETSC_SAFE_CALL(PetscOptionsSetValue(NULL,"-pc_type",type));
        }
    }
 
    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;
    }
 
    Vector<double,PETSC_BASE> 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();
 
       /* // 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_TRUE));
 
        /*PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
        PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));*/
 
        pre_solve_impl(A_,b_,x_);
        solve_simple(A_,b_,x_);
 
        x.update();
 
        return x;
 
        /*pre_solve_impl(A_,b_,x_);
        solve_simple(A_,b_,x_);
        x.update();
 
        return true;*/
    }
 
    Vector<double,PETSC_BASE> solve_successive(const Vector<double,PETSC_BASE> & b,bool initial_guess = false)
    {
        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();
        PETSC_SAFE_CALL(KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED));
        solve_simple(b_,x_);
 
        x.update();
 
        return x;
 
        /*pre_solve_impl(A_,b_,x_);
        solve_simple(A_,b_,x_);
        x.update();
 
        return true;*/
    }
 
    Vector<double,PETSC_BASE> solve_successive(Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
    {
        const Vec & b_ = b.getVec();
        Vec & x_ = x.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_TRUE));
 
        /*PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
        PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));*/
 
        PETSC_SAFE_CALL(KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED));
 
        solve_simple(b_,x_);
 
        x.update();
 
        return x;
 
        /*pre_solve_impl(A_,b_,x_);
        solve_simple(A_,b_,x_);
        x.update();
 
        return true;*/
    }
 
    Vector<double,PETSC_BASE> with_nullspace_solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b, bool initial_guess = false,bool symmetric = 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();
 
        PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
        PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
        PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
 
        Mat      F, work, V;
        PetscInt N, rows;
 
        /* Determine factorability */
        PETSC_SAFE_CALL(MatGetFactor(A_, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
        PETSC_SAFE_CALL(MatGetLocalSize(A_, &rows, NULL));
 
        /* Set MUMPS options, see MUMPS documentation for more information */
        PETSC_SAFE_CALL(MatMumpsSetIcntl(F, 24, 1));
        PETSC_SAFE_CALL(MatMumpsSetIcntl(F, 25, 1));
 
        /* Perform factorization */
        PETSC_SAFE_CALL(MatLUFactorSymbolic(F, A_, NULL, NULL, NULL));
        PETSC_SAFE_CALL(MatLUFactorNumeric(F, A_, NULL));
 
        /* This is the dimension of the null space */
        PETSC_SAFE_CALL(MatMumpsGetInfog(F, 28, &N));
        /* This will contain the null space in the columns */
        PETSC_SAFE_CALL(MatCreateDense(PETSC_COMM_WORLD, rows, N, PETSC_DETERMINE, PETSC_DETERMINE, NULL, &V));
        PETSC_SAFE_CALL(MatDuplicate(V, MAT_DO_NOT_COPY_VALUES, &work));
        PETSC_SAFE_CALL(MatMatSolve(F, work, V));
 
        std::cout<<"Dimension:" << N;
        Vec nvec[N];
        for(int i=0;i<N;i++)
        {
            PETSC_SAFE_CALL(MatGetColumnVector(V,nvec[i],i));
        }
        MatNullSpace  nullspace;
 
        PETSC_SAFE_CALL(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,N,nvec,&nullspace));
        PETSC_SAFE_CALL(MatSetTransposeNullSpace(A_,nullspace));
        PETSC_SAFE_CALL(MatSetNullSpace(A_,nullspace));
        PETSC_SAFE_CALL(MatNullSpaceDestroy(&nullspace));
 
        PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
        PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
        PETSC_SAFE_CALL(KSPSolve(ksp,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);
    }
 
    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_ */