| Re: [eigen] How to get triangular lower part of a matrix with dev branche ? |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>
>>
>>>>
>>
>>
>>