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

11
12
13
14
15
16
17
18
19
20
#define UMFPACK_NONE 0

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

#ifdef HAVE_EIGEN

/////// Compiled with EIGEN support

21
22
#include "Vector/Vector.hpp"
#include "Eigen/UmfPackSupport"
23
24
#include <Eigen/SparseLU>

25
26
27
28
29
30
31
32
33
34
35
36

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

37

38
39
40
template<>
class umfpack_solver<double>
{
41

42
43
public:

44
45
46
47
48
49
50
51
52
	/*! \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
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
53
	static Vector<double,EIGEN_BASE> solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
54
	{
55
		Vcluster & vcl = create_vcluster();
56

57
		Vector<double> x;
58

59
		// only master processor solve
60
		Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
61

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

70
71
72
73
74
75
		// 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);
76

77
78
79
80
81
82
			if(solver.info()!=Eigen::Success)
			{
				// decomposition failed
				std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";
				return x;
			}
83

84
			x_ei = solver.solve(b_ei);
85

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

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

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

99
			x = x_ei;
100
		}
101
102
103
104

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

105
106
107
108
		return x;
	}
};

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#else

/////// Compiled without EIGEN support

#include "Vector/Vector.hpp"

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 << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
	}
};


template<>
class umfpack_solver<double>
{

public:

	template<typename impl> static Vector<double> solve(SparseMatrix<double,int,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;
	}
};

#endif
144
145
146


#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */