interface/zger.c

   1 /*********************************************************************/
   2 /* Copyright 2009, 2010 The University of Texas at Austin.           */
   3 /* All rights reserved.                                              */
   4 /*                                                                   */
   5 /* Redistribution and use in source and binary forms, with or        */
   6 /* without modification, are permitted provided that the following   */
   7 /* conditions are met:                                               */
   8 /*                                                                   */
   9 /*   1. Redistributions of source code must retain the above         */
  10 /*      copyright notice, this list of conditions and the following  */
  11 /*      disclaimer.                                                  */
  12 /*                                                                   */
  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.                              */
  17 /*                                                                   */
  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.                                    */
  32 /*                                                                   */
  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 /*********************************************************************/
  38
  39 #include <stdio.h>
  40 #include "common.h"
  41 #ifdef FUNCTION_PROFILE
  42 #include "functable.h"
  43 #endif
  44
  45 #ifdef SMP
  46 #ifdef __64BIT__
  47 #define SMPTEST 1
  48 #endif
  49 #endif
  50
  51
  52 #ifdef XDOUBLE
  53 #ifndef CONJ
  54 #define ERROR_NAME "XGERU  "
  55 #else
  56 #define ERROR_NAME "XGERC  "
  57 #endif
  58 #elif defined DOUBLE
  59 #ifndef CONJ
  60 #define ERROR_NAME "ZGERU  "
  61 #else
  62 #define ERROR_NAME "ZGERC  "
  63 #endif
  64 #else
  65 #ifndef CONJ
  66 #define ERROR_NAME "CGERU  "
  67 #else
  68 #define ERROR_NAME "CGERC "
  69 #endif
  70 #endif
  71
  72 #if   defined XDOUBLE
  73 #ifndef CONJ
  74 #define GER             GERU_K
  75 #define GER_THREAD      xger_thread_U
  76 #else
  77 #define GER             GERC_K
  78 #define GER_THREAD      xger_thread_C
  79 #define GERV            GERV_K
  80 #define GERV_THREAD     xger_thread_V
  81 #endif
  82 #elif defined DOUBLE
  83 #ifndef CONJ
  84 #define GER             GERU_K
  85 #define GER_THREAD      zger_thread_U
  86 #else
  87 #define GER             GERC_K
  88 #define GER_THREAD      zger_thread_C
  89 #define GERV            GERV_K
  90 #define GERV_THREAD     zger_thread_V
  91 #endif
  92 #else
  93 #ifndef CONJ
  94 #define GER             GERU_K
  95 #define GER_THREAD      cger_thread_U
  96 #else
  97 #define GER             GERC_K
  98 #define GER_THREAD      cger_thread_C
  99 #define GERV            GERV_K
 100 #define GERV_THREAD     cger_thread_V
 101 #endif
 102 #endif
 103
 104 #ifndef CBLAS
 105
 106 void NAME(blasint *M, blasint *N, FLOAT *Alpha,
 107           FLOAT *x, blasint *INCX,
 108           FLOAT *y, blasint *INCY,
 109           FLOAT *a, blasint *LDA){
 110
 111   blasint    m     = *M;
 112   blasint    n     = *N;
 113   FLOAT  alpha_r = Alpha[0];
 114   FLOAT  alpha_i = Alpha[1];
 115   blasint    incx  = *INCX;
 116   blasint    incy  = *INCY;
 117   blasint    lda   = *LDA;
 118   FLOAT *buffer;
 119 #ifdef SMPTEST
 120   int nthreads;
 121 #endif
 122
 123   blasint info;
 124
 125   PRINT_DEBUG_NAME;
 126
 127   info = 0;
 128
 129   if (lda < MAX(1,m)) info = 9;
 130   if (incy == 0)      info = 7;
 131   if (incx == 0)      info = 5;
 132   if (n < 0)          info = 2;
 133   if (m < 0)          info = 1;
 134
 135   if (info){
 136     BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
 137     return;
 138   }
 139
 140 #else
 141
 142 void CNAME(enum CBLAS_ORDER order,
 143            blasint m, blasint n,
 144            void *VAlpha,
 145            void  *vx, blasint incx,
 146            void  *vy, blasint incy,
 147            void  *va, blasint lda) {
 148
 149   FLOAT* Alpha = (FLOAT*) VAlpha;
 150   FLOAT* a = (FLOAT*) va;
 151   FLOAT* x = (FLOAT*) vx;
 152   FLOAT* y = (FLOAT*) vy;
 153
 154   FLOAT  alpha_r = Alpha[0];
 155   FLOAT  alpha_i = Alpha[1];
 156
 157   FLOAT *buffer;
 158   blasint info, t;
 159 #ifdef SMPTEST
 160   int nthreads;
 161 #endif
 162
 163   PRINT_DEBUG_CNAME;
 164
 165   info  =  0;
 166
 167   if (order == CblasColMajor) {
 168     info = -1;
 169
 170     if (lda < MAX(1,m)) info = 9;
 171     if (incy == 0)      info = 7;
 172     if (incx == 0)      info = 5;
 173     if (n < 0)          info = 2;
 174     if (m < 0)          info = 1;
 175   }
 176
 177   if (order == CblasRowMajor) {
 178     info = -1;
 179
 180     t = n;
 181     n = m;
 182     m = t;
 183
 184     t    = incx;
 185     incx = incy;
 186     incy = t;
 187
 188     buffer = x;
 189     x = y;
 190     y = buffer;
 191
 192     if (lda < MAX(1,m)) info = 9;
 193     if (incy == 0)      info = 7;
 194     if (incx == 0)      info = 5;
 195     if (n < 0)          info = 2;
 196     if (m < 0)          info = 1;
 197   }
 198
 199   if (info >= 0) {
 200     BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
 201     return;
 202   }
 203
 204 #endif
 205
 206   /*     Quick return if possible. */
 207   if (m == 0 || n == 0) return;
 208
 209   if ((alpha_r == 0.) && (alpha_i == 0.)) return;
 210
 211   IDEBUG_START;
 212
 213   FUNCTION_PROFILE_START();
 214
 215   if (incy < 0) y -= (n - 1) * incy * 2;
 216   if (incx < 0) x -= (m - 1) * incx * 2;
 217
 218   STACK_ALLOC(2 * m, FLOAT, buffer);
 219
 220 #ifdef SMPTEST
 221   // Threshold chosen so that speed-up is > 1 on a Xeon E5-2630
 222   if(1L * m * n > 36L * sizeof(FLOAT) * sizeof(FLOAT) * GEMM_MULTITHREAD_THRESHOLD)
 223     nthreads = num_cpu_avail(2);
 224   else
 225     nthreads = 1;
 226
 227   if (nthreads == 1) {
 228 #endif
 229
 230 #if !defined(CBLAS) || !defined(CONJ)
 231   GER(m, n, 0, alpha_r, alpha_i, x, incx, y, incy, a, lda, buffer);
 232 #else
 233   if (order == CblasColMajor) {
 234     GER(m, n, 0, alpha_r, alpha_i, x, incx, y, incy, a, lda, buffer);
 235   } else {
 236     GERV(m, n, 0, alpha_r, alpha_i, x, incx, y, incy, a, lda, buffer);
 237   }
 238 #endif
 239
 240 #ifdef SMPTEST
 241
 242   } else {
 243
 244 #if !defined(CBLAS) || !defined(CONJ)
 245       GER_THREAD(m, n, Alpha, x, incx, y, incy, a, lda, buffer, nthreads);
 246 #else
 247       if (order == CblasColMajor) {
 248         GER_THREAD(m, n, Alpha, x, incx, y, incy, a, lda, buffer, nthreads);
 249       } else {
 250         GERV_THREAD(m, n, Alpha, x, incx, y, incy, a, lda, buffer, nthreads);
 251       }
 252 #endif
 253
 254   }
 255 #endif
 256
 257   STACK_FREE(buffer);
 258
 259   FUNCTION_PROFILE_END(4, m * n + m + n, 2 * m * n);
 260
 261   IDEBUG_END;
 262
 263   return;
 264
 265 }