src/gmxlib/gmx_blas/dgemm.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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
  40 #include "gmx_blas.h"
  41
  42 void
  43 F77_FUNC(dgemm,DGEMM)(const char *transa,
  44                       const char *transb,
  45                       int *m__,
  46                       int *n__,
  47                       int *k__,
  48                       double *alpha__,
  49                       double *a,
  50                       int *lda__,
  51                       double *b,
  52                       int *ldb__,
  53                       double *beta__,
  54                       double *c,
  55                       int *ldc__)
  56 {
  57   const char tra=toupper(*transa);
  58   const char trb=toupper(*transb);
  59   double temp;
  60   int i,j,l;
  61   int nrowa,ncola,nrowb;
  62
  63   int m = *m__;
  64   int n = *n__;
  65   int k = *k__;
  66   int lda = *lda__;
  67   int ldb = *ldb__;
  68   int ldc = *ldc__;
  69
  70   double alpha = *alpha__;
  71   double beta  = *beta__;
  72
  73   if(tra=='N') {
  74     nrowa = m;
  75     ncola = k;
  76   } else {
  77     nrowa = k;
  78     ncola = m;
  79   }
  80
  81   if(trb=='N')
  82     nrowb = k;
  83    else
  84     nrowb = n;
  85
  86   if(m==0 || n==0 || (( fabs(alpha)<GMX_DOUBLE_MIN || k==0) && fabs(beta-1.0)<GMX_DOUBLE_EPS))
  87     return;
  88
  89   if(fabs(alpha)<GMX_DOUBLE_MIN) {
  90     if(fabs(beta)<GMX_DOUBLE_MIN) {
  91       for(j=0;j<n;j++)
  92         for(i=0;i<m;i++)
  93           c[j*(ldc)+i] = 0.0;
  94     } else {
  95       /* nonzero beta */
  96       for(j=0;j<n;j++)
  97         for(i=0;i<m;i++)
  98           c[j*(ldc)+i] *= beta;
  99     }
 100     return;
 101   }
 102
 103   if(trb=='N') {
 104     if(tra=='N') {
 105
 106       for(j=0;j<n;j++) {
 107         if(fabs(beta)<GMX_DOUBLE_MIN) {
 108           for(i=0;i<m;i++)
 109             c[j*(ldc)+i] = 0.0;
 110         } else if(fabs(beta-1.0)>GMX_DOUBLE_EPS) {
 111           for(i=0;i<m;i++)
 112             c[j*(ldc)+i] *= beta;
 113         }
 114         for(l=0;l<k;l++) {
 115           if( fabs(b[ j*(ldb) + l ])>GMX_DOUBLE_MIN) {
 116             temp = alpha * b[ j*(ldb) + l ];
 117             for(i=0;i<m;i++)
 118               c[j*(ldc)+i] += temp * a[l*(lda)+i];
 119           }
 120         }
 121       }
 122     } else {
 123       /* transpose A, but not B */
 124       for(j=0;j<n;j++) {
 125         for(i=0;i<m;i++) {
 126           temp = 0.0;
 127           for(l=0;l<k;l++)
 128             temp += a[i*(lda)+l] * b[j*(ldb)+l];
 129           if(fabs(beta)<GMX_DOUBLE_MIN)
 130             c[j*(ldc)+i] = alpha * temp;
 131           else
 132             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
 133         }
 134       }
 135     }
 136   } else {
 137     /* transpose B */
 138     if(tra=='N') {
 139
 140       /* transpose B, but not A */
 141
 142       for(j=0;j<n;j++) {
 143         if(fabs(beta)<GMX_DOUBLE_MIN) {
 144           for(i=0;i<m;i++)
 145             c[j*(ldc)+i] = 0.0;
 146         } else if(fabs(beta-1.0)>GMX_DOUBLE_EPS) {
 147           for(i=0;i<m;i++)
 148             c[j*(ldc)+i] *= beta;
 149         }
 150         for(l=0;l<k;l++) {
 151           if( fabs(b[ l*(ldb) + j ])>GMX_DOUBLE_MIN) {
 152             temp = alpha * b[ l*(ldb) + j ];
 153             for(i=0;i<m;i++)
 154               c[j*(ldc)+i] += temp * a[l*(lda)+i];
 155           }
 156         }
 157       }
 158
 159     } else {
 160       /* Transpose both A and B */
 161        for(j=0;j<n;j++) {
 162         for(i=0;i<m;i++) {
 163           temp = 0.0;
 164           for(l=0;l<k;l++)
 165             temp += a[i*(lda)+l] * b[l*(ldb)+j];
 166           if(fabs(beta)<GMX_DOUBLE_MIN)
 167             c[j*(ldc)+i] = alpha * temp;
 168           else
 169             c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
 170         }
 171        }
 172     }
 173   }
 174 }