Re: [eigen] Blocked QR algorithm - lapack compatible ? |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: <eigen@xxxxxxxxxxxxxxxxxxx>*Subject*: Re: [eigen] Blocked QR algorithm - lapack compatible ?*From*: <vincent.lejeune@xxxxxxxxxx>*Date*: Mon, 07 Jun 2010 17:19:44 +0200

Hi, so there is some error in my code, but I don't see where... I wished to write a blocked-Householder applying function like DLARFB in Lapack. This routine compute C:=C-V*T*V'*C, where the i-th column of V store the vector that defines the i-th householder factor. Indeed H = H1 * H2 * H3... * Hn can be written H = I - V * T * V', where T is a triangular factor, whose elements depends on V and on taus factor. T is computed by DLARFT in Lapack. I've attached a header file that contain my implemantation of such functions. It seems that they return the same results than Lapack routines, but the result of a QR decomposition is not the same than Eigen's one. Maybe you could give it an eye and find where is the issue...Thanks Vincent On Wed, 2 Jun 2010 11:33:13 +0200, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote: > sorry but our householder factor should really be the same than Lapack. Our > makeHouseholder computes the same thing than the Lapack dlarfg routine. > > gael > > On Fri, May 21, 2010 at 6:30 PM, Benoit Jacob > <jacob.benoit.1@xxxxxxxxx>wrote: > >> True, our hCoeffs aren't quite the same as LAPACK's tau. I don't >> remember exactly but I'm sure Gael does, as he coded that, so, pinging >> him. >> >> Benoit >> >> 2010/5/21 <vincent.lejeune@xxxxxxxxxx>: >> > >> > Hi, >> > >> > >> > >> > I'm trying to implement a blocked QR algorithm, using the code from >> > >> > HouseholderQR. >> > >> > I'm using the algorithm used in lapack routine dgeqrf : first I make >> > the >> > >> > reduction of a panel of the input matrix, then I build the "T" >> > triangular >> > >> > factor in H=I-VT'V' that defines the blocked householder >> > transformation, >> > >> > then I apply this transformation to the trailing submatrix by a left >> > >> > multiplication. Then I iterate the process. >> > >> > >> > >> > I think that the current implemantation of householderQR does not >> > comply >> > >> > with lapack routine. The coefficients in hCoeffs() does actually not >> > have >> > >> > the same meaning than the one in the TAU vector found in dgeqrf. When >> > >> > applying dlarft and dlarfb with hCoeffs to a partially reduced matrix, >> > I >> > >> > end with mostly different results than the one furnished by a complete >> > >> > eigen householder decomposition (strangly the first row is correct). >> > >> > >> > >> > I suspect that the applyhousolderfrom the left use the coefficient in a >> > >> > different war than dlarf in lapack does. >> > >> > >> > >> > Is there some doc on how hcoeffs from eigen and tau from lapack are >> > >> > related ? >> > >> > >> > >> > Thanks, Vincent >> > >> > >> > >>

