Loading ...
Sorry, an error occurred while loading the content.

Re: [ublas-dev] matrix_assign_traits & symmetry

Expand Messages
  • Gunter Winkler
    ... This should map to a matrix_expression which makes sense. ... But what should row(S,k) = x; mean? Does it assign x to row k and column k
    Message 1 of 12 , Jul 1, 2004
    • 0 Attachment
      On Wednesday 30 June 2004 20:36, Michael Stevens wrote:
      > > http://www.bauv.unibw-muenchen.de/~winkler/ublas/matrix_assign_traits.htm
      > >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.

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

      > > IMHO it would be much better to only allow write access to elements that
      > > 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.

      But what should row(S,k) = x; mean? Does it assign x to row k and column k or
      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
      > 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.

      I see, C++ differentiates between const and non const references -
      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
      > 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.

      Ok. I am going to reactivate my SF account.

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

      > Sounds good. uBLAS already has the row/column functions which are similar.
      > We could follow this example and provide make_ for all the adaptors.

      I chose make_symmetric and make_symmetric_const as names, because the compiler
      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 ;
      }
    Your message has been successfully submitted and would be delivered to recipients shortly.