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

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

Sorry, TriangularView doesn't allow that. You'll have to do this by hand.

Benoit


2010/5/21  <vincent.lejeune@xxxxxxxxxx>:
>
> Hi,
>
>
>
> I have another problem : I would like to get a block from a
>
> TriangularView<Matrixtype,UnitLower> matrix.
>
> However there is no block() member in a TriangularView. I would like to
>
> get a block below the diagonal, with only the right upper corner on the
>
> diagonal, without making any matrix copy just to change one coefficient....
>
>
>
> Thx,
>
> Vincent
>
>
>
> On Thu, 20 May 2010 09:51:44 -0400, Benoit Jacob
>
> <jacob.benoit.1@xxxxxxxxx>
>
> wrote:
>
>> 2010/5/20  <vincent.lejeune@xxxxxxxxxx>:
>
>>>
>
>>> I have another issue with triangularview :
>
>>>
>
>>>
>
>>>
>
>>> when compiling the following code :
>
>>>
>
>>>
>
>>>
>
>>>  template<typename Derived>
>
>>>
>
>>>  void applyBlockHouseholderOnTheRight(Derived& A,const Derived& V,
>
>>>
>
>>>      const Matrix<typename
>
>>> Derived::Scalar,Derived::RowsAtCompileTime,1>&
>
>>>
>
>>> coeffs)
>
>>>
>
>>>  {
>
>>>
>
>>>    typedef Matrix<typename
>
>>>
>
>>> Derived::Scalar,Derived::RowsAtCompileTime,Derived::RowsAtCompileTime>
>
>>>
>
>>> TriangularDerivedBlock;
>
>>>
>
>>>    const TriangularDerivedBlock&
>
>>>
>
>>> T=BlockHouseholderTriangularFactor(V,coeffs);
>
>>>
>
>>>    A -= A* V.triangularView<UnitLower>() * T.transpose() *
>
>>>
>
>>> V.triangularView<UnitLower>().transpose() ;
>
>>>
>
>>>  }
>
>>>
>
>>>
>
>>>
>
>>> int main()
>
>>>
>
>>> {
>
>>>
>
>>>  MAT A=MAT::Random();
>
>>>
>
>>>  MAT Abis(A);
>
>>>
>
>>>  MAT mat=MAT::Random();
>
>>>
>
>>>  VECT vect=VECT::Random();
>
>>>
>
>>>  applyBlockHouseholderOnTheRight(A,mat,vect);
>
>>>
>
>>>
>
>>>
>
>>>  return 0;
>
>>>
>
>>> }
>
>>>
>
>>>
>
>>>
>
>>> I get this (rather long) error message, telling me that gcc is trying
>
> to
>
>>>
>
>>> consider  V.triangularView<UnitLower>() as (V.triangularView) <
>
>>>
>
>>> (UnitLower>(...something...))  :
>
>>
>
>>
>
>> This is again a general c++ issue, not at all eigen specific: since
>
>> the type of V depends on template params, when you call a templated
>
>> method on it, you must write "template" like this:
>
>>
>
>>     V.template triangularView<UnitLower>()
>
>>
>
>> Benoit
>
>>
>
>>>
>
>>>
>
>>>
>
>>> /home/vlj/Sources/cudalinalg/eigen-fork/src/QR.h: In function ‘void
>
>>>
>
>>> GPU::applyBlockHouseholderOnTheRight(Derived&, const Derived&, const
>
>>>
>
>>> Eigen::Matrix<typename _MatrixType::Scalar,
>
>>> _MatrixType::RowsAtCompileTime,
>
>>>
>
>>> 1, ((Eigen::._89)0u | (((_MatrixType::RowsAtCompileTime == 1) && (1 !=
>
>>> 1))
>
>>>
>
>>> ? (Eigen::._89)1u : (((1 == 1) && (_MatrixType::RowsAtCompileTime !=
>
> 1))
>
>>> ?
>
>>>
>
>>> (Eigen::._89)0u : (Eigen::._89)0u))), _MatrixType::RowsAtCompileTime,
>
>>>
>
>>> 1>&)’:
>
>>>
>
>>> /home/vlj/Sources/cudalinalg/eigen-fork/src/QR.h:45: error: expected
>
>>>
>
>>> primary-expression before ‘)’ token
>
>>>
>
>>> /home/vlj/Sources/cudalinalg/eigen-fork/src/QR.h:45: error: expected
>
>>>
>
>>> primary-expression before ‘)’ token
>
>>>
>
>>> /home/vlj/Sources/cudalinalg/eigen-fork/src/QR.h: In function ‘void
>
>>>
>
>>> GPU::applyBlockHouseholderOnTheRight(Derived&, const Derived&, const
>
>>>
>
>>> Eigen::Matrix<typename _MatrixType::Scalar,
>
>>> _MatrixType::RowsAtCompileTime,
>
>>>
>
>>> 1, ((Eigen::._89)0u | (((_MatrixType::RowsAtCompileTime == 1) && (1 !=
>
>>> 1))
>
>>>
>
>>> ? (Eigen::._89)1u : (((1 == 1) && (_MatrixType::RowsAtCompileTime !=
>
> 1))
>
>>> ?
>
>>>
>
>>> (Eigen::._89)0u : (Eigen::._89)0u))), _MatrixType::RowsAtCompileTime,
>
>>> 1>&)
>
>>>
>
>>> [with Derived = MAT]’:
>
>>>
>
>>> /home/vlj/Sources/cudalinalg/eigen-fork/tests/blockhouseholder.cpp:44:
>
>>>
>
>>> instantiated from here
>
>>>
>
>>> /home/vlj/Sources/cudalinalg/eigen-fork/src/QR.h:45: error: no match
>
> for
>
>>>
>
>>> ‘operator*’ in ‘A * V->Eigen::MatrixBase::triangularView [with unsigned
>
>>> int
>
>>>
>
>>> Mode = Mode, Derived = Eigen::Matrix<double, 4, 4, 0, 4, 4>]’
>
>>>
>
>>>
>
> /home/vlj/Sources/cudalinalg/eigen-fork/eigen-ref/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:64:
>
>>>
>
>>> note: candidates are: const
>
>>>
>
>>> Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar>, Derived>
>
>>>
>
>>> Eigen::MatrixBase<Derived>::operator*(const typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar&) const [with Derived =
>
>>>
>
>>> Eigen::Matrix<double, 4, 4, 0, 4, 4>]
>
>>>
>
>>>
>
> /home/vlj/Sources/cudalinalg/eigen-fork/eigen-ref/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:84:
>
>>>
>
>>> note:                 const
>
>>>
>
>>> Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple2_op<typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar, std::complex<typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar> >, Derived>
>
>>>
>
>>> Eigen::MatrixBase<Derived>::operator*(const std::complex<typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar>&) const [with Derived =
>
>>>
>
>>> Eigen::Matrix<double, 4, 4, 0, 4, 4>]
>
>>>
>
>>>
>
> /home/vlj/Sources/cudalinalg/eigen-fork/eigen-ref/Eigen/src/Geometry/Scaling.h:119:
>
>>>
>
>>> note:                 const
>
>>>
>
>>> Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar>, Derived>
>
>>>
>
>>> Eigen::MatrixBase<Derived>::operator*(const
>
>>> Eigen::UniformScaling<typename
>
>>>
>
>>> Eigen::ei_traits<Derived>::Scalar>&) const [with Derived =
>
>>>
>
>>> Eigen::Matrix<double, 4, 4, 0, 4, 4>]
>
>>>
>
>>> /usr/include/c++/4.4/complex:397: note:
>
>>> std::complex<_Tp>
>
>>>
>
>>> std::operator*(const _Tp&, const std::complex<_Tp>&) [with _Tp =
>
>>>
>
>>> Eigen::Matrix<double, 4, 4, 0, 4, 4>]
>
>>>
>
>>>
>
> /home/vlj/Sources/cudalinalg/eigen-fork/eigen-ref/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:95:
>
>>>
>
>>> note:                 const
>
>>>
>
>>> Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple2_op<double,
>
>>>
>
>>> std::complex<double> >, Eigen::Matrix<double, 4, 4, 0, 4, 4> >
>
>>>
>
>>> Eigen::operator*(const std::complex<double>&, const
>
>>>
>
>>> Eigen::MatrixBase<Eigen::Matrix<double, 4, 4, 0, 4, 4> >&)
>
>>>
>
>>>
>
> /home/vlj/Sources/cudalinalg/eigen-fork/eigen-ref/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:91:
>
>>>
>
>>> note:                 const
>
>>>
>
>>> Eigen::CwiseUnaryOp<Eigen::ei_scalar_multiple_op<double>,
>
>>>
>
>>> Eigen::Matrix<double, 4, 4, 0, 4, 4> > Eigen::operator*(const double&,
>
>>>
>
>>> const Eigen::MatrixBase<Eigen::Matrix<double, 4, 4, 0, 4, 4> >&)
>
>>>
>
>>>
>
>>>
>
>>>
>
>>>
>
>>> I'm a bit stuck in front of this compilation error...
>
>>>
>
>>>
>
>>>
>
>>>
>
>>>
>
>>> On Thu, 20 May 2010 08:08:51 -0400, Benoit Jacob
>
>>>
>
>>> <jacob.benoit.1@xxxxxxxxx>
>
>>>
>
>>> wrote:
>
>>>
>
>>>> 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
>
>>>
>
>>>>>
>
>>>
>
>>>>>>>
>
>>>
>
>>>>>
>
>>>
>
>>>>>>>
>
>>>
>
>>>>>
>
>>>
>
>>>>>>>
>
>>>
>
>>>>>
>
>>>
>
>>>>>
>
>>>
>
>>>>>
>
>>>
>
>>>
>
>>>
>
>
>