1 /*
2 * Copyright (c) 2009, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 #ifndef CXXBLAS_LEVEL1_DOT_TCC
34 #define CXXBLAS_LEVEL1_DOT_TCC 1
35
36 #include "xflens/cxxblas/cxxblas.h"
37
38 namespace cxxblas {
39
40 template <typename IndexType, typename X, typename Y, typename Result>
41 void
dotu_generic(IndexType n,const X * x,IndexType incX,const Y * y,IndexType incY,Result & result)42 dotu_generic(IndexType n,
43 const X *x, IndexType incX, const Y *y, IndexType incY,
44 Result &result)
45 {
46 CXXBLAS_DEBUG_OUT("dotu_generic");
47
48 result = Result(0);
49 for (IndexType i=0, iX=0, iY=0; i<n; ++i, iX+=incX, iY+=incY) {
50 result += Result(x[iX])*Result(y[iY]);
51 }
52 }
53
54
55 template <typename IndexType, typename X, typename Y, typename Result>
56 void
dot_generic(IndexType n,const X * x,IndexType incX,const Y * y,IndexType incY,Result & result)57 dot_generic(IndexType n,
58 const X *x, IndexType incX, const Y *y, IndexType incY,
59 Result &result)
60 {
61 CXXBLAS_DEBUG_OUT("dot_generic");
62
63 result = Result(0);
64 for (IndexType i=0, iX=0, iY=0; i<n; ++i, iX+=incX, iY+=incY) {
65 result += Result(conjugate(x[iX]))*Result(y[iY]);
66 }
67 }
68
69 //------------------------------------------------------------------------------
70
71 template <typename IndexType, typename X, typename Y, typename Result>
72 void
dotu(IndexType n,const X * x,IndexType incX,const Y * y,IndexType incY,Result & result)73 dotu(IndexType n,
74 const X *x, IndexType incX, const Y *y, IndexType incY,
75 Result &result)
76 {
77 if (incX<0) {
78 x -= incX*(n-1);
79 }
80 if (incY<0) {
81 y -= incY*(n-1);
82 }
83 dotu_generic(n, x, incX, y, incY, result);
84 }
85
86 template <typename IndexType, typename X, typename Y, typename Result>
87 void
dot(IndexType n,const X * x,IndexType incX,const Y * y,IndexType incY,Result & result)88 dot(IndexType n,
89 const X *x, IndexType incX, const Y *y, IndexType incY,
90 Result &result)
91 {
92 if (incX<0) {
93 x -= incX*(n-1);
94 }
95 if (incY<0) {
96 y -= incY*(n-1);
97 }
98 dot_generic(n, x, incX, y, incY, result);
99 }
100
101 #ifdef HAVE_CBLAS
102
103 // sdsdot
104 template <typename IndexType>
105 typename If<IndexType>::isBlasCompatibleInteger
sdot(IndexType n,float alpha,const float * x,IndexType incX,const float * y,IndexType incY,float & result)106 sdot(IndexType n, float alpha,
107 const float *x, IndexType incX,
108 const float *y, IndexType incY,
109 float &result)
110 {
111 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_sdsdot");
112
113 result = cblas_sdsdot(n, alpha, x, incX, y, incY);
114 }
115
116 // dsdot
117 template <typename IndexType>
118 typename If<IndexType>::isBlasCompatibleInteger
dot(IndexType n,const float * x,IndexType incX,const float * y,IndexType incY,double & result)119 dot(IndexType n,
120 const float *x, IndexType incX,
121 const float *y, IndexType incY,
122 double &result)
123 {
124 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_dsdot");
125
126 result = cblas_dsdot(n, x, incX, y, incY);
127 }
128
129 // sdot
130 template <typename IndexType>
131 typename If<IndexType>::isBlasCompatibleInteger
dot(IndexType n,const float * x,IndexType incX,const float * y,IndexType incY,float & result)132 dot(IndexType n,
133 const float *x, IndexType incX,
134 const float *y, IndexType incY,
135 float &result)
136 {
137 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_sdot");
138
139 result = cblas_sdot(n, x, incX, y, incY);
140 }
141
142 // ddot
143 template <typename IndexType>
144 typename If<IndexType>::isBlasCompatibleInteger
dot(IndexType n,const double * x,IndexType incX,const double * y,IndexType incY,double & result)145 dot(IndexType n,
146 const double *x, IndexType incX,
147 const double *y, IndexType incY,
148 double &result)
149 {
150 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_ddot");
151
152 result = cblas_ddot(n, x, incX, y, incY);
153 }
154
155 // cdotu_sub
156 template <typename IndexType>
157 typename If<IndexType>::isBlasCompatibleInteger
dotu(IndexType n,const ComplexFloat * x,IndexType incX,const ComplexFloat * y,IndexType incY,ComplexFloat & result)158 dotu(IndexType n,
159 const ComplexFloat *x, IndexType incX,
160 const ComplexFloat *y, IndexType incY,
161 ComplexFloat &result)
162 {
163 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_cdotu_sub");
164
165 cblas_cdotu_sub(n, reinterpret_cast<const float *>(x), incX,
166 reinterpret_cast<const float *>(y), incY,
167 reinterpret_cast<float *>(&result));
168 }
169
170
171 // cdotc_sub
172 template <typename IndexType>
173 typename If<IndexType>::isBlasCompatibleInteger
dot(IndexType n,const ComplexFloat * x,IndexType incX,const ComplexFloat * y,IndexType incY,ComplexFloat & result)174 dot(IndexType n,
175 const ComplexFloat *x, IndexType incX,
176 const ComplexFloat *y, IndexType incY,
177 ComplexFloat &result)
178 {
179 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_cdotc_sub");
180
181 cblas_cdotc_sub(n, reinterpret_cast<const float *>(x), incX,
182 reinterpret_cast<const float *>(y), incY,
183 reinterpret_cast<float *>(&result));
184 }
185
186 // zdotu_sub
187 template <typename IndexType>
188 typename If<IndexType>::isBlasCompatibleInteger
dotu(IndexType n,const ComplexDouble * x,IndexType incX,const ComplexDouble * y,IndexType incY,ComplexDouble & result)189 dotu(IndexType n,
190 const ComplexDouble *x, IndexType incX,
191 const ComplexDouble *y, IndexType incY,
192 ComplexDouble &result)
193 {
194 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_zdotu_sub");
195
196 cblas_zdotu_sub(n, reinterpret_cast<const double *>(x), incX,
197 reinterpret_cast<const double *>(y), incY,
198 reinterpret_cast<double *>(&result));
199 }
200
201 // zdotc_sub
202 template <typename IndexType>
203 typename If<IndexType>::isBlasCompatibleInteger
dot(IndexType n,const ComplexDouble * x,IndexType incX,const ComplexDouble * y,IndexType incY,ComplexDouble & result)204 dot(IndexType n,
205 const ComplexDouble *x, IndexType incX,
206 const ComplexDouble *y, IndexType incY,
207 ComplexDouble &result)
208 {
209 CXXBLAS_DEBUG_OUT("[" BLAS_IMPL "] cblas_zdotc_sub");
210
211 cblas_zdotc_sub(n, reinterpret_cast<const double *>(x), incX,
212 reinterpret_cast<const double *>(y), incY,
213 reinterpret_cast<double *>(&result));
214 }
215
216 #endif // HAVE_CBLAS
217
218 } // namespace cxxblas
219
220 #endif // CXXBLAS_LEVEL1_DOT_TCC
221