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

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


Hi,

Indeed, the way lapack and eigen stores the householder transformations is different.

Here you can see the code i use to go from one to the other:

http://bitbucket.org/eigen/eigen/src/tip/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h#cl-470

only those lines are important:

------------------------------------------------------------------------

ColPivHouseholderQR<JacobianType> qrfac(fjac);

        fjac = qrfac.matrixQR();
        wa1 = fjac.diagonal();
        fjac.diagonal() = qrfac.hCoeffs();
        permutation = qrfac.colsPermutation();
        // TODO : avoid this:
        for(Index ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors
------------------------------------------------------------------------

I use the eigen qr decompozition, but the code later one is raw lapack and relies on the lapack 'packing'.

In both cases, the packing is done within the (strict) lower triangular part of a matrix, a vector stored on the diagonal, and another vector stored "outside" the matrix.

There are two differences, though:

* 'diagonal' and 'external vector' are inversed between eigen and lapack. That is, the lapack diagonal is what eigen calls hCoeffs, and the eigen diagonal is what lapack calls 'wa1'.

* vectors in the lower triangular parts are 'scaled' (check my code) in one implementation and not in the other one.

Hope this helps!

Thomas

--

Thomas Capricelli <orzel@xxxxxxxxxxxxxxx>

http://www.freehackers.org/thomas

In data lunedì 07 giugno 2010 17:35:19, Benoit Jacob ha scritto:

> Actually, pinging Thomas, because I remember he faced the exact same

> issue. The conversion is not hard to do at all, but there is something

> to do.

>

> Benoit

>

> 2010/6/7 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:

> > I'm pretty sure that out householder stuff isn't LAPACK compatible:

> > either our hCoeffs or our housholder vectors are different. Don't

> > remember which. Since Gael says the hCoeffs are compatible, I guess

> > it's the vectors then, that are different.



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