src/gromacs/linearalgebra/gmx_blas/dgemm.c

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