src/gromacs/linearalgebra/gmx_blas/dgemm.cpp

   1 #include <cctype>
   2 #include <cmath>
   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=std::toupper(*transa);
  24   const char trb=std::toupper(*transb);
  25   double temp;
  26   int i,j,l;
  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   double alpha = *alpha__;
  36   double beta  = *beta__;
  37
  38   if(m==0 || n==0 || (( std::abs(alpha)<GMX_DOUBLE_MIN || k==0) && std::abs(beta-1.0)<GMX_DOUBLE_EPS))
  39     return;
  40
  41   if(std::abs(alpha)<GMX_DOUBLE_MIN) {
  42     if(std::abs(beta)<GMX_DOUBLE_MIN) {
  43       for(j=0;j<n;j++)
  44         for(i=0;i<m;i++)
  45           c[j*(ldc)+i] = 0.0;
  46     } else {
  47       /* nonzero beta */
  48       for(j=0;j<n;j++)
  49         for(i=0;i<m;i++)
  50           c[j*(ldc)+i] *= beta;
  51     }
  52     return;
  53   }
  54
  55   if(trb=='N') {
  56     if(tra=='N') {
  57
  58       for(j=0;j<n;j++) {
  59         if(std::abs(beta)<GMX_DOUBLE_MIN) {
  60           for(i=0;i<m;i++)
  61             c[j*(ldc)+i] = 0.0;
  62         } else if(std::abs(beta-1.0)>GMX_DOUBLE_EPS) {
  63           for(i=0;i<m;i++)
  64             c[j*(ldc)+i] *= beta;
  65         }
  66         for(l=0;l<k;l++) {
  67           if( std::abs(b[ j*(ldb) + l ])>GMX_DOUBLE_MIN) {
  68             temp = alpha * b[ j*(ldb) + l ];
  69             for(i=0;i<m;i++)
  70               c[j*(ldc)+i] += temp * a[l*(lda)+i];
  71           }
  72         }
  73       }
  74     } else {
  75       /* transpose A, but not B */
  76       for(j=0;j<n;j++) {
  77         for(i=0;i<m;i++) {
  78           temp = 0.0;
  79           for(l=0;l<k;l++)
  80             temp += a[i*(lda)+l] * b[j*(ldb)+l];
  81           if(std::abs(beta)<GMX_DOUBLE_MIN)
  82             c[j*(ldc)+i] = alpha * temp;
  83           else
  84             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
  85         }
  86       }
  87     }
  88   } else {
  89     /* transpose B */
  90     if(tra=='N') {
  91
  92       /* transpose B, but not A */
  93
  94       for(j=0;j<n;j++) {
  95         if(std::abs(beta)<GMX_DOUBLE_MIN) {
  96           for(i=0;i<m;i++)
  97             c[j*(ldc)+i] = 0.0;
  98         } else if(std::abs(beta-1.0)>GMX_DOUBLE_EPS) {
  99           for(i=0;i<m;i++)
 100             c[j*(ldc)+i] *= beta;
 101         }
 102         for(l=0;l<k;l++) {
 103           if( std::abs(b[ l*(ldb) + j ])>GMX_DOUBLE_MIN) {
 104             temp = alpha * b[ l*(ldb) + j ];
 105             for(i=0;i<m;i++)
 106               c[j*(ldc)+i] += temp * a[l*(lda)+i];
 107           }
 108         }
 109       }
 110
 111     } else {
 112       /* Transpose both A and B */
 113        for(j=0;j<n;j++) {
 114         for(i=0;i<m;i++) {
 115           temp = 0.0;
 116           for(l=0;l<k;l++)
 117             temp += a[i*(lda)+l] * b[l*(ldb)+j];
 118           if(std::abs(beta)<GMX_DOUBLE_MIN)
 119             c[j*(ldc)+i] = alpha * temp;
 120           else
 121             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
 122         }
 123        }
 124     }
 125   }
 126 }