Extend support for decompositions on CUDA devices

Extend support for decompositions on CUDA devices
  - support QR decomposition
  - support LLT and LDLT Cholesky decomposition

Author: Viktor Khristenko

Eigen/src/Cholesky/LDLT.h

@@ -82,6 +82,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * The default constructor is useful in cases in which the user intends to
       * perform decompositions via LDLT::compute(const MatrixType&).
       */
+    EIGEN_DEVICE_FUNC
     LDLT()
       : m_matrix(),
         m_transpositions(),
@@ -95,6 +96,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * according to the specified problem \a size.
       * \sa LDLT()
       */
+    EIGEN_DEVICE_FUNC
     explicit LDLT(Index size)
       : m_matrix(size, size),
         m_transpositions(size),
@@ -110,6 +112,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * \sa LDLT(Index size)
       */
     template<typename InputType>
+    EIGEN_DEVICE_FUNC
     explicit LDLT(const EigenBase<InputType>& matrix)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_transpositions(matrix.rows()),
@@ -127,6 +130,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * \sa LDLT(const EigenBase&)
       */
     template<typename InputType>
+    EIGEN_DEVICE_FUNC
     explicit LDLT(EigenBase<InputType>& matrix)
       : m_matrix(matrix.derived()),
         m_transpositions(matrix.rows()),
@@ -140,19 +144,22 @@ template<typename _MatrixType, int _UpLo> class LDLT
     /** Clear any existing decomposition
      * \sa rankUpdate(w,sigma)
      */
+    EIGEN_DEVICE_FUNC
     void setZero()
     {
       m_isInitialized = false;
     }
 
     /** \returns a view of the upper triangular matrix U */
+    EIGEN_DEVICE_FUNC
     inline typename Traits::MatrixU matrixU() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
       return Traits::getU(m_matrix);
     }
 
     /** \returns a view of the lower triangular matrix L */
+    EIGEN_DEVICE_FUNC
     inline typename Traits::MatrixL matrixL() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
@@ -161,27 +168,31 @@ template<typename _MatrixType, int _UpLo> class LDLT
 
     /** \returns the permutation matrix P as a transposition sequence.
       */
+    EIGEN_DEVICE_FUNC
     inline const TranspositionType& transpositionsP() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
       return m_transpositions;
     }
 
     /** \returns the coefficients of the diagonal matrix D */
+    EIGEN_DEVICE_FUNC
     inline Diagonal<const MatrixType> vectorD() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
       return m_matrix.diagonal();
     }
 
     /** \returns true if the matrix is positive (semidefinite) */
+    EIGEN_DEVICE_FUNC
     inline bool isPositive() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
       return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
     }
 
     /** \returns true if the matrix is negative (semidefinite) */
+    EIGEN_DEVICE_FUNC
     inline bool isNegative(void) const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
@@ -205,45 +216,53 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt()
       */
     template<typename Rhs>
+    EIGEN_DEVICE_FUNC
     inline const Solve<LDLT, Rhs>
     solve(const MatrixBase<Rhs>& b) const;
     #endif
 
     template<typename Derived>
+    EIGEN_DEVICE_FUNC
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
 
     template<typename InputType>
+    EIGEN_DEVICE_FUNC
     LDLT& compute(const EigenBase<InputType>& matrix);
 
     /** \returns an estimate of the reciprocal condition number of the matrix of
      *  which \c *this is the LDLT decomposition.
      */
+    EIGEN_DEVICE_FUNC
     RealScalar rcond() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
       return internal::rcond_estimate_helper(m_l1_norm, *this);
     }
 
     template <typename Derived>
+    EIGEN_DEVICE_FUNC
     LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1);
 
     /** \returns the internal LDLT decomposition matrix
       *
       * TODO: document the storage layout
       */
+    EIGEN_DEVICE_FUNC
     inline const MatrixType& matrixLDLT() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
       return m_matrix;
     }
 
+    EIGEN_DEVICE_FUNC
     MatrixType reconstructedMatrix() const;
 
     /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint.
       *
       * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
       * \code x = decomposition.adjoint().solve(b) \endcode
       */
+    EIGEN_DEVICE_FUNC
     const LDLT& adjoint() const { return *this; };
 
     EIGEN_DEVICE_FUNC inline Index rows() const { return m_matrix.rows(); }
@@ -254,6 +273,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * \returns \c Success if computation was successful,
       *          \c NumericalIssue if the factorization failed because of a zero pivot.
       */
+    EIGEN_DEVICE_FUNC
     ComputationInfo info() const
     {
       eigen_assert(m_isInitialized && "LDLT is not initialized.");
@@ -262,6 +282,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
 
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename RhsType, typename DstType>
+    EIGEN_DEVICE_FUNC
     void _solve_impl(const RhsType &rhs, DstType &dst) const;
 
     template<bool Conjugate, typename RhsType, typename DstType>
@@ -270,6 +291,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
 
   protected:
 
+    EIGEN_DEVICE_FUNC
     static void check_template_parameters()
     {
       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
@@ -297,6 +319,7 @@ template<int UpLo> struct ldlt_inplace;
 template<> struct ldlt_inplace<Lower>
 {
   template<typename MatrixType, typename TranspositionType, typename Workspace>
+  EIGEN_DEVICE_FUNC
   static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
   {
     using std::abs;
@@ -333,7 +356,7 @@ template<> struct ldlt_inplace<Lower>
         Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
         mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
         mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
-        std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
+        numext::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
         for(Index i=k+1;i<index_of_biggest_in_corner;++i)
         {
           Scalar tmp = mat.coeffRef(i,k);
@@ -410,6 +433,7 @@ template<> struct ldlt_inplace<Lower>
   // Here only rank-1 updates are implemented, to reduce the
   // requirement for intermediate storage and improve accuracy
   template<typename MatrixType, typename WDerived>
+  EIGEN_DEVICE_FUNC
   static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1)
   {
     using numext::isfinite;
@@ -448,6 +472,7 @@ template<> struct ldlt_inplace<Lower>
   }
 
   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
+  EIGEN_DEVICE_FUNC
   static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
   {
     // Apply the permutation to the input w
@@ -460,13 +485,15 @@ template<> struct ldlt_inplace<Lower>
 template<> struct ldlt_inplace<Upper>
 {
   template<typename MatrixType, typename TranspositionType, typename Workspace>
+  EIGEN_DEVICE_FUNC
   static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
   {
     Transpose<MatrixType> matt(mat);
     return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
   }
 
   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
+  EIGEN_DEVICE_FUNC
   static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1)
   {
     Transpose<MatrixType> matt(mat);
@@ -478,15 +505,19 @@ template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
 {
   typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
+  EIGEN_DEVICE_FUNC
   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
+  EIGEN_DEVICE_FUNC
   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
 };
 
 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
 {
   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
   typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
+  EIGEN_DEVICE_FUNC
   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
+  EIGEN_DEVICE_FUNC
   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
 };
 
@@ -496,6 +527,7 @@ template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
   */
 template<typename MatrixType, int _UpLo>
 template<typename InputType>
+EIGEN_DEVICE_FUNC
 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
 {
   check_template_parameters();
@@ -536,6 +568,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputTyp
   */
 template<typename MatrixType, int _UpLo>
 template<typename Derived>
+EIGEN_DEVICE_FUNC
 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename LDLT<MatrixType,_UpLo>::RealScalar& sigma)
 {
   typedef typename TranspositionType::StorageIndex IndexType;
@@ -564,6 +597,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Deri
 #ifndef EIGEN_PARSED_BY_DOXYGEN
 template<typename _MatrixType, int _UpLo>
 template<typename RhsType, typename DstType>
+EIGEN_DEVICE_FUNC
 void LDLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
 {
   _solve_impl_transposed<true>(rhs, dst);
@@ -626,6 +660,7 @@ void LDLT<_MatrixType,_UpLo>::_solve_impl_transposed(const RhsType &rhs, DstType
   */
 template<typename MatrixType,int _UpLo>
 template<typename Derived>
+EIGEN_DEVICE_FUNC
 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
 {
   eigen_assert(m_isInitialized && "LDLT is not initialized.");
@@ -640,6 +675,7 @@ bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
  * i.e., it returns the product: P^T L D L^* P.
  * This function is provided for debug purpose. */
 template<typename MatrixType, int _UpLo>
+EIGEN_DEVICE_FUNC
 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
 {
   eigen_assert(m_isInitialized && "LDLT is not initialized.");
@@ -666,6 +702,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
   * \sa MatrixBase::ldlt()
   */
 template<typename MatrixType, unsigned int UpLo>
+EIGEN_DEVICE_FUNC
 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
 SelfAdjointView<MatrixType, UpLo>::ldlt() const
 {
@@ -677,6 +714,7 @@ SelfAdjointView<MatrixType, UpLo>::ldlt() const
   * \sa SelfAdjointView::ldlt()
   */
 template<typename Derived>
+EIGEN_DEVICE_FUNC
 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
 MatrixBase<Derived>::ldlt() const
 {