Chemical Data Processing Library C++ API - Version 1.4.0
LinearSolve.hpp
Go to the documentation of this file.
1 /*
2  * LinearSolve.hpp
3  *
4  * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with this library; see the file COPYING. If not, write to
18  * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19  * Boston, MA 02111-1307, USA.
20  */
21 
27 #ifndef CDPL_MATH_LINEARSOLVE_HPP
28 #define CDPL_MATH_LINEARSOLVE_HPP
29 
30 #include "CDPL/Math/CommonType.hpp"
31 
32 
33 namespace CDPL
34 {
35 
36  namespace Math
37  {
38 
39  template <typename E>
40  class VectorExpression;
41  template <typename E>
42  class MatrixExpression;
43 
52  template <typename E1, typename E2>
53  bool
54  solveLower(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
55  {
56  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
57  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
58 
59  if (e1().getSize1() != e1().getSize2())
60  return false;
61 
62  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
63  return false;
64 
65  SizeType size = e2().getSize();
66 
67  for (SizeType n = 0; n < size; n++) {
68  if (e1()(n, n) == typename E1::ValueType(0))
69  return false;
70 
71  ValueType t = e2()(n) /= e1()(n, n);
72 
73  if (t != ValueType(0))
74  for (SizeType m = n + 1; m < size; m++)
75  e2()(m) -= e1()(m, n) * t;
76  }
77 
78  return true;
79  }
80 
89  template <typename E1, typename E2>
90  bool
91  solveUnitLower(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
92  {
93  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
94  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
95 
96  if (e1().getSize1() != e1().getSize2())
97  return false;
98 
99  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
100  return false;
101 
102  SizeType size = e2().getSize();
103 
104  for (SizeType n = 0; n < size; n++) {
105  ValueType t = e2()(n);
106 
107  if (t != ValueType(0))
108  for (SizeType m = n + 1; m < size; m++)
109  e2()(m) -= e1()(m, n) * t;
110  }
111 
112  return true;
113  }
114 
123  template <typename E1, typename E2>
124  bool
125  solveLower(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
126  {
127  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
128  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
129 
130  if (e1().getSize1() != e1().getSize2())
131  return false;
132 
133  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
134  return false;
135 
136  SizeType size1 = e2().getSize1();
137  SizeType size2 = e2().getSize2();
138 
139  for (SizeType n = 0; n < size1; n++) {
140  if (e1()(n, n) == typename E1::ValueType(0))
141  return false;
142 
143  for (SizeType l = 0; l < size2; l++) {
144  ValueType t = e2()(n, l) /= e1()(n, n);
145 
146  if (t != ValueType(0))
147  for (SizeType m = n + 1; m < size1; m++)
148  e2()(m, l) -= e1()(m, n) * t;
149  }
150  }
151 
152  return true;
153  }
154 
163  template <typename E1, typename E2>
164  bool
165  solveUnitLower(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
166  {
167  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
168  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
169 
170  if (e1().getSize1() != e1().getSize2())
171  return false;
172 
173  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
174  return false;
175 
176  SizeType size1 = e2().getSize1();
177  SizeType size2 = e2().getSize2();
178 
179  for (SizeType n = 0; n < size1; n++) {
180  for (SizeType l = 0; l < size2; l++) {
181  ValueType t = e2()(n, l);
182 
183  if (t != ValueType(0))
184  for (SizeType m = n + 1; m < size1; m++)
185  e2()(m, l) -= e1()(m, n) * t;
186  }
187  }
188 
189  return true;
190  }
191 
200  template <typename E1, typename E2>
201  bool
202  solveUpper(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
203  {
204  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
205  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
206  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
207 
208  if (e1().getSize1() != e1().getSize2())
209  return false;
210 
211  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
212  return false;
213 
214  SizeType size = e2().getSize();
215 
216  for (DifferenceType n = size - 1; n >= 0; n--) {
217  if (e1()(n, n) == typename E1::ValueType(0))
218  return false;
219 
220  ValueType t = e2()(n) /= e1()(n, n);
221 
222  if (t != ValueType(0))
223  for (DifferenceType m = n - 1; m >= 0; m--)
224  e2()(m) -= e1()(m, n) * t;
225  }
226 
227  return true;
228  }
229 
238  template <typename E1, typename E2>
239  bool
240  solveUnitUpper(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
241  {
242  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
243  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
244  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
245 
246  if (e1().getSize1() != e1().getSize2())
247  return false;
248 
249  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
250  return false;
251 
252  SizeType size = e2().getSize();
253 
254  for (DifferenceType n = size - 1; n >= 0; n--) {
255  ValueType t = e2()(n);
256 
257  if (t != ValueType(0))
258  for (DifferenceType m = n - 1; m >= 0; m--)
259  e2()(m) -= e1()(m, n) * t;
260  }
261 
262  return true;
263  }
264 
273  template <typename E1, typename E2>
274  bool
275  solveUpper(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
276  {
277  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
278  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
279  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
280 
281  if (e1().getSize1() != e1().getSize2())
282  return false;
283 
284  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
285  return false;
286 
287  SizeType size1 = e2().getSize1();
288  SizeType size2 = e2().getSize2();
289 
290  for (DifferenceType n = size1 - 1; n >= 0; n--) {
291  if (e1()(n, n) == typename E1::ValueType(0))
292  return false;
293 
294  for (DifferenceType l = size2 - 1; l >= 0; l--) {
295  ValueType t = e2()(n, l) /= e1()(n, n);
296 
297  if (t != ValueType(0))
298  for (DifferenceType m = n - 1; m >= 0; m--)
299  e2()(m, l) -= e1()(m, n) * t;
300  }
301  }
302 
303  return true;
304  }
305 
314  template <typename E1, typename E2>
315  bool
316  solveUnitUpper(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
317  {
318  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
319  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
320  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
321 
322  if (e1().getSize1() != e1().getSize2())
323  return false;
324 
325  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
326  return false;
327 
328  SizeType size1 = e2().getSize1();
329  SizeType size2 = e2().getSize2();
330 
331  for (DifferenceType n = size1 - 1; n >= 0; n--) {
332  for (DifferenceType l = size2 - 1; l >= 0; l--) {
333  ValueType t = e2()(n, l);
334 
335  if (t != ValueType(0))
336  for (DifferenceType m = n - 1; m >= 0; m--)
337  e2()(m, l) -= e1()(m, n) * t;
338  }
339  }
340 
341  return true;
342  }
343  } // namespace Math
344 } // namespace CDPL
345 
346 #endif // CDPL_MATH_LINEARSOLVE_HPP
CommonType.hpp
Common type deduction.
CDPL::Math::MatrixExpression
CRTP base class for all matrix expression types.
Definition: Expression.hpp:104
CDPL::Math::VectorExpression
CRTP base class for all vector expression types.
Definition: Expression.hpp:66
CDPL::Chem::CIPDescriptor::m
constexpr unsigned int m
Specifies that the stereocenter has m configuration.
Definition: CIPDescriptor.hpp:116
CDPL::Math::solveUpper
bool solveUpper(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Solves in place by back-substitution, where e1 is an upper-triangular matrix.
Definition: LinearSolve.hpp:202
CDPL::Math::solveUnitUpper
bool solveUnitUpper(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Solves in place by back-substitution, where e1 is a unit upper-triangular matrix (1 on the diagonal)...
Definition: LinearSolve.hpp:240
CDPL::Math::solveUnitLower
bool solveUnitLower(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Solves in place by forward-substitution, where e1 is a unit lower-triangular matrix (1 on the diagon...
Definition: LinearSolve.hpp:91
CDPL::Math::solveLower
bool solveLower(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Solves in place by forward-substitution, where e1 is a lower-triangular matrix.
Definition: LinearSolve.hpp:54
CDPL
The namespace of the Chemical Data Processing Library.
CDPL::Math::CommonType::Type
std::common_type< T1, T2 >::type Type
The common type.
Definition: CommonType.hpp:49
Generated by doxygen 1.9.1