Eigen  3.3.3
SuiteSparseQRSupport.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
00005 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SUITESPARSEQRSUPPORT_H
00012 #define EIGEN_SUITESPARSEQRSUPPORT_H
00013 
00014 namespace Eigen {
00015   
00016   template<typename MatrixType> class SPQR; 
00017   template<typename SPQRType> struct SPQRMatrixQReturnType; 
00018   template<typename SPQRType> struct SPQRMatrixQTransposeReturnType; 
00019   template <typename SPQRType, typename Derived> struct SPQR_QProduct;
00020   namespace internal {
00021     template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> >
00022     {
00023       typedef typename SPQRType::MatrixType ReturnType;
00024     };
00025     template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
00026     {
00027       typedef typename SPQRType::MatrixType ReturnType;
00028     };
00029     template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> >
00030     {
00031       typedef typename Derived::PlainObject ReturnType;
00032     };
00033   } // End namespace internal
00034   
00059 template<typename _MatrixType>
00060 class SPQR : public SparseSolverBase<SPQR<_MatrixType> >
00061 {
00062   protected:
00063     typedef SparseSolverBase<SPQR<_MatrixType> > Base;
00064     using Base::m_isInitialized;
00065   public:
00066     typedef typename _MatrixType::Scalar Scalar;
00067     typedef typename _MatrixType::RealScalar RealScalar;
00068     typedef SuiteSparse_long StorageIndex ;
00069     typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType;
00070     typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType;
00071     enum {
00072       ColsAtCompileTime = Dynamic,
00073       MaxColsAtCompileTime = Dynamic
00074     };
00075   public:
00076     SPQR() 
00077       : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
00078     { 
00079       cholmod_l_start(&m_cc);
00080     }
00081     
00082     explicit SPQR(const _MatrixType& matrix)
00083     : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
00084     {
00085       cholmod_l_start(&m_cc);
00086       compute(matrix);
00087     }
00088     
00089     ~SPQR()
00090     {
00091       SPQR_free();
00092       cholmod_l_finish(&m_cc);
00093     }
00094     void SPQR_free()
00095     {
00096       cholmod_l_free_sparse(&m_H, &m_cc);
00097       cholmod_l_free_sparse(&m_cR, &m_cc);
00098       cholmod_l_free_dense(&m_HTau, &m_cc);
00099       std::free(m_E);
00100       std::free(m_HPinv);
00101     }
00102 
00103     void compute(const _MatrixType& matrix)
00104     {
00105       if(m_isInitialized) SPQR_free();
00106 
00107       MatrixType mat(matrix);
00108       
00109       /* Compute the default threshold as in MatLab, see:
00110        * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
00111        * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 
00112        */
00113       RealScalar pivotThreshold = m_tolerance;
00114       if(m_useDefaultThreshold) 
00115       {
00116         RealScalar max2Norm = 0.0;
00117         for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
00118         if(max2Norm==RealScalar(0))
00119           max2Norm = RealScalar(1);
00120         pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
00121       }
00122       cholmod_sparse A; 
00123       A = viewAsCholmod(mat);
00124       m_rows = matrix.rows();
00125       Index col = matrix.cols();
00126       m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A, 
00127                              &m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
00128 
00129       if (!m_cR)
00130       {
00131         m_info = NumericalIssue;
00132         m_isInitialized = false;
00133         return;
00134       }
00135       m_info = Success;
00136       m_isInitialized = true;
00137       m_isRUpToDate = false;
00138     }
00142     inline Index rows() const {return m_rows; }
00143     
00147     inline Index cols() const { return m_cR->ncol; }
00148     
00149     template<typename Rhs, typename Dest>
00150     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00151     {
00152       eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
00153       eigen_assert(b.cols()==1 && "This method is for vectors only");
00154 
00155       //Compute Q^T * b
00156       typename Dest::PlainObject y, y2;
00157       y = matrixQ().transpose() * b;
00158       
00159       // Solves with the triangular matrix R
00160       Index rk = this->rank();
00161       y2 = y;
00162       y.resize((std::max)(cols(),Index(y.rows())),y.cols());
00163       y.topRows(rk) = this->matrixR().topLeftCorner(rk, rk).template triangularView<Upper>().solve(y2.topRows(rk));
00164 
00165       // Apply the column permutation 
00166       // colsPermutation() performs a copy of the permutation,
00167       // so let's apply it manually:
00168       for(Index i = 0; i < rk; ++i) dest.row(m_E[i]) = y.row(i);
00169       for(Index i = rk; i < cols(); ++i) dest.row(m_E[i]).setZero();
00170       
00171 //       y.bottomRows(y.rows()-rk).setZero();
00172 //       dest = colsPermutation() * y.topRows(cols());
00173       
00174       m_info = Success;
00175     }
00176     
00179     const MatrixType matrixR() const
00180     {
00181       eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
00182       if(!m_isRUpToDate) {
00183         m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
00184         m_isRUpToDate = true;
00185       }
00186       return m_R;
00187     }
00189     SPQRMatrixQReturnType<SPQR> matrixQ() const
00190     {
00191       return SPQRMatrixQReturnType<SPQR>(*this);
00192     }
00194     PermutationType colsPermutation() const
00195     { 
00196       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00197       return PermutationType(m_E, m_cR->ncol);
00198     }
00203     Index rank() const
00204     {
00205       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00206       return m_cc.SPQR_istat[4];
00207     }
00209     void setSPQROrdering(int ord) { m_ordering = ord;}
00211     void setPivotThreshold(const RealScalar& tol)
00212     {
00213       m_useDefaultThreshold = false;
00214       m_tolerance = tol;
00215     }
00216     
00218     cholmod_common *cholmodCommon() const { return &m_cc; }
00219     
00220     
00226     ComputationInfo info() const
00227     {
00228       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00229       return m_info;
00230     }
00231   protected:
00232     bool m_analysisIsOk;
00233     bool m_factorizationIsOk;
00234     mutable bool m_isRUpToDate;
00235     mutable ComputationInfo m_info;
00236     int m_ordering; // Ordering method to use, see SPQR's manual
00237     int m_allow_tol; // Allow to use some tolerance during numerical factorization.
00238     RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
00239     mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
00240     mutable MatrixType m_R; // The sparse matrix R in Eigen format
00241     mutable StorageIndex *m_E; // The permutation applied to columns
00242     mutable cholmod_sparse *m_H;  //The householder vectors
00243     mutable StorageIndex *m_HPinv; // The row permutation of H
00244     mutable cholmod_dense *m_HTau; // The Householder coefficients
00245     mutable Index m_rank; // The rank of the matrix
00246     mutable cholmod_common m_cc; // Workspace and parameters
00247     bool m_useDefaultThreshold;     // Use default threshold
00248     Index m_rows;
00249     template<typename ,typename > friend struct SPQR_QProduct;
00250 };
00251 
00252 template <typename SPQRType, typename Derived>
00253 struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
00254 {
00255   typedef typename SPQRType::Scalar Scalar;
00256   typedef typename SPQRType::StorageIndex StorageIndex;
00257   //Define the constructor to get reference to argument types
00258   SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
00259   
00260   inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); }
00261   inline Index cols() const { return m_other.cols(); }
00262   // Assign to a vector
00263   template<typename ResType>
00264   void evalTo(ResType& res) const
00265   {
00266     cholmod_dense y_cd;
00267     cholmod_dense *x_cd; 
00268     int method = m_transpose ? SPQR_QTX : SPQR_QX; 
00269     cholmod_common *cc = m_spqr.cholmodCommon();
00270     y_cd = viewAsCholmod(m_other.const_cast_derived());
00271     x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc);
00272     res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
00273     cholmod_l_free_dense(&x_cd, cc);
00274   }
00275   const SPQRType& m_spqr; 
00276   const Derived& m_other; 
00277   bool m_transpose; 
00278   
00279 };
00280 template<typename SPQRType>
00281 struct SPQRMatrixQReturnType{
00282   
00283   SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
00284   template<typename Derived>
00285   SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
00286   {
00287     return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
00288   }
00289   SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const
00290   {
00291     return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
00292   }
00293   // To use for operations with the transpose of Q
00294   SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
00295   {
00296     return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
00297   }
00298   const SPQRType& m_spqr;
00299 };
00300 
00301 template<typename SPQRType>
00302 struct SPQRMatrixQTransposeReturnType{
00303   SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
00304   template<typename Derived>
00305   SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
00306   {
00307     return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
00308   }
00309   const SPQRType& m_spqr;
00310 };
00311 
00312 }// End namespace Eigen
00313 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends