SkylineInplaceLU.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr>
5 //
6 // Eigen 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 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SKYLINEINPLACELU_H
26 #define EIGEN_SKYLINEINPLACELU_H
27 
28 namespace Eigen {
29 
39 template<typename MatrixType>
41 protected:
42  typedef typename MatrixType::Scalar Scalar;
43  typedef typename MatrixType::Index Index;
44 
45  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
46 
47 public:
48 
51  SkylineInplaceLU(MatrixType& matrix, int flags = 0)
52  : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) {
53  m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
54  m_lu.IsRowMajor ? computeRowMajor() : compute();
55  }
56 
67  void setPrecision(RealScalar v) {
68  m_precision = v;
69  }
70 
74  RealScalar precision() const {
75  return m_precision;
76  }
77 
86  void setFlags(int f) {
87  m_flags = f;
88  }
89 
91  int flags() const {
92  return m_flags;
93  }
94 
95  void setOrderingMethod(int m) {
96  m_flags = m;
97  }
98 
99  int orderingMethod() const {
100  return m_flags;
101  }
102 
104  void compute();
105  void computeRowMajor();
106 
108  //inline const MatrixType& matrixL() const { return m_matrixL; }
109 
111  //inline const MatrixType& matrixU() const { return m_matrixU; }
112 
113  template<typename BDerived, typename XDerived>
114  bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
115  const int transposed = 0) const;
116 
118  inline bool succeeded(void) const {
119  return m_succeeded;
120  }
121 
122 protected:
123  RealScalar m_precision;
124  int m_flags;
125  mutable int m_status;
126  bool m_succeeded;
127  MatrixType& m_lu;
128 };
129 
133 template<typename MatrixType>
134 //template<typename _Scalar>
136  const size_t rows = m_lu.rows();
137  const size_t cols = m_lu.cols();
138 
139  eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
140  eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage");
141 
142  for (Index row = 0; row < rows; row++) {
143  const double pivot = m_lu.coeffDiag(row);
144 
145  //Lower matrix Columns update
146  const Index& col = row;
147  for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
148  lIt.valueRef() /= pivot;
149  }
150 
151  //Upper matrix update -> contiguous memory access
152  typename MatrixType::InnerLowerIterator lIt(m_lu, col);
153  for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
154  typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
155  typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
156  const double coef = lIt.value();
157 
158  uItPivot += (rrow - row - 1);
159 
160  //update upper part -> contiguous memory access
161  for (++uItPivot; uIt && uItPivot;) {
162  uIt.valueRef() -= uItPivot.value() * coef;
163 
164  ++uIt;
165  ++uItPivot;
166  }
167  ++lIt;
168  }
169 
170  //Upper matrix update -> non contiguous memory access
171  typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
172  for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
173  typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
174  const double coef = lIt3.value();
175 
176  //update lower part -> non contiguous memory access
177  for (Index i = 0; i < rrow - row - 1; i++) {
178  m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
179  ++uItPivot;
180  }
181  ++lIt3;
182  }
183  //update diag -> contiguous
184  typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
185  for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
186 
187  typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
188  typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
189  const double coef = lIt2.value();
190 
191  uItPivot += (rrow - row - 1);
192  m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
193  ++lIt2;
194  }
195  }
196 }
197 
198 template<typename MatrixType>
200  const size_t rows = m_lu.rows();
201  const size_t cols = m_lu.cols();
202 
203  eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
204  eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !");
205 
206  for (Index row = 0; row < rows; row++) {
207  typename MatrixType::InnerLowerIterator llIt(m_lu, row);
208 
209 
210  for (Index col = llIt.col(); col < row; col++) {
211  if (m_lu.coeffExistLower(row, col)) {
212  const double diag = m_lu.coeffDiag(col);
213 
214  typename MatrixType::InnerLowerIterator lIt(m_lu, row);
215  typename MatrixType::InnerUpperIterator uIt(m_lu, col);
216 
217 
218  const Index offset = lIt.col() - uIt.row();
219 
220 
221  Index stop = offset > 0 ? col - lIt.col() : col - uIt.row();
222 
223  //#define VECTORIZE
224 #ifdef VECTORIZE
225  Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
226  Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
227 
228 
229  Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
230 #else
231  if (offset > 0) //Skip zero value of lIt
232  uIt += offset;
233  else //Skip zero values of uIt
234  lIt += -offset;
235  Scalar newCoeff = m_lu.coeffLower(row, col);
236 
237  for (Index k = 0; k < stop; ++k) {
238  const Scalar tmp = newCoeff;
239  newCoeff = tmp - lIt.value() * uIt.value();
240  ++lIt;
241  ++uIt;
242  }
243 #endif
244 
245  m_lu.coeffRefLower(row, col) = newCoeff / diag;
246  }
247  }
248 
249  //Upper matrix update
250  const Index col = row;
251  typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
252  for (Index rrow = uuIt.row(); rrow < col; rrow++) {
253 
254  typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
255  typename MatrixType::InnerUpperIterator uIt(m_lu, col);
256  const Index offset = lIt.col() - uIt.row();
257 
258  Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
259 
260 #ifdef VECTORIZE
261  Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
262  Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
263 
264  Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
265 #else
266  if (offset > 0) //Skip zero value of lIt
267  uIt += offset;
268  else //Skip zero values of uIt
269  lIt += -offset;
270  Scalar newCoeff = m_lu.coeffUpper(rrow, col);
271  for (Index k = 0; k < stop; ++k) {
272  const Scalar tmp = newCoeff;
273  newCoeff = tmp - lIt.value() * uIt.value();
274 
275  ++lIt;
276  ++uIt;
277  }
278 #endif
279  m_lu.coeffRefUpper(rrow, col) = newCoeff;
280  }
281 
282 
283  //Diag matrix update
284  typename MatrixType::InnerLowerIterator lIt(m_lu, row);
285  typename MatrixType::InnerUpperIterator uIt(m_lu, row);
286 
287  const Index offset = lIt.col() - uIt.row();
288 
289 
290  Index stop = offset > 0 ? lIt.size() : uIt.size();
291 #ifdef VECTORIZE
292  Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
293  Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
294  Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
295 #else
296  if (offset > 0) //Skip zero value of lIt
297  uIt += offset;
298  else //Skip zero values of uIt
299  lIt += -offset;
300  Scalar newCoeff = m_lu.coeffDiag(row);
301  for (Index k = 0; k < stop; ++k) {
302  const Scalar tmp = newCoeff;
303  newCoeff = tmp - lIt.value() * uIt.value();
304  ++lIt;
305  ++uIt;
306  }
307 #endif
308  m_lu.coeffRefDiag(row) = newCoeff;
309  }
310 }
311 
320 template<typename MatrixType>
321 template<typename BDerived, typename XDerived>
322 bool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const {
323  const size_t rows = m_lu.rows();
324  const size_t cols = m_lu.cols();
325 
326 
327  for (Index row = 0; row < rows; row++) {
328  x->coeffRef(row) = b.coeff(row);
329  Scalar newVal = x->coeff(row);
330  typename MatrixType::InnerLowerIterator lIt(m_lu, row);
331 
332  Index col = lIt.col();
333  while (lIt.col() < row) {
334 
335  newVal -= x->coeff(col++) * lIt.value();
336  ++lIt;
337  }
338 
339  x->coeffRef(row) = newVal;
340  }
341 
342 
343  for (Index col = rows - 1; col > 0; col--) {
344  x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
345 
346  const Scalar x_col = x->coeff(col);
347 
348  typename MatrixType::InnerUpperIterator uIt(m_lu, col);
349  uIt += uIt.size()-1;
350 
351 
352  while (uIt) {
353  x->coeffRef(uIt.row()) -= x_col * uIt.value();
354  //TODO : introduce --operator
355  uIt += -1;
356  }
357 
358 
359  }
360  x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
361 
362  return true;
363 }
364 
365 } // end namespace Eigen
366 
367 #endif // EIGEN_SKYLINELU_H