Skip to content

Commit 35c750d

Browse files
PeizeLinPeizeLin
committed
1. update Blas_Interface::dotu() and Blas_Interface::dotc()
1 parent e7e6a4d commit 35c750d

File tree

1 file changed

+26
-44
lines changed

1 file changed

+26
-44
lines changed

‎include/RI/global/Blas_Interface.h‎

Lines changed: 26 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -86,71 +86,53 @@ namespace Blas_Interface
8686
}
8787

8888
// d = Vx * Vy
89-
inline float dotu(const int n, const float*const X, const int incX, const float*const Y, const int incY)
89+
template<typename Tdata,
90+
typename std::enable_if<std::is_floating_point<Tdata>::value,bool>::type=0>
91+
inline Tdata dotu(const int n, const Tdata*const X, const int incX, const Tdata*const Y, const int incY)
9092
{
91-
return sdot_(&n, X, &incX, Y, &incY);
92-
}
93-
inline double dotu(const int n, const double*const X, const int incX, const double*const Y, const int incY)
94-
{
95-
return ddot_(&n, X, &incX, Y, &incY);
93+
return dot(n, X, incX, Y, incY);
9694
}
97-
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)
95+
template<typename Tdata,
96+
typename std::enable_if<std::is_floating_point<Tdata>::value,bool>::type=0>
97+
inline std::complex<Tdata> dotu(const int n, const std::complex<Tdata>*const X, const int incX, const std::complex<Tdata>*const Y, const int incY)
9898
{
9999
//cdotu_(&result, &n, X, &incX, Y, &incY);
100-
const int incX2 = 2 * incX;
101-
const int incY2 = 2 * incY;
102-
auto x = reinterpret_cast<const float*>(X);
103-
auto y = reinterpret_cast<const float*>(Y);
104-
//Re(result)=Re(x)*Re(y)-Im(x)*Im(y)
105-
//Im(result)=Re(x)*Im(y)+Im(x)*Re(y)
106-
return std::complex<float>(sdot_(&n, x, &incX2, y, &incY2) - sdot_(&n, x + 1, &incX2, y + 1, &incY2),
107-
sdot_(&n, x, &incX2, y + 1, &incY2) + sdot_(&n, x + 1, &incX2, y, &incY2));
108-
}
109-
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)
110-
{
111100
//zdotu_(&result, &n, X, &incX, Y, &incY);
101+
112102
const int incX2 = 2 * incX;
113103
const int incY2 = 2 * incY;
114-
auto x = reinterpret_cast<const double*>(X);
115-
auto y = reinterpret_cast<const double*>(Y);
104+
const Tdata*const x = reinterpret_cast<const Tdata*const>(X);
105+
const Tdata*const y = reinterpret_cast<const Tdata*const>(Y);
116106
//Re(result)=Re(x)*Re(y)-Im(x)*Im(y)
117107
//Im(result)=Re(x)*Im(y)+Im(x)*Re(y)
118-
return std::complex<double>(ddot_(&n, x, &incX2, y, &incY2) - ddot_(&n, x + 1, &incX2, y + 1, &incY2),
119-
ddot_(&n, x, &incX2, y + 1, &incY2) + ddot_(&n, x + 1, &incX2, y, &incY2));
108+
return std::complex<Tdata>(
109+
dot(n, x, incX2, y, incY2) - dot(n, x+1, incX2, y+1, incY2),
110+
dot(n, x, incX2, y+1, incY2) + dot(n, x+1, incX2, y, incY2));
120111
}
121112

122113
// d = Vx.conj() * Vy
123-
inline float dotc(const int n, const float*const X, const int incX, const float*const Y, const int incY)
114+
template<typename Tdata,
115+
typename std::enable_if<std::is_floating_point<Tdata>::value,bool>::type=0>
116+
inline Tdata dotc(const int n, const Tdata*const X, const int incX, const Tdata*const Y, const int incY)
124117
{
125-
return sdot_(&n, X, &incX, Y, &incY);
126-
}
127-
inline double dotc(const int n, const double*const X, const int incX, const double*const Y, const int incY)
128-
{
129-
return ddot_(&n, X, &incX, Y, &incY);
118+
return dot(n, X, incX, Y, incY);
130119
}
131-
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)
120+
template<typename Tdata,
121+
typename std::enable_if<std::is_floating_point<Tdata>::value,bool>::type=0>
122+
inline std::complex<Tdata> dotc(const int n, const std::complex<Tdata>*const X, const int incX, const std::complex<Tdata>*const Y, const int incY)
132123
{
133124
//cdotc_(&result, &n, X, &incX, Y, &incY);
134-
const int incX2 = 2 * incX;
135-
const int incY2 = 2 * incY;
136-
auto x = reinterpret_cast<const float*>(X);
137-
auto y = reinterpret_cast<const float*>(Y);
138-
// Re(result)=Re(X)*Re(Y)+Im(X)*Im(Y)
139-
// Im(result)=Re(X)*Im(Y)-Im(X)*Re(Y)
140-
return std::complex<float>(sdot_(&n, x, &incX2, y, &incY2) + sdot_(&n, x + 1, &incX2, y + 1, &incY2),
141-
sdot_(&n, x, &incX2, y + 1, &incY2) - sdot_(&n, x + 1, &incX2, y, &incY2));
142-
}
143-
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)
144-
{
145125
//zdotc_(&result, &n, X, &incX, Y, &incY);
126+
146127
const int incX2 = 2 * incX;
147128
const int incY2 = 2 * incY;
148-
auto x = reinterpret_cast<const double*>(X);
149-
auto y = reinterpret_cast<const double*>(Y);
129+
const Tdata*const x = reinterpret_cast<const Tdata*const>(X);
130+
const Tdata*const y = reinterpret_cast<const Tdata*const>(Y);
150131
// Re(result)=Re(X)*Re(Y)+Im(X)*Im(Y)
151132
// Im(result)=Re(X)*Im(Y)-Im(X)*Re(Y)
152-
return std::complex<double>(ddot_(&n, x, &incX2, y, &incY2) + ddot_(&n, x + 1, &incX2, y + 1, &incY2),
153-
ddot_(&n, x, &incX2, y+1, &incY2) - ddot_(&n, x + 1, &incX2, y, &incY2));
133+
return std::complex<Tdata>(
134+
dot(n, x, incX2, y, incY2) + dot(n, x+1, incX2, y+1, incY2),
135+
dot(n, x, incX2, y+1, incY2) - dot(n, x+1, incX2, y, incY2));
154136
}
155137

156138
// Vy = alpha * Ma.? * Vx + beta * Vy

0 commit comments

Comments
 (0)