1 /*********************************************************************/
2 /* Copyright 2009, 2010 The University of Texas at Austin. */
3 /* All rights reserved. */
5 /* Redistribution and use in source and binary forms, with or */
6 /* without modification, are permitted provided that the following */
7 /* conditions are met: */
9 /* 1. Redistributions of source code must retain the above */
10 /* copyright notice, this list of conditions and the following */
13 /* 2. Redistributions in binary form must reproduce the above */
14 /* copyright notice, this list of conditions and the following */
15 /* disclaimer in the documentation and/or other materials */
16 /* provided with the distribution. */
18 /* THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF TEXAS AT */
19 /* AUSTIN ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, */
20 /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
21 /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
22 /* DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF TEXAS AT */
23 /* AUSTIN OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
24 /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES */
25 /* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE */
26 /* GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR */
27 /* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF */
28 /* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */
29 /* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT */
30 /* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
31 /* POSSIBILITY OF SUCH DAMAGE. */
33 /* The views and conclusions contained in the software and */
34 /* documentation are those of the authors and should not be */
35 /* interpreted as representing official policies, either expressed */
36 /* or implied, of The University of Texas at Austin. */
37 /*********************************************************************/
41 #ifdef FUNCTION_PROFILE
42 #include "functable.h"
46 #define ERROR_NAME "XGEMV "
48 #define ERROR_NAME "ZGEMV "
50 #define ERROR_NAME "CGEMV "
54 static int (*gemv_thread[])(BLASLONG, BLASLONG, FLOAT *, FLOAT *, BLASLONG, FLOAT * , BLASLONG, FLOAT *, BLASLONG, FLOAT *, int) = {
56 xgemv_thread_n, xgemv_thread_t, xgemv_thread_r, xgemv_thread_c, xgemv_thread_o, xgemv_thread_u, xgemv_thread_s, xgemv_thread_d,
58 zgemv_thread_n, zgemv_thread_t, zgemv_thread_r, zgemv_thread_c, zgemv_thread_o, zgemv_thread_u, zgemv_thread_s, zgemv_thread_d,
60 cgemv_thread_n, cgemv_thread_t, cgemv_thread_r, cgemv_thread_c, cgemv_thread_o, cgemv_thread_u, cgemv_thread_s, cgemv_thread_d,
67 void NAME(char *TRANS, blasint *M, blasint *N,
68 FLOAT *ALPHA, FLOAT *a, blasint *LDA,
69 FLOAT *x, blasint *INCX,
70 FLOAT *BETA, FLOAT *y, blasint *INCY){
85 int (*gemv[])(BLASLONG, BLASLONG, BLASLONG, FLOAT, FLOAT, FLOAT *, BLASLONG,
86 FLOAT * , BLASLONG, FLOAT *, BLASLONG, FLOAT *) = {
87 GEMV_N, GEMV_T, GEMV_R, GEMV_C,
88 GEMV_O, GEMV_U, GEMV_S, GEMV_D,
95 FLOAT alpha_r = *(ALPHA + 0);
96 FLOAT alpha_i = *(ALPHA + 1);
98 FLOAT beta_r = *(BETA + 0);
99 FLOAT beta_i = *(BETA + 1);
109 if (trans == 'N') i = 0;
110 if (trans == 'T') i = 1;
111 if (trans == 'R') i = 2;
112 if (trans == 'C') i = 3;
113 if (trans == 'O') i = 4;
114 if (trans == 'U') i = 5;
115 if (trans == 'S') i = 6;
116 if (trans == 'D') i = 7;
118 if (incy == 0) info = 11;
119 if (incx == 0) info = 8;
120 if (lda < MAX(1,m)) info = 6;
128 BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
134 void CNAME(enum CBLAS_ORDER order,
135 enum CBLAS_TRANSPOSE TransA,
136 blasint m, blasint n,
138 void *va, blasint lda,
139 void *vx, blasint incx,
141 void *vy, blasint incy){
143 FLOAT *ALPHA = (FLOAT*) VALPHA;
144 FLOAT *a = (FLOAT*) va;
145 FLOAT *x = (FLOAT*) vx;
146 FLOAT *BETA = (FLOAT*) VBETA;
147 FLOAT *y = (FLOAT*) vy;
150 int trans, buffer_size;
156 int (*gemv[])(BLASLONG, BLASLONG, BLASLONG, FLOAT, FLOAT, FLOAT *, BLASLONG,
157 FLOAT * , BLASLONG, FLOAT *, BLASLONG, FLOAT *) = {
158 GEMV_N, GEMV_T, GEMV_R, GEMV_C,
159 GEMV_O, GEMV_U, GEMV_S, GEMV_D,
162 FLOAT alpha_r = *(ALPHA + 0);
163 FLOAT alpha_i = *(ALPHA + 1);
165 FLOAT beta_r = *(BETA + 0);
166 FLOAT beta_i = *(BETA + 1);
173 if (order == CblasColMajor) {
174 if (TransA == CblasNoTrans) trans = 0;
175 if (TransA == CblasTrans) trans = 1;
176 if (TransA == CblasConjNoTrans) trans = 2;
177 if (TransA == CblasConjTrans) trans = 3;
181 if (incy == 0) info = 11;
182 if (incx == 0) info = 8;
183 if (lda < MAX(1, m)) info = 6;
186 if (trans < 0) info = 1;
190 if (order == CblasRowMajor) {
191 if (TransA == CblasNoTrans) trans = 1;
192 if (TransA == CblasTrans) trans = 0;
193 if (TransA == CblasConjNoTrans) trans = 3;
194 if (TransA == CblasConjTrans) trans = 2;
202 if (incy == 0) info = 11;
203 if (incx == 0) info = 8;
204 if (lda < MAX(1, m)) info = 6;
207 if (trans < 0) info = 1;
212 BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
218 /* Quick return if possible. */
220 if (m == 0 || n == 0) return;
225 if (trans & 1) lenx = m;
226 if (trans & 1) leny = n;
228 if (beta_r != ONE || beta_i != ZERO) SCAL_K(leny, 0, 0, beta_r, beta_i, y, blasabs(incy), NULL, 0, NULL, 0);
230 if (alpha_r == ZERO && alpha_i == ZERO) return;
234 FUNCTION_PROFILE_START();
236 if (incx < 0) x -= (lenx - 1) * incx * 2;
237 if (incy < 0) y -= (leny - 1) * incy * 2;
239 buffer_size = 2 * (m + n) + 128 / sizeof(FLOAT);
241 buffer_size += 160 / sizeof(FLOAT) ;
244 buffer_size = (buffer_size + 3) & ~3;
245 STACK_ALLOC(buffer_size, FLOAT, buffer);
247 #if defined(ARCH_X86_64) && defined(MAX_STACK_ALLOC) && MAX_STACK_ALLOC > 0
248 // cgemv_t.S return NaN if there are NaN or Inf in the buffer (see bug #746)
249 if(trans && stack_alloc_size)
250 memset(buffer, 0, MIN(BUFFER_SIZE, sizeof(FLOAT) * buffer_size));
255 if ( 1L * m * n < 1024L * GEMM_MULTITHREAD_THRESHOLD )
258 nthreads = num_cpu_avail(2);
263 (gemv[(int)trans])(m, n, 0, alpha_r, alpha_i, a, lda, x, incx, y, incy, buffer);
269 (gemv_thread[(int)trans])(m, n, ALPHA, a, lda, x, incx, y, incy, buffer, nthreads);
276 FUNCTION_PROFILE_END(4, m * n + m + n, 2 * m * n);