Re: [ublas] Element-wise operations are really slow - if done wrong

Subject: Re: [ublas] Element-wise operations are really slow - if done wrong
From: Gunter Winkler (guwi17_at_[hidden])
Date: 2010-04-25 15:23:49

Hello x,

Am Monday 19 April 2010 schrieb xiss burg:

>                 for(unsigned int r=0; r<3; ++r)
>                     for(unsigned int s=0; s<3; ++s)
>                     {
>                         m_RKR_1(3*jj+r, 3*kk+s) +=
> t->getCorotatedStiffness0(3*j+r, 3*k+s);
>                         m_RK(3*jj+r, 3*kk+s) +=
> t->getCorotatedStiffness1(3*j+r, 3*k+s);
>                     }
>             }

I think the biggest performance loss is that you access the matrix 
elements by index. This is the slowest possible access for sparse 
matrices. If you want best speed then you must use a coordinate matrix 
and the push_back operation.

Why? Because append_element(i,j,value) is a O(1) Operation für 
coordinate matrix. Additionally coordinate_matrix::append_element(int, 
int, T) behaves like the += operator - which is exactly what you need 
here. (This is a special feature of COO.) The only cost is an 
additional sort and compress step before you can (efficiently) use the 
matrix in computations. (This is done explicitly by a call to sort() or 
implicitly by trying to modify the matrix with any method except 
append_element).

I used the attached code in my FE program. The flag USE_INSERT is set 
for coordinate matrix.

My favorite matrix types are:

// type of stiffnes matrix
typedef boost::numeric::ublas::compressed_matrix<
  double, 
  boost::numeric::ublas::basic_row_major<unsigned int, int>,
  0,
  boost::numeric::ublas::unbounded_array<unsigned int>,
  boost::numeric::ublas::unbounded_array<double>
>  MY_STIMA;

typedef MY_STIMA MY_STIMA_B;

and for assembly

  // type of temporary stiffnes matrix
# ifdef USE_INSERT
  typedef boost::numeric::ublas::coordinate_matrix<
    double, 
    boost::numeric::ublas::row_major,
    0,
    boost::numeric::ublas::unbounded_array<std::size_t>,
    boost::numeric::ublas::unbounded_array<double>
    >  MY_STIMA_ASS;
# else
    typedef boost::numeric::ublas::generalized_vector_of_vector<
      double, 
      boost::numeric::ublas::row_major,
      boost::numeric::ublas::vector< 
boost::numeric::ublas::compressed_vector<double> >
    >  MY_STIMA_ASS;
# endif

HTH
Gunter