src/gromacs/linearalgebra/gmx_lapack/dlarft.c

   1 #include <math.h>
   2 #include "types/simple.h"
   3
   4 #include "../gmx_blas.h"
   5 #include "../gmx_lapack.h"
   6
   7 void
   8 F77_FUNC(dlarft,DLARFT)(const char *direct,
   9         const char *storev,
  10         int *n,
  11         int *k,
  12         double *v,
  13         int *ldv,
  14         double *tau,
  15         double *t,
  16         int *ldt)
  17 {
  18     /* System generated locals */
  19     int t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3;
  20     double d__1;
  21
  22     /* Local variables */
  23     int i__, j;
  24     double vii;
  25     int c__1 = 1;
  26     double zero = 0.0;
  27
  28     v_dim1 = *ldv;
  29     v_offset = 1 + v_dim1;
  30     v -= v_offset;
  31     --tau;
  32     t_dim1 = *ldt;
  33     t_offset = 1 + t_dim1;
  34     t -= t_offset;
  35
  36     if (*n == 0) {
  37         return;
  38     }
  39
  40     if (*direct=='F' || *direct=='f') {
  41         i__1 = *k;
  42         for (i__ = 1; i__ <= i__1; ++i__) {
  43             if (fabs(tau[i__])<GMX_DOUBLE_MIN) {
  44
  45                 i__2 = i__;
  46                 for (j = 1; j <= i__2; ++j) {
  47                     t[j + i__ * t_dim1] = 0.;
  48                 }
  49             } else {
  50
  51                 vii = v[i__ + i__ * v_dim1];
  52                 v[i__ + i__ * v_dim1] = 1.;
  53                 if (*storev=='C' || *storev=='c') {
  54
  55                     i__2 = *n - i__ + 1;
  56                     i__3 = i__ - 1;
  57                     d__1 = -tau[i__];
  58                     F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &d__1, &v[i__ + v_dim1],
  59                              ldv, &v[i__ + i__ * v_dim1], &c__1, &zero, &t[
  60                             i__ * t_dim1 + 1], &c__1);
  61                 } else {
  62
  63                     i__2 = i__ - 1;
  64                     i__3 = *n - i__ + 1;
  65                     d__1 = -tau[i__];
  66                     F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &d__1, &v[i__ *
  67                             v_dim1 + 1], ldv, &v[i__ + i__ * v_dim1], ldv, &
  68                             zero, &t[i__ * t_dim1 + 1], &c__1);
  69                 }
  70                 v[i__ + i__ * v_dim1] = vii;
  71
  72
  73                 i__2 = i__ - 1;
  74                 F77_FUNC(dtrmv,DTRMV)("Upper", "No transpose", "Non-unit", &i__2, &t[
  75                         t_offset], ldt, &t[i__ * t_dim1 + 1], &c__1);
  76                 t[i__ + i__ * t_dim1] = tau[i__];
  77             }
  78         }
  79     } else {
  80         for (i__ = *k; i__ >= 1; --i__) {
  81             if (fabs(tau[i__])<GMX_DOUBLE_MIN) {
  82
  83                 i__1 = *k;
  84                 for (j = i__; j <= i__1; ++j) {
  85                     t[j + i__ * t_dim1] = 0.;
  86                 }
  87             } else {
  88
  89                 if (i__ < *k) {
  90                     if (*storev=='C' || *storev=='c') {
  91                         vii = v[*n - *k + i__ + i__ * v_dim1];
  92                         v[*n - *k + i__ + i__ * v_dim1] = 1.;
  93
  94                         i__1 = *n - *k + i__;
  95                         i__2 = *k - i__;
  96                         d__1 = -tau[i__];
  97                         F77_FUNC(dgemv,DGEMV)("Transpose", &i__1, &i__2, &d__1, &v[(i__ + 1)
  98                                 * v_dim1 + 1], ldv, &v[i__ * v_dim1 + 1], &
  99                                 c__1, &zero, &t[i__ + 1 + i__ * t_dim1], &
 100                                 c__1);
 101                         v[*n - *k + i__ + i__ * v_dim1] = vii;
 102                     } else {
 103                         vii = v[i__ + (*n - *k + i__) * v_dim1];
 104                         v[i__ + (*n - *k + i__) * v_dim1] = 1.;
 105
 106                         i__1 = *k - i__;
 107                         i__2 = *n - *k + i__;
 108                         d__1 = -tau[i__];
 109                         F77_FUNC(dgemv,DGEMV)("No transpose", &i__1, &i__2, &d__1, &v[i__ +
 110                                 1 + v_dim1], ldv, &v[i__ + v_dim1], ldv, &
 111                                 zero, &t[i__ + 1 + i__ * t_dim1], &c__1);
 112                         v[i__ + (*n - *k + i__) * v_dim1] = vii;
 113                     }
 114
 115                     i__1 = *k - i__;
 116                     F77_FUNC(dtrmv,DTRMV)("Lower", "No transpose", "Non-unit", &i__1, &t[i__
 117                             + 1 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ *
 118                              t_dim1], &c__1)
 119                             ;
 120                 }
 121                 t[i__ + i__ * t_dim1] = tau[i__];
 122             }
 123         }
 124     }
 125     return;
 126
 127
 128 }