| Re: [eigen] Blocked QR algorithm - lapack compatible ? |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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;
}
};