Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions include/RI/global/Blas-Fortran.h
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,12 @@ extern "C"
// d = Vx * Vy
// reason for passing results as argument instead of returning it:
// https://www.numbercrunch.de/blog/2014/07/lost-in-translation/
void cdotu_(std::complex<float>*const result, const int*const n, const std::complex<float>*const X, const int*const incX, const std::complex<float>*const Y, const int*const incY);
void zdotu_(std::complex<double>*const result, const int*const n, const std::complex<double>*const X, const int*const incX, const std::complex<double>*const Y, const int*const incY);
// void cdotu_(std::complex<float>*const result, const int*const n, const std::complex<float>*const X, const int*const incX, const std::complex<float>*const Y, const int*const incY);
// void zdotu_(std::complex<double>*const result, const int*const n, const std::complex<double>*const X, const int*const incX, const std::complex<double>*const Y, const int*const incY);

// d = Vx * Vy
void cdotc_(std::complex<float>*const result, const int*const n, const std::complex<float>*const X, const int*const incX, const std::complex<float>*const Y, const int*const incY);
void zdotc_(std::complex<double>*const result, const int*const n, const std::complex<double>*const X, const int*const incX, const std::complex<double>*const Y, const int*const incY);
// void cdotc_(std::complex<float>*const result, const int*const n, const std::complex<float>*const X, const int*const incX, const std::complex<float>*const Y, const int*const incY);
// void zdotc_(std::complex<double>*const result, const int*const n, const std::complex<double>*const X, const int*const incX, const std::complex<double>*const Y, const int*const incY);

// Vy = alpha * Ma.? * Vx + beta * Vy
void sgemv_(const char*const transA, const int*const m, const int*const n,
Expand Down
48 changes: 36 additions & 12 deletions include/RI/global/Blas_Interface.h
Original file line number Diff line number Diff line change
Expand Up @@ -96,15 +96,27 @@ namespace Blas_Interface
}
inline std::complex<float> dotu(const int n, const std::complex<float>*const X, const int incX, const std::complex<float>*const Y, const int incY)
{
std::complex<float> result;
cdotu_(&result, &n, X, &incX, Y, &incY);
return result;
//cdotu_(&result, &n, X, &incX, Y, &incY);
const int incX2 = 2 * incX;
const int incY2 = 2 * incY;
auto x = reinterpret_cast<const float*>(X);
auto y = reinterpret_cast<const float*>(Y);
//Re(result)=Re(x)*Re(y)-Im(x)*Im(y)
//Im(result)=Re(x)*Im(y)+Im(x)*Re(y)
return std::complex<float>(sdot_(&n, x, &incX2, y, &incY2) - sdot_(&n, x + 1, &incX2, y + 1, &incY2),
sdot_(&n, x, &incX2, y + 1, &incY2) + sdot_(&n, x + 1, &incX2, y, &incY2));
}
inline std::complex<double> dotu(const int n, const std::complex<double>*const X, const int incX, const std::complex<double>*const Y, const int incY)
{
std::complex<double> result;
zdotu_(&result, &n, X, &incX, Y, &incY);
return result;
//zdotu_(&result, &n, X, &incX, Y, &incY);
const int incX2 = 2 * incX;
const int incY2 = 2 * incY;
auto x = reinterpret_cast<const double*>(X);
auto y = reinterpret_cast<const double*>(Y);
//Re(result)=Re(x)*Re(y)-Im(x)*Im(y)
//Im(result)=Re(x)*Im(y)+Im(x)*Re(y)
return std::complex<double>(ddot_(&n, x, &incX2, y, &incY2) - ddot_(&n, x + 1, &incX2, y + 1, &incY2),
ddot_(&n, x, &incX2, y + 1, &incY2) + ddot_(&n, x + 1, &incX2, y, &incY2));
}

// d = Vx.conj() * Vy
Expand All @@ -118,15 +130,27 @@ namespace Blas_Interface
}
inline std::complex<float> dotc(const int n, const std::complex<float>*const X, const int incX, const std::complex<float>*const Y, const int incY)
{
std::complex<float> result;
cdotc_(&result, &n, X, &incX, Y, &incY);
return result;
//cdotc_(&result, &n, X, &incX, Y, &incY);
const int incX2 = 2 * incX;
const int incY2 = 2 * incY;
auto x = reinterpret_cast<const float*>(X);
auto y = reinterpret_cast<const float*>(Y);
// Re(result)=Re(X)*Re(Y)+Im(X)*Im(Y)
// Im(result)=Re(X)*Im(Y)-Im(X)*Re(Y)
return std::complex<float>(sdot_(&n, x, &incX2, y, &incY2) + sdot_(&n, x + 1, &incX2, y + 1, &incY2),
sdot_(&n, x, &incX2, y + 1, &incY2) - sdot_(&n, x + 1, &incX2, y, &incY2));
}
inline std::complex<double> dotc(const int n, const std::complex<double>*const X, const int incX, const std::complex<double>*const Y, const int incY)
{
std::complex<double> result;
zdotc_(&result, &n, X, &incX, Y, &incY);
return result;
//zdotc_(&result, &n, X, &incX, Y, &incY);
const int incX2 = 2 * incX;
const int incY2 = 2 * incY;
auto x = reinterpret_cast<const double*>(X);
auto y = reinterpret_cast<const double*>(Y);
// Re(result)=Re(X)*Re(Y)+Im(X)*Im(Y)
// Im(result)=Re(X)*Im(Y)-Im(X)*Re(Y)
return std::complex<double>(ddot_(&n, x, &incX2, y, &incY2) + ddot_(&n, x + 1, &incX2, y + 1, &incY2),
ddot_(&n, x, &incX2, y+1, &incY2) - ddot_(&n, x + 1, &incX2, y, &incY2));
}

// Vy = alpha * Ma.? * Vx + beta * Vy
Expand Down
8 changes: 8 additions & 0 deletions unittests/global/Blas-test.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,14 @@ namespace Blas_Test
/* -18+68i */
std::cout<<RI::Blas_Interface::dotc(a.size(), a.data(), 1, b.data(), 1)<<std::endl;
/* 70-8i */

// test for different incX and incY
const std::vector<Tdata> a3 = { {1,2}, {1,4}, {4,2}, {3,4}, {2,8}, {1,8}};
const std::vector<Tdata> b2 = { {5,6}, {5,7}, {7,8}, {7,5}};
std::cout << RI::Blas_Interface::dotu(a.size(), a3.data(), 3, b2.data(), 2) << std::endl;
/* -18+68i */
std::cout << RI::Blas_Interface::dotc(a.size(), a3.data(), 3, b2.data(), 2) << std::endl;
/* 70-8i */
}

template<typename Tdata>
Expand Down