Re: [eigen] How to get triangular lower part of a matrix with dev branche ?

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Fun coincidence: the problem you were hitting was exaclty issue A2 in
the link that Hauke just sent to this list (see thread "very
interesting...").

You can't bind a ref-to-non-const to a temporary.

The solution in your case is to make that temporary a named variable.
Not any more inefficient than what you were trying to do.

See attached file.

Benoit

2010/5/20  <vincent.lejeune@xxxxxxxxxx>:
>
> Of course :
>
>
>
> #include <QR.h>
>
> using namespace Eigen;
>
>
>
> typedef  Matrix<double,4,4> MAT;
>
> typedef  Matrix<double,4,1> VECT;
>
>
>
>  template<typename MatrixType>
>
>  MatrixType triangfactor(MatrixType& V,
>
>      Matrix<typename MatrixType::Scalar,MatrixType::RowsAtCompileTime,1>&
>
> coeffs)
>
>  {
>
>    typedef typename MatrixType::Scalar Scalar;
>
>    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");
>
>
>
>
>
>    MatrixType ret=MatrixType::Zero(k,k);
>
>    const int j=k;
>
>
>
>    for(int i=0;i<k;i++)
>
>    {
>
>      Scalar Vii=V(i,i);
>
>      V(i,i)=1;
>
>      ret.block(0,i,i,1) = - coeffs(i,0) * V.block(i,0,j-i,i).transpose()
>
> * V.block(i,i,j-i,1);
>
>      V(i,i)=Vii;
>
>      ret.block(0,i,i,1)=ret.block(0,0,i,i)*ret.block(0,i,i,1);
>
>      ret(i,i)=coeffs(i,0);
>
>    }
>
>    return ret;
>
>  }
>
>
>
> int main()
>
> {
>
>  MAT mat=MAT::Random();
>
>  VECT vect=VECT::Random();
>
>  MAT tst=triangfactor(MAT(mat.triangularView<Lower>()),vect);
>
>  cout<<tst<<endl<<endl;
>
>
>
>  return 0;
>
> }
>
>
>
>
>
>
>
>
>
>
>
> On Thu, 20 May 2010 07:16:26 -0400, Benoit Jacob
>
> <jacob.benoit.1@xxxxxxxxx>
>
> wrote:
>
>> 2010/5/20  <vincent.lejeune@xxxxxxxxxx>:
>
>>>
>
>>> Hi,
>
>>>
>
>>>
>
>>>
>
>>> I'm trying to get the lower triangular part of a (non-square) matrix,
>
> but
>
>>>
>
>>> it seems that the part() function is deprecated in Eigen 2.91.
>
>>>
>
>>> Using mat.triangularView<Lower>(), I meet the following compile error :
>
>>
>
>> Sounds like  a bug. Can you send us a test program reproducing it?
>
>>
>
>> Benoit
>
>>
>
>>
>
>>>
>
>>>
>
>>>
>
>>> error: invalid initialization of non-const reference of type
>
>>>
>
>>> ‘Eigen::TriangularView<Eigen::Matrix<double, 4, 4, 0, 4, 4>, 1u>&’ from
>
>>> a
>
>>>
>
>>> temporary of type ‘Eigen::TriangularView<Eigen::Matrix<double, 4, 4, 0,
>
>>> 4,
>
>>>
>
>>> 4>, 1u>’
>
>>>
>
>>>
>
>>>
>
>>> That I do not really understand...
>
>>>
>
>>>
>
>>>
>
>>> Thx
>
>>>
>
>>>
>
>>>
>
>
>

#include<iostream>
#include <Eigen/QR>

using namespace std;
using namespace Eigen;



typedef  Matrix<double,4,4> MAT;

typedef  Matrix<double,4,1> VECT;



 template<typename MatrixType>

 MatrixType triangfactor(MatrixType& V,

     Matrix<typename MatrixType::Scalar,MatrixType::RowsAtCompileTime,1>&

coeffs)

 {

   typedef typename MatrixType::Scalar Scalar;

   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");





   MatrixType ret=MatrixType::Zero(k,k);

   const int j=k;



   for(int i=0;i<k;i++)

   {

     Scalar Vii=V(i,i);

     V(i,i)=1;

     ret.block(0,i,i,1) = - coeffs(i,0) * V.block(i,0,j-i,i).transpose()

* V.block(i,i,j-i,1);

     V(i,i)=Vii;

     ret.block(0,i,i,1)=ret.block(0,0,i,i)*ret.block(0,i,i,1);

     ret(i,i)=coeffs(i,0);

   }

   return ret;

 }



int main()

{

 MAT mat=MAT::Random();

 VECT vect=VECT::Random();

 MAT tmp(mat.triangularView<Lower>());
 MAT tst=triangfactor(tmp,vect);

 cout<<tst<<endl<<endl;



 return 0;

}