src/gromacs/linearalgebra/gmx_blas/dsymv.c

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