umfpack_solver.hpp 2.34 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
/*
 * 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_

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

15
#define UMFPACK_NONE 0
16 17 18 19 20 21 22 23 24 25 26 27

template<typename T>
class umfpack_solver
{
public:

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

28 29 30 31
#define SOLVER_NOOPTION 0
#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
#define SOLVER_PRINT_DETERMINANT 2

32 33 34
template<>
class umfpack_solver<double>
{
35

36 37
public:

38 39 40 41 42 43 44 45 46
	/*! \brief Here we invert the matrix and solve the system
	 *
	 *  \warning umfpack is not a parallel solver, this function work only with one processor
	 *
	 *  \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet
	 *
	 *	\tparam impl Implementation of the SparseMatrix
	 *
	 */
47
	template<typename impl> static Vector<double> solve(SparseMatrix<double,int,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE)
48
	{
49
		Vcluster & vcl = create_vcluster();
50

51
		Vector<double> x;
52

53
		// only master processor solve
54
		Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
55

56 57 58 59 60 61 62
		// Collect the matrix on master
		auto mat_ei = A.getMat();

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

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

64 65 66 67 68 69
		// 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);
70

71 72 73 74 75 76
			if(solver.info()!=Eigen::Success)
			{
				// decomposition failed
				std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";
				return x;
			}
77

78
			x_ei = solver.solve(b_ei);
79

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

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

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

93
			x = x_ei;
94
		}
95 96 97 98

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

99 100 101 102 103 104 105
		return x;
	}
};



#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */