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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static doublereal c_b5 = -1.;
   6 static doublereal c_b6 = 1.;
   7 static integer c__1 = 1;
   8 static doublereal c_b16 = 0.;
   9
  10 /* Subroutine */ int dlatrd_(char *uplo, integer *n, integer *nb, doublereal *
  11         a, integer *lda, doublereal *e, doublereal *tau, doublereal *w,
  12         integer *ldw)
  13 {
  14     /* System generated locals */
  15     integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3;
  16
  17     /* Local variables */
  18     integer i__, iw;
  19     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  20             integer *);
  21     doublereal alpha;
  22     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  23             integer *);
  24     extern logical lsame_(char *, char *);
  25     extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
  26             doublereal *, doublereal *, integer *, doublereal *, integer *,
  27             doublereal *, doublereal *, integer *), daxpy_(integer *,
  28             doublereal *, doublereal *, integer *, doublereal *, integer *),
  29             dsymv_(char *, integer *, doublereal *, doublereal *, integer *,
  30             doublereal *, integer *, doublereal *, doublereal *, integer *), dlarfg_(integer *, doublereal *, doublereal *, integer *,
  31              doublereal *);
  32
  33
  34 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  35 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  36 /*     November 2006 */
  37
  38 /*     .. Scalar Arguments .. */
  39 /*     .. */
  40 /*     .. Array Arguments .. */
  41 /*     .. */
  42
  43 /*  Purpose */
  44 /*  ======= */
  45
  46 /*  DLATRD reduces NB rows and columns of a real symmetric matrix A to */
  47 /*  symmetric tridiagonal form by an orthogonal similarity */
  48 /*  transformation Q' * A * Q, and returns the matrices V and W which are */
  49 /*  needed to apply the transformation to the unreduced part of A. */
  50
  51 /*  If UPLO = 'U', DLATRD reduces the last NB rows and columns of a */
  52 /*  matrix, of which the upper triangle is supplied; */
  53 /*  if UPLO = 'L', DLATRD reduces the first NB rows and columns of a */
  54 /*  matrix, of which the lower triangle is supplied. */
  55
  56 /*  This is an auxiliary routine called by DSYTRD. */
  57
  58 /*  Arguments */
  59 /*  ========= */
  60
  61 /*  UPLO    (input) CHARACTER*1 */
  62 /*          Specifies whether the upper or lower triangular part of the */
  63 /*          symmetric matrix A is stored: */
  64 /*          = 'U': Upper triangular */
  65 /*          = 'L': Lower triangular */
  66
  67 /*  N       (input) INTEGER */
  68 /*          The order of the matrix A. */
  69
  70 /*  NB      (input) INTEGER */
  71 /*          The number of rows and columns to be reduced. */
  72
  73 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
  74 /*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading */
  75 /*          n-by-n upper triangular part of A contains the upper */
  76 /*          triangular part of the matrix A, and the strictly lower */
  77 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
  78 /*          leading n-by-n lower triangular part of A contains the lower */
  79 /*          triangular part of the matrix A, and the strictly upper */
  80 /*          triangular part of A is not referenced. */
  81 /*          On exit: */
  82 /*          if UPLO = 'U', the last NB columns have been reduced to */
  83 /*            tridiagonal form, with the diagonal elements overwriting */
  84 /*            the diagonal elements of A; the elements above the diagonal */
  85 /*            with the array TAU, represent the orthogonal matrix Q as a */
  86 /*            product of elementary reflectors; */
  87 /*          if UPLO = 'L', the first NB columns have been reduced to */
  88 /*            tridiagonal form, with the diagonal elements overwriting */
  89 /*            the diagonal elements of A; the elements below the diagonal */
  90 /*            with the array TAU, represent the  orthogonal matrix Q as a */
  91 /*            product of elementary reflectors. */
  92 /*          See Further Details. */
  93
  94 /*  LDA     (input) INTEGER */
  95 /*          The leading dimension of the array A.  LDA >= (1,N). */
  96
  97 /*  E       (output) DOUBLE PRECISION array, dimension (N-1) */
  98 /*          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal */
  99 /*          elements of the last NB columns of the reduced matrix; */
 100 /*          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of */
 101 /*          the first NB columns of the reduced matrix. */
 102
 103 /*  TAU     (output) DOUBLE PRECISION array, dimension (N-1) */
 104 /*          The scalar factors of the elementary reflectors, stored in */
 105 /*          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. */
 106 /*          See Further Details. */
 107
 108 /*  W       (output) DOUBLE PRECISION array, dimension (LDW,NB) */
 109 /*          The n-by-nb matrix W required to update the unreduced part */
 110 /*          of A. */
 111
 112 /*  LDW     (input) INTEGER */
 113 /*          The leading dimension of the array W. LDW >= max(1,N). */
 114
 115 /*  Further Details */
 116 /*  =============== */
 117
 118 /*  If UPLO = 'U', the matrix Q is represented as a product of elementary */
 119 /*  reflectors */
 120
 121 /*     Q = H(n) H(n-1) . . . H(n-nb+1). */
 122
 123 /*  Each H(i) has the form */
 124
 125 /*     H(i) = I - tau * v * v' */
 126
 127 /*  where tau is a real scalar, and v is a real vector with */
 128 /*  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), */
 129 /*  and tau in TAU(i-1). */
 130
 131 /*  If UPLO = 'L', the matrix Q is represented as a product of elementary */
 132 /*  reflectors */
 133
 134 /*     Q = H(1) H(2) . . . H(nb). */
 135
 136 /*  Each H(i) has the form */
 137
 138 /*     H(i) = I - tau * v * v' */
 139
 140 /*  where tau is a real scalar, and v is a real vector with */
 141 /*  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), */
 142 /*  and tau in TAU(i). */
 143
 144 /*  The elements of the vectors v together form the n-by-nb matrix V */
 145 /*  which is needed, with W, to apply the transformation to the unreduced */
 146 /*  part of the matrix, using a symmetric rank-2k update of the form: */
 147 /*  A := A - V*W' - W*V'. */
 148
 149 /*  The contents of A on exit are illustrated by the following examples */
 150 /*  with n = 5 and nb = 2: */
 151
 152 /*  if UPLO = 'U':                       if UPLO = 'L': */
 153
 154 /*    (  a   a   a   v4  v5 )              (  d                  ) */
 155 /*    (      a   a   v4  v5 )              (  1   d              ) */
 156 /*    (          a   1   v5 )              (  v1  1   a          ) */
 157 /*    (              d   1  )              (  v1  v2  a   a      ) */
 158 /*    (                  d  )              (  v1  v2  a   a   a  ) */
 159
 160 /*  where d denotes a diagonal element of the reduced matrix, a denotes */
 161 /*  an element of the original matrix that is unchanged, and vi denotes */
 162 /*  an element of the vector defining H(i). */
 163
 164 /*  ===================================================================== */
 165
 166 /*     .. Parameters .. */
 167 /*     .. */
 168 /*     .. Local Scalars .. */
 169 /*     .. */
 170 /*     .. External Subroutines .. */
 171 /*     .. */
 172 /*     .. External Functions .. */
 173 /*     .. */
 174 /*     .. Intrinsic Functions .. */
 175 /*     .. */
 176 /*     .. Executable Statements .. */
 177
 178 /*     Quick return if possible */
 179
 180     /* Parameter adjustments */
 181     a_dim1 = *lda;
 182     a_offset = 1 + a_dim1;
 183     a -= a_offset;
 184     --e;
 185     --tau;
 186     w_dim1 = *ldw;
 187     w_offset = 1 + w_dim1;
 188     w -= w_offset;
 189
 190     /* Function Body */
 191     if (*n <= 0) {
 192         return 0;
 193     }
 194
 195     if (lsame_(uplo, "U")) {
 196
 197 /*        Reduce last NB columns of upper triangle */
 198
 199         i__1 = *n - *nb + 1;
 200         for (i__ = *n; i__ >= i__1; --i__) {
 201             iw = i__ - *n + *nb;
 202             if (i__ < *n) {
 203
 204 /*              Update A(1:i,i) */
 205
 206                 i__2 = *n - i__;
 207                 dgemv_("No transpose", &i__, &i__2, &c_b5, &a[(i__ + 1) *
 208                         a_dim1 + 1], lda, &w[i__ + (iw + 1) * w_dim1], ldw, &
 209                         c_b6, &a[i__ * a_dim1 + 1], &c__1);
 210                 i__2 = *n - i__;
 211                 dgemv_("No transpose", &i__, &i__2, &c_b5, &w[(iw + 1) *
 212                         w_dim1 + 1], ldw, &a[i__ + (i__ + 1) * a_dim1], lda, &
 213                         c_b6, &a[i__ * a_dim1 + 1], &c__1);
 214             }
 215             if (i__ > 1) {
 216
 217 /*              Generate elementary reflector H(i) to annihilate */
 218 /*              A(1:i-2,i) */
 219
 220                 i__2 = i__ - 1;
 221                 dlarfg_(&i__2, &a[i__ - 1 + i__ * a_dim1], &a[i__ * a_dim1 +
 222                         1], &c__1, &tau[i__ - 1]);
 223                 e[i__ - 1] = a[i__ - 1 + i__ * a_dim1];
 224                 a[i__ - 1 + i__ * a_dim1] = 1.;
 225
 226 /*              Compute W(1:i-1,i) */
 227
 228                 i__2 = i__ - 1;
 229                 dsymv_("Upper", &i__2, &c_b6, &a[a_offset], lda, &a[i__ *
 230                         a_dim1 + 1], &c__1, &c_b16, &w[iw * w_dim1 + 1], &
 231                         c__1);
 232                 if (i__ < *n) {
 233                     i__2 = i__ - 1;
 234                     i__3 = *n - i__;
 235                     dgemv_("Transpose", &i__2, &i__3, &c_b6, &w[(iw + 1) *
 236                             w_dim1 + 1], ldw, &a[i__ * a_dim1 + 1], &c__1, &
 237                             c_b16, &w[i__ + 1 + iw * w_dim1], &c__1);
 238                     i__2 = i__ - 1;
 239                     i__3 = *n - i__;
 240                     dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[(i__ + 1) *
 241                              a_dim1 + 1], lda, &w[i__ + 1 + iw * w_dim1], &
 242                             c__1, &c_b6, &w[iw * w_dim1 + 1], &c__1);
 243                     i__2 = i__ - 1;
 244                     i__3 = *n - i__;
 245                     dgemv_("Transpose", &i__2, &i__3, &c_b6, &a[(i__ + 1) *
 246                             a_dim1 + 1], lda, &a[i__ * a_dim1 + 1], &c__1, &
 247                             c_b16, &w[i__ + 1 + iw * w_dim1], &c__1);
 248                     i__2 = i__ - 1;
 249                     i__3 = *n - i__;
 250                     dgemv_("No transpose", &i__2, &i__3, &c_b5, &w[(iw + 1) *
 251                             w_dim1 + 1], ldw, &w[i__ + 1 + iw * w_dim1], &
 252                             c__1, &c_b6, &w[iw * w_dim1 + 1], &c__1);
 253                 }
 254                 i__2 = i__ - 1;
 255                 dscal_(&i__2, &tau[i__ - 1], &w[iw * w_dim1 + 1], &c__1);
 256                 i__2 = i__ - 1;
 257                 alpha = tau[i__ - 1] * -.5 * ddot_(&i__2, &w[iw * w_dim1 + 1],
 258                          &c__1, &a[i__ * a_dim1 + 1], &c__1);
 259                 i__2 = i__ - 1;
 260                 daxpy_(&i__2, &alpha, &a[i__ * a_dim1 + 1], &c__1, &w[iw *
 261                         w_dim1 + 1], &c__1);
 262             }
 263
 264 /* L10: */
 265         }
 266     } else {
 267
 268 /*        Reduce first NB columns of lower triangle */
 269
 270         i__1 = *nb;
 271         for (i__ = 1; i__ <= i__1; ++i__) {
 272
 273 /*           Update A(i:n,i) */
 274
 275             i__2 = *n - i__ + 1;
 276             i__3 = i__ - 1;
 277             dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + a_dim1], lda,
 278                      &w[i__ + w_dim1], ldw, &c_b6, &a[i__ + i__ * a_dim1], &
 279                     c__1);
 280             i__2 = *n - i__ + 1;
 281             i__3 = i__ - 1;
 282             dgemv_("No transpose", &i__2, &i__3, &c_b5, &w[i__ + w_dim1], ldw,
 283                      &a[i__ + a_dim1], lda, &c_b6, &a[i__ + i__ * a_dim1], &
 284                     c__1);
 285             if (i__ < *n) {
 286
 287 /*              Generate elementary reflector H(i) to annihilate */
 288 /*              A(i+2:n,i) */
 289
 290                 i__2 = *n - i__;
 291 /* Computing MIN */
 292                 i__3 = i__ + 2;
 293                 dlarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *n)+
 294                         i__ * a_dim1], &c__1, &tau[i__]);
 295                 e[i__] = a[i__ + 1 + i__ * a_dim1];
 296                 a[i__ + 1 + i__ * a_dim1] = 1.;
 297
 298 /*              Compute W(i+1:n,i) */
 299
 300                 i__2 = *n - i__;
 301                 dsymv_("Lower", &i__2, &c_b6, &a[i__ + 1 + (i__ + 1) * a_dim1]
 302 , lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &w[
 303                         i__ + 1 + i__ * w_dim1], &c__1);
 304                 i__2 = *n - i__;
 305                 i__3 = i__ - 1;
 306                 dgemv_("Transpose", &i__2, &i__3, &c_b6, &w[i__ + 1 + w_dim1],
 307                          ldw, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &w[
 308                         i__ * w_dim1 + 1], &c__1);
 309                 i__2 = *n - i__;
 310                 i__3 = i__ - 1;
 311                 dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 +
 312                         a_dim1], lda, &w[i__ * w_dim1 + 1], &c__1, &c_b6, &w[
 313                         i__ + 1 + i__ * w_dim1], &c__1);
 314                 i__2 = *n - i__;
 315                 i__3 = i__ - 1;
 316                 dgemv_("Transpose", &i__2, &i__3, &c_b6, &a[i__ + 1 + a_dim1],
 317                          lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &w[
 318                         i__ * w_dim1 + 1], &c__1);
 319                 i__2 = *n - i__;
 320                 i__3 = i__ - 1;
 321                 dgemv_("No transpose", &i__2, &i__3, &c_b5, &w[i__ + 1 +
 322                         w_dim1], ldw, &w[i__ * w_dim1 + 1], &c__1, &c_b6, &w[
 323                         i__ + 1 + i__ * w_dim1], &c__1);
 324                 i__2 = *n - i__;
 325                 dscal_(&i__2, &tau[i__], &w[i__ + 1 + i__ * w_dim1], &c__1);
 326                 i__2 = *n - i__;
 327                 alpha = tau[i__] * -.5 * ddot_(&i__2, &w[i__ + 1 + i__ *
 328                         w_dim1], &c__1, &a[i__ + 1 + i__ * a_dim1], &c__1);
 329                 i__2 = *n - i__;
 330                 daxpy_(&i__2, &alpha, &a[i__ + 1 + i__ * a_dim1], &c__1, &w[
 331                         i__ + 1 + i__ * w_dim1], &c__1);
 332             }
 333
 334 /* L20: */
 335         }
 336     }
 337
 338     return 0;
 339
 340 /*     End of DLATRD */
 341
 342 } /* dlatrd_ */