umfpack_solver.hpp 4.87 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

template<typename T>
class umfpack_solver
{
public:

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

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

42

43
44
45
template<>
class umfpack_solver<double>
{
46

47
48
public:

49
50
51
52
53
	/*! \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
Pietro Incardona's avatar
Pietro Incardona committed
54
	 *
55
	 *	\tparam impl Implementation of the SparseMatrix
Pietro Incardona's avatar
Pietro Incardona committed
56
57
	 *
	 */
58
59
60
61
	static 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);
	}
Pietro Incardona's avatar
Pietro Incardona committed
62

63
64
65
66
67
68
69
70
71
	/*! \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
72
	static Vector<double,EIGEN_BASE> solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
73
	{
incardon's avatar
incardon committed
74
		Vcluster<> & vcl = create_vcluster();
75

76
		Vector<double> x;
77

78
		// only master processor solve
79
		Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
80

81
82
83
84
85
86
87
		// 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();
88

89
90
91
92
93
94
		// 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);
95

96
97
98
99
			if(solver.info()!=Eigen::Success)
			{
				// decomposition failed
				std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";
100
101
102

				x.scatter();

103
104
				return x;
			}
105

106
			x_ei = solver.solve(b_ei);
107

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

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

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

121
			x = x_ei;
122
		}
123
124
125
126

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

127
128
129
130
		return x;
	}
};

131
132
133
134
135
136
#else

/////// Compiled without EIGEN support

#include "Vector/Vector.hpp"

137
//! stub when library compiled without eigen
138
139
140
141
142
template<typename T>
class umfpack_solver
{
public:

143
	//! stub solve
incardon's avatar
incardon committed
144
	template<typename impl, typename id_type> static Vector<T> solve(const SparseMatrix<T,id_type,impl> & A, const Vector<T,impl> & b)
145
146
147
	{
		std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
	}
148

149
	//! stub solve
150
151
152
153
	void best_solve()
	{
		std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
	}
154
155

	//! stub solve
incardon's avatar
incardon committed
156
	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)
157
158
159
	{
		std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
	}
160
161
};

162
//! stub when library compiled without eigen
163
164
165
166
167
168
template<>
class umfpack_solver<double>
{

public:

169
	//! stub solve
170
	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)
171
172
173
174
175
176
177
	{
		std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";

		Vector<double> x;

		return x;
	}
178

179
	//! stub solve
180
181
182
183
	void best_solve()
	{
		std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
	}
184
185
186
187
188
189

	//! stub solve
	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";
	}
190
191
192
};

#endif
193
194
195


#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */