Chemical Data Processing Library C++ API - Version 1.1.1
LUDecomposition.hpp
Go to the documentation of this file.
1 /*
2  * LUDecomposition.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_LUDECOMPOSITION_HPP
28 #define CDPL_MATH_LUDECOMPOSITION_HPP
29 
30 #include <algorithm>
31 
35 
36 
37 namespace CDPL
38 {
39 
40  namespace Math
41  {
42 
43  template <typename E>
44  typename E::SizeType
46  {
47  typedef typename E::SizeType SizeType;
48  typedef typename E::ValueType ValueType;
49  typedef typename MatrixRange<E>::RangeType MatrixRangeType;
50  typedef typename VectorRange<MatrixRow<E> >::RangeType RowRangeType;
51  typedef typename VectorRange<MatrixColumn<E> >::RangeType ColumnRangeType;
52 
53  SizeType size1 = e().getSize1();
54  SizeType size2 = e().getSize2();
55  SizeType size = std::min(size1, size2);
56 
57  if (size == 0)
58  return 0;
59 
60  SizeType singular = 0;
61 
62  for (SizeType i = 0; i < size; i++) {
63  MatrixColumn<E> col_i(column(e, i));
64  MatrixRow<E> row_i(row(e, i));
65 
66  if (e()(i, i) != ValueType(0)) {
67  ValueType m_inv = ValueType(1) / e()(i, i);
68  range(col_i, ColumnRangeType(i + 1, size1)) *= m_inv;
69 
70  } else if (singular == 0)
71  singular = i + 1;
72 
73  range(e, MatrixRangeType(i + 1, size1), MatrixRangeType(i + 1, size2))
74  .minusAssign(outerProd(range(col_i, ColumnRangeType(i + 1, size1)),
75  range(row_i, RowRangeType(i + 1, size2))));
76  }
77 
78  return singular;
79  }
80 
81  template <typename E, typename PV, typename T>
82  typename E::SizeType
83  luDecompose(MatrixExpression<E>& e, PV& pv, T& num_row_swaps)
84  {
85  typedef typename E::SizeType SizeType;
86  typedef typename E::ValueType ValueType;
87  typedef typename MatrixRange<E>::RangeType MatrixRangeType;
88  typedef typename VectorRange<MatrixRow<E> >::RangeType RowRangeType;
89  typedef typename VectorRange<MatrixColumn<E> >::RangeType ColumnRangeType;
90 
91  SizeType size1 = e().getSize1();
92  SizeType size2 = e().getSize2();
93  SizeType size = std::min(size1, size2);
94 
95  num_row_swaps = 0;
96 
97  if (size == 0)
98  return 0;
99 
100  SizeType singular = 0;
101 
102  for (SizeType i = 0; i < size; i++) {
103  MatrixColumn<E> col_i(column(e, i));
104  MatrixRow<E> row_i(row(e, i));
105  SizeType norm_inf_idx = i + normInfIndex(range(col_i, ColumnRangeType(i, size1)));
106 
107  if (e()(norm_inf_idx, i) != ValueType(0)) {
108  pv[i] = norm_inf_idx;
109 
110  if (norm_inf_idx != i) {
111  row(e, norm_inf_idx).swap(row_i);
112  num_row_swaps++;
113  }
114 
115  ValueType m_inv = ValueType(1) / e()(i, i);
116  range(col_i, ColumnRangeType(i + 1, size1)) *= m_inv;
117 
118  } else if (singular == 0)
119  singular = i + 1;
120 
121  range(e, MatrixRangeType(i + 1, size1), MatrixRangeType(i + 1, size2))
122  .minusAssign(outerProd(range(col_i, ColumnRangeType(i + 1, size1)),
123  range(row_i, RowRangeType(i + 1, size2))));
124  }
125 
126  return singular;
127  }
128 
129  template <typename E, typename PV>
130  void
131  swapRows(VectorExpression<E>& e, const PV& pv)
132  {
133  typedef typename E::SizeType SizeType;
134 
135  for (SizeType i = 0, size = e().getSize(); i < size; i++) {
136  if (i != SizeType(pv[i]))
137  std::swap(e()(i), e()(pv[i]));
138  }
139  }
140 
141  template <typename E, typename PV>
142  void
143  swapRows(MatrixExpression<E>& e, const PV& pv)
144  {
145  typedef typename E::SizeType SizeType;
146 
147  for (SizeType i = 0, size = e().getSize1(); i < size; i++) {
148  if (i != SizeType(pv[i])) {
149  MatrixRow<E> other_row(e(), pv[i]);
150  row(e, i).swap(other_row);
151  }
152  }
153  }
154 
155  template <typename E1, typename E2>
156  bool
158  {
159  if (!solveUnitLower(lu, b))
160  return false;
161 
162  return solveUpper(lu, b);
163  }
164 
165  template <typename E1, typename E2, typename PV>
166  bool
168  {
169  swapRows(b, pv);
170 
171  return luSubstitute(lu, b);
172  }
173 
174  template <typename E1, typename E2>
175  bool
177  {
178  if (!solveUnitLower(lu, b))
179  return false;
180 
181  return solveUpper(lu, b);
182  }
183 
184  template <typename E1, typename E2, typename PV>
185  bool
187  {
188  swapRows(b, pv);
189 
190  return luSubstitute(lu, b);
191  }
192  } // namespace Math
193 } // namespace CDPL
194 
195 #endif // CDPL_MATH_LUDECOMPOSITION_HPP
CDPL::Math::outerProd
VectorMatrixBinaryTraits< E1, E2, ScalarMultiplication< typename E1::ValueType, typename E2::ValueType > >::ResultType outerProd(const VectorExpression< E1 > &e1, const VectorExpression< E2 > &e2)
Definition: MatrixExpression.hpp:794
CDPL::Math::column
MatrixColumn< M > column(MatrixExpression< M > &e, typename MatrixColumn< M >::SizeType j)
Definition: MatrixProxy.hpp:730
CDPL::Math::VectorExpression
Definition: Expression.hpp:54
CDPL::Math::solveUpper
bool solveUpper(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Definition: LinearSolve.hpp:162
CDPL::Math::Range< SizeType >
CDPL::Math::MatrixExpression
Definition: Expression.hpp:76
CDPL::Math::luDecompose
E::SizeType luDecompose(MatrixExpression< E > &e)
Definition: LUDecomposition.hpp:45
VectorProxy.hpp
Definition of vector proxy types.
CDPL::Math::MatrixRow
Definition: MatrixProxy.hpp:49
CDPL::Math::MatrixColumn
Definition: MatrixProxy.hpp:196
CDPL::Math::swapRows
void swapRows(VectorExpression< E > &e, const PV &pv)
Definition: LUDecomposition.hpp:131
CDPL::Chem::AtomType::T
const unsigned int T
Specifies Hydrogen (Tritium).
Definition: AtomType.hpp:67
CDPL::Math::solveUnitLower
bool solveUnitLower(const MatrixExpression< E1 > &e1, VectorExpression< E2 > &e2)
Definition: LinearSolve.hpp:75
CDPL
The namespace of the Chemical Data Processing Library.
LinearSolve.hpp
Functions for solving linear equations.
MatrixProxy.hpp
Definition of matrix proxy types.
CDPL::Math::row
MatrixRow< M > row(MatrixExpression< M > &e, typename MatrixRow< M >::SizeType i)
Definition: MatrixProxy.hpp:716
CDPL::Math::VectorRange
Definition: VectorProxy.hpp:48
CDPL::Math::range
MatrixRange< E > range(MatrixExpression< E > &e, const typename MatrixRange< E >::RangeType &r1, const typename MatrixRange< E >::RangeType &r2)
Definition: MatrixProxy.hpp:744
CDPL::Math::luSubstitute
bool luSubstitute(const MatrixExpression< E1 > &lu, VectorExpression< E2 > &b)
Definition: LUDecomposition.hpp:157
CDPL::Math::normInfIndex
VectorNormInfinityIndex< E >::ResultType normInfIndex(const VectorExpression< E > &e)
Definition: VectorExpression.hpp:546