1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
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 EIGEN2_LU_H
11#define EIGEN2_LU_H
12
13namespace Eigen {
14
15template<typename MatrixType>
16class LU : public FullPivLU<MatrixType>
17{
18  public:
19
20    typedef typename MatrixType::Scalar Scalar;
21    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
22    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> IntRowVectorType;
23    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> IntColVectorType;
24    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> RowVectorType;
25    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> ColVectorType;
26
27    typedef Matrix<typename MatrixType::Scalar,
28                  MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix" is the number of cols of the original matrix
29                                                 // so that the product "matrix * kernel = zero" makes sense
30                  Dynamic,                       // we don't know at compile-time the dimension of the kernel
31                  MatrixType::Options,
32                  MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter
33                  MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, whose dimension is the number
34                                                   // of columns of the original matrix
35    > KernelResultType;
36
37    typedef Matrix<typename MatrixType::Scalar,
38                   MatrixType::RowsAtCompileTime, // the image is a subspace of the destination space, whose dimension is the number
39                                                  // of rows of the original matrix
40                   Dynamic,                       // we don't know at compile time the dimension of the image (the rank)
41                   MatrixType::Options,
42                   MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix,
43                   MatrixType::MaxColsAtCompileTime  // so it has the same number of rows and at most as many columns.
44    > ImageResultType;
45
46    typedef FullPivLU<MatrixType> Base;
47
48    template<typename T>
49    explicit LU(const T& t) : Base(t), m_originalMatrix(t) {}
50
51    template<typename OtherDerived, typename ResultType>
52    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
53    {
54      *result = static_cast<const Base*>(this)->solve(b);
55      return true;
56    }
57
58    template<typename ResultType>
59    inline void computeInverse(ResultType *result) const
60    {
61      solve(MatrixType::Identity(this->rows(), this->cols()), result);
62    }
63
64    template<typename KernelMatrixType>
65    void computeKernel(KernelMatrixType *result) const
66    {
67      *result = static_cast<const Base*>(this)->kernel();
68    }
69
70    template<typename ImageMatrixType>
71    void computeImage(ImageMatrixType *result) const
72    {
73      *result = static_cast<const Base*>(this)->image(m_originalMatrix);
74    }
75
76    const ImageResultType image() const
77    {
78      return static_cast<const Base*>(this)->image(m_originalMatrix);
79    }
80
81    const MatrixType& m_originalMatrix;
82};
83
84#if EIGEN2_SUPPORT_STAGE < STAGE20_RESOLVE_API_CONFLICTS
85/** \lu_module
86  *
87  * Synonym of partialPivLu().
88  *
89  * \return the partial-pivoting LU decomposition of \c *this.
90  *
91  * \sa class PartialPivLU
92  */
93template<typename Derived>
94inline const LU<typename MatrixBase<Derived>::PlainObject>
95MatrixBase<Derived>::lu() const
96{
97  return LU<PlainObject>(eval());
98}
99#endif
100
101#ifdef EIGEN2_SUPPORT
102/** \lu_module
103  *
104  * Synonym of partialPivLu().
105  *
106  * \return the partial-pivoting LU decomposition of \c *this.
107  *
108  * \sa class PartialPivLU
109  */
110template<typename Derived>
111inline const LU<typename MatrixBase<Derived>::PlainObject>
112MatrixBase<Derived>::eigen2_lu() const
113{
114  return LU<PlainObject>(eval());
115}
116#endif
117
118} // end namespace Eigen
119
120#endif // EIGEN2_LU_H
121