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

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