> > http://www.bauv.unibw-muenchen.de/~winkler/ublas/matrix_assign_traits.htm

This should map to a matrix_expression<const_self_type> which makes sense.

> >l created by tidy and the perl script:

> > http://www.bauv.unibw-muenchen.de/~winkler/ublas/show_traits.pl

>

> Oh this is really useful. This code is horribly hard to understand. On

> additional thing that often gets me into trouble is that the storage_tags

> are in fact derived one from the other. They are so ordered that the most

> general algorithm will match if there is nothing more specific. This caused

> the recent bug as the sparse_proxy_tag was able match as I has messed up

> the definition of the sparse_tag version.

>

> > The question is: Should uBlas allow such an assignement? I can not

> > imagine a usesul application and it would be very easy to drop the

> > function symmetric_matrix::operator= (matrix_expression ...).

>

> It would be nice to still support things like

> S = S * 1.2;

> without a matrix_expression assignment would loose a lot of uBLAS nice

> syntax.

> > IMHO it would be much better to only allow write access to elements that

But what should row(S,k) = x; mean? Does it assign x to row k and column k or

> > are really stored: s(1,0) = 1.0 is ok, but s(0,1) = 1.0 gives a runtime

> > error. The assignement of a symmetric matrix would look like this:

>

> uBLAS used to have this as the default setting. But it turned out to be

> very bad for two reasons.

> a) Many algorithms are hard to write if they have to precheck the indexes

> to avoid these errors. A simple example is

> row(S,1) = vector;

> If we know S is symmetric we can split the expression. But this is slow. If

> we are writting generic code then expression will fail for some types but

> not others.

could it mean to assign only part of x to row k and column k. The docs does

not mention anything about these features. If we virtually mirror the storage

we have to write a precise documentation!

IMHO, If I write symmetric algorithms I should take care of the symmetries

myself. The above example could look like:

// if I have to take care of the structure

symmetric_matrix<double, lower> S(size,size);

row (triangular_adaptor<lower>(S), k) = project( x, range(1, k) );

column (triangular_adaptor<lower>(S), k) = project( x, range(k, size) );

I do not think that it is slower then checking on every element if it is in

the upper or lower part of S.

// if writes to mirrored elements are silently ignored I could write

symmetric_matrix<double, lower> S(size,size);

row (S, k) = x;

column (S, k) = x;

this one even works if the type of S is not symmetric, but I don't know if a

compiler can detect the needless operations. However, it is very clear what I

intend to do. (Note: this is more complicated for hermitian_matrices)

> b) C++ does not differentiate between an s(0,1) on the LHS or RHS of an

I see, C++ differentiates between const and non const references -

> expression. In either case the element access function simple returns a

> reference to the actual storage location. The only way to enforce this is

> to use the element proxies in the same way that sparse matrices use them.

> These are complex and have a runtime cost.

unfortunatly it only looks at the "constsness" of the intance. The expression

S(2,0) = S(0,1);

should call

const_ref operator()(0,1) const

and

ref operator()(2,0)

if I have operator=(const_ref rhs). However, since S is not const, it really

calls both times

ref operator()(...) ;

The only (ugly) workaround is:

S(2,0) = const_cast<const MATRIX&>(S)(0,1);

or

S(2,0) = make_const(S)(0,1);

While thinking about structured matrices I wondered if we really need the

symmetric_matrix class. The full functionality should still be possible with

symmetric_adaptor< triangular_matrix > and similar. This way we could drop to

large classes which are very similar. Moreover we could join

symmetric_adaptor and hermitian_adaptor by introducing a

symmetry_operation_functor:

typedef lower F;

symmetric_operation_functor t;

t(F::other(1, 0), 2.0) gives 2.0

t(F::other(0, 1), 2.0) gives 2.0

hermitian_operation_functor h;

h(F::other(1, 0), (2.0, 1.0)) gives (2.0, 1.0)

h(F::other(0, 1), (2.0, 1.0)) gives (2.0, -1.0)

which could be extended to skew symmetric/hermitian and more ...

> I would think that would be a very good idea. If you want to take on the

Ok. I am going to reactivate my SF account.

> initial work yourself that would be very good. I could help out proof

> reading and checking consistency with the code. If you would like me to

> give you access to the SF project I think this would be appropriate.

> > PS: How about some functions make_symmetric, make_hermitian, ... ?

I chose make_symmetric and make_symmetric_const as names, because the compiler

> Sounds good. uBLAS already has the row/column functions which are similar.

> We could follow this example and provide make_ for all the adaptors.

