diff --git a/test/triqs/arrays/dot.cpp b/test/triqs/arrays/dot.cpp index b61704b3..740d55da 100644 --- a/test/triqs/arrays/dot.cpp +++ b/test/triqs/arrays/dot.cpp @@ -33,7 +33,7 @@ template void test() { Vd a(2),aa(2),c(2) ;a()=2.0; c() = 1; Vi b(2);b()=3; - std::cout << blas::dot(a,b) << std::endl; + std::cerr << blas::dot(a,b) << std::endl; aa = 2*a; @@ -51,7 +51,19 @@ void test() { int main(int argc, char **argv) { test , vector >(); - + + /// Added by I. Krivenko, #122 + /// test the complex version, specially with the zdotu workaround on Os X. + vector> v(2); + v(0) = 0; + v(1) = {0, 1}; + + std::cerr << v << std::endl; + std::cerr << blas::dot(v, v) << std::endl; + std::cerr << blas::dot(v, v) << std::endl; + assert_close( dot(v,v), -1); + assert_close( dotc(v,v), 1); + // does not work for array just because of .size() vs .shape(0)... //test , array >(); //test , vector >(); diff --git a/triqs/arrays/blas_lapack/dot.hpp b/triqs/arrays/blas_lapack/dot.hpp index 49eb4724..3011fc33 100644 --- a/triqs/arrays/blas_lapack/dot.hpp +++ b/triqs/arrays/blas_lapack/dot.hpp @@ -67,13 +67,10 @@ namespace triqs { namespace arrays { namespace blas { //return f77::dot(X.size(), Cx().data_start(), Cx().stride(), Cy().data_start(), Cy().stride()); } - template< bool Star, typename T> - typename std::enable_if::value && Star,T>::type - _conj(T && x) { return conj(std::forward(x));} - - template< bool Star, typename T> - typename std::enable_if::value && Star),T>::type - _conj(T && x) { return std::forward(x);} + template std::complex _conj(std::complex const& x); + template <> std::complex _conj(std::complex const& x) { return conj(x); } + template <> std::complex _conj(std::complex const& x) { return x;} + template double _conj(double x) { return x; } /** * Calls dot product of 2 vectors. @@ -94,7 +91,7 @@ namespace triqs { namespace arrays { namespace blas { auto * restrict X_ = X.data_start(); auto * restrict Y_ = Y.data_start(); if ((incx==1) && (incy==1)) { - for (size_t i=0; i(X_[i]) * Y_[i]; + for (size_t i = 0; i < N; ++i) res += _conj(X_[i]) * Y_[i]; } else { // code for unequal increments or equal increments not equal to 1 for (size_t i=0, ix=0, iy=0; i(X_[ix]) * Y_[iy]; }