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 FLENS_BLAS_LEVEL2_MV_TCC
34 #define FLENS_BLAS_LEVEL2_MV_TCC 1
35
36 #include <flens/blas/closures/debugclosure.h>
37 #include <flens/typedefs.h>
38
39 namespace flens { namespace blas {
40
41 //== GeneralMatrix - Vector products ===========================================
42
43 //-- gemv
44 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
45 void
46 mv(Transpose trans,
47 const ALPHA &alpha, const GeMatrix<MA> &A, const DenseVector<VX> &x,
48 const BETA &beta, DenseVector<VY> &y)
49 {
50 ASSERT(x.length()==((trans==NoTrans) ? A.numCols()
51 : A.numRows()));
52
53 typedef typename GeMatrix<MA>::IndexType IndexType;
54 IndexType yLength = (trans==NoTrans) ? A.numRows()
55 : A.numCols();
56
57 ASSERT(beta==BETA(0) || y.length()==yLength || y.length()==0);
58
59 if (y.length()!=yLength) {
60 typedef typename DenseVector<VY>::ElementType T;
61 const T Zero(0);
62
63 FLENS_BLASLOG_RESIZE_VECTOR(y, yLength);
64 y.resize(yLength, y.firstIndex(), Zero);
65 }
66
67
68 # ifndef FLENS_DEBUG_CLOSURES
69 ASSERT(!DEBUGCLOSURE::identical(x, y));
70 # else
71 //
72 // If x and y are identical an temporary is needed if we want to use mv
73 //
74 if (DEBUGCLOSURE::identical(x, y)) {
75 typename Result<DenseVector<VX> >::Type _x;
76 FLENS_BLASLOG_TMP_ADD(_x);
77 _x = x;
78
79 mv(trans, alpha, A, _x, beta, y);
80
81 FLENS_BLASLOG_TMP_REMOVE(_x, x);
82 return;
83 }
84 # endif
85
86 FLENS_BLASLOG_SETTAG("--> ");
87 FLENS_BLASLOG_BEGIN_GEMV(trans, alpha, A, x, beta, y);
88
89 # ifdef HAVE_CXXBLAS_GEMV
90 cxxblas::gemv(MA::order,
91 trans,
92 A.numRows(), A.numCols(),
93 alpha,
94 A.data(), A.leadingDimension(),
95 x.data(), x.stride(),
96 beta,
97 y.data(), y.stride());
98 # else
99 ASSERT(0);
100 # endif
101
102 FLENS_BLASLOG_END;
103 FLENS_BLASLOG_UNSETTAG;
104 }
105
106
107 //== TriangularMatrix - Vector products ========================================
108
109 //-- trmv
110 template <typename MA, typename VX>
111 void
112 mv(Transpose trans, const TrMatrix<MA> &A, DenseVector<VX> &x)
113 {
114 FLENS_BLASLOG_SETTAG("--> ");
115 FLENS_BLASLOG_BEGIN_TRMV(trans, A, x);
116
117 ASSERT(x.length()==A.dim());
118
119 # ifdef HAVE_CXXBLAS_TRMV
120 cxxblas::trmv(MA::order, A.upLo(),
121 trans, A.diag(),
122 A.dim(),
123 A.data(), A.leadingDimension(),
124 x.data(), x.stride());
125 # else
126 ASSERT(0);
127 # endif
128
129 FLENS_BLASLOG_END;
130 FLENS_BLASLOG_UNSETTAG;
131 }
132
133
134 //== SymmetricMatrix - Vector products =========================================
135
136 //-- symv
137 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
138 void
139 mv(const ALPHA &alpha, const SyMatrix<MA> &A, const DenseVector<VX> &x,
140 const BETA &beta, DenseVector<VY> &y)
141 {
142 ASSERT(ADDRESS(y)!=ADDRESS(x));
143 ASSERT(x.length()==A.dim());
144 ASSERT((beta==BETA(0)) || (y.length()==A.dim()));
145
146 if (y.length()!=A.dim()) {
147 y.resize(A.dim(), 0);
148 }
149
150 # ifdef HAVE_CXXBLAS_SYMV
151 cxxblas::symv(MA::order, A.upLo(),
152 A.dim(),
153 alpha,
154 A.data(), A.leadingDimension(),
155 x.data(), x.stride(),
156 beta,
157 y.data(), y.stride());
158 # else
159 ASSERT(0);
160 # endif
161 }
162
163
164 //== HermitianMatrix - Vector products =========================================
165
166 //-- hemv
167 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
168 void
169 mv(const ALPHA &alpha, const HeMatrix<MA> &A, const DenseVector<VX> &x,
170 const BETA &beta, DenseVector<VY> &y)
171 {
172 ASSERT(ADDRESS(y)!=ADDRESS(x));
173 ASSERT(x.length()==A.dim());
174 ASSERT((beta==BETA(0)) || (y.length()==A.dim()));
175
176 if (y.length()!=A.dim()) {
177 y.resize(A.dim(), 0);
178 }
179
180 # ifdef HAVE_CXXBLAS_HEMV
181 cxxblas::hemv(MA::order, A.upLo(),
182 A.dim(),
183 alpha,
184 A.data(), A.leadingDimension(),
185 x.data(), x.stride(),
186 beta,
187 y.data(), y.stride());
188 # else
189 ASSERT(0);
190 # endif
191 }
192
193
194 //== forwarding ================================================================
195
196 //-- GeneralMatrix - Vector products
197 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
198 typename RestrictTo<IsGeneralMatrix<MA>::value &&
199 !IsClosure<MA>::value && !IsClosure<VX>::value &&
200 IsSame<VY, typename VY::Impl>::value,
201 void>::Type
202 mv(Transpose trans,
203 const ALPHA &alpha, const MA &A, const VX &x, const BETA &beta, VY &&y)
204 {
205 CHECKPOINT_ENTER;
206 mv(trans, alpha, A, x, beta, y);
207 CHECKPOINT_LEAVE;
208 }
209
210 //-- TriangularMatrix - Vector products
211 template <typename MA, typename VX>
212 typename RestrictTo<IsTriangularMatrix<MA>::value &&
213 !IsClosure<MA>::value &&
214 IsSame<VX, typename VX::Impl>::value,
215 void>::Type
216 mv(Transpose trans, const MA &A, VX &&x)
217 {
218 CHECKPOINT_ENTER;
219 mv(trans, A, x);
220 CHECKPOINT_LEAVE;
221 }
222
223 //-- Symmetric Matrix - Vector products
224 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
225 typename RestrictTo<IsSymmetricMatrix<MA>::value &&
226 !IsClosure<MA>::value && !IsClosure<VX>::value &&
227 IsSame<VY, typename VY::Impl>::value,
228 void>::Type
229 mv(const ALPHA &alpha, const MA &A, const VX &x, const BETA &beta, VY &&y)
230 {
231 CHECKPOINT_ENTER;
232 mv(alpha, A, x, beta, y);
233 CHECKPOINT_LEAVE;
234 }
235
236 //-- Hermitian Matrix - Vector products
237 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
238 typename RestrictTo<IsHermitianMatrix<MA>::value &&
239 !IsClosure<MA>::value && !IsClosure<VX>::value &&
240 IsSame<VY, typename VY::Impl>::value,
241 void>::Type
242 mv(const ALPHA &alpha, const MA &A, const VX &x, const BETA &beta, VY &&y)
243 {
244 CHECKPOINT_ENTER;
245 mv(alpha, A, x, beta, y);
246 CHECKPOINT_LEAVE;
247 }
248
249
250 } } // namespace blas, flens
251
252 #endif // FLENS_BLAS_LEVEL2_MV_TCC
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 FLENS_BLAS_LEVEL2_MV_TCC
34 #define FLENS_BLAS_LEVEL2_MV_TCC 1
35
36 #include <flens/blas/closures/debugclosure.h>
37 #include <flens/typedefs.h>
38
39 namespace flens { namespace blas {
40
41 //== GeneralMatrix - Vector products ===========================================
42
43 //-- gemv
44 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
45 void
46 mv(Transpose trans,
47 const ALPHA &alpha, const GeMatrix<MA> &A, const DenseVector<VX> &x,
48 const BETA &beta, DenseVector<VY> &y)
49 {
50 ASSERT(x.length()==((trans==NoTrans) ? A.numCols()
51 : A.numRows()));
52
53 typedef typename GeMatrix<MA>::IndexType IndexType;
54 IndexType yLength = (trans==NoTrans) ? A.numRows()
55 : A.numCols();
56
57 ASSERT(beta==BETA(0) || y.length()==yLength || y.length()==0);
58
59 if (y.length()!=yLength) {
60 typedef typename DenseVector<VY>::ElementType T;
61 const T Zero(0);
62
63 FLENS_BLASLOG_RESIZE_VECTOR(y, yLength);
64 y.resize(yLength, y.firstIndex(), Zero);
65 }
66
67
68 # ifndef FLENS_DEBUG_CLOSURES
69 ASSERT(!DEBUGCLOSURE::identical(x, y));
70 # else
71 //
72 // If x and y are identical an temporary is needed if we want to use mv
73 //
74 if (DEBUGCLOSURE::identical(x, y)) {
75 typename Result<DenseVector<VX> >::Type _x;
76 FLENS_BLASLOG_TMP_ADD(_x);
77 _x = x;
78
79 mv(trans, alpha, A, _x, beta, y);
80
81 FLENS_BLASLOG_TMP_REMOVE(_x, x);
82 return;
83 }
84 # endif
85
86 FLENS_BLASLOG_SETTAG("--> ");
87 FLENS_BLASLOG_BEGIN_GEMV(trans, alpha, A, x, beta, y);
88
89 # ifdef HAVE_CXXBLAS_GEMV
90 cxxblas::gemv(MA::order,
91 trans,
92 A.numRows(), A.numCols(),
93 alpha,
94 A.data(), A.leadingDimension(),
95 x.data(), x.stride(),
96 beta,
97 y.data(), y.stride());
98 # else
99 ASSERT(0);
100 # endif
101
102 FLENS_BLASLOG_END;
103 FLENS_BLASLOG_UNSETTAG;
104 }
105
106
107 //== TriangularMatrix - Vector products ========================================
108
109 //-- trmv
110 template <typename MA, typename VX>
111 void
112 mv(Transpose trans, const TrMatrix<MA> &A, DenseVector<VX> &x)
113 {
114 FLENS_BLASLOG_SETTAG("--> ");
115 FLENS_BLASLOG_BEGIN_TRMV(trans, A, x);
116
117 ASSERT(x.length()==A.dim());
118
119 # ifdef HAVE_CXXBLAS_TRMV
120 cxxblas::trmv(MA::order, A.upLo(),
121 trans, A.diag(),
122 A.dim(),
123 A.data(), A.leadingDimension(),
124 x.data(), x.stride());
125 # else
126 ASSERT(0);
127 # endif
128
129 FLENS_BLASLOG_END;
130 FLENS_BLASLOG_UNSETTAG;
131 }
132
133
134 //== SymmetricMatrix - Vector products =========================================
135
136 //-- symv
137 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
138 void
139 mv(const ALPHA &alpha, const SyMatrix<MA> &A, const DenseVector<VX> &x,
140 const BETA &beta, DenseVector<VY> &y)
141 {
142 ASSERT(ADDRESS(y)!=ADDRESS(x));
143 ASSERT(x.length()==A.dim());
144 ASSERT((beta==BETA(0)) || (y.length()==A.dim()));
145
146 if (y.length()!=A.dim()) {
147 y.resize(A.dim(), 0);
148 }
149
150 # ifdef HAVE_CXXBLAS_SYMV
151 cxxblas::symv(MA::order, A.upLo(),
152 A.dim(),
153 alpha,
154 A.data(), A.leadingDimension(),
155 x.data(), x.stride(),
156 beta,
157 y.data(), y.stride());
158 # else
159 ASSERT(0);
160 # endif
161 }
162
163
164 //== HermitianMatrix - Vector products =========================================
165
166 //-- hemv
167 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
168 void
169 mv(const ALPHA &alpha, const HeMatrix<MA> &A, const DenseVector<VX> &x,
170 const BETA &beta, DenseVector<VY> &y)
171 {
172 ASSERT(ADDRESS(y)!=ADDRESS(x));
173 ASSERT(x.length()==A.dim());
174 ASSERT((beta==BETA(0)) || (y.length()==A.dim()));
175
176 if (y.length()!=A.dim()) {
177 y.resize(A.dim(), 0);
178 }
179
180 # ifdef HAVE_CXXBLAS_HEMV
181 cxxblas::hemv(MA::order, A.upLo(),
182 A.dim(),
183 alpha,
184 A.data(), A.leadingDimension(),
185 x.data(), x.stride(),
186 beta,
187 y.data(), y.stride());
188 # else
189 ASSERT(0);
190 # endif
191 }
192
193
194 //== forwarding ================================================================
195
196 //-- GeneralMatrix - Vector products
197 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
198 typename RestrictTo<IsGeneralMatrix<MA>::value &&
199 !IsClosure<MA>::value && !IsClosure<VX>::value &&
200 IsSame<VY, typename VY::Impl>::value,
201 void>::Type
202 mv(Transpose trans,
203 const ALPHA &alpha, const MA &A, const VX &x, const BETA &beta, VY &&y)
204 {
205 CHECKPOINT_ENTER;
206 mv(trans, alpha, A, x, beta, y);
207 CHECKPOINT_LEAVE;
208 }
209
210 //-- TriangularMatrix - Vector products
211 template <typename MA, typename VX>
212 typename RestrictTo<IsTriangularMatrix<MA>::value &&
213 !IsClosure<MA>::value &&
214 IsSame<VX, typename VX::Impl>::value,
215 void>::Type
216 mv(Transpose trans, const MA &A, VX &&x)
217 {
218 CHECKPOINT_ENTER;
219 mv(trans, A, x);
220 CHECKPOINT_LEAVE;
221 }
222
223 //-- Symmetric Matrix - Vector products
224 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
225 typename RestrictTo<IsSymmetricMatrix<MA>::value &&
226 !IsClosure<MA>::value && !IsClosure<VX>::value &&
227 IsSame<VY, typename VY::Impl>::value,
228 void>::Type
229 mv(const ALPHA &alpha, const MA &A, const VX &x, const BETA &beta, VY &&y)
230 {
231 CHECKPOINT_ENTER;
232 mv(alpha, A, x, beta, y);
233 CHECKPOINT_LEAVE;
234 }
235
236 //-- Hermitian Matrix - Vector products
237 template <typename ALPHA, typename MA, typename VX, typename BETA, typename VY>
238 typename RestrictTo<IsHermitianMatrix<MA>::value &&
239 !IsClosure<MA>::value && !IsClosure<VX>::value &&
240 IsSame<VY, typename VY::Impl>::value,
241 void>::Type
242 mv(const ALPHA &alpha, const MA &A, const VX &x, const BETA &beta, VY &&y)
243 {
244 CHECKPOINT_ENTER;
245 mv(alpha, A, x, beta, y);
246 CHECKPOINT_LEAVE;
247 }
248
249
250 } } // namespace blas, flens
251
252 #endif // FLENS_BLAS_LEVEL2_MV_TCC