src/gmxlib/gmx_blas/sgemm.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <ctype.h>
  36 #include <math.h>
  37
  38 #include <types/simple.h>
  39 #include "gmx_blas.h"
  40
  41 void
  42 F77_FUNC(sgemm,SGEMM)(const char *transa,
  43        const char *transb,
  44        int *m__,
  45        int *n__,
  46        int *k__,
  47        float *alpha__,
  48        float *a,
  49        int *lda__,
  50        float *b,
  51        int *ldb__,
  52        float *beta__,
  53        float *c,
  54        int *ldc__)
  55 {
  56   const char tra=toupper(*transa);
  57   const char trb=toupper(*transb);
  58   float temp;
  59   int i,j,l;
  60   int nrowa,ncola,nrowb;
  61
  62   int m   = *m__;
  63   int n   = *n__;
  64   int k   = *k__;
  65   int lda = *lda__;
  66   int ldb = *ldb__;
  67   int ldc = *ldc__;
  68
  69   float alpha = *alpha__;
  70   float beta  = *beta__;
  71
  72   if(tra=='N') {
  73     nrowa = m;
  74     ncola = k;
  75   } else {
  76     nrowa = k;
  77     ncola = m;
  78   }
  79
  80   if(trb=='N')
  81     nrowb = k;
  82    else
  83     nrowb = n;
  84
  85   if(m==0 || n==0 || (( fabs(alpha)<GMX_FLOAT_MIN || k==0) && fabs(beta-1.0)<GMX_FLOAT_EPS))
  86     return;
  87
  88   if(fabs(alpha)<GMX_FLOAT_MIN) {
  89     if(fabs(beta)<GMX_FLOAT_MIN) {
  90       for(j=0;j<n;j++)
  91         for(i=0;i<m;i++)
  92           c[j*(ldc)+i] = 0.0;
  93     } else {
  94       /* nonzero beta */
  95       for(j=0;j<n;j++)
  96         for(i=0;i<m;i++)
  97           c[j*(ldc)+i] *= beta;
  98     }
  99     return;
 100   }
 101
 102   if(trb=='N') {
 103     if(tra=='N') {
 104
 105       for(j=0;j<n;j++) {
 106         if(fabs(beta)<GMX_FLOAT_MIN) {
 107           for(i=0;i<m;i++)
 108             c[j*(ldc)+i] = 0.0;
 109         } else if(fabs(beta-1.0)>GMX_FLOAT_EPS) {
 110           for(i=0;i<m;i++)
 111             c[j*(ldc)+i] *= beta;
 112         }
 113         for(l=0;l<k;l++) {
 114           if( fabs(b[ j*(ldb) + l ])>GMX_FLOAT_MIN) {
 115             temp = alpha * b[ j*(ldb) + l ];
 116             for(i=0;i<m;i++)
 117               c[j*(ldc)+i] += temp * a[l*(lda)+i];
 118           }
 119         }
 120       }
 121     } else {
 122       /* transpose A, but not B */
 123       for(j=0;j<n;j++) {
 124         for(i=0;i<m;i++) {
 125           temp = 0.0;
 126           for(l=0;l<k;l++)
 127             temp += a[i*(lda)+l] * b[j*(ldb)+l];
 128           if(fabs(beta)<GMX_FLOAT_MIN)
 129             c[j*(ldc)+i] = alpha * temp;
 130           else
 131             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
 132         }
 133       }
 134     }
 135   } else {
 136     /* transpose B */
 137     if(tra=='N') {
 138
 139       /* transpose B, but not A */
 140
 141       for(j=0;j<n;j++) {
 142         if(fabs(beta)<GMX_FLOAT_MIN) {
 143           for(i=0;i<m;i++)
 144             c[j*(ldc)+i] = 0.0;
 145         } else if(fabs(beta-1.0)>GMX_FLOAT_EPS) {
 146           for(i=0;i<m;i++)
 147             c[j*(ldc)+i] *= beta;
 148         }
 149         for(l=0;l<k;l++) {
 150           if( fabs(b[ l*(ldb) + j ])>GMX_FLOAT_MIN) {
 151             temp = alpha * b[ l*(ldb) + j ];
 152             for(i=0;i<m;i++)
 153               c[j*(ldc)+i] += temp * a[l*(lda)+i];
 154           }
 155         }
 156       }
 157
 158     } else {
 159       /* Transpose both A and B */
 160        for(j=0;j<n;j++) {
 161         for(i=0;i<m;i++) {
 162           temp = 0.0;
 163           for(l=0;l<k;l++)
 164             temp += a[i*(lda)+l] * b[l*(ldb)+j];
 165           if(fabs(beta)<GMX_FLOAT_MIN)
 166             c[j*(ldc)+i] = alpha * temp;
 167           else
 168             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
 169         }
 170        }
 171     }
 172   }
 173 }