src/gromacs/linearalgebra/gmx_blas/sgemm.c

   1 #include <ctype.h>
   2 #include <math.h>
   3
   4 #include "types/simple.h"
   5 #include "../gmx_blas.h"
   6
   7 void
   8 F77_FUNC(sgemm,SGEMM)(const char *transa,
   9        const char *transb,
  10        int *m__,
  11        int *n__,
  12        int *k__,
  13        float *alpha__,
  14        float *a,
  15        int *lda__,
  16        float *b,
  17        int *ldb__,
  18        float *beta__,
  19        float *c,
  20        int *ldc__)
  21 {
  22   const char tra=toupper(*transa);
  23   const char trb=toupper(*transb);
  24   float temp;
  25   int i,j,l;
  26   int nrowa,ncola,nrowb;
  27
  28   int m   = *m__;
  29   int n   = *n__;
  30   int k   = *k__;
  31   int lda = *lda__;
  32   int ldb = *ldb__;
  33   int ldc = *ldc__;
  34
  35   float alpha = *alpha__;
  36   float beta  = *beta__;
  37
  38   if(tra=='N') {
  39     nrowa = m;
  40     ncola = k;
  41   } else {
  42     nrowa = k;
  43     ncola = m;
  44   }
  45
  46   if(trb=='N')
  47     nrowb = k;
  48    else
  49     nrowb = n;
  50
  51   if(m==0 || n==0 || (( fabs(alpha)<GMX_FLOAT_MIN || k==0) && fabs(beta-1.0)<GMX_FLOAT_EPS))
  52     return;
  53
  54   if(fabs(alpha)<GMX_FLOAT_MIN) {
  55     if(fabs(beta)<GMX_FLOAT_MIN) {
  56       for(j=0;j<n;j++)
  57         for(i=0;i<m;i++)
  58           c[j*(ldc)+i] = 0.0;
  59     } else {
  60       /* nonzero beta */
  61       for(j=0;j<n;j++)
  62         for(i=0;i<m;i++)
  63           c[j*(ldc)+i] *= beta;
  64     }
  65     return;
  66   }
  67
  68   if(trb=='N') {
  69     if(tra=='N') {
  70
  71       for(j=0;j<n;j++) {
  72         if(fabs(beta)<GMX_FLOAT_MIN) {
  73           for(i=0;i<m;i++)
  74             c[j*(ldc)+i] = 0.0;
  75         } else if(fabs(beta-1.0)>GMX_FLOAT_EPS) {
  76           for(i=0;i<m;i++)
  77             c[j*(ldc)+i] *= beta;
  78         }
  79         for(l=0;l<k;l++) {
  80           if( fabs(b[ j*(ldb) + l ])>GMX_FLOAT_MIN) {
  81             temp = alpha * b[ j*(ldb) + l ];
  82             for(i=0;i<m;i++)
  83               c[j*(ldc)+i] += temp * a[l*(lda)+i];
  84           }
  85         }
  86       }
  87     } else {
  88       /* transpose A, but not B */
  89       for(j=0;j<n;j++) {
  90         for(i=0;i<m;i++) {
  91           temp = 0.0;
  92           for(l=0;l<k;l++)
  93             temp += a[i*(lda)+l] * b[j*(ldb)+l];
  94           if(fabs(beta)<GMX_FLOAT_MIN)
  95             c[j*(ldc)+i] = alpha * temp;
  96           else
  97             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
  98         }
  99       }
 100     }
 101   } else {
 102     /* transpose B */
 103     if(tra=='N') {
 104
 105       /* transpose B, but not A */
 106
 107       for(j=0;j<n;j++) {
 108         if(fabs(beta)<GMX_FLOAT_MIN) {
 109           for(i=0;i<m;i++)
 110             c[j*(ldc)+i] = 0.0;
 111         } else if(fabs(beta-1.0)>GMX_FLOAT_EPS) {
 112           for(i=0;i<m;i++)
 113             c[j*(ldc)+i] *= beta;
 114         }
 115         for(l=0;l<k;l++) {
 116           if( fabs(b[ l*(ldb) + j ])>GMX_FLOAT_MIN) {
 117             temp = alpha * b[ l*(ldb) + j ];
 118             for(i=0;i<m;i++)
 119               c[j*(ldc)+i] += temp * a[l*(lda)+i];
 120           }
 121         }
 122       }
 123
 124     } else {
 125       /* Transpose both A and B */
 126        for(j=0;j<n;j++) {
 127         for(i=0;i<m;i++) {
 128           temp = 0.0;
 129           for(l=0;l<k;l++)
 130             temp += a[i*(lda)+l] * b[l*(ldb)+j];
 131           if(fabs(beta)<GMX_FLOAT_MIN)
 132             c[j*(ldc)+i] = alpha * temp;
 133           else
 134             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
 135         }
 136        }
 137     }
 138   }
 139 }