1/*
2 * Umfpack_solver.hpp
3 *
4 * Created on: Nov 27, 2015
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_
9#define OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_
10
11#define UMFPACK_NONE 0
12
13#define SOLVER_NOOPTION 0
14#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
15#define SOLVER_PRINT_DETERMINANT 2
16
17#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
18
19/////// Compiled with EIGEN support
20
21#include "Vector/Vector.hpp"
22#include "Eigen/UmfPackSupport"
23#include <Eigen/SparseLU>
24
25
26template<typename T>
27class umfpack_solver
28{
29public:
30
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 void best_solve()
37 {
38 std::cerr << "Error Umfpack only support double precision, and int ad id type" << "/n";
39 }
40};
41
42
43template<>
44class umfpack_solver<double>
45{
46
47public:
48
49 /*! \brief Here we invert the matrix and solve the system
50 *
51 * \warning umfpack is not a parallel solver, this function work only with one processor
52 *
53 * \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet
54 *
55 * \tparam impl Implementation of the SparseMatrix
56 *
57 */
58 static Vector<double,EIGEN_BASE> try_solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
59 {
60 return solve(A,b,opt);
61 }
62
63 /*! \brief Here we invert the matrix and solve the system
64 *
65 * \warning umfpack is not a parallel solver, this function work only with one processor
66 *
67 * \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet
68 *
69 * \tparam impl Implementation of the SparseMatrix
70 *
71 */
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 {
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 // Collect the matrix on master
82 auto mat_ei = A.getMat();
83
84 Eigen::Matrix<double, Eigen::Dynamic, 1> x_ei;
85
86 // Collect the vector on master
87 auto b_ei = b.getVec();
88
89 // Copy b into x, this also copy the information on how to scatter back the information on x
90 x = b;
91
92 if (vcl.getProcessUnitID() == 0)
93 {
94 solver.compute(mat_ei);
95
96 if(solver.info()!=Eigen::Success)
97 {
98 // decomposition failed
99 std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n";
100
101 x.scatter();
102
103 return x;
104 }
105
106 x_ei = solver.solve(b_ei);
107
108 if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY)
109 {
110 Eigen::Matrix<double, Eigen::Dynamic, 1> res;
111 res = mat_ei * x_ei - b_ei;
112
113 std::cout << "Infinity norm: " << res.lpNorm<Eigen::Infinity>() << "\n";
114 }
115
116 if (opt & SOLVER_PRINT_DETERMINANT)
117 {
118 std::cout << " Determinant: " << solver.determinant() << "\n";
119 }
120
121 x = x_ei;
122 }
123
124 // Vector is only on master, scatter back the information
125 x.scatter();
126
127 return x;
128 }
129};
130
131#else
132
133/////// Compiled without EIGEN support
134
135#include "Vector/Vector.hpp"
136
137//! stub when library compiled without eigen
138template<typename T>
139class umfpack_solver
140{
141public:
142
143 //! stub solve
144 template<unsigned int impl, typename id_type> static Vector<T> solve(const SparseMatrix<T,id_type,impl> & A, const Vector<T,impl> & b)
145 {
146 std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
147 }
148
149 //! stub solve
150 void best_solve()
151 {
152 std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
153 }
154
155 //! stub solve
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 std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
159 }
160};
161
162//! stub when library compiled without eigen
163template<>
164class umfpack_solver<double>
165{
166
167public:
168
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 std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
173
174 Vector<double> x;
175
176 return x;
177 }
178
179 //! stub solve
180 void best_solve()
181 {
182 std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
183 }
184
185 //! stub solve
186 static Vector<double,EIGEN_BASE> try_solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
187 {
188 std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
189 return Vector<double,EIGEN_BASE>();
190 }
191};
192
193#endif
194
195
196#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */
197