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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static doublereal c_b8 = 0.;
   7 static doublereal c_b14 = -1.;
   8
   9 /* Subroutine */ int dsytd2_(char *uplo, integer *n, doublereal *a, integer *
  10         lda, doublereal *d__, doublereal *e, doublereal *tau, integer *info)
  11 {
  12     /* System generated locals */
  13     integer a_dim1, a_offset, i__1, i__2, i__3;
  14
  15     /* Local variables */
  16     integer i__;
  17     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  18             integer *);
  19     doublereal taui;
  20     extern /* Subroutine */ int dsyr2_(char *, integer *, doublereal *,
  21             doublereal *, integer *, doublereal *, integer *, doublereal *,
  22             integer *);
  23     doublereal alpha;
  24     extern logical lsame_(char *, char *);
  25     extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
  26             integer *, doublereal *, integer *);
  27     logical upper;
  28     extern /* Subroutine */ int dsymv_(char *, integer *, doublereal *,
  29             doublereal *, integer *, doublereal *, integer *, doublereal *,
  30             doublereal *, integer *), dlarfg_(integer *, doublereal *,
  31              doublereal *, integer *, doublereal *), xerbla_(char *, integer *
  32 );
  33
  34
  35 /*  -- LAPACK routine (version 3.1) -- */
  36 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  37 /*     November 2006 */
  38
  39 /*     .. Scalar Arguments .. */
  40 /*     .. */
  41 /*     .. Array Arguments .. */
  42 /*     .. */
  43
  44 /*  Purpose */
  45 /*  ======= */
  46
  47 /*  DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal */
  48 /*  form T by an orthogonal similarity transformation: Q' * A * Q = T. */
  49
  50 /*  Arguments */
  51 /*  ========= */
  52
  53 /*  UPLO    (input) CHARACTER*1 */
  54 /*          Specifies whether the upper or lower triangular part of the */
  55 /*          symmetric matrix A is stored: */
  56 /*          = 'U':  Upper triangular */
  57 /*          = 'L':  Lower triangular */
  58
  59 /*  N       (input) INTEGER */
  60 /*          The order of the matrix A.  N >= 0. */
  61
  62 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
  63 /*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading */
  64 /*          n-by-n upper triangular part of A contains the upper */
  65 /*          triangular part of the matrix A, and the strictly lower */
  66 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
  67 /*          leading n-by-n lower triangular part of A contains the lower */
  68 /*          triangular part of the matrix A, and the strictly upper */
  69 /*          triangular part of A is not referenced. */
  70 /*          On exit, if UPLO = 'U', the diagonal and first superdiagonal */
  71 /*          of A are overwritten by the corresponding elements of the */
  72 /*          tridiagonal matrix T, and the elements above the first */
  73 /*          superdiagonal, with the array TAU, represent the orthogonal */
  74 /*          matrix Q as a product of elementary reflectors; if UPLO */
  75 /*          = 'L', the diagonal and first subdiagonal of A are over- */
  76 /*          written by the corresponding elements of the tridiagonal */
  77 /*          matrix T, and the elements below the first subdiagonal, with */
  78 /*          the array TAU, represent the orthogonal matrix Q as a product */
  79 /*          of elementary reflectors. See Further Details. */
  80
  81 /*  LDA     (input) INTEGER */
  82 /*          The leading dimension of the array A.  LDA >= max(1,N). */
  83
  84 /*  D       (output) DOUBLE PRECISION array, dimension (N) */
  85 /*          The diagonal elements of the tridiagonal matrix T: */
  86 /*          D(i) = A(i,i). */
  87
  88 /*  E       (output) DOUBLE PRECISION array, dimension (N-1) */
  89 /*          The off-diagonal elements of the tridiagonal matrix T: */
  90 /*          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. */
  91
  92 /*  TAU     (output) DOUBLE PRECISION array, dimension (N-1) */
  93 /*          The scalar factors of the elementary reflectors (see Further */
  94 /*          Details). */
  95
  96 /*  INFO    (output) INTEGER */
  97 /*          = 0:  successful exit */
  98 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
  99
 100 /*  Further Details */
 101 /*  =============== */
 102
 103 /*  If UPLO = 'U', the matrix Q is represented as a product of elementary */
 104 /*  reflectors */
 105
 106 /*     Q = H(n-1) . . . H(2) H(1). */
 107
 108 /*  Each H(i) has the form */
 109
 110 /*     H(i) = I - tau * v * v' */
 111
 112 /*  where tau is a real scalar, and v is a real vector with */
 113 /*  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in */
 114 /*  A(1:i-1,i+1), and tau in TAU(i). */
 115
 116 /*  If UPLO = 'L', the matrix Q is represented as a product of elementary */
 117 /*  reflectors */
 118
 119 /*     Q = H(1) H(2) . . . H(n-1). */
 120
 121 /*  Each H(i) has the form */
 122
 123 /*     H(i) = I - tau * v * v' */
 124
 125 /*  where tau is a real scalar, and v is a real vector with */
 126 /*  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), */
 127 /*  and tau in TAU(i). */
 128
 129 /*  The contents of A on exit are illustrated by the following examples */
 130 /*  with n = 5: */
 131
 132 /*  if UPLO = 'U':                       if UPLO = 'L': */
 133
 134 /*    (  d   e   v2  v3  v4 )              (  d                  ) */
 135 /*    (      d   e   v3  v4 )              (  e   d              ) */
 136 /*    (          d   e   v4 )              (  v1  e   d          ) */
 137 /*    (              d   e  )              (  v1  v2  e   d      ) */
 138 /*    (                  d  )              (  v1  v2  v3  e   d  ) */
 139
 140 /*  where d and e denote diagonal and off-diagonal elements of T, and vi */
 141 /*  denotes an element of the vector defining H(i). */
 142
 143 /*  ===================================================================== */
 144
 145 /*     .. Parameters .. */
 146 /*     .. */
 147 /*     .. Local Scalars .. */
 148 /*     .. */
 149 /*     .. External Subroutines .. */
 150 /*     .. */
 151 /*     .. External Functions .. */
 152 /*     .. */
 153 /*     .. Intrinsic Functions .. */
 154 /*     .. */
 155 /*     .. Executable Statements .. */
 156
 157 /*     Test the input parameters */
 158
 159     /* Parameter adjustments */
 160     a_dim1 = *lda;
 161     a_offset = 1 + a_dim1;
 162     a -= a_offset;
 163     --d__;
 164     --e;
 165     --tau;
 166
 167     /* Function Body */
 168     *info = 0;
 169     upper = lsame_(uplo, "U");
 170     if (! upper && ! lsame_(uplo, "L")) {
 171         *info = -1;
 172     } else if (*n < 0) {
 173         *info = -2;
 174     } else if (*lda < max(1,*n)) {
 175         *info = -4;
 176     }
 177     if (*info != 0) {
 178         i__1 = -(*info);
 179         xerbla_("DSYTD2", &i__1);
 180         return 0;
 181     }
 182
 183 /*     Quick return if possible */
 184
 185     if (*n <= 0) {
 186         return 0;
 187     }
 188
 189     if (upper) {
 190
 191 /*        Reduce the upper triangle of A */
 192
 193         for (i__ = *n - 1; i__ >= 1; --i__) {
 194
 195 /*           Generate elementary reflector H(i) = I - tau * v * v' */
 196 /*           to annihilate A(1:i-1,i+1) */
 197
 198             dlarfg_(&i__, &a[i__ + (i__ + 1) * a_dim1], &a[(i__ + 1) * a_dim1
 199                     + 1], &c__1, &taui);
 200             e[i__] = a[i__ + (i__ + 1) * a_dim1];
 201
 202             if (taui != 0.) {
 203
 204 /*              Apply H(i) from both sides to A(1:i,1:i) */
 205
 206                 a[i__ + (i__ + 1) * a_dim1] = 1.;
 207
 208 /*              Compute  x := tau * A * v  storing x in TAU(1:i) */
 209
 210                 dsymv_(uplo, &i__, &taui, &a[a_offset], lda, &a[(i__ + 1) *
 211                         a_dim1 + 1], &c__1, &c_b8, &tau[1], &c__1);
 212
 213 /*              Compute  w := x - 1/2 * tau * (x'*v) * v */
 214
 215                 alpha = taui * -.5 * ddot_(&i__, &tau[1], &c__1, &a[(i__ + 1)
 216                         * a_dim1 + 1], &c__1);
 217                 daxpy_(&i__, &alpha, &a[(i__ + 1) * a_dim1 + 1], &c__1, &tau[
 218                         1], &c__1);
 219
 220 /*              Apply the transformation as a rank-2 update: */
 221 /*                 A := A - v * w' - w * v' */
 222
 223                 dsyr2_(uplo, &i__, &c_b14, &a[(i__ + 1) * a_dim1 + 1], &c__1,
 224                         &tau[1], &c__1, &a[a_offset], lda);
 225
 226                 a[i__ + (i__ + 1) * a_dim1] = e[i__];
 227             }
 228             d__[i__ + 1] = a[i__ + 1 + (i__ + 1) * a_dim1];
 229             tau[i__] = taui;
 230 /* L10: */
 231         }
 232         d__[1] = a[a_dim1 + 1];
 233     } else {
 234
 235 /*        Reduce the lower triangle of A */
 236
 237         i__1 = *n - 1;
 238         for (i__ = 1; i__ <= i__1; ++i__) {
 239
 240 /*           Generate elementary reflector H(i) = I - tau * v * v' */
 241 /*           to annihilate A(i+2:n,i) */
 242
 243             i__2 = *n - i__;
 244 /* Computing MIN */
 245             i__3 = i__ + 2;
 246             dlarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *n)+ i__ *
 247                      a_dim1], &c__1, &taui);
 248             e[i__] = a[i__ + 1 + i__ * a_dim1];
 249
 250             if (taui != 0.) {
 251
 252 /*              Apply H(i) from both sides to A(i+1:n,i+1:n) */
 253
 254                 a[i__ + 1 + i__ * a_dim1] = 1.;
 255
 256 /*              Compute  x := tau * A * v  storing y in TAU(i:n-1) */
 257
 258                 i__2 = *n - i__;
 259                 dsymv_(uplo, &i__2, &taui, &a[i__ + 1 + (i__ + 1) * a_dim1],
 260                         lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b8, &tau[
 261                         i__], &c__1);
 262
 263 /*              Compute  w := x - 1/2 * tau * (x'*v) * v */
 264
 265                 i__2 = *n - i__;
 266                 alpha = taui * -.5 * ddot_(&i__2, &tau[i__], &c__1, &a[i__ +
 267                         1 + i__ * a_dim1], &c__1);
 268                 i__2 = *n - i__;
 269                 daxpy_(&i__2, &alpha, &a[i__ + 1 + i__ * a_dim1], &c__1, &tau[
 270                         i__], &c__1);
 271
 272 /*              Apply the transformation as a rank-2 update: */
 273 /*                 A := A - v * w' - w * v' */
 274
 275                 i__2 = *n - i__;
 276                 dsyr2_(uplo, &i__2, &c_b14, &a[i__ + 1 + i__ * a_dim1], &c__1,
 277                          &tau[i__], &c__1, &a[i__ + 1 + (i__ + 1) * a_dim1],
 278                         lda);
 279
 280                 a[i__ + 1 + i__ * a_dim1] = e[i__];
 281             }
 282             d__[i__] = a[i__ + i__ * a_dim1];
 283             tau[i__] = taui;
 284 /* L20: */
 285         }
 286         d__[*n] = a[*n + *n * a_dim1];
 287     }
 288
 289     return 0;
 290
 291 /*     End of DSYTD2 */
 292
 293 } /* dsytd2_ */