src/gmxlib/gmx_lapack/dlabrd.c

   1 #include <math.h>
   2 #include "gmx_blas.h"
   3 #include "gmx_lapack.h"
   4
   5
   6 void
   7 F77_FUNC(dlabrd,DLABRD)(int *m,
   8         int *n,
   9         int *nb,
  10         double *a,
  11         int *lda,
  12         double *d__,
  13         double *e,
  14         double *tauq,
  15         double *taup,
  16         double *x,
  17         int *ldx,
  18         double *y,
  19         int *ldy)
  20 {
  21     int a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset;
  22     int i__1, i__2, i__3;
  23     double one = 1.0;
  24     double minusone = -1.0;
  25     double zero = 0.0;
  26     int c__1 = 1;
  27     int i__;
  28
  29     a_dim1 = *lda;
  30     a_offset = 1 + a_dim1;
  31     a -= a_offset;
  32     --d__;
  33     --e;
  34     --tauq;
  35     --taup;
  36     x_dim1 = *ldx;
  37     x_offset = 1 + x_dim1;
  38     x -= x_offset;
  39     y_dim1 = *ldy;
  40     y_offset = 1 + y_dim1;
  41     y -= y_offset;
  42
  43     if (*m <= 0 || *n <= 0) {
  44         return;
  45     }
  46
  47     if (*m >= *n) {
  48
  49         i__1 = *nb;
  50         for (i__ = 1; i__ <= i__1; ++i__) {
  51
  52             i__2 = *m - i__ + 1;
  53             i__3 = i__ - 1;
  54             F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &a[i__ + a_dim1], lda,
  55                      &y[i__ + y_dim1], ldy, &one, &a[i__ + i__ * a_dim1], &c__1);
  56             i__2 = *m - i__ + 1;
  57             i__3 = i__ - 1;
  58             F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &x[i__ + x_dim1], ldx,
  59                    &a[i__*a_dim1+1],&c__1,&one,&a[i__+i__*a_dim1],&c__1);
  60
  61             i__2 = *m - i__ + 1;
  62             i__3 = i__ + 1;
  63             if(*m<i__3)
  64               i__3 = *m;
  65             F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + i__ * a_dim1], &a[i__3 + i__ * a_dim1],
  66                     &c__1, &tauq[i__]);
  67             d__[i__] = a[i__ + i__ * a_dim1];
  68             if (i__ < *n) {
  69                 a[i__ + i__ * a_dim1] = 1.;
  70
  71                 i__2 = *m - i__ + 1;
  72                 i__3 = *n - i__;
  73                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + (i__ + 1) *
  74                         a_dim1], lda, &a[i__ + i__ * a_dim1], &c__1, &zero, &
  75                         y[i__ + 1 + i__ * y_dim1], &c__1);
  76                 i__2 = *m - i__ + 1;
  77                 i__3 = i__ - 1;
  78                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + a_dim1],
  79                         lda, &a[i__ + i__ * a_dim1], &c__1, &zero, &y[i__ *
  80                         y_dim1 + 1], &c__1);
  81                 i__2 = *n - i__;
  82                 i__3 = i__ - 1;
  83                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &y[i__ + 1 +
  84                         y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &one, &y[
  85                         i__ + 1 + i__ * y_dim1], &c__1);
  86                 i__2 = *m - i__ + 1;
  87                 i__3 = i__ - 1;
  88                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &x[i__ + x_dim1],
  89                         ldx, &a[i__ + i__ * a_dim1], &c__1, &zero, &y[i__ *
  90                         y_dim1 + 1], &c__1);
  91                 i__2 = i__ - 1;
  92                 i__3 = *n - i__;
  93                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &minusone, &a[(i__ + 1) *
  94                         a_dim1 + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &one,
  95                         &y[i__ + 1 + i__ * y_dim1], &c__1);
  96                 i__2 = *n - i__;
  97                 F77_FUNC(dscal,DSCAL)(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
  98
  99                 i__2 = *n - i__;
 100                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__, &minusone, &y[i__ + 1 +
 101                         y_dim1], ldy, &a[i__ + a_dim1], lda, &one, &a[i__ + (
 102                         i__ + 1) * a_dim1], lda);
 103                 i__2 = i__ - 1;
 104                 i__3 = *n - i__;
 105                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &minusone, &a[(i__ + 1) *
 106                         a_dim1 + 1], lda, &x[i__ + x_dim1], ldx, &one, &a[
 107                         i__ + (i__ + 1) * a_dim1], lda);
 108
 109                 i__2 = *n - i__;
 110                 i__3 = i__ + 2;
 111                 if(*n<i__3)
 112                   i__3 = *n;
 113                 F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + (i__ + 1) * a_dim1],
 114                         &a[i__ + i__3 * a_dim1], lda, &taup[i__]);
 115                 e[i__] = a[i__ + (i__ + 1) * a_dim1];
 116                 a[i__ + (i__ + 1) * a_dim1] = 1.;
 117
 118                 i__2 = *m - i__;
 119                 i__3 = *n - i__;
 120                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &one, &a[i__ + 1 + (i__
 121                         + 1) * a_dim1], lda, &a[i__ + (i__ + 1) * a_dim1],
 122                         lda, &zero, &x[i__ + 1 + i__ * x_dim1], &c__1);
 123                 i__2 = *n - i__;
 124                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__, &one, &y[i__ + 1 + y_dim1],
 125                         ldy, &a[i__ + (i__ + 1) * a_dim1], lda, &zero, &x[
 126                         i__ * x_dim1 + 1], &c__1);
 127                 i__2 = *m - i__;
 128                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__, &minusone, &a[i__ + 1 +
 129                         a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 130                         i__ + 1 + i__ * x_dim1], &c__1);
 131                 i__2 = i__ - 1;
 132                 i__3 = *n - i__;
 133                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &one, &a[(i__ + 1) *
 134                         a_dim1 + 1], lda, &a[i__ + (i__ + 1) * a_dim1], lda, &
 135                         zero, &x[i__ * x_dim1 + 1], &c__1);
 136                 i__2 = *m - i__;
 137                 i__3 = i__ - 1;
 138                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &x[i__ + 1 +
 139                         x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 140                         i__ + 1 + i__ * x_dim1], &c__1);
 141                 i__2 = *m - i__;
 142                 F77_FUNC(dscal,DSCAL)(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
 143             }
 144         }
 145     } else {
 146
 147         i__1 = *nb;
 148         for (i__ = 1; i__ <= i__1; ++i__) {
 149
 150             i__2 = *n - i__ + 1;
 151             i__3 = i__ - 1;
 152             F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &y[i__ + y_dim1], ldy,
 153                      &a[i__ + a_dim1], lda, &one, &a[i__ + i__ * a_dim1],lda);
 154             i__2 = i__ - 1;
 155             i__3 = *n - i__ + 1;
 156             F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &minusone, &a[i__ * a_dim1 + 1],
 157                     lda, &x[i__ + x_dim1], ldx, &one,&a[i__+i__*a_dim1],lda);
 158
 159             i__2 = *n - i__ + 1;
 160             i__3 = i__ + 1;
 161             if(*n<i__3)
 162               i__3 = *n;
 163             F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + i__ * a_dim1],
 164                     &a[i__ + i__3 * a_dim1], lda, &taup[i__]);
 165             d__[i__] = a[i__ + i__ * a_dim1];
 166             if (i__ < *m) {
 167                 a[i__ + i__ * a_dim1] = 1.;
 168
 169                 i__2 = *m - i__;
 170                 i__3 = *n - i__ + 1;
 171                 F77_FUNC(dgemv,DGEMV)("No transpose",&i__2,&i__3,&one,&a[i__+1+i__*a_dim1],
 172                        lda, &a[i__ + i__ * a_dim1], lda, &zero,
 173                        &x[i__ + 1 + i__ * x_dim1], &c__1);
 174                 i__2 = *n - i__ + 1;
 175                 i__3 = i__ - 1;
 176                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &y[i__ + y_dim1],
 177                         ldy, &a[i__ + i__ * a_dim1], lda, &zero, &x[i__ *
 178                         x_dim1 + 1], &c__1);
 179                 i__2 = *m - i__;
 180                 i__3 = i__ - 1;
 181                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &a[i__ + 1 +
 182                         a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 183                         i__ + 1 + i__ * x_dim1], &c__1);
 184                 i__2 = i__ - 1;
 185                 i__3 = *n - i__ + 1;
 186                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &one, &a[i__ * a_dim1 +
 187                         1], lda, &a[i__ + i__ * a_dim1], lda, &zero, &x[i__ *
 188                          x_dim1 + 1], &c__1);
 189                 i__2 = *m - i__;
 190                 i__3 = i__ - 1;
 191                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &x[i__ + 1 +
 192                         x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &one, &x[
 193                         i__ + 1 + i__ * x_dim1], &c__1);
 194                 i__2 = *m - i__;
 195                 F77_FUNC(dscal,DSCAL)(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
 196
 197                 i__2 = *m - i__;
 198                 i__3 = i__ - 1;
 199                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &a[i__ + 1 +
 200                         a_dim1], lda, &y[i__ + y_dim1], ldy, &one, &a[i__ +
 201                         1 + i__ * a_dim1], &c__1);
 202                 i__2 = *m - i__;
 203                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__, &minusone, &x[i__ + 1 +
 204                         x_dim1], ldx, &a[i__ * a_dim1 + 1], &c__1, &one, &a[
 205                         i__ + 1 + i__ * a_dim1], &c__1);
 206
 207                 i__2 = *m - i__;
 208                 i__3 = i__ + 2;
 209                 if(*m<i__3)
 210                   i__3 = *m;
 211                 F77_FUNC(dlarfg,DLARFG)(&i__2, &a[i__ + 1 + i__ * a_dim1],
 212                         &a[i__3 + i__ * a_dim1], &c__1, &tauq[i__]);
 213                 e[i__] = a[i__ + 1 + i__ * a_dim1];
 214                 a[i__ + 1 + i__ * a_dim1] = 1.;
 215
 216                 i__2 = *m - i__;
 217                 i__3 = *n - i__;
 218                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + 1 + (i__ +
 219                         1) * a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1,
 220                         &zero, &y[i__ + 1 + i__ * y_dim1], &c__1);
 221                 i__2 = *m - i__;
 222                 i__3 = i__ - 1;
 223                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__3, &one, &a[i__ + 1 + a_dim1],
 224                          lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &zero, &y[
 225                         i__ * y_dim1 + 1], &c__1);
 226                 i__2 = *n - i__;
 227                 i__3 = i__ - 1;
 228                 F77_FUNC(dgemv,DGEMV)("No transpose", &i__2, &i__3, &minusone, &y[i__ + 1 +
 229                         y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &one, &y[
 230                         i__ + 1 + i__ * y_dim1], &c__1);
 231                 i__2 = *m - i__;
 232                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__2, &i__, &one, &x[i__ + 1 + x_dim1],
 233                         ldx, &a[i__ + 1 + i__ * a_dim1], &c__1, &zero, &y[
 234                         i__ * y_dim1 + 1], &c__1);
 235                 i__2 = *n - i__;
 236                 F77_FUNC(dgemv,DGEMV)("Transpose", &i__, &i__2, &minusone, &a[(i__ + 1) * a_dim1
 237                         + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &one, &y[i__
 238                         + 1 + i__ * y_dim1], &c__1);
 239                 i__2 = *n - i__;
 240                 F77_FUNC(dscal,DSCAL)(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
 241             }
 242         }
 243     }
 244     return;
 245 }
 246