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// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_SKYLINEINPLACELU_H
11#define EIGEN_SKYLINEINPLACELU_H
12
13namespace Eigen {
14
15/** \ingroup Skyline_Module
16 *
17 * \class SkylineInplaceLU
18 *
19 * \brief Inplace LU decomposition of a skyline matrix and associated features
20 *
21 * \param MatrixType the type of the matrix of which we are computing the LU factorization
22 *
23 */
24template<typename MatrixType>
25class SkylineInplaceLU {
26protected:
27    typedef typename MatrixType::Scalar Scalar;
28    typedef typename MatrixType::Index Index;
29
30    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
31
32public:
33
34    /** Creates a LU object and compute the respective factorization of \a matrix using
35     * flags \a flags. */
36    SkylineInplaceLU(MatrixType& matrix, int flags = 0)
37    : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) {
38        m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
39        m_lu.IsRowMajor ? computeRowMajor() : compute();
40    }
41
42    /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
43     *
44     * Setting a value greater than zero speeds up computation, and yields to an imcomplete
45     * factorization with fewer non zero coefficients. Such approximate factors are especially
46     * useful to initialize an iterative solver.
47     *
48     * Note that the exact meaning of this parameter might depends on the actual
49     * backend. Moreover, not all backends support this feature.
50     *
51     * \sa precision() */
52    void setPrecision(RealScalar v) {
53        m_precision = v;
54    }
55
56    /** \returns the current precision.
57     *
58     * \sa setPrecision() */
59    RealScalar precision() const {
60        return m_precision;
61    }
62
63    /** Sets the flags. Possible values are:
64     *  - CompleteFactorization
65     *  - IncompleteFactorization
66     *  - MemoryEfficient
67     *  - one of the ordering methods
68     *  - etc...
69     *
70     * \sa flags() */
71    void setFlags(int f) {
72        m_flags = f;
73    }
74
75    /** \returns the current flags */
76    int flags() const {
77        return m_flags;
78    }
79
80    void setOrderingMethod(int m) {
81        m_flags = m;
82    }
83
84    int orderingMethod() const {
85        return m_flags;
86    }
87
88    /** Computes/re-computes the LU factorization */
89    void compute();
90    void computeRowMajor();
91
92    /** \returns the lower triangular matrix L */
93    //inline const MatrixType& matrixL() const { return m_matrixL; }
94
95    /** \returns the upper triangular matrix U */
96    //inline const MatrixType& matrixU() const { return m_matrixU; }
97
98    template<typename BDerived, typename XDerived>
99    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
100            const int transposed = 0) const;
101
102    /** \returns true if the factorization succeeded */
103    inline bool succeeded(void) const {
104        return m_succeeded;
105    }
106
107protected:
108    RealScalar m_precision;
109    int m_flags;
110    mutable int m_status;
111    bool m_succeeded;
112    MatrixType& m_lu;
113};
114
115/** Computes / recomputes the in place LU decomposition of the SkylineInplaceLU.
116 * using the default algorithm.
117 */
118template<typename MatrixType>
119//template<typename _Scalar>
120void SkylineInplaceLU<MatrixType>::compute() {
121    const size_t rows = m_lu.rows();
122    const size_t cols = m_lu.cols();
123
124    eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
125    eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage");
126
127    for (Index row = 0; row < rows; row++) {
128        const double pivot = m_lu.coeffDiag(row);
129
130        //Lower matrix Columns update
131        const Index& col = row;
132        for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
133            lIt.valueRef() /= pivot;
134        }
135
136        //Upper matrix update -> contiguous memory access
137        typename MatrixType::InnerLowerIterator lIt(m_lu, col);
138        for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
139            typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
140            typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
141            const double coef = lIt.value();
142
143            uItPivot += (rrow - row - 1);
144
145            //update upper part  -> contiguous memory access
146            for (++uItPivot; uIt && uItPivot;) {
147                uIt.valueRef() -= uItPivot.value() * coef;
148
149                ++uIt;
150                ++uItPivot;
151            }
152            ++lIt;
153        }
154
155        //Upper matrix update -> non contiguous memory access
156        typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
157        for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
158            typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
159            const double coef = lIt3.value();
160
161            //update lower part ->  non contiguous memory access
162            for (Index i = 0; i < rrow - row - 1; i++) {
163                m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
164                ++uItPivot;
165            }
166            ++lIt3;
167        }
168        //update diag -> contiguous
169        typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
170        for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
171
172            typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
173            typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
174            const double coef = lIt2.value();
175
176            uItPivot += (rrow - row - 1);
177            m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
178            ++lIt2;
179        }
180    }
181}
182
183template<typename MatrixType>
184void SkylineInplaceLU<MatrixType>::computeRowMajor() {
185    const size_t rows = m_lu.rows();
186    const size_t cols = m_lu.cols();
187
188    eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
189    eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !");
190
191    for (Index row = 0; row < rows; row++) {
192        typename MatrixType::InnerLowerIterator llIt(m_lu, row);
193
194
195        for (Index col = llIt.col(); col < row; col++) {
196            if (m_lu.coeffExistLower(row, col)) {
197                const double diag = m_lu.coeffDiag(col);
198
199                typename MatrixType::InnerLowerIterator lIt(m_lu, row);
200                typename MatrixType::InnerUpperIterator uIt(m_lu, col);
201
202
203                const Index offset = lIt.col() - uIt.row();
204
205
206                Index stop = offset > 0 ? col - lIt.col() : col - uIt.row();
207
208                //#define VECTORIZE
209#ifdef VECTORIZE
210                Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
211                Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
212
213
214                Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
215#else
216                if (offset > 0) //Skip zero value of lIt
217                    uIt += offset;
218                else //Skip zero values of uIt
219                    lIt += -offset;
220                Scalar newCoeff = m_lu.coeffLower(row, col);
221
222                for (Index k = 0; k < stop; ++k) {
223                    const Scalar tmp = newCoeff;
224                    newCoeff = tmp - lIt.value() * uIt.value();
225                    ++lIt;
226                    ++uIt;
227                }
228#endif
229
230                m_lu.coeffRefLower(row, col) = newCoeff / diag;
231            }
232        }
233
234        //Upper matrix update
235        const Index col = row;
236        typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
237        for (Index rrow = uuIt.row(); rrow < col; rrow++) {
238
239            typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
240            typename MatrixType::InnerUpperIterator uIt(m_lu, col);
241            const Index offset = lIt.col() - uIt.row();
242
243            Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
244
245#ifdef VECTORIZE
246            Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
247            Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
248
249            Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
250#else
251            if (offset > 0) //Skip zero value of lIt
252                uIt += offset;
253            else //Skip zero values of uIt
254                lIt += -offset;
255            Scalar newCoeff = m_lu.coeffUpper(rrow, col);
256            for (Index k = 0; k < stop; ++k) {
257                const Scalar tmp = newCoeff;
258                newCoeff = tmp - lIt.value() * uIt.value();
259
260                ++lIt;
261                ++uIt;
262            }
263#endif
264            m_lu.coeffRefUpper(rrow, col) = newCoeff;
265        }
266
267
268        //Diag matrix update
269        typename MatrixType::InnerLowerIterator lIt(m_lu, row);
270        typename MatrixType::InnerUpperIterator uIt(m_lu, row);
271
272        const Index offset = lIt.col() - uIt.row();
273
274
275        Index stop = offset > 0 ? lIt.size() : uIt.size();
276#ifdef VECTORIZE
277        Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
278        Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
279        Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
280#else
281        if (offset > 0) //Skip zero value of lIt
282            uIt += offset;
283        else //Skip zero values of uIt
284            lIt += -offset;
285        Scalar newCoeff = m_lu.coeffDiag(row);
286        for (Index k = 0; k < stop; ++k) {
287            const Scalar tmp = newCoeff;
288            newCoeff = tmp - lIt.value() * uIt.value();
289            ++lIt;
290            ++uIt;
291        }
292#endif
293        m_lu.coeffRefDiag(row) = newCoeff;
294    }
295}
296
297/** Computes *x = U^-1 L^-1 b
298 *
299 * If \a transpose is set to SvTranspose or SvAdjoint, the solution
300 * of the transposed/adjoint system is computed instead.
301 *
302 * Not all backends implement the solution of the transposed or
303 * adjoint system.
304 */
305template<typename MatrixType>
306template<typename BDerived, typename XDerived>
307bool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const {
308    const size_t rows = m_lu.rows();
309    const size_t cols = m_lu.cols();
310
311
312    for (Index row = 0; row < rows; row++) {
313        x->coeffRef(row) = b.coeff(row);
314        Scalar newVal = x->coeff(row);
315        typename MatrixType::InnerLowerIterator lIt(m_lu, row);
316
317        Index col = lIt.col();
318        while (lIt.col() < row) {
319
320            newVal -= x->coeff(col++) * lIt.value();
321            ++lIt;
322        }
323
324        x->coeffRef(row) = newVal;
325    }
326
327
328    for (Index col = rows - 1; col > 0; col--) {
329        x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
330
331        const Scalar x_col = x->coeff(col);
332
333        typename MatrixType::InnerUpperIterator uIt(m_lu, col);
334        uIt += uIt.size()-1;
335
336
337        while (uIt) {
338            x->coeffRef(uIt.row()) -= x_col * uIt.value();
339            //TODO : introduce --operator
340            uIt += -1;
341        }
342
343
344    }
345    x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
346
347    return true;
348}
349
350} // end namespace Eigen
351
352#endif // EIGEN_SKYLINELU_H
353