From 900ac486283bee64c2562933d1a0f229c25aae8e Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 20 Apr 2023 00:32:09 +0800 Subject: [PATCH 1/3] implement dotc and dotu by sdot_ and ddot_ --- include/RI/global/Blas-Fortran.h | 8 ++--- include/RI/global/Blas_Interface.h | 48 ++++++++++++++++++++++-------- 2 files changed, 40 insertions(+), 16 deletions(-) diff --git a/include/RI/global/Blas-Fortran.h b/include/RI/global/Blas-Fortran.h index 181981f..37c4af3 100644 --- a/include/RI/global/Blas-Fortran.h +++ b/include/RI/global/Blas-Fortran.h @@ -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*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); - void zdotu_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); + // void cdotu_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); + // void zdotu_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); // d = Vx * Vy - void cdotc_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); - void zdotc_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); + // void cdotc_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*const Y, const int*const incY); + // void zdotc_(std::complex*const result, const int*const n, const std::complex*const X, const int*const incX, const std::complex*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, diff --git a/include/RI/global/Blas_Interface.h b/include/RI/global/Blas_Interface.h index d15e382..1cf68ff 100644 --- a/include/RI/global/Blas_Interface.h +++ b/include/RI/global/Blas_Interface.h @@ -96,15 +96,27 @@ namespace Blas_Interface } inline std::complex dotu(const int n, const std::complex*const X, const int incX, const std::complex*const Y, const int incY) { - std::complex 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(X); + auto y = reinterpret_cast(Y); + //Re(resut)=Re(x)*Re(y)-Im(x)*Im(y) + //Im(resut)=Re(x)*Im(y)+Im(x)*Re(y) + return std::complex(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 dotu(const int n, const std::complex*const X, const int incX, const std::complex*const Y, const int incY) { - std::complex 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(X); + auto y = reinterpret_cast(Y); + //Re(resut)=Re(x)*Re(y)-Im(x)*Im(y) + //Im(resut)=Re(x)*Im(y)+Im(x)*Re(y) + return std::complex(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 @@ -118,15 +130,27 @@ namespace Blas_Interface } inline std::complex dotc(const int n, const std::complex*const X, const int incX, const std::complex*const Y, const int incY) { - std::complex 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(X); + auto y = reinterpret_cast(Y); + // Re(result)=Re(X)*Re(Y)+Im(X)*Im(Y) + // Im(result)=Re(X)*Im(Y)-Im(X)*Re(Y) + return std::complex(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 dotc(const int n, const std::complex*const X, const int incX, const std::complex*const Y, const int incY) { - std::complex 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(X); + auto y = reinterpret_cast(Y); + // Re(result)=Re(X)*Re(Y)+Im(X)*Im(Y) + // Im(result)=Re(X)*Im(Y)-Im(X)*Re(Y) + return std::complex(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 From 6b699ad3f0905c246d4e98d5631a600efbec0d4d Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 20 Apr 2023 01:31:17 +0800 Subject: [PATCH 2/3] UT for dotc and dotu with different incX and incY --- unittests/global/Blas-test.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/unittests/global/Blas-test.hpp b/unittests/global/Blas-test.hpp index 2430a46..0557a6a 100644 --- a/unittests/global/Blas-test.hpp +++ b/unittests/global/Blas-test.hpp @@ -72,6 +72,14 @@ namespace Blas_Test /* -18+68i */ std::cout< a3 = { {1,2}, {1,4}, {4,2}, {3,4}, {2,8}, {1,8}}; + const std::vector 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 From b321b71a8677a88a42bbb78e6d31c10073454e14 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 20 Apr 2023 01:58:12 +0800 Subject: [PATCH 3/3] fix a typo --- include/RI/global/Blas_Interface.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/RI/global/Blas_Interface.h b/include/RI/global/Blas_Interface.h index 1cf68ff..c729038 100644 --- a/include/RI/global/Blas_Interface.h +++ b/include/RI/global/Blas_Interface.h @@ -101,8 +101,8 @@ namespace Blas_Interface const int incY2 = 2 * incY; auto x = reinterpret_cast(X); auto y = reinterpret_cast(Y); - //Re(resut)=Re(x)*Re(y)-Im(x)*Im(y) - //Im(resut)=Re(x)*Im(y)+Im(x)*Re(y) + //Re(result)=Re(x)*Re(y)-Im(x)*Im(y) + //Im(result)=Re(x)*Im(y)+Im(x)*Re(y) return std::complex(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)); } @@ -113,8 +113,8 @@ namespace Blas_Interface const int incY2 = 2 * incY; auto x = reinterpret_cast(X); auto y = reinterpret_cast(Y); - //Re(resut)=Re(x)*Re(y)-Im(x)*Im(y) - //Im(resut)=Re(x)*Im(y)+Im(x)*Re(y) + //Re(result)=Re(x)*Re(y)-Im(x)*Im(y) + //Im(result)=Re(x)*Im(y)+Im(x)*Re(y) return std::complex(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)); }