src/gmxlib/gmx_blas/dgemv.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 <math.h>
  36 #include <ctype.h>
  37
  38 #include <types/simple.h>
  39
  40 #include "gmx_blas.h"
  41
  42 void
  43 F77_FUNC(dgemv,DGEMV)(const char *trans,
  44        int *m__,
  45        int *n__,
  46        double *alpha__,
  47        double *a,
  48        int *lda__,
  49        double *x,
  50        int *incx__,
  51        double *beta__,
  52        double *y,
  53        int *incy__)
  54 {
  55   const char ch=toupper(*trans);
  56   int lenx,leny,kx,ky;
  57   int i,j,jx,jy,ix,iy;
  58   double temp;
  59
  60   int m = *m__;
  61   int n = *n__;
  62   double alpha = *alpha__;
  63   double beta = *beta__;
  64   int incx = *incx__;
  65   int incy = *incy__;
  66   int lda = *lda__;
  67
  68   if(n<=0 || m<=0 || (fabs(alpha)<GMX_DOUBLE_MIN && fabs(beta-1.0)<GMX_DOUBLE_EPS))
  69     return;
  70
  71   if(ch=='N') {
  72     lenx = n;
  73     leny = m;
  74   } else {
  75     lenx = m;
  76     leny = n;
  77   }
  78
  79    if(incx>0)
  80     kx = 1;
  81   else
  82     kx = 1 - (lenx -1)*(incx);
  83
  84   if(incy>0)
  85     ky = 1;
  86   else
  87     ky = 1 - (leny -1)*(incy);
  88
  89   if(fabs(beta-1.0)>GMX_DOUBLE_EPS) {
  90     if(incy==1) {
  91       if(fabs(beta)<GMX_DOUBLE_MIN)
  92         for(i=0;i<leny;i++)
  93           y[i] = 0.0;
  94       else
  95         for(i=0;i<leny;i++)
  96           y[i] *= beta;
  97     } else {
  98       /* non-unit incr. */
  99       iy = ky;
 100       if(fabs(beta)<GMX_DOUBLE_MIN)
 101         for(i=0;i<leny;i++,iy+=incy)
 102           y[iy] = 0.0;
 103       else
 104         for(i=0;i<leny;i++,iy+=incy)
 105           y[iy] *= beta;
 106     }
 107   }
 108
 109   if(fabs(alpha)<GMX_DOUBLE_MIN)
 110     return;
 111
 112   if(ch=='N') {
 113     jx = kx;
 114     if(incy==1) {
 115       for(j=1;j<=n;j++,jx+=incx)
 116         if(fabs(x[jx-1])>GMX_DOUBLE_MIN) {
 117           temp = alpha * x[jx-1];
 118           for(i=1;i<=m;i++)
 119             y[i-1] += temp * a[(j-1)*(lda)+(i-1)];
 120         }
 121     } else {
 122       /* non-unit y incr. */
 123       for(j=1;j<=n;j++,jx+=incx)
 124         if(fabs(x[jx-1])>GMX_DOUBLE_MIN) {
 125           temp = alpha * x[jx-1];
 126           iy = ky;
 127           for(i=1;i<=m;i++,iy+=incy)
 128             y[iy-1] += temp * a[(j-1)*(lda)+(i-1)];
 129         }
 130     }
 131   } else {
 132     /* transpose */
 133     jy = ky;
 134     if(incx==1) {
 135       for(j=1;j<=n;j++,jy+=incy) {
 136         temp = 0.0;
 137         for(i=1;i<=m;i++)
 138           temp += a[(j-1)*(lda)+(i-1)] * x[i-1];
 139         y[jy-1] += alpha * temp;
 140       }
 141     } else {
 142       /* non-unit y incr. */
 143       for(j=1;j<=n;j++,jy+=incy) {
 144         temp = 0.0;
 145         ix = kx;
 146         for(i=1;i<=m;i++,ix+=incx)
 147           temp += a[(j-1)*(lda)+(i-1)] * x[ix-1];
 148         y[jy-1] += alpha * temp;
 149       }
 150     }
 151   }
 152 }
 153