petsc_solver_AMG_report.hpp

/*
 * petsc_solver_report.hpp
 *
 *  Created on: Jun 22, 2017
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_AMG_REPORT_HPP_
#define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_AMG_REPORT_HPP_

#include "config.h"

#ifdef HAVE_PETSC

#include <fstream>
#include "Solvers/petsc_solver.hpp"
#include "Vector/map_vector.hpp"

struct AMG_time_err_coars
{
    double time_setup;

    double time_solve;

    double residual;

    double score;
};

class petsc_AMG_report
{
    size_t num_sweep_test = 10;

    double target_accuracy = 1;

    openfpm::vector<AMG_time_err_coars> perf_amg_coars;

    openfpm::vector<AMG_time_err_coars> perf_amg_sweep_sym;

    openfpm::vector<AMG_time_err_coars> perf_amg_sweep_asym;

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

    openfpm::vector<size_t> v_cycle_tested_sym;

    openfpm::vector<std::pair<size_t,size_t>> v_cycle_tested_asym;

    openfpm::vector<size_t> asym_tests;

    void benchmark(SparseMatrix<double,int,PETSC_BASE> & A,
                   const Vector<double,PETSC_BASE> & b,
                   petsc_solver<double> & solver,
                   openfpm::vector<AMG_time_err_coars> & perf_amg)
    {
        timer tm_solve;
        tm_solve.start();
        auto x_ = solver.solve(A,b);
        tm_solve.stop();

        solError serr = solver.get_residual_error(A,x_,b);

        // In order to measure the time to solve we solve again the system

        timer tm_solve2;
        tm_solve2.start();
        solver.solve(x_,b);
        tm_solve2.stop();

        double time1 = tm_solve.getwct();
        double time2 = tm_solve2.getwct();

        Vcluster & v_cl = create_vcluster();
        v_cl.max(time1);
        v_cl.max(time2);
        v_cl.execute();

        // Save the result
        AMG_time_err_coars tmp;
        tmp.time_setup = time1 - time2;
        tmp.time_solve = time2;
        tmp.residual = serr.err_inf;

        perf_amg.add(tmp);

        if (v_cl.getProcessUnitID() == 0)
            std::cout << "Time1: " << time1 << "    Time2: " << time2 << std::endl;
    }

    void test_coarsener(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b)
    {
        Vcluster & v_cl = create_vcluster();

        petsc_solver<double> pts;

        petsc_solver<double> solver;
        solver.setSolver(KSPRICHARDSON);
        solver.setAbsTol(0.01);
        solver.setMaxIter(1);


        solver.setPreconditioner(PCHYPRE_BOOMERAMG);
        solver.setPreconditionerAMG_nl(6);
        solver.setPreconditionerAMG_maxit(3);
        solver.setPreconditionerAMG_relax("SOR/Jacobi");
        solver.setPreconditionerAMG_cycleType("V",6,6);
        solver.setPreconditionerAMG_coarsenNodalType(0);
        solver.log_monitor();


        // Reset the Matrix A
        A.getMatrixTriplets();
        if (v_cl.getProcessUnitID() == 0)   {std::cout << "Benchmarking HMIS coarsener" << std::endl;}
        solver.setPreconditionerAMG_coarsen("HMIS");
        benchmark(A,b,solver,perf_amg_coars);
        method.add("HMIS");

        // Reset the Matrix A
        A.getMatrixTriplets();
        if (v_cl.getProcessUnitID() == 0)   {std::cout << "Benchmarking PMIS coarsener" << std::endl;}
        solver.setPreconditionerAMG_coarsen("PMIS");
        benchmark(A,b,solver,perf_amg_coars);
        method.add("PMIS");

        // Reset the Matrix A
        A.getMatrixTriplets();
        if (v_cl.getProcessUnitID() == 0)   {std::cout << "Benchmarking Falgout coarsener" << std::endl;}
        solver.setPreconditionerAMG_coarsen("Falgout");
        benchmark(A,b,solver,perf_amg_coars);
        method.add("Falgout");

        // Reset the Matrix A
        A.getMatrixTriplets();
        if (v_cl.getProcessUnitID() == 0)   {std::cout << "Benchmarking modifiedRuge-Stueben coarsener" << std::endl;}
        solver.setPreconditionerAMG_coarsen("modifiedRuge-Stueben");
        benchmark(A,b,solver,perf_amg_coars);
        method.add("modifiedRuge-Stueben");

        // Reset the Matrix A
        A.getMatrixTriplets();
        if (v_cl.getProcessUnitID() == 0)   {std::cout << "Benchmarking Ruge-Stueben coarsener" << std::endl;}
        solver.setPreconditionerAMG_coarsen("Ruge-Stueben");
        benchmark(A,b,solver,perf_amg_coars);
        method.add("Ruge-Stueben");

        // Reset the Matrix A
        A.getMatrixTriplets();
        if (v_cl.getProcessUnitID() == 0)   {std::cout << "Benchmarking CLJP coarsener" << std::endl;}
        solver.setPreconditionerAMG_coarsen("CLJP");
        benchmark(A,b,solver,perf_amg_coars);
        method.add("CLJP");

    }

    double score_solver(double t_solve, double t_m)
    {
        return 1.0/t_solve * std::min(1.0,target_accuracy/t_m);
    }

    void write_report_coars(GoogleChart & cg)
    {
        std::cout << "Composing coarsening report" << std::endl;

        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 < perf_amg_coars.size() ; i++)
            x.add(method.get(i));

        // Each colum can have multiple data set (in this case 4 dataset)
        // Each dataset can have a name
        yn.add("Norm infinity");
        yn.add("time to setup");
        yn.add("time to solve");

        // Each colums can have multiple data-set
        for (size_t i = 0 ; i < perf_amg_coars.size() ; i++)
            y.add({perf_amg_coars.get(i).residual,perf_amg_coars.get(i).time_setup,perf_amg_coars.get(i).time_solve});

        // Google charts options
        GCoptions options;

        options.title = std::string("Overview of different coarsener on a V cycle with 6 relaxation SOR/Jacobi sweeps for each level up and down");
        options.yAxis = std::string("Error/time");
        options.xAxis = std::string("Method");
        options.stype = std::string("bars");
        options.more = std::string("vAxis: {scaleType: 'log'}");

        std::stringstream ss;

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

    void write_report_cycle(GoogleChart & cg, openfpm::vector<std::string> & coars)
    {
        cg.addHTML("<h2>V cycle sweep</h2>");
        for(size_t k = 0 ; k < perf_amg_sweep_sym.size() / num_sweep_test ; k++)
        {
            openfpm::vector<double> x;
            openfpm::vector<openfpm::vector<double>> y;
            openfpm::vector<std::string> yn;
            openfpm::vector<double> xd;

            x.clear();
            for (size_t i = 0 ; i < num_sweep_test ; i++)
                x.add(v_cycle_tested_sym.get(i));

            // Each colum can have multiple data set (in this case 4 dataset)
            // Each dataset can have a name
            yn.add("error");
            yn.add("time");
            yn.add("score");

            // Each colums can have multiple data-set
            for (size_t i = 0 ; i < num_sweep_test ; i++)
            {
                double score = score_solver(perf_amg_sweep_sym.get(k*num_sweep_test+i).time_solve,perf_amg_sweep_sym.get(k*num_sweep_test+i).residual);
                y.add({perf_amg_sweep_sym.get(k*num_sweep_test+i).residual,perf_amg_sweep_sym.get(k*num_sweep_test+i).time_solve,score});
            }

            // Google charts options
            GCoptions options;

            options.title = std::string("Coarsener: " + coars.get(k)) + std::string(": Overview of V cycle with n sweeps symmetrically up and down");
            options.yAxis = std::string("error/time/score");
            options.xAxis = std::string("sweeps");
            options.stype = std::string("bars");
            options.more = std::string("vAxis: {scaleType: 'log'}");

            std::stringstream ss;

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

    void write_report_cycle_asym(GoogleChart & cg, openfpm::vector<std::string> & coars)
    {
        cg.addHTML("<h2>V cycle sweep asymmetric test</h2>");
        for(size_t k = 0 ; k < coars.size() ; k++)
        {
            openfpm::vector<std::string> x;
            openfpm::vector<openfpm::vector<double>> y;
            openfpm::vector<std::string> yn;
            openfpm::vector<double> xd;

            size_t num_sweep_test = asym_tests.get(k+1) - asym_tests.get(k);
            size_t sweep_start = asym_tests.get(k);

            x.clear();
            for (size_t i = 0 ; i < num_sweep_test ; i++)
            {x.add(std::string("(" + std::to_string(v_cycle_tested_asym.get(i+sweep_start).first) + "," + std::to_string(v_cycle_tested_asym.get(i+sweep_start).second) + ")"));}

            // Each colum can have multiple data set (in this case 4 dataset)
            // Each dataset can have a name
            yn.add("error");
            yn.add("time");
            yn.add("score");

            // Each colums can have multiple data-set
            for (size_t i = 0 ; i < num_sweep_test ; i++)
            {
                double score = score_solver(perf_amg_sweep_asym.get(sweep_start+i).time_solve,perf_amg_sweep_asym.get(sweep_start+i).residual);
                y.add({perf_amg_sweep_asym.get(sweep_start+i).residual,perf_amg_sweep_asym.get(sweep_start+i).time_solve,score});
            }

            // Google charts options
            GCoptions options;

            options.title = std::string("Coarsener: " + coars.get(k)) + std::string(": Overview of V cycle with n sweeps asymmetrically distributed");
            options.yAxis = std::string("error/time/score");
            options.xAxis = std::string("sweeps");
            options.stype = std::string("bars");
            options.more = std::string("vAxis: {scaleType: 'log'}");

            std::stringstream ss;

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

    void test_cycle_type(SparseMatrix<double,int,PETSC_BASE> & A,
                         const Vector<double,PETSC_BASE> & b,
                         openfpm::vector<std::string> & coarsener_to_test)
    {
        Vcluster & v_cl = create_vcluster();

        petsc_solver<double> pts;

        petsc_solver<double> solver;
        solver.setSolver(KSPRICHARDSON);
        solver.setAbsTol(0.01);
        solver.setMaxIter(1);


        for (size_t k = 0 ; k < coarsener_to_test.size(); k++)
        {
            for (size_t i = 1 ; i <= num_sweep_test ; i++)
            {

                solver.setPreconditioner(PCHYPRE_BOOMERAMG);
                solver.setPreconditionerAMG_nl(6);
                solver.setPreconditionerAMG_maxit(3);
                solver.setPreconditionerAMG_relax("SOR/Jacobi");
                solver.setPreconditionerAMG_cycleType("V",i,i);
                solver.setPreconditionerAMG_coarsenNodalType(0);
                solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k));
                solver.log_monitor();

                A.getMatrixTriplets();
                solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k).c_str());
                benchmark(A,b,solver,perf_amg_sweep_sym);

                double score = score_solver(perf_amg_sweep_sym.get(k*num_sweep_test+(i-1)).time_solve,perf_amg_sweep_sym.get(k*num_sweep_test+(i-1)).residual);
                perf_amg_sweep_sym.last().score = score;

                if (v_cl.getProcessUnitID() == 0)
                    std::cout << "Coarsener: " << coarsener_to_test.get(k) << "   Sweep: " << i << "  Score: " << score << std::endl;

                v_cycle_tested_sym.add(i);
            }
        }

    }

    void test_cycle_asym(SparseMatrix<double,int,PETSC_BASE> & A,
                         const Vector<double,PETSC_BASE> & b,
                         openfpm::vector<size_t> & ids_ts,
                         openfpm::vector<std::string> & coarsener_to_test)
    {
        Vcluster & v_cl = create_vcluster();

        petsc_solver<double> pts;

        petsc_solver<double> solver;
        solver.setSolver(KSPRICHARDSON);
        solver.setAbsTol(0.01);
        solver.setMaxIter(1);


        for (size_t k = 0 ; k < coarsener_to_test.size(); k++)
        {
            asym_tests.add(perf_amg_sweep_asym.size());
            size_t num_tot_sweep = 2*ids_ts.get(k);
            for (size_t i = 0 ; i <= num_tot_sweep ; i++)
            {

                solver.setPreconditioner(PCHYPRE_BOOMERAMG);
                solver.setPreconditionerAMG_nl(6);
                solver.setPreconditionerAMG_maxit(3);
                solver.setPreconditionerAMG_relax("SOR/Jacobi");
                solver.setPreconditionerAMG_cycleType("V",i,num_tot_sweep-i);
                solver.setPreconditionerAMG_coarsenNodalType(0);
                solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k));
                solver.log_monitor();

                A.getMatrixTriplets();
                solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k).c_str());

                benchmark(A,b,solver,perf_amg_sweep_asym);

                double score = score_solver(perf_amg_sweep_asym.last().time_solve,perf_amg_sweep_asym.last().residual);

                if (v_cl.getProcessUnitID() == 0)
                    std::cout << "Coarsener: " << coarsener_to_test.get(k) << "   Sweep: (" << i << "," << num_tot_sweep-i << ")  Score: " << score << std::endl;

                v_cycle_tested_asym.add(std::pair<size_t,size_t>(i,num_tot_sweep-i));
            }

        }

        asym_tests.add(perf_amg_sweep_asym.size());
    }

    void write_report(openfpm::vector<std::string> & coars)
    {
        if (create_vcluster().getProcessUnitID() != 0)
            return;

        GoogleChart cg;
        write_report_coars(cg);

        // Write V cycles performances
        write_report_cycle(cg,coars);

        write_report_cycle_asym(cg,coars);

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

    void best_coarsener(openfpm::vector<std::string> & coars)
    {
        int fastest = 0;
        double best_time = perf_amg_coars.get(0).time_solve;
        int accurate = 0;
        double best_accuracy = perf_amg_coars.get(0).residual;

        for (size_t i = 1 ; i < perf_amg_coars.size() ; i++)
        {
            if (perf_amg_coars.get(i).time_solve < best_time)
            {
                fastest = i;
                best_time = perf_amg_coars.get(i).time_solve;
            }

            if (perf_amg_coars.get(i).residual < best_accuracy)
            {
                accurate = i;
                best_accuracy = perf_amg_coars.get(i).residual;
            }
        }

        coars.add(method.get(fastest));
        coars.add(method.get(accurate));
    }

    void best_score(openfpm::vector<AMG_time_err_coars> & perf,
                    openfpm::vector<size_t> & sw_optimal,
                    size_t nm)
    {
        size_t page = perf.size() / nm;

        for (size_t k = 0 ; k < nm ; k++)
        {
            int score_id = page*k;
            double best_score = perf.get(score_id).score;

            for (size_t i = page*k ; i < page*k+page ; i++)
            {
                if (perf.get(i).score > best_score)
                {
                    score_id = i;
                    best_score = perf.get(i).score;
                }
            }

            sw_optimal.add(v_cycle_tested_sym.get(score_id));
        }
    }

public:

    void setTergetAMGAccuracy(double t_a)
    {
        target_accuracy = t_a;
    }

    void setMaxSweep(int max_sw)
    {
        num_sweep_test = max_sw;
    }

    void try_solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b)
    {
        test_coarsener(A,b);

        openfpm::vector<std::string> coa_methods;
        best_coarsener(coa_methods);

        test_cycle_type(A,b,coa_methods);
        openfpm::vector<size_t> sweep_optimal;
        best_score(perf_amg_sweep_sym,sweep_optimal,coa_methods.size());

        test_cycle_asym(A,b,sweep_optimal,coa_methods);

        write_report(coa_methods);
    }
};

#endif

#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_AMG_REPORT_HPP_ */