// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN2_LU_H
#define EIGEN2_LU_H

namespace Eigen { 

template<typename MatrixType>
class LU : public FullPivLU<MatrixType>
{
  public:

    typedef typename MatrixType::Scalar Scalar;
    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> IntRowVectorType;
    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> IntColVectorType;
    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> RowVectorType;
    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> ColVectorType;

    typedef Matrix<typename MatrixType::Scalar,
                  MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix" is the number of cols of the original matrix
                                                 // so that the product "matrix * kernel = zero" makes sense
                  Dynamic,                       // we don't know at compile-time the dimension of the kernel
                  MatrixType::Options,
                  MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter
                  MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, whose dimension is the number
                                                   // of columns of the original matrix
    > KernelResultType;

    typedef Matrix<typename MatrixType::Scalar,
                   MatrixType::RowsAtCompileTime, // the image is a subspace of the destination space, whose dimension is the number
                                                  // of rows of the original matrix
                   Dynamic,                       // we don't know at compile time the dimension of the image (the rank)
                   MatrixType::Options,
                   MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix,
                   MatrixType::MaxColsAtCompileTime  // so it has the same number of rows and at most as many columns.
    > ImageResultType;

    typedef FullPivLU<MatrixType> Base;

    template<typename T>
    explicit LU(const T& t) : Base(t), m_originalMatrix(t) {}

    template<typename OtherDerived, typename ResultType>
    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
    {
      *result = static_cast<const Base*>(this)->solve(b);
      return true;
    }

    template<typename ResultType>
    inline void computeInverse(ResultType *result) const
    {
      solve(MatrixType::Identity(this->rows(), this->cols()), result);
    }
    
    template<typename KernelMatrixType>
    void computeKernel(KernelMatrixType *result) const
    {
      *result = static_cast<const Base*>(this)->kernel();
    }
    
    template<typename ImageMatrixType>
    void computeImage(ImageMatrixType *result) const
    {
      *result = static_cast<const Base*>(this)->image(m_originalMatrix);
    }
    
    const ImageResultType image() const
    {
      return static_cast<const Base*>(this)->image(m_originalMatrix);
    }
    
    const MatrixType& m_originalMatrix;
};

#if EIGEN2_SUPPORT_STAGE < STAGE20_RESOLVE_API_CONFLICTS
/** \lu_module
  *
  * Synonym of partialPivLu().
  *
  * \return the partial-pivoting LU decomposition of \c *this.
  *
  * \sa class PartialPivLU
  */
template<typename Derived>
inline const LU<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::lu() const
{
  return LU<PlainObject>(eval());
}
#endif

#ifdef EIGEN2_SUPPORT
/** \lu_module
  *
  * Synonym of partialPivLu().
  *
  * \return the partial-pivoting LU decomposition of \c *this.
  *
  * \sa class PartialPivLU
  */
template<typename Derived>
inline const LU<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::eigen2_lu() const
{
  return LU<PlainObject>(eval());
}
#endif

} // end namespace Eigen

#endif // EIGEN2_LU_H