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

```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;
}

};
```

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