src/gromacs/linearalgebra/gmx_blas/dgemv.c

   1 #include <math.h>
   2 #include <ctype.h>
   3
   4 #include "gromacs/utility/real.h"
   5
   6 #include "../gmx_blas.h"
   7
   8 void
   9 F77_FUNC(dgemv,DGEMV)(const char *trans,
  10        int *m__,
  11        int *n__,
  12        double *alpha__,
  13        double *a,
  14        int *lda__,
  15        double *x,
  16        int *incx__,
  17        double *beta__,
  18        double *y,
  19        int *incy__)
  20 {
  21   const char ch=toupper(*trans);
  22   int lenx,leny,kx,ky;
  23   int i,j,jx,jy,ix,iy;
  24   double temp;
  25
  26   int m = *m__;
  27   int n = *n__;
  28   double alpha = *alpha__;
  29   double beta = *beta__;
  30   int incx = *incx__;
  31   int incy = *incy__;
  32   int lda = *lda__;
  33
  34   if(n<=0 || m<=0 || (fabs(alpha)<GMX_DOUBLE_MIN && fabs(beta-1.0)<GMX_DOUBLE_EPS))
  35     return;
  36
  37   if(ch=='N') {
  38     lenx = n;
  39     leny = m;
  40   } else {
  41     lenx = m;
  42     leny = n;
  43   }
  44
  45    if(incx>0)
  46     kx = 1;
  47   else
  48     kx = 1 - (lenx -1)*(incx);
  49
  50   if(incy>0)
  51     ky = 1;
  52   else
  53     ky = 1 - (leny -1)*(incy);
  54
  55   if(fabs(beta-1.0)>GMX_DOUBLE_EPS) {
  56     if(incy==1) {
  57       if(fabs(beta)<GMX_DOUBLE_MIN)
  58         for(i=0;i<leny;i++)
  59           y[i] = 0.0;
  60       else
  61         for(i=0;i<leny;i++)
  62           y[i] *= beta;
  63     } else {
  64       /* non-unit incr. */
  65       iy = ky;
  66       if(fabs(beta)<GMX_DOUBLE_MIN)
  67         for(i=0;i<leny;i++,iy+=incy)
  68           y[iy] = 0.0;
  69       else
  70         for(i=0;i<leny;i++,iy+=incy)
  71           y[iy] *= beta;
  72     }
  73   }
  74
  75   if(fabs(alpha)<GMX_DOUBLE_MIN)
  76     return;
  77
  78   if(ch=='N') {
  79     jx = kx;
  80     if(incy==1) {
  81       for(j=1;j<=n;j++,jx+=incx)
  82         if(fabs(x[jx-1])>GMX_DOUBLE_MIN) {
  83           temp = alpha * x[jx-1];
  84           for(i=1;i<=m;i++)
  85             y[i-1] += temp * a[(j-1)*(lda)+(i-1)];
  86         }
  87     } else {
  88       /* non-unit y incr. */
  89       for(j=1;j<=n;j++,jx+=incx)
  90         if(fabs(x[jx-1])>GMX_DOUBLE_MIN) {
  91           temp = alpha * x[jx-1];
  92           iy = ky;
  93           for(i=1;i<=m;i++,iy+=incy)
  94             y[iy-1] += temp * a[(j-1)*(lda)+(i-1)];
  95         }
  96     }
  97   } else {
  98     /* transpose */
  99     jy = ky;
 100     if(incx==1) {
 101       for(j=1;j<=n;j++,jy+=incy) {
 102         temp = 0.0;
 103         for(i=1;i<=m;i++)
 104           temp += a[(j-1)*(lda)+(i-1)] * x[i-1];
 105         y[jy-1] += alpha * temp;
 106       }
 107     } else {
 108       /* non-unit y incr. */
 109       for(j=1;j<=n;j++,jy+=incy) {
 110         temp = 0.0;
 111         ix = kx;
 112         for(i=1;i<=m;i++,ix+=incx)
 113           temp += a[(j-1)*(lda)+(i-1)] * x[ix-1];
 114         y[jy-1] += alpha * temp;
 115       }
 116     }
 117   }
 118 }
 119