romeric / Fastor

A lightweight high performance tensor algebra framework for modern C++
MIT License
720 stars 68 forks source link

Issue with einsum on four tensors #176

Open ThomasYangth opened 1 week ago

ThomasYangth commented 1 week ago

I am using einsum on four tensors, in a line like einsum<Index<1,2,3>,Index<6,1,4,5>,Index<7,4,2>,Index<8,5,3>>(A, B, C, D) the shapes of A,B,C,D are, respectively: {3, 4, 4}, {3, 3, 2, 2}, {4, 2, 4}, {4, 2, 4}. This lines raises an error upon compilation, the error message is:

/home/ty1475/CPPLib/Fastor/meta/einsum_meta.h: In instantiation of ‘constexpr const int Fastor::contraction_impl<Fastor::Index<7, 4, 2, 6, 4, 2, 8>, Fastor::Tensor<std::complex, 4, 2, 4, 4, 3, 2, 4>, Fastor::std_ext::index_sequence<0, 1, 2, 3, 4, 5, 6> >::result [7]’:

/home/ty1475/CPPLib/Fastor/meta/einsum_meta.h:80:74: required from ‘struct Fastor::contraction_impl<Fastor::Index<7, 4, 2, 6, 4, 2, 8>, Fastor::Tensor<std::complex, 4, 2, 4, 4, 3, 2, 4>, Fastor::std_ext::index_sequence<0, 1, 2, 3, 4, 5, 6> >’

/home/ty1475/CPPLib/Fastor/tensor_algebra/contraction.h:288:7: required by substitution of ‘template<class T, long unsigned int ...Rest0, long unsigned int ...Rest1> static typename Fastor::contraction_impl<Fastor::Index<7, 4, 2, 6, 4, 2, 8>, Fastor::Tensor<T, Rest0 ..., Rest1 ...>, typename Fastor::std_ext::make_index_sequence<(sizeof... (Rest0) + sizeof... (Rest1))>::type>::type Fastor::extractor_contract_2<Fastor::Index<7, 4, 2>, Fastor::Index<6, 4, 2, 8>, void>::contract_impl<T, Rest0 ..., Rest1 ...>(const Fastor::Tensor<T, Rest ...>&, const Fastor::Tensor<T, Rest ...>&) [with T = std::complex; long unsigned int ...Rest0 = {4, 2, 4}; long unsigned int ...Rest1 = {4, 3, 2, 4}]’

/home/ty1475/CPPLib/Fastor/tensor_algebra/einsum.h:57:65: required by substitution of ‘template<class Index_I, class Index_J, class T, long unsigned int ...Rest0, long unsigned int ...Rest1, typename std::enable_if<((((! Fastor::is_pair_reduction<Idx0, Idx1>::value) && (! Fastor::internal::is_generalised_matrix_vector<Idx0, Idx1>::value)) && (! Fastor::internal::is_generalised_vector_matrix<Idx0, Idx1>::value)) && (! Fastor::internal::is_generalised_matrix_matrix<Idx0, Idx1>::value)), bool>::type > decltype (Fastor::extractor_contract_2<Index_I, Index_J>::contract_impl(a, b)) Fastor::einsum(const Fastor::Tensor<T, Rest ...>&, const Fastor::Tensor<T, Rest1 ...>&) [with Index_I = Fastor::Index<7, 4, 2>; Index_J = Fastor::Index<6, 4, 2, 8>; T = std::complex; long unsigned int ...Rest0 = {4, 2, 4}; long unsigned int ...Rest1 = {4, 3, 2, 4}; typename std::enable_if<((((! Fastor::is_pair_reduction<Idx0, Idx1>::value) && (! Fastor::internal::is_generalised_matrix_vector<Idx0, Idx1>::value)) && (! Fastor::internal::is_generalised_vector_matrix<Idx0, Idx1>::value)) && (! Fastor::internal::is_generalised_matrix_matrix<Idx0, Idx1>::value)), bool>::type = 0]’

