17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library
27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra.
37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez//
47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez//
67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla
77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed
87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_SUITESPARSEQRSUPPORT_H
117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_SUITESPARSEQRSUPPORT_H
127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen {
147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename MatrixType> class SPQR;
167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename SPQRType> struct SPQRMatrixQReturnType;
177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename SPQRType> struct SPQRMatrixQTransposeReturnType;
187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template <typename SPQRType, typename Derived> struct SPQR_QProduct;
197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  namespace internal {
207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> >
217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      typedef typename SPQRType::MatrixType ReturnType;
237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    };
247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      typedef typename SPQRType::MatrixType ReturnType;
277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    };
287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> >
297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      typedef typename Derived::PlainObject ReturnType;
317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    };
327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  } // End namespace internal
337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/**
357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \ingroup SPQRSupport_Module
367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class SPQR
377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Sparse QR factorization based on SuiteSparseQR library
387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition
407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * of sparse matrices. The result is then used to solve linear leasts_square systems.
417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Clearly, a QR factorization is returned such that A*P = Q*R where :
427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * P is the column permutation. Use colsPermutation() to get it.
447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Q is the orthogonal matrix represented as Householder reflectors.
467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * You can then apply it to a vector.
487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index
517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE
547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */
567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType>
577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass SPQR
587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  public:
607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef typename _MatrixType::Scalar Scalar;
617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef typename _MatrixType::RealScalar RealScalar;
627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef UF_long Index ;
637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType;
647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef PermutationMatrix<Dynamic, Dynamic> PermutationType;
657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  public:
667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    SPQR()
677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      : m_isInitialized(false),
687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_ordering(SPQR_ORDERING_DEFAULT),
697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_allow_tol(SPQR_DEFAULT_TOL),
707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_tolerance (NumTraits<Scalar>::epsilon())
717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_l_start(&m_cc);
737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    SPQR(const _MatrixType& matrix)
767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    : m_isInitialized(false),
777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_ordering(SPQR_ORDERING_DEFAULT),
787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_allow_tol(SPQR_DEFAULT_TOL),
797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_tolerance (NumTraits<Scalar>::epsilon())
807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_l_start(&m_cc);
827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      compute(matrix);
837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    ~SPQR()
867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      SPQR_free();
887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_l_finish(&m_cc);
897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void SPQR_free()
917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_l_free_sparse(&m_H, &m_cc);
937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_l_free_sparse(&m_cR, &m_cc);
947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_l_free_dense(&m_HTau, &m_cc);
957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      std::free(m_E);
967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      std::free(m_HPinv);
977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void compute(const _MatrixType& matrix)
1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if(m_isInitialized) SPQR_free();
1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      MatrixType mat(matrix);
1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      cholmod_sparse A;
1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      A = viewAsCholmod(mat);
1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index col = matrix.cols();
1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A,
1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez                             &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if (!m_cR)
1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      {
1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_info = NumericalIssue;
1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_isInitialized = false;
1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        return;
1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      }
1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_info = Success;
1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_isInitialized = true;
1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_isRUpToDate = false;
1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Get the number of rows of the input matrix and the Q matrix
1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Index rows() const {return m_H->nrow; }
1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Get the number of columns of the input matrix.
1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Index cols() const { return m_cR->ncol; }
1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      *
1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      * \sa compute()
1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      */
1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template<typename Rhs>
1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const
1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(this->rows()==B.rows()
1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez                    && "SPQR::solve(): invalid number of rows of the right hand side matrix B");
1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          return internal::solve_retval<SPQR, Rhs>(*this, B.derived());
1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template<typename Rhs, typename Dest>
1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(b.cols()==1 && "This method is for vectors only");
1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      //Compute Q^T * b
1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      typename Dest::PlainObject y;
1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      y = matrixQ().transpose() * b;
1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // Solves with the triangular matrix R
1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index rk = this->rank();
1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y.topRows(rk));
1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      y.bottomRows(cols()-rk).setZero();
1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      // Apply the column permutation
1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      dest.topRows(cols()) = colsPermutation() * y.topRows(cols());
1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_info = Success;
1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /** \returns the sparse triangular factor R. It is a sparse matrix
1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const MatrixType matrixR() const
1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if(!m_isRUpToDate) {
1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_isRUpToDate = true;
1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      }
1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_R;
1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /// Get an expression of the matrix Q
1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    SPQRMatrixQReturnType<SPQR> matrixQ() const
1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return SPQRMatrixQReturnType<SPQR>(*this);
1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /// Get the permutation that was applied to columns of A
1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    PermutationType colsPermutation() const
1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index n = m_cR->ncol;
1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      PermutationType colsPerm(n);
1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j];
1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return colsPerm;
1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Gets the rank of the matrix.
1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * It should be equal to matrixQR().cols if the matrix is full-rank
1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index rank() const
1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_cc.SPQR_istat[4];
1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /// Set the fill-reducing ordering method to be used
1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void setSPQROrdering(int ord) { m_ordering = ord;}
1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /// Set the tolerance tol to treat columns with 2-norm < =tol as zero
2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void setPivotThreshold(const RealScalar& tol) { m_tolerance = tol; }
2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /** \returns a pointer to the SPQR workspace */
2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    cholmod_common *cholmodCommon() const { return &m_cc; }
2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /** \brief Reports whether previous computation was successful.
2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      *
2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      * \returns \c Success if computation was succesful,
2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      *          \c NumericalIssue if the sparse QR can not be computed
2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      */
2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    ComputationInfo info() const
2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_info;
2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  protected:
2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    bool m_isInitialized;
2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    bool m_analysisIsOk;
2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    bool m_factorizationIsOk;
2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable bool m_isRUpToDate;
2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable ComputationInfo m_info;
2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    int m_ordering; // Ordering method to use, see SPQR's manual
2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    int m_allow_tol; // Allow to use some tolerance during numerical factorization.
2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable MatrixType m_R; // The sparse matrix R in Eigen format
2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable Index *m_E; // The permutation applied to columns
2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable cholmod_sparse *m_H;  //The householder vectors
2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable Index *m_HPinv; // The row permutation of H
2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable cholmod_dense *m_HTau; // The Householder coefficients
2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable Index m_rank; // The rank of the matrix
2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    mutable cholmod_common m_cc; // Workspace and parameters
2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template<typename ,typename > friend struct SPQR_QProduct;
2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename SPQRType, typename Derived>
2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  typedef typename SPQRType::Scalar Scalar;
2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  typedef typename SPQRType::Index Index;
2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  //Define the constructor to get reference to argument types
2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  inline Index cols() const { return m_other.cols(); }
2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  // Assign to a vector
2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename ResType>
2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  void evalTo(ResType& res) const
2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    cholmod_dense y_cd;
2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    cholmod_dense *x_cd;
2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    int method = m_transpose ? SPQR_QTX : SPQR_QX;
2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    cholmod_common *cc = m_spqr.cholmodCommon();
2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    y_cd = viewAsCholmod(m_other.const_cast_derived());
2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    cholmod_l_free_dense(&x_cd, cc);
2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  const SPQRType& m_spqr;
2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  const Derived& m_other;
2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  bool m_transpose;
2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename SPQRType>
2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SPQRMatrixQReturnType{
2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename Derived>
2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  // To use for operations with the transpose of Q
2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  const SPQRType& m_spqr;
2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename SPQRType>
2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SPQRMatrixQTransposeReturnType{
2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename Derived>
2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  const SPQRType& m_spqr;
2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal {
2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, typename Rhs>
2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct solve_retval<SPQR<_MatrixType>, Rhs>
3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  : solve_retval_base<SPQR<_MatrixType>, Rhs>
3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  typedef SPQR<_MatrixType> Dec;
3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  template<typename Dest> void evalTo(Dest& dst) const
3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    dec()._solve(rhs(),dst);
3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal
3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}// End namespace Eigen
3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif
315