|
|
Ublas : |
From: Michael Stevens (mail_at_[hidden])
Date: 2005-08-14 11:50:26
Hi again Darin,
On Friday 12 August 2005 01:03, Darin DeForest wrote:
> I have some issues with triangular matrices.
>
> First issue seems that the triangle matrix allocates too much memory. Right
> now I non-square triangular matrices, where size1() != size2();
>
> size is being computed by this method:
>
> static
> BOOST_UBLAS_INLINE
> size_type triangular_size (size_type size1, size_type size2) {
> size_type size = (std::max) (size1, size2);
> // Guard against size_type overflow - simplified
> BOOST_UBLAS_CHECK (size == 0 || size / 2 <
> (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
> return ((size + 1) * size) / 2;
> }
>
> so for a 2x5 triangle matrix, this will take 6*5/2 = 15 elements. This is
> more than a regular 2x5 matrix which should be 10 elements large.
It to me the 'max' should simply be replaced by a 'min'. We only need storage
for the triangular element which may be non zero.
> Another issue is with assignment, I have
>
> matrix<float> m( M, N);
> triangular_matrix<float, lower> lower (M, N);
>
> // fill m with some values
>
> lower = m;
>
> In method:
>
> // Packed (proxy) row major case
> template<template <class T1, class T2> class F, class R, class M, class
> E> // BOOST_UBLAS_INLINE This function seems to be big. So we do not let
> the compiler inline it. void matrix_assign (M &m, const
> matrix_expression<E> &e, packed_proxy_tag, row_major_tag) {
>
> I get an error from this test failing:
>
> #if BOOST_UBLAS_TYPE_CHECK
> if (! disable_type_check<bool>::value)
> BOOST_UBLAS_CHECK (detail::expression_type_check (m, cm),
> external_logic ()); #endif
The check will be correctly triggered if m is not numerically triangular. That
is if 'm' has non-zero elements where 'lower' must be zero then uBLAS tells
you so.
> when I proceed the above assignment statement with
>
> disable_type_check<bool>::value = true;
>
> then the assignment works without an exception being generated, and some
> hand verification shows that the assignment worked.
>
> The question is that the check seems useful, after all it is in the code.
> However I don't like surrounding assignment statements with check-off or
> check-on. The other option is to disable the BOOST_UBLAS_TYPE_CHECK define.
>
> So given that, what is the best choice that I can make?
Neither. You should make the RHS conform to the numeric 'type' requirements of
the LHS. You can do this with a triangular adaptor.
>
> The last issue is a problem with indexing into the triangular matrix,
> because in some cases an exception can be raised as the following
> triangular_matrix method illustrates:
>
> BOOST_UBLAS_INLINE
> reference operator () (size_type i, size_type j) {
> BOOST_UBLAS_CHECK (i < size1_, bad_index ());
> BOOST_UBLAS_CHECK (j < size2_, bad_index ());
> if (triangular_type::other (i, j))
> return data () [triangular_type::element (layout_type (),
> i, size1_, j, size2_)]; else {
> bad_index ().raise ();
> // arbitary return value
> return const_cast<reference>(zero_);
> }
> }
>
> In this case, I would expect the rest of the matrix to have some default
> value, like zero. This makes the caller if using this api to be aware that
> this matrix has a different accessing strategy.
This is a problem with C++. You cannot distinguish between operator() being
called where the 'reference' is used as an l-value or an r-value. So the
operator() must fail if an attempt is made to access one of the zero values.
If you wish to index elements of a triangular matrix as r-values (for printing
etc) then you must make sure you access the matrix as a constant. In that
case the 'operator () const' is used which returns the values 0 or 1 for
appropriate elements.
One way to avoid the above would be to have specific member functions which
access the elements either as l or r values. But this seems no better then
the const non-const differentiation.
> Also I would point out that
> this would also make it work similar to unit_vector, where that only
> contains one value for a specified index, and all of it's other indexes
> have zero.
unit_vector etc have all constant element so they don't suffer from the above
problems.
> Right I consider this to be a documentation defect, since the documentation
> doesn't indicate that this exception would be raised.
I think the documentation requires some explanation of how elements of packed
type behave when index. This can get quite complex as you have to take into
account the special element behaviour of symmetric and hermitian matrices!
All the best,
Michael
--
___________________________________
Michael Stevens Systems Engineering
34128 Kassel, Germany
Phone/Fax: +49 561 5218038
Navigation Systems, Estimation and
Bayesian Filtering
http://bayesclasses.sf.net
___________________________________