umfpack_solver.hpp 1.52 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
16
17
18
19
20
21
22
23
24
25
26

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";
	}
};

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

31
32
33
34
35
template<>
class umfpack_solver<double>
{
public:

36
	template<typename impl> static Vector<double> solve(const SparseMatrix<double,int,impl> & A, const Vector<double> & b, size_t opt = NONE)
37
38
39
40
41
42
43
44
45
46
47
48
49
	{
		Vector<double> x;

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

		solver.compute(A.getMat());

		if(solver.info()!=Eigen::Success)
		{
			// decomposition failed
			std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";
			return x;
		}
50

51
52
		x.getVec() = solver.solve(b.getVec());

53
54
55
56
57
58
59
60
61
62
63
64
65
		if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY)
		{
			Vector<double> res;
			res.getVec() = A.getMat() * x.getVec() - b.getVec();

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

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

66
67
68
69
70
71
72
		return x;
	}
};



#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */