src/gromacs/linearalgebra/gmx_blas/dtrmv.c

   1 #include <math.h>
   2
   3 #include "gromacs/utility/real.h"
   4 #include "../gmx_blas.h"
   5
   6 void
   7 F77_FUNC(dtrmv,DTRMV)(const char *uplo,
   8        const char *trans,
   9        const char *diag,
  10        int *n__,
  11        double *a,
  12        int *lda__,
  13        double *x,
  14        int *incx__)
  15 {
  16     int a_dim1, a_offset, i__1, i__2;
  17
  18     int i__, j, ix, jx, kx, info;
  19     double temp;
  20     int nounit;
  21
  22     int n = *n__;
  23     int lda = *lda__;
  24     int incx = *incx__;
  25
  26     a_dim1 = lda;
  27     a_offset = 1 + a_dim1;
  28     a -= a_offset;
  29     --x;
  30
  31     info = 0;
  32
  33     if (n == 0) {
  34         return;
  35     }
  36
  37     nounit = (*diag=='n' || *diag=='N');
  38
  39     if (incx <= 0) {
  40         kx = 1 - (n - 1) * incx;
  41     } else {
  42         kx = 1;
  43     }
  44
  45     if (*trans=='N' || *trans=='n') {
  46
  47         if (*uplo=='U' || *uplo=='u') {
  48             if (incx == 1) {
  49                 i__1 = n;
  50                 for (j = 1; j <= i__1; ++j) {
  51                     if (fabs(x[j])>GMX_DOUBLE_MIN) {
  52                         temp = x[j];
  53                         i__2 = j - 1;
  54                         for (i__ = 1; i__ <= i__2; ++i__) {
  55                             x[i__] += temp * a[i__ + j * a_dim1];
  56                         }
  57                         if (nounit) {
  58                             x[j] *= a[j + j * a_dim1];
  59                         }
  60                     }
  61                 }
  62             } else {
  63                 jx = kx;
  64                 i__1 = n;
  65                 for (j = 1; j <= i__1; ++j) {
  66                     if (fabs(x[jx])>GMX_DOUBLE_MIN) {
  67                         temp = x[jx];
  68                         ix = kx;
  69                         i__2 = j - 1;
  70                         for (i__ = 1; i__ <= i__2; ++i__) {
  71                             x[ix] += temp * a[i__ + j * a_dim1];
  72                             ix += incx;
  73                         }
  74                         if (nounit) {
  75                             x[jx] *= a[j + j * a_dim1];
  76                         }
  77                     }
  78                     jx += incx;
  79                 }
  80             }
  81         } else {
  82             if (incx == 1) {
  83                 for (j = n; j >= 1; --j) {
  84                     if (fabs(x[j])>GMX_DOUBLE_MIN) {
  85                         temp = x[j];
  86                         i__1 = j + 1;
  87                         for (i__ = n; i__ >= i__1; --i__) {
  88                             x[i__] += temp * a[i__ + j * a_dim1];
  89                         }
  90                         if (nounit) {
  91                             x[j] *= a[j + j * a_dim1];
  92                         }
  93                     }
  94                 }
  95             } else {
  96                 kx += (n - 1) * incx;
  97                 jx = kx;
  98                 for (j = n; j >= 1; --j) {
  99                     if (fabs(x[jx])>GMX_DOUBLE_MIN) {
 100                         temp = x[jx];
 101                         ix = kx;
 102                         i__1 = j + 1;
 103                         for (i__ = n; i__ >= i__1; --i__) {
 104                             x[ix] += temp * a[i__ + j * a_dim1];
 105                             ix -= incx;
 106                         }
 107                         if (nounit) {
 108                             x[jx] *= a[j + j * a_dim1];
 109                         }
 110                     }
 111                     jx -= incx;
 112                 }
 113             }
 114         }
 115     } else {
 116
 117         if (*uplo=='U' || *uplo=='u') {
 118             if (incx == 1) {
 119                 for (j = n; j >= 1; --j) {
 120                     temp = x[j];
 121                     if (nounit) {
 122                         temp *= a[j + j * a_dim1];
 123                     }
 124                     for (i__ = j - 1; i__ >= 1; --i__) {
 125                         temp += a[i__ + j * a_dim1] * x[i__];
 126                     }
 127                     x[j] = temp;
 128                 }
 129             } else {
 130                 jx = kx + (n - 1) * incx;
 131                 for (j = n; j >= 1; --j) {
 132                     temp = x[jx];
 133                     ix = jx;
 134                     if (nounit) {
 135                         temp *= a[j + j * a_dim1];
 136                     }
 137                     for (i__ = j - 1; i__ >= 1; --i__) {
 138                         ix -= incx;
 139                         temp += a[i__ + j * a_dim1] * x[ix];
 140                     }
 141                     x[jx] = temp;
 142                     jx -= incx;
 143                 }
 144             }
 145         } else {
 146             if (incx == 1) {
 147                 i__1 = n;
 148                 for (j = 1; j <= i__1; ++j) {
 149                     temp = x[j];
 150                     if (nounit) {
 151                         temp *= a[j + j * a_dim1];
 152                     }
 153                     i__2 = n;
 154                     for (i__ = j + 1; i__ <= i__2; ++i__) {
 155                         temp += a[i__ + j * a_dim1] * x[i__];
 156                     }
 157                     x[j] = temp;
 158                 }
 159             } else {
 160                 jx = kx;
 161                 i__1 = n;
 162                 for (j = 1; j <= i__1; ++j) {
 163                     temp = x[jx];
 164                     ix = jx;
 165                     if (nounit) {
 166                         temp *= a[j + j * a_dim1];
 167                     }
 168                     i__2 = n;
 169                     for (i__ = j + 1; i__ <= i__2; ++i__) {
 170                         ix += incx;
 171                         temp += a[i__ + j * a_dim1] * x[ix];
 172                     }
 173                     x[jx] = temp;
 174                     jx += incx;
 175                 }
 176             }
 177         }
 178     }
 179
 180     return;
 181
 182 }
 183
 184