linear_solver.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, 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
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
A generic solver for the linear system [abstract base].
Definition: linear_solver.h:53
Eigen linear solver using QR decomposition with column pivoting.
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.
Eigen linear solver using LU decomposition with full pivoting.
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.
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.
A class for representing permutations.
Definition: permutation.h:70
virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x)
Solve square linear system of size n.
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
Generic Householder linear solver.
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
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.
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.
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.
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.
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
Eigen linear solver using QR decomposition with column pivoting.
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
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
Armadillo linear solver.
Eigen linear solver using LU decomposition with partial pivoting.
Eigen linear solver using LDLT decomposition with full pivoting.
Eigen linear solver using LLT decomposition with full pivoting.
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:248
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
virtual void solve(size_t n, mat_t &a, vec_t &b, vec_t &x)=0
Solve square linear system of size n.
Generic linear solver using QR decomposition.
Definition: linear_solver.h:89
Eigen linear solver using QR decomposition with full pivoting.
Generic linear solver using LU decomposition.
Definition: linear_solver.h:67
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.

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