#include <GPUmat.h> #include <Eigen/Eigen> using namespace Eigen; #define NB 8 namespace GPU { template<typename Derived, typename OtherDerived = Matrix< typename Derived::Scalar, Derived::ColsAtCompileTime, Dynamic> > class HHutilities { public: typedef Matrix<typename Derived::Scalar, Derived::ColsAtCompileTime, Derived::ColsAtCompileTime> TriangularFactorType; typedef typename ei_plain_diag_type<Derived>::type CoeffVect; TriangularFactorType static BlockHouseholderTriangularFactor(const Derived& V, const CoeffVect& coeffs) { const int k = coeffs.size(); ei_assert(V.cols()==k && V.rows()>=k && "Vector storing Matrix must have same columns than coeffs size and more rows than columns"); TriangularFactorType ret = TriangularFactorType::Zero(k, k); const int j = k; for (int i = 0; i < k; i++) { ret.block(0, i, i, 1) = -coeffs(i) * V.block(i, 0, j - i, i).transpose() * V.block(i, i, j - i, 1); ret.block(0, i, i, 1) = ret.block(0, 0, i, i) * ret.block(0, i, i, 1); ret(i, i) = coeffs(i); } return ret; } static void applyBlockHouseholderOnTheLeft(OtherDerived& A, const Derived& V, const CoeffVect& coeffs) { const Derived& Vtmp(V.template triangularView<UnitLower> ()); const TriangularFactorType& T = BlockHouseholderTriangularFactor(Vtmp, coeffs); cout<<T<<endl<<endl; A -= Vtmp * T.transpose() * Vtmp.transpose() * A; } static void applyBlockHouseholderOnTheRight(OtherDerived& A, const Derived& V, const CoeffVect& coeffs) { const Derived& Vtmp(V.template triangularView<UnitLower> ()); const TriangularFactorType& T = BlockHouseholderTriangularFactor(Vtmp, coeffs); A -= A * Vtmp * T * Vtmp.transpose(); } }; template<typename MatrixType> class HouseholderQR2 : public HouseholderQR<MatrixType> { public: typedef HouseholderQR<MatrixType> Base; typedef Block<MatrixType, Dynamic, Dynamic> BlockType; typedef typename MatrixType::Scalar RealScalar; typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrix; typedef typename HHutilities<DynamicMatrix,DynamicMatrix>::TriangularFactorType TriangularFactorType; typedef typename HHutilities<DynamicMatrix,DynamicMatrix>::CoeffVect HCoeffsType; void unblockedCompute(Block<MatrixType>& m_qr, HCoeffsType& m_hCoeffs) { int rows = m_qr.rows(); int cols = m_qr.cols(); int size = std::min(rows, cols); m_hCoeffs.resize(size); //Base::m_temp.resize(cols); for (int k = 0; k < size; ++k) { int remainingRows = rows - k; int remainingCols = cols - k - 1; RealScalar beta; m_qr.col(k).tail(remainingRows).makeHouseholderInPlace( m_hCoeffs.coeffRef(k), beta); m_qr.coeffRef(k, k) = beta; // apply H to remaining part of m_qr from the left m_qr.bottomRightCorner(remainingRows, remainingCols) .applyHouseholderOnTheLeft( m_qr.col(k).tail(remainingRows - 1), m_hCoeffs.coeffRef(k), &(Base::m_temp.coeffRef(k + 1))); } } HouseholderQR2<MatrixType>& compute(const MatrixType& matrix) { int rows = matrix.rows(); int cols = matrix.cols(); int size = std::min(rows, cols); HCoeffsType tmpcoeffs; Base::m_qr = matrix; int BLOCKSIZE = 2; for (int k = 0; k < size; k += BLOCKSIZE) { int remainingRows = rows - k; int remainingCols = cols - k; BlockType tmp=Base::m_qr.block(k,k,rows-k,BLOCKSIZE); unblockedCompute(tmp,tmpcoeffs); tmpcoeffs(1)/=tmpcoeffs(0); DynamicMatrix dtmp(tmp); DynamicMatrix dtmp2(Base::m_qr.block(k,k+BLOCKSIZE,rows-k,cols-k-BLOCKSIZE)); HHutilities<DynamicMatrix,DynamicMatrix>::applyBlockHouseholderOnTheLeft(dtmp2,dtmp,tmpcoeffs); Base::m_qr.block(k,k+BLOCKSIZE,rows-k,cols-k-BLOCKSIZE)=dtmp2; cout<<tmpcoeffs<<endl<<endl; } Base::m_isInitialized = true; return *this; } };

**Follow-Ups**:**Re: [eigen] Blocked QR algorithm - lapack compatible ?***From:*Benoit Jacob

**References**:**Re: [eigen] Blocked QR algorithm - lapack compatible ?***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] runtime API to configure the cache/blocking parameters** - Next by Date:
**Re: [eigen] Blocked QR algorithm - lapack compatible ?** - Previous by thread:
**Re: [eigen] Blocked QR algorithm - lapack compatible ?** - Next by thread:
**Re: [eigen] Blocked QR algorithm - lapack compatible ?**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |