interface/symm.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 <ctype.h>
  41 #include "common.h"
  42 #ifdef FUNCTION_PROFILE
  43 #include "functable.h"
  44 #endif
  45
  46 #ifndef COMPLEX
  47 #ifdef XDOUBLE
  48 #define ERROR_NAME "QSYMM "
  49 #elif defined(DOUBLE)
  50 #define ERROR_NAME "DSYMM "
  51 #else
  52 #define ERROR_NAME "SSYMM "
  53 #endif
  54 #else
  55 #ifndef GEMM3M
  56 #ifndef HEMM
  57 #ifdef XDOUBLE
  58 #define ERROR_NAME "XSYMM "
  59 #elif defined(DOUBLE)
  60 #define ERROR_NAME "ZSYMM "
  61 #else
  62 #define ERROR_NAME "CSYMM "
  63 #endif
  64 #else
  65 #ifdef XDOUBLE
  66 #define ERROR_NAME "XHEMM "
  67 #elif defined(DOUBLE)
  68 #define ERROR_NAME "ZHEMM "
  69 #else
  70 #define ERROR_NAME "CHEMM "
  71 #endif
  72 #endif
  73 #else
  74 #ifndef HEMM
  75 #ifdef XDOUBLE
  76 #define ERROR_NAME "XSYMM3M "
  77 #elif defined(DOUBLE)
  78 #define ERROR_NAME "ZSYMM3M "
  79 #else
  80 #define ERROR_NAME "CSYMM3M "
  81 #endif
  82 #else
  83 #ifdef XDOUBLE
  84 #define ERROR_NAME "XHEMM3M "
  85 #elif defined(DOUBLE)
  86 #define ERROR_NAME "ZHEMM3M "
  87 #else
  88 #define ERROR_NAME "CHEMM3M "
  89 #endif
  90 #endif
  91 #endif
  92 #endif
  93
  94
  95 #ifdef SMP
  96 #ifndef COMPLEX
  97 #ifdef XDOUBLE
  98 #define MODE    (BLAS_XDOUBLE | BLAS_REAL)
  99 #elif defined(DOUBLE)
 100 #define MODE    (BLAS_DOUBLE  | BLAS_REAL)
 101 #else
 102 #define MODE    (BLAS_SINGLE  | BLAS_REAL)
 103 #endif
 104 #else
 105 #ifdef XDOUBLE
 106 #define MODE    (BLAS_XDOUBLE | BLAS_COMPLEX)
 107 #elif defined(DOUBLE)
 108 #define MODE    (BLAS_DOUBLE  | BLAS_COMPLEX)
 109 #else
 110 #define MODE    (BLAS_SINGLE  | BLAS_COMPLEX)
 111 #endif
 112 #endif
 113 #endif
 114
 115 static int (*symm[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
 116 #ifndef GEMM3M
 117 #ifndef HEMM
 118   SYMM_LU, SYMM_LL, SYMM_RU, SYMM_RL,
 119 #if defined(SMP) && !defined(USE_SIMPLE_THREADED_LEVEL3)
 120   SYMM_THREAD_LU, SYMM_THREAD_LL, SYMM_THREAD_RU, SYMM_THREAD_RL,
 121 #endif
 122 #else
 123   HEMM_LU, HEMM_LL, HEMM_RU, HEMM_RL,
 124 #if defined(SMP) && !defined(USE_SIMPLE_THREADED_LEVEL3)
 125   HEMM_THREAD_LU, HEMM_THREAD_LL, HEMM_THREAD_RU, HEMM_THREAD_RL,
 126 #endif
 127 #endif
 128 #else
 129 #ifndef HEMM
 130   SYMM3M_LU, SYMM3M_LL, SYMM3M_RU, SYMM3M_RL,
 131 #if defined(SMP) && !defined(USE_SIMPLE_THREADED_LEVEL3)
 132   SYMM3M_THREAD_LU, SYMM3M_THREAD_LL, SYMM3M_THREAD_RU, SYMM3M_THREAD_RL,
 133 #endif
 134 #else
 135   HEMM3M_LU, HEMM3M_LL, HEMM3M_RU, HEMM3M_RL,
 136 #if defined(SMP) && !defined(USE_SIMPLE_THREADED_LEVEL3)
 137   HEMM3M_THREAD_LU, HEMM3M_THREAD_LL, HEMM3M_THREAD_RU, HEMM3M_THREAD_RL,
 138 #endif
 139 #endif
 140 #endif
 141 };
 142
 143 #ifndef CBLAS
 144
 145 void NAME(char *SIDE, char *UPLO,
 146          blasint *M, blasint *N,
 147          FLOAT *alpha, FLOAT *a, blasint *ldA,
 148          FLOAT *b, blasint *ldB,
 149          FLOAT *beta,  FLOAT *c, blasint *ldC){
 150
 151   char side_arg  = *SIDE;
 152   char uplo_arg  = *UPLO;
 153
 154   blas_arg_t args;
 155
 156   FLOAT *buffer;
 157   FLOAT *sa, *sb;
 158
 159 #if defined(SMP) && !defined(NO_AFFINITY)
 160   int nodes;
 161 #endif
 162
 163   blasint info;
 164   int side;
 165   int uplo;
 166
 167   PRINT_DEBUG_NAME;
 168
 169   args.alpha = (void *)alpha;
 170   args.beta  = (void *)beta;
 171
 172   TOUPPER(side_arg);
 173   TOUPPER(uplo_arg);
 174
 175   side  = -1;
 176   uplo  = -1;
 177
 178   if (side_arg  == 'L') side  = 0;
 179   if (side_arg  == 'R') side  = 1;
 180
 181   if (uplo_arg  == 'U') uplo  = 0;
 182   if (uplo_arg  == 'L') uplo  = 1;
 183
 184   args.m = *M;
 185   args.n = *N;
 186
 187   args.c = (void *)c;
 188   args.ldc = *ldC;
 189
 190   info = 0;
 191
 192   if (args.ldc < MAX(1, args.m)) info = 12;
 193
 194   if (!side) {
 195     args.a = (void *)a;
 196     args.b = (void *)b;
 197
 198     args.lda = *ldA;
 199     args.ldb = *ldB;
 200
 201     if (args.ldb < MAX(1, args.m)) info =  9;
 202     if (args.lda < MAX(1, args.m)) info =  7;
 203
 204   } else {
 205     args.a = (void *)b;
 206     args.b = (void *)a;
 207
 208     args.lda = *ldB;
 209     args.ldb = *ldA;
 210
 211   if (args.lda < MAX(1, args.m)) info =  9;
 212   if (args.ldb < MAX(1, args.n)) info =  7;
 213   }
 214
 215   if (args.n   < 0)              info =  4;
 216   if (args.m   < 0)              info =  3;
 217   if (uplo     < 0)              info =  2;
 218   if (side     < 0)              info =  1;
 219
 220   if (info != 0) {
 221     BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
 222     return;
 223   }
 224
 225 #else
 226
 227 void CNAME(enum CBLAS_ORDER order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo,
 228            blasint m, blasint n,
 229 #ifndef COMPLEX
 230            FLOAT alpha,
 231            FLOAT *a, blasint lda,
 232            FLOAT *b, blasint ldb,
 233            FLOAT beta,
 234            FLOAT *c, blasint ldc) {
 235 #else
 236            void *valpha,
 237            void *va, blasint lda,
 238            void *vb, blasint ldb,
 239            void *vbeta,
 240            void *vc, blasint ldc) {
 241   FLOAT *alpha = (FLOAT*) valpha;
 242   FLOAT *beta  = (FLOAT*) vbeta;
 243   FLOAT *a = (FLOAT*) va;
 244   FLOAT *b = (FLOAT*) vb;
 245   FLOAT *c = (FLOAT*) vc;
 246 #endif
 247
 248   blas_arg_t args;
 249   int side, uplo;
 250   blasint info;
 251
 252   FLOAT *buffer;
 253   FLOAT *sa, *sb;
 254
 255 #if defined(SMP) && !defined(NO_AFFINITY)
 256   int nodes;
 257 #endif
 258
 259   PRINT_DEBUG_CNAME;
 260
 261 #ifndef COMPLEX
 262   args.alpha = (void *)&alpha;
 263   args.beta  = (void *)&beta;
 264 #else
 265   args.alpha = (void *)alpha;
 266   args.beta  = (void *)beta;
 267 #endif
 268
 269   args.c = (void *)c;
 270   args.ldc = ldc;
 271
 272   side  = -1;
 273   uplo  = -1;
 274   info  =  0;
 275
 276   if (order == CblasColMajor) {
 277     if (Side == CblasLeft)  side = 0;
 278     if (Side == CblasRight) side = 1;
 279
 280     if (Uplo == CblasUpper) uplo  = 0;
 281     if (Uplo == CblasLower) uplo  = 1;
 282
 283     info = -1;
 284
 285     args.m = m;
 286     args.n = n;
 287
 288     if (args.ldc < MAX(1, args.m)) info = 12;
 289
 290     if (!side) {
 291       args.a = (void *)a;
 292       args.b = (void *)b;
 293
 294       args.lda = lda;
 295       args.ldb = ldb;
 296
 297       if (args.ldb < MAX(1, args.m)) info =  9;
 298       if (args.lda < MAX(1, args.m)) info =  7;
 299
 300     } else {
 301       args.a = (void *)b;
 302       args.b = (void *)a;
 303
 304       args.lda = ldb;
 305       args.ldb = lda;
 306
 307       if (args.lda < MAX(1, args.m)) info =  9;
 308       if (args.ldb < MAX(1, args.n)) info =  7;
 309     }
 310
 311     if (args.n   < 0)              info =  4;
 312     if (args.m   < 0)              info =  3;
 313     if (uplo     < 0)              info =  2;
 314     if (side     < 0)              info =  1;
 315   }
 316
 317   if (order == CblasRowMajor) {
 318     if (Side == CblasLeft)  side = 1;
 319     if (Side == CblasRight) side = 0;
 320
 321     if (Uplo == CblasUpper) uplo  = 1;
 322     if (Uplo == CblasLower) uplo  = 0;
 323
 324     info = -1;
 325
 326     args.m = n;
 327     args.n = m;
 328
 329     if (args.ldc < MAX(1, args.m)) info = 12;
 330
 331     if (!side) {
 332       args.a = (void *)a;
 333       args.b = (void *)b;
 334
 335       args.lda = lda;
 336       args.ldb = ldb;
 337
 338       if (args.ldb < MAX(1, args.m)) info =  9;
 339       if (args.lda < MAX(1, args.m)) info =  7;
 340
 341     } else {
 342       args.a = (void *)b;
 343       args.b = (void *)a;
 344
 345       args.lda = ldb;
 346       args.ldb = lda;
 347
 348       if (args.lda < MAX(1, args.m)) info =  9;
 349       if (args.ldb < MAX(1, args.n)) info =  7;
 350     }
 351
 352     if (args.n   < 0)              info =  4;
 353     if (args.m   < 0)              info =  3;
 354     if (uplo     < 0)              info =  2;
 355     if (side     < 0)              info =  1;
 356   }
 357
 358   if (info >= 0) {
 359     BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
 360     return;
 361   }
 362
 363 #endif
 364
 365   if (args.m == 0 || args.n == 0) return;
 366
 367   IDEBUG_START;
 368
 369   FUNCTION_PROFILE_START();
 370
 371   buffer = (FLOAT *)blas_memory_alloc(0);
 372
 373   sa = (FLOAT *)((BLASLONG)buffer + GEMM_OFFSET_A);
 374   sb = (FLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
 375
 376 #ifdef SMP
 377   args.common = NULL;
 378   args.nthreads = num_cpu_avail(3);
 379
 380   if (args.nthreads == 1) {
 381 #endif
 382
 383     (symm[(side << 1) | uplo ])(&args, NULL, NULL, sa, sb, 0);
 384
 385 #ifdef SMP
 386
 387   } else {
 388
 389 #ifndef NO_AFFINITY
 390     nodes = get_num_nodes();
 391
 392     if (nodes > 1) {
 393
 394       args.nthreads /= nodes;
 395
 396       gemm_thread_mn(MODE, &args, NULL, NULL,
 397                      symm[4 | (side << 1) | uplo ], sa, sb, nodes);
 398
 399     } else {
 400 #endif
 401
 402 #ifndef USE_SIMPLE_THREADED_LEVEL3
 403
 404       (symm[4 | (side << 1) | uplo ])(&args, NULL, NULL, sa, sb, 0);
 405
 406 #else
 407
 408       GEMM_THREAD(MODE, &args, NULL, NULL, symm[(side << 1) | uplo ], sa, sb, args.nthreads);
 409
 410 #endif
 411
 412 #ifndef NO_AFFINITY
 413     }
 414 #endif
 415
 416   }
 417 #endif
 418
 419  blas_memory_free(buffer);
 420
 421   FUNCTION_PROFILE_END(COMPSIZE * COMPSIZE,
 422                        (!side)? args.m * (args.m / 2 + args.n) : args.n * (args.m + args.n / 2),
 423                        (!side)? 2 * args.m * args.m * args.n : 2 * args.m * args.n * args.n);
 424
 425   IDEBUG_END;
 426
 427   return;
 428 }