interface/trsm.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 TRMM
  47 #ifndef COMPLEX
  48 #ifdef XDOUBLE
  49 #define ERROR_NAME "QTRSM "
  50 #elif defined(DOUBLE)
  51 #define ERROR_NAME "DTRSM "
  52 #else
  53 #define ERROR_NAME "STRSM "
  54 #endif
  55 #else
  56 #ifdef XDOUBLE
  57 #define ERROR_NAME "XTRSM "
  58 #elif defined(DOUBLE)
  59 #define ERROR_NAME "ZTRSM "
  60 #else
  61 #define ERROR_NAME "CTRSM "
  62 #endif
  63 #endif
  64 #else
  65 #ifndef COMPLEX
  66 #ifdef XDOUBLE
  67 #define ERROR_NAME "QTRMM "
  68 #elif defined(DOUBLE)
  69 #define ERROR_NAME "DTRMM "
  70 #else
  71 #define ERROR_NAME "STRMM "
  72 #endif
  73 #else
  74 #ifdef XDOUBLE
  75 #define ERROR_NAME "XTRMM "
  76 #elif defined(DOUBLE)
  77 #define ERROR_NAME "ZTRMM "
  78 #else
  79 #define ERROR_NAME "CTRMM "
  80 #endif
  81 #endif
  82 #endif
  83
  84 static int (*trsm[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
  85 #ifndef TRMM
  86   TRSM_LNUU, TRSM_LNUN, TRSM_LNLU, TRSM_LNLN,
  87   TRSM_LTUU, TRSM_LTUN, TRSM_LTLU, TRSM_LTLN,
  88   TRSM_LRUU, TRSM_LRUN, TRSM_LRLU, TRSM_LRLN,
  89   TRSM_LCUU, TRSM_LCUN, TRSM_LCLU, TRSM_LCLN,
  90   TRSM_RNUU, TRSM_RNUN, TRSM_RNLU, TRSM_RNLN,
  91   TRSM_RTUU, TRSM_RTUN, TRSM_RTLU, TRSM_RTLN,
  92   TRSM_RRUU, TRSM_RRUN, TRSM_RRLU, TRSM_RRLN,
  93   TRSM_RCUU, TRSM_RCUN, TRSM_RCLU, TRSM_RCLN,
  94 #else
  95   TRMM_LNUU, TRMM_LNUN, TRMM_LNLU, TRMM_LNLN,
  96   TRMM_LTUU, TRMM_LTUN, TRMM_LTLU, TRMM_LTLN,
  97   TRMM_LRUU, TRMM_LRUN, TRMM_LRLU, TRMM_LRLN,
  98   TRMM_LCUU, TRMM_LCUN, TRMM_LCLU, TRMM_LCLN,
  99   TRMM_RNUU, TRMM_RNUN, TRMM_RNLU, TRMM_RNLN,
 100   TRMM_RTUU, TRMM_RTUN, TRMM_RTLU, TRMM_RTLN,
 101   TRMM_RRUU, TRMM_RRUN, TRMM_RRLU, TRMM_RRLN,
 102   TRMM_RCUU, TRMM_RCUN, TRMM_RCLU, TRMM_RCLN,
 103 #endif
 104 };
 105
 106 #ifndef CBLAS
 107
 108 void NAME(char *SIDE, char *UPLO, char *TRANS, char *DIAG,
 109            blasint *M, blasint *N, FLOAT *alpha,
 110            FLOAT *a, blasint *ldA, FLOAT *b, blasint *ldB){
 111
 112   char side_arg  = *SIDE;
 113   char uplo_arg  = *UPLO;
 114   char trans_arg = *TRANS;
 115   char diag_arg  = *DIAG;
 116
 117   blas_arg_t args;
 118
 119   FLOAT *buffer;
 120   FLOAT *sa, *sb;
 121
 122 #ifdef SMP
 123 #ifndef COMPLEX
 124 #ifdef XDOUBLE
 125   int mode  =  BLAS_XDOUBLE | BLAS_REAL;
 126 #elif defined(DOUBLE)
 127   int mode  =  BLAS_DOUBLE  | BLAS_REAL;
 128 #else
 129   int mode  =  BLAS_SINGLE  | BLAS_REAL;
 130 #endif
 131 #else
 132 #ifdef XDOUBLE
 133   int mode  =  BLAS_XDOUBLE | BLAS_COMPLEX;
 134 #elif defined(DOUBLE)
 135   int mode  =  BLAS_DOUBLE  | BLAS_COMPLEX;
 136 #else
 137   int mode  =  BLAS_SINGLE  | BLAS_COMPLEX;
 138 #endif
 139 #endif
 140 #endif
 141
 142   blasint info;
 143   int side;
 144   int uplo;
 145   int unit;
 146   int trans;
 147   int nrowa;
 148
 149   PRINT_DEBUG_NAME;
 150
 151   args.m = *M;
 152   args.n = *N;
 153
 154   args.a = (void *)a;
 155   args.b = (void *)b;
 156
 157   args.lda = *ldA;
 158   args.ldb = *ldB;
 159
 160   args.beta = (void *)alpha;
 161
 162   TOUPPER(side_arg);
 163   TOUPPER(uplo_arg);
 164   TOUPPER(trans_arg);
 165   TOUPPER(diag_arg);
 166
 167   side  = -1;
 168   trans = -1;
 169   unit  = -1;
 170   uplo  = -1;
 171
 172   if (side_arg  == 'L') side  = 0;
 173   if (side_arg  == 'R') side  = 1;
 174
 175   if (trans_arg == 'N') trans = 0;
 176   if (trans_arg == 'T') trans = 1;
 177   if (trans_arg == 'R') trans = 2;
 178   if (trans_arg == 'C') trans = 3;
 179
 180   if (diag_arg  == 'U') unit  = 0;
 181   if (diag_arg  == 'N') unit  = 1;
 182
 183   if (uplo_arg  == 'U') uplo  = 0;
 184   if (uplo_arg  == 'L') uplo  = 1;
 185
 186   nrowa = args.m;
 187   if (side & 1) nrowa = args.n;
 188
 189   info = 0;
 190
 191   if (args.ldb < MAX(1,args.m)) info = 11;
 192   if (args.lda < MAX(1,nrowa))  info =  9;
 193   if (args.n < 0)               info =  6;
 194   if (args.m < 0)               info =  5;
 195   if (unit < 0)                 info =  4;
 196   if (trans < 0)                info =  3;
 197   if (uplo  < 0)                info =  2;
 198   if (side  < 0)                info =  1;
 199
 200   if (info != 0) {
 201     BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
 202     return;
 203   }
 204
 205 #else
 206
 207 void CNAME(enum CBLAS_ORDER order,
 208            enum CBLAS_SIDE Side,  enum CBLAS_UPLO Uplo,
 209            enum CBLAS_TRANSPOSE Trans, enum CBLAS_DIAG Diag,
 210            blasint m, blasint n,
 211 #ifndef COMPLEX
 212            FLOAT alpha,
 213            FLOAT *a, blasint lda,
 214            FLOAT *b, blasint ldb) {
 215 #else
 216            void *valpha,
 217            void *va, blasint lda,
 218            void *vb, blasint ldb) {
 219   FLOAT *alpha = (FLOAT*) valpha;
 220   FLOAT *a = (FLOAT*) va;
 221   FLOAT *b = (FLOAT*) vb;
 222 #endif
 223
 224   blas_arg_t args;
 225   int side, uplo, trans, unit;
 226   blasint info, nrowa;
 227
 228   XFLOAT *buffer;
 229   XFLOAT *sa, *sb;
 230
 231 #ifdef SMP
 232 #ifndef COMPLEX
 233 #ifdef XDOUBLE
 234   int mode  =  BLAS_XDOUBLE | BLAS_REAL;
 235 #elif defined(DOUBLE)
 236   int mode  =  BLAS_DOUBLE  | BLAS_REAL;
 237 #else
 238   int mode  =  BLAS_SINGLE  | BLAS_REAL;
 239 #endif
 240 #else
 241 #ifdef XDOUBLE
 242   int mode  =  BLAS_XDOUBLE | BLAS_COMPLEX;
 243 #elif defined(DOUBLE)
 244   int mode  =  BLAS_DOUBLE  | BLAS_COMPLEX;
 245 #else
 246   int mode  =  BLAS_SINGLE  | BLAS_COMPLEX;
 247 #endif
 248 #endif
 249 #endif
 250
 251   PRINT_DEBUG_CNAME;
 252
 253   args.a = (void *)a;
 254   args.b = (void *)b;
 255
 256   args.lda = lda;
 257   args.ldb = ldb;
 258
 259 #ifndef COMPLEX
 260   args.beta = (void *)&alpha;
 261 #else
 262   args.beta = (void *)alpha;
 263 #endif
 264
 265   side   = -1;
 266   uplo   = -1;
 267   trans  = -1;
 268   unit   = -1;
 269   info   =  0;
 270
 271   if (order == CblasColMajor) {
 272     args.m = m;
 273     args.n = n;
 274
 275     if (Side == CblasLeft)         side  = 0;
 276     if (Side == CblasRight)        side  = 1;
 277
 278     if (Uplo == CblasUpper)        uplo  = 0;
 279     if (Uplo == CblasLower)        uplo  = 1;
 280
 281     if (Trans == CblasNoTrans)     trans = 0;
 282     if (Trans == CblasTrans)       trans = 1;
 283 #ifndef COMPLEX
 284     if (Trans == CblasConjNoTrans) trans = 0;
 285     if (Trans == CblasConjTrans)   trans = 1;
 286 #else
 287     if (Trans == CblasConjNoTrans) trans = 2;
 288     if (Trans == CblasConjTrans)   trans = 3;
 289 #endif
 290
 291     if (Diag == CblasUnit)          unit  = 0;
 292     if (Diag == CblasNonUnit)       unit  = 1;
 293
 294     info = -1;
 295
 296     nrowa = args.m;
 297     if (side & 1) nrowa = args.n;
 298
 299     if (args.ldb < MAX(1,args.m)) info = 11;
 300     if (args.lda < MAX(1,nrowa))  info =  9;
 301     if (args.n < 0)               info =  6;
 302     if (args.m < 0)               info =  5;
 303     if (unit < 0)                 info =  4;
 304     if (trans < 0)                info =  3;
 305     if (uplo  < 0)                info =  2;
 306     if (side  < 0)                info =  1;
 307   }
 308
 309   if (order == CblasRowMajor) {
 310     args.m = n;
 311     args.n = m;
 312
 313     if (Side == CblasLeft)         side  = 1;
 314     if (Side == CblasRight)        side  = 0;
 315
 316     if (Uplo == CblasUpper)        uplo  = 1;
 317     if (Uplo == CblasLower)        uplo  = 0;
 318
 319     if (Trans == CblasNoTrans)     trans = 0;
 320     if (Trans == CblasTrans)       trans = 1;
 321 #ifndef COMPLEX
 322     if (Trans == CblasConjNoTrans) trans = 0;
 323     if (Trans == CblasConjTrans)   trans = 1;
 324 #else
 325     if (Trans == CblasConjNoTrans) trans = 2;
 326     if (Trans == CblasConjTrans)   trans = 3;
 327 #endif
 328
 329     if (Diag == CblasUnit)         unit  = 0;
 330     if (Diag == CblasNonUnit)      unit  = 1;
 331
 332     info = -1;
 333
 334     nrowa = args.m;
 335     if (side & 1) nrowa = args.n;
 336
 337     if (args.ldb < MAX(1,args.m)) info = 11;
 338     if (args.lda < MAX(1,nrowa))  info =  9;
 339     if (args.n < 0)               info =  6;
 340     if (args.m < 0)               info =  5;
 341     if (unit < 0)                 info =  4;
 342     if (trans < 0)                info =  3;
 343     if (uplo  < 0)                info =  2;
 344     if (side  < 0)                info =  1;
 345   }
 346
 347   if (info >= 0) {
 348     BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
 349     return;
 350   }
 351
 352 #endif
 353
 354   if ((args.m == 0) || (args.n == 0)) return;
 355
 356   IDEBUG_START;
 357
 358   FUNCTION_PROFILE_START();
 359
 360   buffer = (FLOAT *)blas_memory_alloc(0);
 361
 362   sa = (FLOAT *)((BLASLONG)buffer + GEMM_OFFSET_A);
 363   sb = (FLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
 364
 365 #ifdef SMP
 366   mode |= (trans << BLAS_TRANSA_SHIFT);
 367   mode |= (side  << BLAS_RSIDE_SHIFT);
 368
 369   if ( args.m < 2*GEMM_MULTITHREAD_THRESHOLD )
 370         args.nthreads = 1;
 371   else
 372         if ( args.n < 2*GEMM_MULTITHREAD_THRESHOLD )
 373                 args.nthreads = 1;
 374   else
 375         args.nthreads = num_cpu_avail(3);
 376
 377
 378   if (args.nthreads == 1) {
 379 #endif
 380
 381     (trsm[(side<<4) | (trans<<2) | (uplo<<1) | unit])(&args, NULL, NULL, sa, sb, 0);
 382
 383 #ifdef SMP
 384   } else {
 385     if (!side) {
 386       gemm_thread_n(mode, &args, NULL, NULL, trsm[(side<<4) | (trans<<2) | (uplo<<1) | unit], sa, sb, args.nthreads);
 387     } else {
 388       gemm_thread_m(mode, &args, NULL, NULL, trsm[(side<<4) | (trans<<2) | (uplo<<1) | unit], sa, sb, args.nthreads);
 389     }
 390   }
 391 #endif
 392
 393   blas_memory_free(buffer);
 394
 395   FUNCTION_PROFILE_END(COMPSIZE * COMPSIZE,
 396                        (!side) ? args.m * (args.m + args.n) : args.n * (args.m + args.n),
 397                        (!side) ? args.m * args.m * args.n : args.m * args.n * args.n);
 398
 399   IDEBUG_END;
 400
 401   return;
 402 }
 403