git.maemo.org Git - opencv/blob - 3rdparty/lapack/sgemm.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int sgemm_(char *transa, char *transb, integer *m, integer *
   4         n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *
   5         ldb, real *beta, real *c__, integer *ldc)
   6 {
   7     /* System generated locals */
   8     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
   9             i__3;
  10
  11     /* Local variables */
  12     integer i__, j, l, info;
  13     logical nota, notb;
  14     real temp;
  15     integer ncola;
  16     extern logical lsame_(char *, char *);
  17     integer nrowa, nrowb;
  18     extern /* Subroutine */ int xerbla_(char *, integer *);
  19
  20 /*     .. Scalar Arguments .. */
  21 /*     .. */
  22 /*     .. Array Arguments .. */
  23 /*     .. */
  24
  25 /*  Purpose */
  26 /*  ======= */
  27
  28 /*  SGEMM  performs one of the matrix-matrix operations */
  29
  30 /*     C := alpha*op( A )*op( B ) + beta*C, */
  31
  32 /*  where  op( X ) is one of */
  33
  34 /*     op( X ) = X   or   op( X ) = X', */
  35
  36 /*  alpha and beta are scalars, and A, B and C are matrices, with op( A ) */
  37 /*  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix. */
  38
  39 /*  Arguments */
  40 /*  ========== */
  41
  42 /*  TRANSA - CHARACTER*1. */
  43 /*           On entry, TRANSA specifies the form of op( A ) to be used in */
  44 /*           the matrix multiplication as follows: */
  45
  46 /*              TRANSA = 'N' or 'n',  op( A ) = A. */
  47
  48 /*              TRANSA = 'T' or 't',  op( A ) = A'. */
  49
  50 /*              TRANSA = 'C' or 'c',  op( A ) = A'. */
  51
  52 /*           Unchanged on exit. */
  53
  54 /*  TRANSB - CHARACTER*1. */
  55 /*           On entry, TRANSB specifies the form of op( B ) to be used in */
  56 /*           the matrix multiplication as follows: */
  57
  58 /*              TRANSB = 'N' or 'n',  op( B ) = B. */
  59
  60 /*              TRANSB = 'T' or 't',  op( B ) = B'. */
  61
  62 /*              TRANSB = 'C' or 'c',  op( B ) = B'. */
  63
  64 /*           Unchanged on exit. */
  65
  66 /*  M      - INTEGER. */
  67 /*           On entry,  M  specifies  the number  of rows  of the  matrix */
  68 /*           op( A )  and of the  matrix  C.  M  must  be at least  zero. */
  69 /*           Unchanged on exit. */
  70
  71 /*  N      - INTEGER. */
  72 /*           On entry,  N  specifies the number  of columns of the matrix */
  73 /*           op( B ) and the number of columns of the matrix C. N must be */
  74 /*           at least zero. */
  75 /*           Unchanged on exit. */
  76
  77 /*  K      - INTEGER. */
  78 /*           On entry,  K  specifies  the number of columns of the matrix */
  79 /*           op( A ) and the number of rows of the matrix op( B ). K must */
  80 /*           be at least  zero. */
  81 /*           Unchanged on exit. */
  82
  83 /*  ALPHA  - REAL            . */
  84 /*           On entry, ALPHA specifies the scalar alpha. */
  85 /*           Unchanged on exit. */
  86
  87 /*  A      - REAL             array of DIMENSION ( LDA, ka ), where ka is */
  88 /*           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise. */
  89 /*           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k */
  90 /*           part of the array  A  must contain the matrix  A,  otherwise */
  91 /*           the leading  k by m  part of the array  A  must contain  the */
  92 /*           matrix A. */
  93 /*           Unchanged on exit. */
  94
  95 /*  LDA    - INTEGER. */
  96 /*           On entry, LDA specifies the first dimension of A as declared */
  97 /*           in the calling (sub) program. When  TRANSA = 'N' or 'n' then */
  98 /*           LDA must be at least  max( 1, m ), otherwise  LDA must be at */
  99 /*           least  max( 1, k ). */
 100 /*           Unchanged on exit. */
 101
 102 /*  B      - REAL             array of DIMENSION ( LDB, kb ), where kb is */
 103 /*           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise. */
 104 /*           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n */
 105 /*           part of the array  B  must contain the matrix  B,  otherwise */
 106 /*           the leading  n by k  part of the array  B  must contain  the */
 107 /*           matrix B. */
 108 /*           Unchanged on exit. */
 109
 110 /*  LDB    - INTEGER. */
 111 /*           On entry, LDB specifies the first dimension of B as declared */
 112 /*           in the calling (sub) program. When  TRANSB = 'N' or 'n' then */
 113 /*           LDB must be at least  max( 1, k ), otherwise  LDB must be at */
 114 /*           least  max( 1, n ). */
 115 /*           Unchanged on exit. */
 116
 117 /*  BETA   - REAL            . */
 118 /*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is */
 119 /*           supplied as zero then C need not be set on input. */
 120 /*           Unchanged on exit. */
 121
 122 /*  C      - REAL             array of DIMENSION ( LDC, n ). */
 123 /*           Before entry, the leading  m by n  part of the array  C must */
 124 /*           contain the matrix  C,  except when  beta  is zero, in which */
 125 /*           case C need not be set on entry. */
 126 /*           On exit, the array  C  is overwritten by the  m by n  matrix */
 127 /*           ( alpha*op( A )*op( B ) + beta*C ). */
 128
 129 /*  LDC    - INTEGER. */
 130 /*           On entry, LDC specifies the first dimension of C as declared */
 131 /*           in  the  calling  (sub)  program.   LDC  must  be  at  least */
 132 /*           max( 1, m ). */
 133 /*           Unchanged on exit. */
 134
 135
 136 /*  Level 3 Blas routine. */
 137
 138 /*  -- Written on 8-February-1989. */
 139 /*     Jack Dongarra, Argonne National Laboratory. */
 140 /*     Iain Duff, AERE Harwell. */
 141 /*     Jeremy Du Croz, Numerical Algorithms Group Ltd. */
 142 /*     Sven Hammarling, Numerical Algorithms Group Ltd. */
 143
 144
 145 /*     .. External Functions .. */
 146 /*     .. */
 147 /*     .. External Subroutines .. */
 148 /*     .. */
 149 /*     .. Intrinsic Functions .. */
 150 /*     .. */
 151 /*     .. Local Scalars .. */
 152 /*     .. */
 153 /*     .. Parameters .. */
 154 /*     .. */
 155
 156 /*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not */
 157 /*     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows */
 158 /*     and  columns of  A  and the  number of  rows  of  B  respectively. */
 159
 160     /* Parameter adjustments */
 161     a_dim1 = *lda;
 162     a_offset = 1 + a_dim1;
 163     a -= a_offset;
 164     b_dim1 = *ldb;
 165     b_offset = 1 + b_dim1;
 166     b -= b_offset;
 167     c_dim1 = *ldc;
 168     c_offset = 1 + c_dim1;
 169     c__ -= c_offset;
 170
 171     /* Function Body */
 172     nota = lsame_(transa, "N");
 173     notb = lsame_(transb, "N");
 174     if (nota) {
 175         nrowa = *m;
 176         ncola = *k;
 177     } else {
 178         nrowa = *k;
 179         ncola = *m;
 180     }
 181     if (notb) {
 182         nrowb = *k;
 183     } else {
 184         nrowb = *n;
 185     }
 186
 187 /*     Test the input parameters. */
 188
 189     info = 0;
 190     if (! nota && ! lsame_(transa, "C") && ! lsame_(
 191             transa, "T")) {
 192         info = 1;
 193     } else if (! notb && ! lsame_(transb, "C") && !
 194             lsame_(transb, "T")) {
 195         info = 2;
 196     } else if (*m < 0) {
 197         info = 3;
 198     } else if (*n < 0) {
 199         info = 4;
 200     } else if (*k < 0) {
 201         info = 5;
 202     } else if (*lda < max(1,nrowa)) {
 203         info = 8;
 204     } else if (*ldb < max(1,nrowb)) {
 205         info = 10;
 206     } else if (*ldc < max(1,*m)) {
 207         info = 13;
 208     }
 209     if (info != 0) {
 210         xerbla_("SGEMM ", &info);
 211         return 0;
 212     }
 213
 214 /*     Quick return if possible. */
 215
 216     if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
 217         return 0;
 218     }
 219
 220 /*     And if  alpha.eq.zero. */
 221
 222     if (*alpha == 0.f) {
 223         if (*beta == 0.f) {
 224             i__1 = *n;
 225             for (j = 1; j <= i__1; ++j) {
 226                 i__2 = *m;
 227                 for (i__ = 1; i__ <= i__2; ++i__) {
 228                     c__[i__ + j * c_dim1] = 0.f;
 229 /* L10: */
 230                 }
 231 /* L20: */
 232             }
 233         } else {
 234             i__1 = *n;
 235             for (j = 1; j <= i__1; ++j) {
 236                 i__2 = *m;
 237                 for (i__ = 1; i__ <= i__2; ++i__) {
 238                     c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
 239 /* L30: */
 240                 }
 241 /* L40: */
 242             }
 243         }
 244         return 0;
 245     }
 246
 247 /*     Start the operations. */
 248
 249     if (notb) {
 250         if (nota) {
 251
 252 /*           Form  C := alpha*A*B + beta*C. */
 253
 254             i__1 = *n;
 255             for (j = 1; j <= i__1; ++j) {
 256                 if (*beta == 0.f) {
 257                     i__2 = *m;
 258                     for (i__ = 1; i__ <= i__2; ++i__) {
 259                         c__[i__ + j * c_dim1] = 0.f;
 260 /* L50: */
 261                     }
 262                 } else if (*beta != 1.f) {
 263                     i__2 = *m;
 264                     for (i__ = 1; i__ <= i__2; ++i__) {
 265                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
 266 /* L60: */
 267                     }
 268                 }
 269                 i__2 = *k;
 270                 for (l = 1; l <= i__2; ++l) {
 271                     if (b[l + j * b_dim1] != 0.f) {
 272                         temp = *alpha * b[l + j * b_dim1];
 273                         i__3 = *m;
 274                         for (i__ = 1; i__ <= i__3; ++i__) {
 275                             c__[i__ + j * c_dim1] += temp * a[i__ + l *
 276                                     a_dim1];
 277 /* L70: */
 278                         }
 279                     }
 280 /* L80: */
 281                 }
 282 /* L90: */
 283             }
 284         } else {
 285
 286 /*           Form  C := alpha*A'*B + beta*C */
 287
 288             i__1 = *n;
 289             for (j = 1; j <= i__1; ++j) {
 290                 i__2 = *m;
 291                 for (i__ = 1; i__ <= i__2; ++i__) {
 292                     temp = 0.f;
 293                     i__3 = *k;
 294                     for (l = 1; l <= i__3; ++l) {
 295                         temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
 296 /* L100: */
 297                     }
 298                     if (*beta == 0.f) {
 299                         c__[i__ + j * c_dim1] = *alpha * temp;
 300                     } else {
 301                         c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
 302                                 i__ + j * c_dim1];
 303                     }
 304 /* L110: */
 305                 }
 306 /* L120: */
 307             }
 308         }
 309     } else {
 310         if (nota) {
 311
 312 /*           Form  C := alpha*A*B' + beta*C */
 313
 314             i__1 = *n;
 315             for (j = 1; j <= i__1; ++j) {
 316                 if (*beta == 0.f) {
 317                     i__2 = *m;
 318                     for (i__ = 1; i__ <= i__2; ++i__) {
 319                         c__[i__ + j * c_dim1] = 0.f;
 320 /* L130: */
 321                     }
 322                 } else if (*beta != 1.f) {
 323                     i__2 = *m;
 324                     for (i__ = 1; i__ <= i__2; ++i__) {
 325                         c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
 326 /* L140: */
 327                     }
 328                 }
 329                 i__2 = *k;
 330                 for (l = 1; l <= i__2; ++l) {
 331                     if (b[j + l * b_dim1] != 0.f) {
 332                         temp = *alpha * b[j + l * b_dim1];
 333                         i__3 = *m;
 334                         for (i__ = 1; i__ <= i__3; ++i__) {
 335                             c__[i__ + j * c_dim1] += temp * a[i__ + l *
 336                                     a_dim1];
 337 /* L150: */
 338                         }
 339                     }
 340 /* L160: */
 341                 }
 342 /* L170: */
 343             }
 344         } else {
 345
 346 /*           Form  C := alpha*A'*B' + beta*C */
 347
 348             i__1 = *n;
 349             for (j = 1; j <= i__1; ++j) {
 350                 i__2 = *m;
 351                 for (i__ = 1; i__ <= i__2; ++i__) {
 352                     temp = 0.f;
 353                     i__3 = *k;
 354                     for (l = 1; l <= i__3; ++l) {
 355                         temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
 356 /* L180: */
 357                     }
 358                     if (*beta == 0.f) {
 359                         c__[i__ + j * c_dim1] = *alpha * temp;
 360                     } else {
 361                         c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
 362                                 i__ + j * c_dim1];
 363                     }
 364 /* L190: */
 365                 }
 366 /* L200: */
 367             }
 368         }
 369     }
 370
 371     return 0;
 372
 373 /*     End of SGEMM . */
 374
 375 } /* sgemm_ */