cannot distinguish a const and a mutable version, when I pass the result to a

(template-) function.

btw. the documentation of symmetric_matrix does only mention lower and upper

as type. What about unit_lower, unit_upper, strict_lower and strict_upper? At

first glance they should work, too.

I wrote a short test program for symmetric matrices, but it aborts at

make_hermitian<lower>(a)(0,1) = 6;

with SIGABORT. (0,1) is not stored but should be accessible.

mfg

Gunter

------8<-------------------------------

// undocumented define: #define BOOST_UBLAS_STRICT_HERMITIAN

#include <boost/numeric/ublas/matrix_sparse.hpp>

#include <boost/numeric/ublas/symmetric.hpp>

#include <boost/numeric/ublas/hermitian.hpp>

#include <boost/numeric/ublas/triangular.hpp>

#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric::ublas;

template<class MATRIX>

void print_const(const MATRIX & M)

{

for (size_t i = 0; i < M.size1(); ++i) {

matrix_row<const MATRIX> r(M, i);

typename matrix_row<const MATRIX>::const_iterator ri=r.begin();

typename matrix_row<const MATRIX>::const_iterator rn=r.end();

size_t n = 0;

for (; ri != rn; ++ri, ++n) {

std::cout << "(" << i << "," << ri.index()

<< ":" << *ri << ") ";

}

if (n>0) { std::cout << "\n"; }

}

std::cout << std::endl;

return;

}

template <class TYPE, class MATRIX>

const symmetric_adaptor<const MATRIX, TYPE> make_symmetric_const(const MATRIX

& mat) {

return symmetric_adaptor<const MATRIX, TYPE>(mat);

}

template <class TYPE, class MATRIX>

symmetric_adaptor<MATRIX, TYPE> make_symmetric(MATRIX & mat) {

return symmetric_adaptor<MATRIX, TYPE>(mat);

}

template <class TYPE, class MATRIX>

const hermitian_adaptor<const MATRIX, TYPE> make_hermitian_const(const MATRIX

& mat) {

return hermitian_adaptor<const MATRIX, TYPE>(mat);

}

template <class TYPE, class MATRIX>

hermitian_adaptor<MATRIX, TYPE> make_hermitian(MATRIX & mat) {

return hermitian_adaptor<MATRIX, TYPE>(mat);

}

template <class TYPE, class MATRIX>

const triangular_adaptor<const MATRIX, TYPE> make_triangular_const(const

MATRIX & mat) {

return triangular_adaptor<const MATRIX, TYPE>(mat);

}

template <class TYPE, class MATRIX>

triangular_adaptor<MATRIX, TYPE> make_triangular(MATRIX & mat) {

return triangular_adaptor<MATRIX, TYPE>(mat);

}

template <class TYPE, class MATRIX>

void test_make(const MATRIX& m)

{

{

std::cout << "symmetric:\n";

print_const( make_symmetric_const<TYPE>(m) );

MATRIX a(m);

make_symmetric<TYPE>(a)(1,1) = 5;

make_symmetric<TYPE>(a)(0,1) = 6;

make_symmetric<TYPE>(a)(1,0) = 7;

print_const( a );

}

{

std::cout << "hermitian:\n";

print_const( make_hermitian_const<TYPE>(m) );

MATRIX a(m);

make_hermitian<TYPE>(a)(1,1) = 5;

make_hermitian<TYPE>(a)(0,1) = 6;

make_hermitian<TYPE>(a)(1,0) = 7;

print_const( a );

}

{

std::cout << "triangular:\n";

print_const( make_hermitian_const<TYPE>(m) );

MATRIX a(m);

make_triangular<TYPE>(a)(1,1) = 5;

make_triangular<TYPE>(a)(0,1) = 6;

make_triangular<TYPE>(a)(1,0) = 7;

print_const( a );

}

}

int main (int argc, char * argv[]) {

int size = 5;

// if (argc > 1)

// size = ::atoi (argv [1]);

sparse_matrix<double> src(size,size,2*size);

src(0,0) = 10;

src(0,2) = 2;

src(2,0) = -2;

src(size-1,size-1) = double(size);

std::cout << "lower:\n";

test_make<lower>(src);

std::cout << "unit_lower:\n";

test_make<unit_lower>(src);

std::cout << "strict_lower:\n";

test_make<strict_lower>(src);

std::cout << "upper:\n";

test_make<upper>(src);

std::cout << "unit_upper:\n";

test_make<unit_upper>(src);

std::cout << "strict_upper:\n";

test_make<strict_upper>(src);

return 0 ;

}