umfpack_solver.hpp

/*
 * Umfpack_solver.hpp
 *
 *  Created on: Nov 27, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_
#define OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_

#define UMFPACK_NONE 0

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


#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)


#include "Vector/Vector.hpp"
#include "Eigen/UmfPackSupport"
#include <Eigen/SparseLU>


template<typename T>
class umfpack_solver
{
public:

    template<unsigned int impl, typename id_type> static Vector<T> solve(const SparseMatrix<T,id_type,impl> & A, const Vector<T> & b)
    {
        std::cerr << "Error Umfpack only support double precision, and int ad id type" << "/n";
    }

    void best_solve()
    {
        std::cerr << "Error Umfpack only support double precision, and int ad id type" << "/n";
    }
};


template<>
class umfpack_solver<double>
{

    Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;

    Eigen::SparseMatrix<double,0,int> mat_ei;

public:

    Vector<double,EIGEN_BASE> try_solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
    {
        return solve(A,b,opt);
    }

    Vector<double,EIGEN_BASE> solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
    {
        Vcluster<> & vcl = create_vcluster();

        Vector<double> x;

        // Collect the matrix on master
        mat_ei = A.getMat();

        Eigen::Matrix<double, Eigen::Dynamic, 1> x_ei;

        // Collect the vector on master
        auto b_ei = b.getVec();

        // Copy b into x, this also copy the information on how to scatter back the information on x
        x = b;

        if (vcl.getProcessUnitID() == 0)
        {
            solver.compute(mat_ei);

            if(solver.info()!=Eigen::Success)
            {
                // Linear solver failed
                std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";

                x.scatter();

                return x;
            }

            x_ei = solver.solve(b_ei);

            //if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY)
            //{
                Eigen::Matrix<double, Eigen::Dynamic, 1> res;
                res = mat_ei * x_ei - b_ei;

                std::cout << "Infinity norm: " << res.lpNorm<Eigen::Infinity>() << "\n";
            //}

            if (opt & SOLVER_PRINT_DETERMINANT)
            {
                std::cout << " Determinant: " << solver.determinant() << "\n";
            }
            if (opt & SOLVER_PRINT)
            {
                std::cout << mat_ei << "\n";
                std::cout << b_ei << "\n";
            }

            x = x_ei;
        }

        // Vector is only on master, scatter back the information
        x.scatter();

        return x;
    }

    Vector<double,EIGEN_BASE> solve(const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
    {
        Vcluster<> & vcl = create_vcluster();

        Vector<double> x;

        Eigen::Matrix<double, Eigen::Dynamic, 1> x_ei;

        // Collect the vector on master
        auto b_ei = b.getVec();

        // Copy b into x, this also copy the information on how to scatter back the information on x
        x = b;

        if (vcl.getProcessUnitID() == 0)
        {
            x_ei = solver.solve(b_ei);

            if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY)
            {
                std::cout << "umfpack_solver: unsupported you have to pass the matrix for the option SOLVER_PRINT_RESIDUAL_NORM_INFINITY "  << "\n";
            }

            if (opt & SOLVER_PRINT_DETERMINANT)
            {
                std::cout << " Determinant: " << solver.determinant() << "\n";
            }

            x = x_ei;
        }

        // Vector is only on master, scatter back the information
        x.scatter();

        return x;
    }
};

#else


#include "Vector/Vector.hpp"

template<typename T>
class umfpack_solver
{
public:

    template<unsigned int impl, typename id_type> static Vector<T> solve(const SparseMatrix<T,id_type,impl> & A, const Vector<T,impl> & b)
    {
        std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
    }

    void best_solve()
    {
        std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
    }

    template<unsigned int impl, typename id_type> static Vector<T,impl> try_solve(SparseMatrix<T,id_type,impl> & A, const Vector<T,impl> & b, size_t opt = UMFPACK_NONE)
    {
        std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
    }
};

template<>
class umfpack_solver<double>
{

public:

    template<unsigned int impl, typename id_type> static Vector<double> solve(SparseMatrix<double,id_type,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE)
    {
        std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";

        Vector<double> x;

        return x;
    }

    void best_solve()
    {
        std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
    }

    static Vector<double,EIGEN_BASE> try_solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
    {
        std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
        return Vector<double,EIGEN_BASE>();
    }
};

#endif


#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */