src/gmxlib/gmx_lapack/ssytrd.c

   1 #include "gmx_blas.h"
   2 #include "gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 void
   6 F77_FUNC(ssytrd,SSYTRD)(const char *uplo, int *n, float *a, int *
   7         lda, float *d__, float *e, float *tau, float *
   8         work, int *lwork, int *info)
   9 {
  10     /* System generated locals */
  11     int a_dim1, a_offset, i__1, i__2, i__3;
  12
  13     /* Local variables */
  14     int i__, j, nb, kk, nx, iws;
  15     int nbmin, iinfo;
  16     int upper;
  17     int ldwork, lwkopt;
  18     int lquery;
  19     float c_b22 = -1.;
  20     float c_b23 = 1.;
  21
  22
  23     /* Parameter adjustments */
  24     a_dim1 = *lda;
  25     a_offset = 1 + a_dim1;
  26     a -= a_offset;
  27     --d__;
  28     --e;
  29     --tau;
  30     --work;
  31
  32     /* Function Body */
  33     *info = 0;
  34     upper = (*uplo=='U' || *uplo=='u');
  35     lquery = (*lwork == -1);
  36
  37     if (! upper && ! (*uplo=='L' || *uplo=='l')) {
  38         *info = -1;
  39     } else if (*n < 0) {
  40         *info = -2;
  41     } else if (*lda < ((1>*n) ? 1 : *n)) {
  42         *info = -4;
  43     } else if (*lwork < 1 && ! lquery) {
  44         *info = -9;
  45     }
  46
  47     if (*info == 0) {
  48
  49       nb = DSYTRD_BLOCKSIZE;
  50       lwkopt = *n * nb;
  51       work[1] = (float) lwkopt;
  52     } else
  53       return;
  54
  55     if (lquery)
  56       return;
  57
  58     if (*n == 0) {
  59         work[1] = 1.;
  60         return;
  61     }
  62
  63     nx = *n;
  64     iws = 1;
  65     if (nb > 1 && nb < *n) {
  66
  67         nx = DSYTRD_CROSSOVER;
  68         if (nx < *n) {
  69
  70             ldwork = *n;
  71             iws = ldwork * nb;
  72             if (*lwork < iws) {
  73
  74                 i__1 = *lwork / ldwork;
  75                 nb = (i__1>1) ? i__1 : 1;
  76                 nbmin = DSYTRD_MINBLOCKSIZE;
  77                 if (nb < nbmin) {
  78                     nx = *n;
  79                 }
  80             }
  81         } else {
  82             nx = *n;
  83         }
  84     } else {
  85         nb = 1;
  86     }
  87
  88     if (upper) {
  89
  90         kk = *n - (*n - nx + nb - 1) / nb * nb;
  91         i__1 = kk + 1;
  92         i__2 = -nb;
  93         for (i__ = *n - nb + 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
  94                 i__2) {
  95
  96             i__3 = i__ + nb - 1;
  97             F77_FUNC(slatrd,SLATRD)(uplo, &i__3, &nb, &a[a_offset], lda, &e[1], &tau[1], &
  98                     work[1], &ldwork);
  99
 100             i__3 = i__ - 1;
 101             F77_FUNC(ssyr2k,SSYR2K)(uplo, "No transpose", &i__3, &nb, &c_b22, &a[i__ * a_dim1
 102                     + 1], lda, &work[1], &ldwork, &c_b23, &a[a_offset], lda);
 103
 104             i__3 = i__ + nb - 1;
 105             for (j = i__; j <= i__3; ++j) {
 106                 a[j - 1 + j * a_dim1] = e[j - 1];
 107                 d__[j] = a[j + j * a_dim1];
 108
 109             }
 110
 111         }
 112
 113         F77_FUNC(ssytd2,SSYTD2)(uplo, &kk, &a[a_offset], lda, &d__[1], &e[1], &tau[1], &iinfo);
 114     } else {
 115
 116         i__2 = *n - nx;
 117         i__1 = nb;
 118         for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
 119
 120
 121             i__3 = *n - i__ + 1;
 122             F77_FUNC(slatrd,SLATRD)(uplo, &i__3, &nb, &a[i__ + i__ * a_dim1], lda, &e[i__], &
 123                     tau[i__], &work[1], &ldwork);
 124
 125             i__3 = *n - i__ - nb + 1;
 126             F77_FUNC(ssyr2k,SSYR2K)(uplo, "No transpose", &i__3, &nb, &c_b22, &a[i__ + nb +
 127                     i__ * a_dim1], lda, &work[nb + 1], &ldwork, &c_b23, &a[
 128                     i__ + nb + (i__ + nb) * a_dim1], lda);
 129
 130
 131             i__3 = i__ + nb - 1;
 132             for (j = i__; j <= i__3; ++j) {
 133                 a[j + 1 + j * a_dim1] = e[j];
 134                 d__[j] = a[j + j * a_dim1];
 135
 136             }
 137
 138         }
 139
 140
 141         i__1 = *n - i__ + 1;
 142         F77_FUNC(ssytd2,SSYTD2)(uplo, &i__1, &a[i__ + i__ * a_dim1], lda, &d__[i__], &e[i__],
 143                 &tau[i__], &iinfo);
 144     }
 145
 146     work[1] = (float) lwkopt;
 147     return;
 148
 149 }
 150
 151