
Ublas : 
From: Michael Stevens (mail_at_[hidden])
Date: 20050814 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 nonsquare 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 nonzero 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 checkoff or
> checkon. 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 lvalue or an rvalue. 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 rvalues (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 nonconst 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 ___________________________________