/home/ty1475/CPPLib/Fastor/tensor_algebra/network_contraction.h:108:60: required from ‘static typename Fastor::contraction_impl<Fastor::Index<Idx0 ..., Idx1 ..., Idx2 ..., Idx3 ...>, Fastor::Tensor<T, Rest0 ..., Rest1 ..., Rest2 ..., Rest3 ...>, typename Fastor::std_ext::make_index_sequence<(((sizeof... (Rest0) + sizeof... (Rest1)) + sizeof... (Rest2)) + sizeof... (Rest3))>::type>::type Fastor::extractor_contract_4<Fastor::Index<I1 ...>, Fastor::Index<I2 ...>, Fastor::Index<Ind ...>, Fastor::Index<Rest1 ...> >::contract_impl(const Fastor::Tensor<Int, IterSizes ...>&, const Fastor::Tensor<T, Rest1 ...>&, const Fastor::Tensor<T, Rest2 ...>&, const Fastor::Tensor<T, Rest3 ...>&) [with T = std::complex; long unsigned int ...Rest0 = {3, 4, 4}; long unsigned int ...Rest1 = {3, 3, 2, 2}; long unsigned int ...Rest2 = {4, 2, 4}; long unsigned int ...Rest3 = {4, 2, 4}; long unsigned int ...Idx0 = {1, 2, 3}; long unsigned int ...Idx1 = {6, 1, 4, 5}; long unsigned int ...Idx2 = {7, 4, 2}; long unsigned int ...Idx3 = {8, 5, 3}; typename Fastor::contraction_impl<Fastor::Index<Idx0 ..., Idx1 ..., Idx2 ..., Idx3 ...>, Fastor::Tensor<T, Rest0 ..., Rest1 ..., Rest2 ..., Rest3 ...>, typename Fastor::std_ext::make_index_sequence<(((sizeof... (Rest0) + sizeof... (Rest1)) + sizeof... (Rest2)) + sizeof... (Rest3))>::type>::type = Fastor::Tensor<std::complex, 3, 4, 4>]’

/home/ty1475/CPPLib/Fastor/tensor_algebra/network_einsum.h:33:80: required from ‘decltype (Fastor::extractor_contract_4<Index_I, Index_J, Index_K, Index_L>::contract_impl(a, b, c, d)) Fastor::einsum(const Fastor::Tensor<T, Rest0 ...>&, const Fastor::Tensor<T, Rest1 ...>&, const Fastor::Tensor<T, Rest2 ...>&, const Fastor::Tensor<T, Rest3 ...>&) [with Index_I = Fastor::Index<1, 2, 3>; Index_J = Fastor::Index<6, 1, 4, 5>; Index_K = Fastor::Index<7, 4, 2>; Index_L = Fastor::Index<8, 5, 3>; T = std::complex; long unsigned int ...Rest0 = {3, 4, 4}; long unsigned int ...Rest1 = {3, 3, 2, 2}; long unsigned int ...Rest2 = {4, 2, 4}; long unsigned int ...Rest3 = {4, 2, 4}; decltype (Fastor::extractor_contract_4<Index_I, Index_J, Index_K, Index_L>::contract_impl(a, b, c, d)) = Fastor::Tensor<std::complex, 3, 4, 4>]’

/home/ty1475/CPPLib/Fastor/meta/einsum_meta.h:75:56: in ‘constexpr’ expansion of ‘Fastor::calc<7>(1, Fastor::contraction_impl<Fastor::Index<7, 4, 2, 6, 4, 2, 8>, Fastor::Tensor<std::complex, 4, 2, 4, 4, 3, 2, 4>, Fastor::std_ext::index_sequence<0, 1, 2, 3, 4, 5, 6> >::ind, Fastor::contraction_impl<Fastor::Index<7, 4, 2, 6, 4, 2, 8>, Fastor::Tensor<std::complex, 4, 2, 4, 4, 3, 2, 4>, Fastor::std_ext::index_sequence<0, 1, 2, 3, 4, 5, 6> >::dim)’ /home/ty1475/CPPLib/Fastor/meta/einsum_meta.h:41:69: error: expression ‘’ is not a constant expression check_all_eq(ind[i], dim[i], ind, dim) ? 1001001 : throw "dimension mismatch";

As I see it, einsum first contracts A, B, and D. Yet there seems to be a mismatch in the ordering of indices and tensor dimensions after this contraction. The indices are ordered as <6,4,2,8>, which seems to use the order {B,A,D}; however, the resulting tensor dimension is {4,3,2,4}, which is ordered like {A,B,D}, and should actually correspond to indices <2,6,4,8>.

Is this an issue with the implementation of einsum, or a misuse of the function?

romeric commented 1 week ago

Are these with complex numbers? Tensor contraction on complex numbers have some known bugs

ThomasYangth commented 1 week ago

Are these with complex numbers? Tensor contraction on complex numbers have some known bugs

Yes, they are with complex numbers. I managed to overcome the problem by doing the contractions pair by pair.