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

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