linear_solver.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file linear_solver.h
24  \brief File for linear solvers
25 */
26 #ifndef O2SCL_LINEAR_SOLVER_H
27 #define O2SCL_LINEAR_SOLVER_H
28 
29 #include <gsl/gsl_linalg.h>
30 
31 #include <o2scl/permutation.h>
32 #include <o2scl/lu.h>
33 #include <o2scl/qr.h>
34 #include <o2scl/hh.h>
35 
36 namespace o2scl_linalg {
37 
38  /** \brief A generic solver for the linear system \f$ A x = b \f$
39  [abstract base]
40 
41  A generic solver for dense linear systems.
42 
43  Those writing production level code should consider calling
44  LAPACK directly using \o2 objects as described in the \ref
45  linalg_section section of the User's Guide.
46 
47  \future The test code uses a Hilbert matrix, which is known
48  to be ill-conditioned, especially for the larger sizes. This
49  should probably be changed.
50  */
51  template<class vec_t=boost::numeric::ublas::vector<double>,
52  class mat_t=boost::numeric::ublas::matrix<double> >
53  class linear_solver {
54 
55  public:
56 
57  virtual ~linear_solver() {}
58 
59  /// Solve square linear system \f$ A x = b \f$ of size \c n
60  virtual void solve(size_t n, mat_t &a, vec_t &b, vec_t &x)=0;
61  };
62 
63  /** \brief Generic linear solver using LU decomposition
64  */
65  template <class vec_t=boost::numeric::ublas::vector<double>,
66  class mat_t=boost::numeric::ublas::matrix<double> >
67  class linear_solver_LU : public linear_solver<vec_t, mat_t> {
68 
69  public:
70 
71  /// Solve square linear system \f$ A x = b \f$ of size \c n
72  virtual void solve(size_t n, mat_t &A,
73  vec_t &b, vec_t &x) {
74  int sig;
75  o2scl::permutation p(n);
76  LU_decomp(n,A,p,sig);
77  LU_solve(n,A,p,b,x);
78  return;
79  };
80 
81  virtual ~linear_solver_LU() {}
82 
83  };
84 
85  /** \brief Generic linear solver using QR decomposition
86  */
87  template <class vec_t=boost::numeric::ublas::vector<double>,
88  class mat_t=boost::numeric::ublas::matrix<double> >
89  class linear_solver_QR : public linear_solver<vec_t, mat_t> {
90 
91  public:
92 
93  /// Solve square linear system \f$ A x = b \f$ of size \c n
94  virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x) {
96  QR_decomp(n,n,A,tau);
97  QR_solve(n,A,tau,b,x);
98  return;
99  };
100 
101  virtual ~linear_solver_QR() {}
102 
103  };
104 
105  /** \brief Generic Householder linear solver
106  */
107  template <class vec_t=boost::numeric::ublas::vector<double>,
108  class mat_t=boost::numeric::ublas::matrix<double> >
110  public linear_solver<vec_t, mat_t> {
111 
112  public:
113 
114  /// Solve square linear system \f$ A x = b \f$ of size \c n
115  virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x) {
116  HH_solve(n,A,b,x);
117  return;
118  };
119 
120  virtual ~linear_solver_HH() {}
121 
122  };
123 
124  // End of namespace o2scl_linalg
125 }
126 
127 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
128 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
129 #include <armadillo>
130 namespace o2scl_linalg {
131  /** \brief Armadillo linear solver
132 
133  This class is only defined if Armadillo support was enabled
134  during installation
135  */
136  template<class arma_vec_t, class arma_mat_t> class linear_solver_arma :
137  public linear_solver<arma_vec_t,arma_mat_t> {
138  virtual void solve(size_t n, arma_mat_t &A, arma_vec_t &b,
139  arma_vec_t &x) {
140  x=arma::solve(A,b);
141  return;
142  }
143  };
144 }
145 #endif
146 
147 #if defined (O2SCL_EIGEN) || defined (DOXYGEN)
148 #include <eigen3/Eigen/Dense>
149 namespace o2scl_linalg {
150  /** \brief Eigen linear solver using QR decomposition with
151  column pivoting
152 
153  This class is only defined if Eigen support was enabled during
154  installation.
155 
156  */
157  template<class eigen_vec_t, class eigen_mat_t>
159  public linear_solver<eigen_vec_t,eigen_mat_t> {
160  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
161  eigen_vec_t &x) {
162  x=A.householderQr().solve(b);
163  return;
164  }
165  };
166 
167  /** \brief Eigen linear solver using QR decomposition with
168  column pivoting
169 
170  This class is only defined if Eigen support was enabled during
171  installation.
172 
173  */
174  template<class eigen_vec_t, class eigen_mat_t>
176  public linear_solver<eigen_vec_t,eigen_mat_t> {
177  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
178  eigen_vec_t &x) {
179  x=A.colPivHouseholderQr().solve(b);
180  return;
181  }
182  };
183 
184  /** \brief Eigen linear solver using QR decomposition with
185  full pivoting
186 
187  This class is only defined if Eigen support was enabled during
188  installation.
189 
190  */
191  template<class eigen_vec_t, class eigen_mat_t>
193  public linear_solver<eigen_vec_t,eigen_mat_t> {
194  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
195  eigen_vec_t &x) {
196  x=A.fullPivHouseholderQr().solve(b);
197  return;
198  }
199  };
200 
201  /** \brief Eigen linear solver using LU decomposition with
202  partial pivoting
203 
204  This requires the matrix \c A to be invertible.
205 
206  This class is only defined if Eigen support was enabled during
207  installation.
208 
209  */
210  template<class eigen_vec_t, class eigen_mat_t>
212  public linear_solver<eigen_vec_t,eigen_mat_t> {
213  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
214  eigen_vec_t &x) {
215  x=A.partialPivLu().solve(b);
216  return;
217  }
218  };
219 
220  /** \brief Eigen linear solver using LU decomposition with
221  full pivoting
222 
223  This class is only defined if Eigen support was enabled during
224  installation.
225 
226  */
227  template<class eigen_vec_t, class eigen_mat_t>
229  public linear_solver<eigen_vec_t,eigen_mat_t> {
230  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
231  eigen_vec_t &x) {
232  x=A.fullPivLu().solve(b);
233  return;
234  }
235  };
236 
237  /** \brief Eigen linear solver using LLT decomposition with
238  full pivoting
239 
240  This requires the matrix \c A to be positive definite.
241 
242  This class is only defined if Eigen support was enabled during
243  installation.
244 
245  */
246  template<class eigen_vec_t, class eigen_mat_t>
248  public linear_solver<eigen_vec_t,eigen_mat_t> {
249  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
250  eigen_vec_t &x) {
251  x=A.llt().solve(b);
252  return;
253  }
254  };
255 
256  /** \brief Eigen linear solver using LDLT decomposition with
257  full pivoting
258 
259  This requires the matrix \c A to be positive or negative
260  semidefinite.
261 
262  This class is only defined if Eigen support was enabled during
263  installation.
264  */
265  template<class eigen_vec_t, class eigen_mat_t>
267  public linear_solver<eigen_vec_t,eigen_mat_t> {
268  virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b,
269  eigen_vec_t &x) {
270  x=A.ldlt().solve(b);
271  return;
272  }
273  };
274 }
275 #endif
276 
277 #else
278 #include <o2scl/linear_special.h>
279 #endif
280 
281 #endif
o2scl_linalg::linear_solver_LU::solve
virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:72
o2scl_linalg::linear_solver_eigen_colQR::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:177
o2scl_linalg::linear_solver_eigen_houseQR::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:160
o2scl_linalg::linear_solver_QR::solve
virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:94
o2scl_linalg
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
o2scl::permutation
A class for representing permutations.
Definition: permutation.h:70
o2scl_linalg::linear_solver_eigen_houseQR
Eigen linear solver using QR decomposition with column pivoting.
Definition: linear_solver.h:158
o2scl_linalg::linear_solver
A generic solver for the linear system [abstract base].
Definition: linear_solver.h:53
boost::numeric::ublas::vector< double >
o2scl_linalg::linear_solver_QR
Generic linear solver using QR decomposition.
Definition: linear_solver.h:89
o2scl_linalg::linear_solver_eigen_fullQR::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:194
o2scl_linalg::linear_solver_arma
Armadillo linear solver.
Definition: linear_solver.h:136
o2scl_linalg::LU_decomp
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Definition: lu_base.h:86
o2scl_linalg::linear_solver_eigen_LLT::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:249
o2scl_linalg::linear_solver_eigen_LLT
Eigen linear solver using LLT decomposition with full pivoting.
Definition: linear_solver.h:247
o2scl_linalg::linear_solver_eigen_fullLU::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:230
o2scl_linalg::QR_decomp
void QR_decomp(size_t M, size_t N, mat_t &A, vec_t &tau)
Compute the QR decomposition of matrix A.
Definition: qr_base.h:54
o2scl_linalg::linear_solver::solve
virtual void solve(size_t n, mat_t &a, vec_t &b, vec_t &x)=0
Solve square linear system of size n.
o2scl_linalg::linear_solver_eigen_LDLT
Eigen linear solver using LDLT decomposition with full pivoting.
Definition: linear_solver.h:266
o2scl_linalg::linear_solver_HH
Generic Householder linear solver.
Definition: linear_solver.h:109
o2scl_linalg::linear_solver_eigen_partLU::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:213
o2scl_linalg::HH_solve
int HH_solve(size_t n, mat_t &A, const vec_t &b, vec2_t &x)
Solve linear system after Householder decomposition.
Definition: hh_base.h:151
o2scl_linalg::linear_solver_eigen_LDLT::solve
virtual void solve(size_t n, eigen_mat_t &A, eigen_vec_t &b, eigen_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:268
o2scl_linalg::linear_solver_eigen_fullQR
Eigen linear solver using QR decomposition with full pivoting.
Definition: linear_solver.h:192
o2scl_linalg::linear_solver_eigen_colQR
Eigen linear solver using QR decomposition with column pivoting.
Definition: linear_solver.h:175
o2scl_linalg::linear_solver_LU
Generic linear solver using LU decomposition.
Definition: linear_solver.h:67
o2scl_linalg::linear_solver_arma::solve
virtual void solve(size_t n, arma_mat_t &A, arma_vec_t &b, arma_vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:138
o2scl_linalg::LU_solve
int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x)
Solve a linear system after LU decomposition.
Definition: lu_base.h:335
o2scl_linalg::linear_solver_HH::solve
virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:115
o2scl_linalg::QR_solve
void QR_solve(size_t N, const mat_t &QR, const vec_t &tau, const vec2_t &b, vec3_t &x)
Solve the system A x = b using the QR factorization.
Definition: qr_base.h:141
o2scl_linalg::linear_solver_eigen_fullLU
Eigen linear solver using LU decomposition with full pivoting.
Definition: linear_solver.h:228
o2scl_linalg::linear_solver_eigen_partLU
Eigen linear solver using LU decomposition with partial pivoting.
Definition: linear_solver.h:211

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).