Einstein summation using uBLAS/tensor

[gear-head]

Hi there,

I am currently trying to use the Tensor class to compute tensor expressions similar to NumPy's einsum function. I get this to work with second and higher order tensors, however, a simple matrix vector multiplication I only got to work using the code shown below:
`
#include <boost/numeric/ublas/tensor.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include

int main() {

using namespace boost::numeric::ublas::index;
using tensor  = boost::numeric::ublas::tensor<double>;
typedef boost::numeric::ublas::vector<double> vec;

tensor A{2,2};
A.at(0,0) = 1.0;
A.at(0,1) = 3.0;
A.at(1,0) = 4.0;
A.at(1,1) = 2.0;

vec tmp(2);

tensor b(tmp);
b.at(0) = 5.0;
b.at(1) = 1.0;

tensor C = A(_i,_j)*b(_j,_k);

std::cout << C << std::endl;

return 0;

}
`

I find this notation quite cumbersome. I would expect a proper matrix vector multiplication to be expressed as A_ij b_j = C_i, however, the code only works if I add the dummy index _k. Am I missing something here?

Thanks for any suggestions in advance!

Regards

[bassoy]

Hi @gear-head,

That is true, you need an extra dummy index k although b is only a vector. This is because valid ublas tensors always have at least two indices in order distinguish between row and column vectors.

Note that you can use the placeholder _ for the dummy index:

tensor A{2,2};
tensor b{2,1};
tensor C = A(_i,_j)*b(_j,_);

or you could use the ublas' prod function if b is a ublas::vector:

tensor C = ublas::prod( A, b, 2 ) ;

See also: https://godbolt.org/z/9Yd7fqj73