git.maemo.org Git - opencv/blob - 3rdparty/lapack/ssytrd.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c_n1 = -1;
   7 static integer c__3 = 3;
   8 static integer c__2 = 2;
   9 static real c_b22 = -1.f;
  10 static real c_b23 = 1.f;
  11
  12 /* Subroutine */ int ssytrd_(char *uplo, integer *n, real *a, integer *lda,
  13         real *d__, real *e, real *tau, real *work, integer *lwork, integer *
  14         info)
  15 {
  16     /* System generated locals */
  17     integer a_dim1, a_offset, i__1, i__2, i__3;
  18
  19     /* Local variables */
  20     integer i__, j, nb, kk, nx, iws;
  21     extern logical lsame_(char *, char *);
  22     integer nbmin, iinfo;
  23     logical upper;
  24     extern /* Subroutine */ int ssytd2_(char *, integer *, real *, integer *,
  25             real *, real *, real *, integer *), ssyr2k_(char *, char *
  26 , integer *, integer *, real *, real *, integer *, real *,
  27             integer *, real *, real *, integer *), xerbla_(
  28             char *, integer *);
  29     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  30             integer *, integer *);
  31     extern /* Subroutine */ int slatrd_(char *, integer *, integer *, real *,
  32             integer *, real *, real *, real *, integer *);
  33     integer ldwork, lwkopt;
  34     logical lquery;
  35
  36
  37 /*  -- LAPACK routine (version 3.1) -- */
  38 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  39 /*     November 2006 */
  40
  41 /*     .. Scalar Arguments .. */
  42 /*     .. */
  43 /*     .. Array Arguments .. */
  44 /*     .. */
  45
  46 /*  Purpose */
  47 /*  ======= */
  48
  49 /*  SSYTRD reduces a real symmetric matrix A to real symmetric */
  50 /*  tridiagonal form T by an orthogonal similarity transformation: */
  51 /*  Q**T * A * Q = T. */
  52
  53 /*  Arguments */
  54 /*  ========= */
  55
  56 /*  UPLO    (input) CHARACTER*1 */
  57 /*          = 'U':  Upper triangle of A is stored; */
  58 /*          = 'L':  Lower triangle of A is stored. */
  59
  60 /*  N       (input) INTEGER */
  61 /*          The order of the matrix A.  N >= 0. */
  62
  63 /*  A       (input/output) REAL array, dimension (LDA,N) */
  64 /*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading */
  65 /*          N-by-N upper triangular part of A contains the upper */
  66 /*          triangular part of the matrix A, and the strictly lower */
  67 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
  68 /*          leading N-by-N lower triangular part of A contains the lower */
  69 /*          triangular part of the matrix A, and the strictly upper */
  70 /*          triangular part of A is not referenced. */
  71 /*          On exit, if UPLO = 'U', the diagonal and first superdiagonal */
  72 /*          of A are overwritten by the corresponding elements of the */
  73 /*          tridiagonal matrix T, and the elements above the first */
  74 /*          superdiagonal, with the array TAU, represent the orthogonal */
  75 /*          matrix Q as a product of elementary reflectors; if UPLO */
  76 /*          = 'L', the diagonal and first subdiagonal of A are over- */
  77 /*          written by the corresponding elements of the tridiagonal */
  78 /*          matrix T, and the elements below the first subdiagonal, with */
  79 /*          the array TAU, represent the orthogonal matrix Q as a product */
  80 /*          of elementary reflectors. See Further Details. */
  81
  82 /*  LDA     (input) INTEGER */
  83 /*          The leading dimension of the array A.  LDA >= max(1,N). */
  84
  85 /*  D       (output) REAL array, dimension (N) */
  86 /*          The diagonal elements of the tridiagonal matrix T: */
  87 /*          D(i) = A(i,i). */
  88
  89 /*  E       (output) REAL array, dimension (N-1) */
  90 /*          The off-diagonal elements of the tridiagonal matrix T: */
  91 /*          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. */
  92
  93 /*  TAU     (output) REAL array, dimension (N-1) */
  94 /*          The scalar factors of the elementary reflectors (see Further */
  95 /*          Details). */
  96
  97 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
  98 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  99
 100 /*  LWORK   (input) INTEGER */
 101 /*          The dimension of the array WORK.  LWORK >= 1. */
 102 /*          For optimum performance LWORK >= N*NB, where NB is the */
 103 /*          optimal blocksize. */
 104
 105 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 106 /*          only calculates the optimal size of the WORK array, returns */
 107 /*          this value as the first entry of the WORK array, and no error */
 108 /*          message related to LWORK is issued by XERBLA. */
 109
 110 /*  INFO    (output) INTEGER */
 111 /*          = 0:  successful exit */
 112 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 113
 114 /*  Further Details */
 115 /*  =============== */
 116
 117 /*  If UPLO = 'U', the matrix Q is represented as a product of elementary */
 118 /*  reflectors */
 119
 120 /*     Q = H(n-1) . . . H(2) H(1). */
 121
 122 /*  Each H(i) has the form */
 123
 124 /*     H(i) = I - tau * v * v' */
 125
 126 /*  where tau is a real scalar, and v is a real vector with */
 127 /*  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in */
 128 /*  A(1:i-1,i+1), and tau in TAU(i). */
 129
 130 /*  If UPLO = 'L', the matrix Q is represented as a product of elementary */
 131 /*  reflectors */
 132
 133 /*     Q = H(1) H(2) . . . H(n-1). */
 134
 135 /*  Each H(i) has the form */
 136
 137 /*     H(i) = I - tau * v * v' */
 138
 139 /*  where tau is a real scalar, and v is a real vector with */
 140 /*  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), */
 141 /*  and tau in TAU(i). */
 142
 143 /*  The contents of A on exit are illustrated by the following examples */
 144 /*  with n = 5: */
 145
 146 /*  if UPLO = 'U':                       if UPLO = 'L': */
 147
 148 /*    (  d   e   v2  v3  v4 )              (  d                  ) */
 149 /*    (      d   e   v3  v4 )              (  e   d              ) */
 150 /*    (          d   e   v4 )              (  v1  e   d          ) */
 151 /*    (              d   e  )              (  v1  v2  e   d      ) */
 152 /*    (                  d  )              (  v1  v2  v3  e   d  ) */
 153
 154 /*  where d and e denote diagonal and off-diagonal elements of T, and vi */
 155 /*  denotes an element of the vector defining H(i). */
 156
 157 /*  ===================================================================== */
 158
 159 /*     .. Parameters .. */
 160 /*     .. */
 161 /*     .. Local Scalars .. */
 162 /*     .. */
 163 /*     .. External Subroutines .. */
 164 /*     .. */
 165 /*     .. Intrinsic Functions .. */
 166 /*     .. */
 167 /*     .. External Functions .. */
 168 /*     .. */
 169 /*     .. Executable Statements .. */
 170
 171 /*     Test the input parameters */
 172
 173     /* Parameter adjustments */
 174     a_dim1 = *lda;
 175     a_offset = 1 + a_dim1;
 176     a -= a_offset;
 177     --d__;
 178     --e;
 179     --tau;
 180     --work;
 181
 182     /* Function Body */
 183     *info = 0;
 184     upper = lsame_(uplo, "U");
 185     lquery = *lwork == -1;
 186     if (! upper && ! lsame_(uplo, "L")) {
 187         *info = -1;
 188     } else if (*n < 0) {
 189         *info = -2;
 190     } else if (*lda < max(1,*n)) {
 191         *info = -4;
 192     } else if (*lwork < 1 && ! lquery) {
 193         *info = -9;
 194     }
 195
 196     if (*info == 0) {
 197
 198 /*        Determine the block size. */
 199
 200         nb = ilaenv_(&c__1, "SSYTRD", uplo, n, &c_n1, &c_n1, &c_n1);
 201         lwkopt = *n * nb;
 202         work[1] = (real) lwkopt;
 203     }
 204
 205     if (*info != 0) {
 206         i__1 = -(*info);
 207         xerbla_("SSYTRD", &i__1);
 208         return 0;
 209     } else if (lquery) {
 210         return 0;
 211     }
 212
 213 /*     Quick return if possible */
 214
 215     if (*n == 0) {
 216         work[1] = 1.f;
 217         return 0;
 218     }
 219
 220     nx = *n;
 221     iws = 1;
 222     if (nb > 1 && nb < *n) {
 223
 224 /*        Determine when to cross over from blocked to unblocked code */
 225 /*        (last block is always handled by unblocked code). */
 226
 227 /* Computing MAX */
 228         i__1 = nb, i__2 = ilaenv_(&c__3, "SSYTRD", uplo, n, &c_n1, &c_n1, &
 229                 c_n1);
 230         nx = max(i__1,i__2);
 231         if (nx < *n) {
 232
 233 /*           Determine if workspace is large enough for blocked code. */
 234
 235             ldwork = *n;
 236             iws = ldwork * nb;
 237             if (*lwork < iws) {
 238
 239 /*              Not enough workspace to use optimal NB:  determine the */
 240 /*              minimum value of NB, and reduce NB or force use of */
 241 /*              unblocked code by setting NX = N. */
 242
 243 /* Computing MAX */
 244                 i__1 = *lwork / ldwork;
 245                 nb = max(i__1,1);
 246                 nbmin = ilaenv_(&c__2, "SSYTRD", uplo, n, &c_n1, &c_n1, &c_n1);
 247                 if (nb < nbmin) {
 248                     nx = *n;
 249                 }
 250             }
 251         } else {
 252             nx = *n;
 253         }
 254     } else {
 255         nb = 1;
 256     }
 257
 258     if (upper) {
 259
 260 /*        Reduce the upper triangle of A. */
 261 /*        Columns 1:kk are handled by the unblocked method. */
 262
 263         kk = *n - (*n - nx + nb - 1) / nb * nb;
 264         i__1 = kk + 1;
 265         i__2 = -nb;
 266         for (i__ = *n - nb + 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
 267                 i__2) {
 268
 269 /*           Reduce columns i:i+nb-1 to tridiagonal form and form the */
 270 /*           matrix W which is needed to update the unreduced part of */
 271 /*           the matrix */
 272
 273             i__3 = i__ + nb - 1;
 274             slatrd_(uplo, &i__3, &nb, &a[a_offset], lda, &e[1], &tau[1], &
 275                     work[1], &ldwork);
 276
 277 /*           Update the unreduced submatrix A(1:i-1,1:i-1), using an */
 278 /*           update of the form:  A := A - V*W' - W*V' */
 279
 280             i__3 = i__ - 1;
 281             ssyr2k_(uplo, "No transpose", &i__3, &nb, &c_b22, &a[i__ * a_dim1
 282                     + 1], lda, &work[1], &ldwork, &c_b23, &a[a_offset], lda);
 283
 284 /*           Copy superdiagonal elements back into A, and diagonal */
 285 /*           elements into D */
 286
 287             i__3 = i__ + nb - 1;
 288             for (j = i__; j <= i__3; ++j) {
 289                 a[j - 1 + j * a_dim1] = e[j - 1];
 290                 d__[j] = a[j + j * a_dim1];
 291 /* L10: */
 292             }
 293 /* L20: */
 294         }
 295
 296 /*        Use unblocked code to reduce the last or only block */
 297
 298         ssytd2_(uplo, &kk, &a[a_offset], lda, &d__[1], &e[1], &tau[1], &iinfo);
 299     } else {
 300
 301 /*        Reduce the lower triangle of A */
 302
 303         i__2 = *n - nx;
 304         i__1 = nb;
 305         for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
 306
 307 /*           Reduce columns i:i+nb-1 to tridiagonal form and form the */
 308 /*           matrix W which is needed to update the unreduced part of */
 309 /*           the matrix */
 310
 311             i__3 = *n - i__ + 1;
 312             slatrd_(uplo, &i__3, &nb, &a[i__ + i__ * a_dim1], lda, &e[i__], &
 313                     tau[i__], &work[1], &ldwork);
 314
 315 /*           Update the unreduced submatrix A(i+ib:n,i+ib:n), using */
 316 /*           an update of the form:  A := A - V*W' - W*V' */
 317
 318             i__3 = *n - i__ - nb + 1;
 319             ssyr2k_(uplo, "No transpose", &i__3, &nb, &c_b22, &a[i__ + nb +
 320                     i__ * a_dim1], lda, &work[nb + 1], &ldwork, &c_b23, &a[
 321                     i__ + nb + (i__ + nb) * a_dim1], lda);
 322
 323 /*           Copy subdiagonal elements back into A, and diagonal */
 324 /*           elements into D */
 325
 326             i__3 = i__ + nb - 1;
 327             for (j = i__; j <= i__3; ++j) {
 328                 a[j + 1 + j * a_dim1] = e[j];
 329                 d__[j] = a[j + j * a_dim1];
 330 /* L30: */
 331             }
 332 /* L40: */
 333         }
 334
 335 /*        Use unblocked code to reduce the last or only block */
 336
 337         i__1 = *n - i__ + 1;
 338         ssytd2_(uplo, &i__1, &a[i__ + i__ * a_dim1], lda, &d__[i__], &e[i__],
 339                 &tau[i__], &iinfo);
 340     }
 341
 342     work[1] = (real) lwkopt;
 343     return 0;
 344
 345 /*     End of SSYTRD */
 346
 347 } /* ssytrd_ */