Chemical Data Processing Library C++ API - Version 1.2.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 
44  template <typename E1, typename E2>
45  bool
46  solveLower(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
47  {
48  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
49  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
50 
51  if (e1().getSize1() != e1().getSize2())
52  return false;
53 
54  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
55  return false;
56 
57  SizeType size = e2().getSize();
58 
59  for (SizeType n = 0; n < size; n++) {
60  if (e1()(n, n) == typename E1::ValueType(0))
61  return false;
62 
63  ValueType t = e2()(n) /= e1()(n, n);
64 
65  if (t != ValueType(0))
66  for (SizeType m = n + 1; m < size; m++)
67  e2()(m) -= e1()(m, n) * t;
68  }
69 
70  return true;
71  }
72 
73  template <typename E1, typename E2>
74  bool
75  solveUnitLower(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
76  {
77  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
78  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
79 
80  if (e1().getSize1() != e1().getSize2())
81  return false;
82 
83  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
84  return false;
85 
86  SizeType size = e2().getSize();
87 
88  for (SizeType n = 0; n < size; n++) {
89  ValueType t = e2()(n);
90 
91  if (t != ValueType(0))
92  for (SizeType m = n + 1; m < size; m++)
93  e2()(m) -= e1()(m, n) * t;
94  }
95 
96  return true;
97  }
98 
99  template <typename E1, typename E2>
100  bool
101  solveLower(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
102  {
103  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
104  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
105 
106  if (e1().getSize1() != e1().getSize2())
107  return false;
108 
109  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
110  return false;
111 
112  SizeType size1 = e2().getSize1();
113  SizeType size2 = e2().getSize2();
114 
115  for (SizeType n = 0; n < size1; n++) {
116  if (e1()(n, n) == typename E1::ValueType(0))
117  return false;
118 
119  for (SizeType l = 0; l < size2; l++) {
120  ValueType t = e2()(n, l) /= e1()(n, n);
121 
122  if (t != ValueType(0))
123  for (SizeType m = n + 1; m < size1; m++)
124  e2()(m, l) -= e1()(m, n) * t;
125  }
126  }
127 
128  return true;
129  }
130 
131  template <typename E1, typename E2>
132  bool
133  solveUnitLower(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
134  {
135  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
136  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
137 
138  if (e1().getSize1() != e1().getSize2())
139  return false;
140 
141  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
142  return false;
143 
144  SizeType size1 = e2().getSize1();
145  SizeType size2 = e2().getSize2();
146 
147  for (SizeType n = 0; n < size1; n++) {
148  for (SizeType l = 0; l < size2; l++) {
149  ValueType t = e2()(n, l);
150 
151  if (t != ValueType(0))
152  for (SizeType m = n + 1; m < size1; m++)
153  e2()(m, l) -= e1()(m, n) * t;
154  }
155  }
156 
157  return true;
158  }
159 
160  template <typename E1, typename E2>
161  bool
162  solveUpper(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
163  {
164  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
165  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
166  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
167 
168  if (e1().getSize1() != e1().getSize2())
169  return false;
170 
171  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
172  return false;
173 
174  SizeType size = e2().getSize();
175 
176  for (DifferenceType n = size - 1; n >= 0; n--) {
177  if (e1()(n, n) == typename E1::ValueType(0))
178  return false;
179 
180  ValueType t = e2()(n) /= e1()(n, n);
181 
182  if (t != ValueType(0))
183  for (DifferenceType m = n - 1; m >= 0; m--)
184  e2()(m) -= e1()(m, n) * t;
185  }
186 
187  return true;
188  }
189 
190  template <typename E1, typename E2>
191  bool
192  solveUnitUpper(const MatrixExpression<E1>& e1, VectorExpression<E2>& e2)
193  {
194  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
195  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
196  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
197 
198  if (e1().getSize1() != e1().getSize2())
199  return false;
200 
201  if (SizeType(e1().getSize2()) != SizeType(e2().getSize()))
202  return false;
203 
204  SizeType size = e2().getSize();
205 
206  for (DifferenceType n = size - 1; n >= 0; n--) {
207  ValueType t = e2()(n);
208 
209  if (t != ValueType(0))
210  for (DifferenceType m = n - 1; m >= 0; m--)
211  e2()(m) -= e1()(m, n) * t;
212  }
213 
214  return true;
215  }
216 
217  template <typename E1, typename E2>
218  bool
219  solveUpper(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
220  {
221  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
222  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
223  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
224 
225  if (e1().getSize1() != e1().getSize2())
226  return false;
227 
228  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
229  return false;
230 
231  SizeType size1 = e2().getSize1();
232  SizeType size2 = e2().getSize2();
233 
234  for (DifferenceType n = size1 - 1; n >= 0; n--) {
235  if (e1()(n, n) == typename E1::ValueType(0))
236  return false;
237 
238  for (DifferenceType l = size2 - 1; l >= 0; l--) {
239  ValueType t = e2()(n, l) /= e1()(n, n);
240 
241  if (t != ValueType(0))
242  for (DifferenceType m = n - 1; m >= 0; m--)
243  e2()(m, l) -= e1()(m, n) * t;
244  }
245  }
246 
247  return true;
248  }
249 
250  template <typename E1, typename E2>
251  bool
252  solveUnitUpper(const MatrixExpression<E1>& e1, MatrixExpression<E2>& e2)
253  {
254  typedef typename CommonType<typename E1::SizeType, typename E2::SizeType>::Type SizeType;
255  typedef typename CommonType<typename E1::DifferenceType, typename E2::DifferenceType>::Type DifferenceType;
256  typedef typename CommonType<typename E1::ValueType, typename E2::ValueType>::Type ValueType;
257 
258  if (e1().getSize1() != e1().getSize2())
259  return false;
260 
261  if (SizeType(e1().getSize2()) != SizeType(e2().getSize1()))
262  return false;
263 
264  SizeType size1 = e2().getSize1();
265  SizeType size2 = e2().getSize2();
266 
267  for (DifferenceType n = size1 - 1; n >= 0; n--) {
268  for (DifferenceType l = size2 - 1; l >= 0; l--) {
269  ValueType t = e2()(n, l);
270 
271  if (t != ValueType(0))
272  for (DifferenceType m = n - 1; m >= 0; m--)
273  e2()(m, l) -= e1()(m, n) * t;
274  }
275  }
276 
277  return true;
278  }
279  } // namespace Math
280 } // namespace CDPL
281 
282 #endif // CDPL_MATH_LINEARSOLVE_HPP
CommonType.hpp
Common type deduction.
CDPL::Math::MatrixExpression
Definition: Expression.hpp:76
CDPL::Math::VectorExpression
Definition: Expression.hpp:54
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)
Definition: LinearSolve.hpp:162
CDPL::Math::solveUnitUpper
bool solveUnitUpper(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Definition: LinearSolve.hpp:192
CDPL::Math::solveUnitLower
bool solveUnitLower(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Definition: LinearSolve.hpp:75
CDPL::Math::solveLower
bool solveLower(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Definition: LinearSolve.hpp:46
CDPL
The namespace of the Chemical Data Processing Library.
CDPL::Math::CommonType::Type
std::common_type< T1, T2 >::type Type
Definition: CommonType.hpp:43
Generated by doxygen 1.9.1