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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static real c_b3 = -1.f;
   6 static integer c__1 = 1;
   7
   8 /* Subroutine */ int slaed8_(integer *icompq, integer *k, integer *n, integer
   9         *qsiz, real *d__, real *q, integer *ldq, integer *indxq, real *rho,
  10         integer *cutpnt, real *z__, real *dlamda, real *q2, integer *ldq2,
  11         real *w, integer *perm, integer *givptr, integer *givcol, real *
  12         givnum, integer *indxp, integer *indx, integer *info)
  13 {
  14     /* System generated locals */
  15     integer q_dim1, q_offset, q2_dim1, q2_offset, i__1;
  16     real r__1;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal);
  20
  21     /* Local variables */
  22     real c__;
  23     integer i__, j;
  24     real s, t;
  25     integer k2, n1, n2, jp, n1p1;
  26     real eps, tau, tol;
  27     integer jlam, imax, jmax;
  28     extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
  29             integer *, real *, real *), sscal_(integer *, real *, real *,
  30             integer *), scopy_(integer *, real *, integer *, real *, integer *
  31 );
  32     extern doublereal slapy2_(real *, real *), slamch_(char *);
  33     extern /* Subroutine */ int xerbla_(char *, integer *);
  34     extern integer isamax_(integer *, real *, integer *);
  35     extern /* Subroutine */ int slamrg_(integer *, integer *, real *, integer
  36             *, integer *, integer *), slacpy_(char *, integer *, integer *,
  37             real *, integer *, real *, integer *);
  38
  39
  40 /*  -- LAPACK routine (version 3.1) -- */
  41 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  42 /*     November 2006 */
  43
  44 /*     .. Scalar Arguments .. */
  45 /*     .. */
  46 /*     .. Array Arguments .. */
  47 /*     .. */
  48
  49 /*  Purpose */
  50 /*  ======= */
  51
  52 /*  SLAED8 merges the two sets of eigenvalues together into a single */
  53 /*  sorted set.  Then it tries to deflate the size of the problem. */
  54 /*  There are two ways in which deflation can occur:  when two or more */
  55 /*  eigenvalues are close together or if there is a tiny element in the */
  56 /*  Z vector.  For each such occurrence the order of the related secular */
  57 /*  equation problem is reduced by one. */
  58
  59 /*  Arguments */
  60 /*  ========= */
  61
  62 /*  ICOMPQ  (input) INTEGER */
  63 /*          = 0:  Compute eigenvalues only. */
  64 /*          = 1:  Compute eigenvectors of original dense symmetric matrix */
  65 /*                also.  On entry, Q contains the orthogonal matrix used */
  66 /*                to reduce the original matrix to tridiagonal form. */
  67
  68 /*  K      (output) INTEGER */
  69 /*         The number of non-deflated eigenvalues, and the order of the */
  70 /*         related secular equation. */
  71
  72 /*  N      (input) INTEGER */
  73 /*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */
  74
  75 /*  QSIZ   (input) INTEGER */
  76 /*         The dimension of the orthogonal matrix used to reduce */
  77 /*         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1. */
  78
  79 /*  D      (input/output) REAL array, dimension (N) */
  80 /*         On entry, the eigenvalues of the two submatrices to be */
  81 /*         combined.  On exit, the trailing (N-K) updated eigenvalues */
  82 /*         (those which were deflated) sorted into increasing order. */
  83
  84 /*  Q      (input/output) REAL array, dimension (LDQ,N) */
  85 /*         If ICOMPQ = 0, Q is not referenced.  Otherwise, */
  86 /*         on entry, Q contains the eigenvectors of the partially solved */
  87 /*         system which has been previously updated in matrix */
  88 /*         multiplies with other partially solved eigensystems. */
  89 /*         On exit, Q contains the trailing (N-K) updated eigenvectors */
  90 /*         (those which were deflated) in its last N-K columns. */
  91
  92 /*  LDQ    (input) INTEGER */
  93 /*         The leading dimension of the array Q.  LDQ >= max(1,N). */
  94
  95 /*  INDXQ  (input) INTEGER array, dimension (N) */
  96 /*         The permutation which separately sorts the two sub-problems */
  97 /*         in D into ascending order.  Note that elements in the second */
  98 /*         half of this permutation must first have CUTPNT added to */
  99 /*         their values in order to be accurate. */
 100
 101 /*  RHO    (input/output) REAL */
 102 /*         On entry, the off-diagonal element associated with the rank-1 */
 103 /*         cut which originally split the two submatrices which are now */
 104 /*         being recombined. */
 105 /*         On exit, RHO has been modified to the value required by */
 106 /*         SLAED3. */
 107
 108 /*  CUTPNT (input) INTEGER */
 109 /*         The location of the last eigenvalue in the leading */
 110 /*         sub-matrix.  min(1,N) <= CUTPNT <= N. */
 111
 112 /*  Z      (input) REAL array, dimension (N) */
 113 /*         On entry, Z contains the updating vector (the last row of */
 114 /*         the first sub-eigenvector matrix and the first row of the */
 115 /*         second sub-eigenvector matrix). */
 116 /*         On exit, the contents of Z are destroyed by the updating */
 117 /*         process. */
 118
 119 /*  DLAMDA (output) REAL array, dimension (N) */
 120 /*         A copy of the first K eigenvalues which will be used by */
 121 /*         SLAED3 to form the secular equation. */
 122
 123 /*  Q2     (output) REAL array, dimension (LDQ2,N) */
 124 /*         If ICOMPQ = 0, Q2 is not referenced.  Otherwise, */
 125 /*         a copy of the first K eigenvectors which will be used by */
 126 /*         SLAED7 in a matrix multiply (SGEMM) to update the new */
 127 /*         eigenvectors. */
 128
 129 /*  LDQ2   (input) INTEGER */
 130 /*         The leading dimension of the array Q2.  LDQ2 >= max(1,N). */
 131
 132 /*  W      (output) REAL array, dimension (N) */
 133 /*         The first k values of the final deflation-altered z-vector and */
 134 /*         will be passed to SLAED3. */
 135
 136 /*  PERM   (output) INTEGER array, dimension (N) */
 137 /*         The permutations (from deflation and sorting) to be applied */
 138 /*         to each eigenblock. */
 139
 140 /*  GIVPTR (output) INTEGER */
 141 /*         The number of Givens rotations which took place in this */
 142 /*         subproblem. */
 143
 144 /*  GIVCOL (output) INTEGER array, dimension (2, N) */
 145 /*         Each pair of numbers indicates a pair of columns to take place */
 146 /*         in a Givens rotation. */
 147
 148 /*  GIVNUM (output) REAL array, dimension (2, N) */
 149 /*         Each number indicates the S value to be used in the */
 150 /*         corresponding Givens rotation. */
 151
 152 /*  INDXP  (workspace) INTEGER array, dimension (N) */
 153 /*         The permutation used to place deflated values of D at the end */
 154 /*         of the array.  INDXP(1:K) points to the nondeflated D-values */
 155 /*         and INDXP(K+1:N) points to the deflated eigenvalues. */
 156
 157 /*  INDX   (workspace) INTEGER array, dimension (N) */
 158 /*         The permutation used to sort the contents of D into ascending */
 159 /*         order. */
 160
 161 /*  INFO   (output) INTEGER */
 162 /*          = 0:  successful exit. */
 163 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 164
 165 /*  Further Details */
 166 /*  =============== */
 167
 168 /*  Based on contributions by */
 169 /*     Jeff Rutter, Computer Science Division, University of California */
 170 /*     at Berkeley, USA */
 171
 172 /*  ===================================================================== */
 173
 174 /*     .. Parameters .. */
 175 /*     .. */
 176 /*     .. Local Scalars .. */
 177
 178 /*     .. */
 179 /*     .. External Functions .. */
 180 /*     .. */
 181 /*     .. External Subroutines .. */
 182 /*     .. */
 183 /*     .. Intrinsic Functions .. */
 184 /*     .. */
 185 /*     .. Executable Statements .. */
 186
 187 /*     Test the input parameters. */
 188
 189     /* Parameter adjustments */
 190     --d__;
 191     q_dim1 = *ldq;
 192     q_offset = 1 + q_dim1;
 193     q -= q_offset;
 194     --indxq;
 195     --z__;
 196     --dlamda;
 197     q2_dim1 = *ldq2;
 198     q2_offset = 1 + q2_dim1;
 199     q2 -= q2_offset;
 200     --w;
 201     --perm;
 202     givcol -= 3;
 203     givnum -= 3;
 204     --indxp;
 205     --indx;
 206
 207     /* Function Body */
 208     *info = 0;
 209
 210     if (*icompq < 0 || *icompq > 1) {
 211         *info = -1;
 212     } else if (*n < 0) {
 213         *info = -3;
 214     } else if (*icompq == 1 && *qsiz < *n) {
 215         *info = -4;
 216     } else if (*ldq < max(1,*n)) {
 217         *info = -7;
 218     } else if (*cutpnt < min(1,*n) || *cutpnt > *n) {
 219         *info = -10;
 220     } else if (*ldq2 < max(1,*n)) {
 221         *info = -14;
 222     }
 223     if (*info != 0) {
 224         i__1 = -(*info);
 225         xerbla_("SLAED8", &i__1);
 226         return 0;
 227     }
 228
 229 /*     Quick return if possible */
 230
 231     if (*n == 0) {
 232         return 0;
 233     }
 234
 235     n1 = *cutpnt;
 236     n2 = *n - n1;
 237     n1p1 = n1 + 1;
 238
 239     if (*rho < 0.f) {
 240         sscal_(&n2, &c_b3, &z__[n1p1], &c__1);
 241     }
 242
 243 /*     Normalize z so that norm(z) = 1 */
 244
 245     t = 1.f / sqrt(2.f);
 246     i__1 = *n;
 247     for (j = 1; j <= i__1; ++j) {
 248         indx[j] = j;
 249 /* L10: */
 250     }
 251     sscal_(n, &t, &z__[1], &c__1);
 252     *rho = (r__1 = *rho * 2.f, dabs(r__1));
 253
 254 /*     Sort the eigenvalues into increasing order */
 255
 256     i__1 = *n;
 257     for (i__ = *cutpnt + 1; i__ <= i__1; ++i__) {
 258         indxq[i__] += *cutpnt;
 259 /* L20: */
 260     }
 261     i__1 = *n;
 262     for (i__ = 1; i__ <= i__1; ++i__) {
 263         dlamda[i__] = d__[indxq[i__]];
 264         w[i__] = z__[indxq[i__]];
 265 /* L30: */
 266     }
 267     i__ = 1;
 268     j = *cutpnt + 1;
 269     slamrg_(&n1, &n2, &dlamda[1], &c__1, &c__1, &indx[1]);
 270     i__1 = *n;
 271     for (i__ = 1; i__ <= i__1; ++i__) {
 272         d__[i__] = dlamda[indx[i__]];
 273         z__[i__] = w[indx[i__]];
 274 /* L40: */
 275     }
 276
 277 /*     Calculate the allowable deflation tolerence */
 278
 279     imax = isamax_(n, &z__[1], &c__1);
 280     jmax = isamax_(n, &d__[1], &c__1);
 281     eps = slamch_("Epsilon");
 282     tol = eps * 8.f * (r__1 = d__[jmax], dabs(r__1));
 283
 284 /*     If the rank-1 modifier is small enough, no more needs to be done */
 285 /*     except to reorganize Q so that its columns correspond with the */
 286 /*     elements in D. */
 287
 288     if (*rho * (r__1 = z__[imax], dabs(r__1)) <= tol) {
 289         *k = 0;
 290         if (*icompq == 0) {
 291             i__1 = *n;
 292             for (j = 1; j <= i__1; ++j) {
 293                 perm[j] = indxq[indx[j]];
 294 /* L50: */
 295             }
 296         } else {
 297             i__1 = *n;
 298             for (j = 1; j <= i__1; ++j) {
 299                 perm[j] = indxq[indx[j]];
 300                 scopy_(qsiz, &q[perm[j] * q_dim1 + 1], &c__1, &q2[j * q2_dim1
 301                         + 1], &c__1);
 302 /* L60: */
 303             }
 304             slacpy_("A", qsiz, n, &q2[q2_dim1 + 1], ldq2, &q[q_dim1 + 1], ldq);
 305         }
 306         return 0;
 307     }
 308
 309 /*     If there are multiple eigenvalues then the problem deflates.  Here */
 310 /*     the number of equal eigenvalues are found.  As each equal */
 311 /*     eigenvalue is found, an elementary reflector is computed to rotate */
 312 /*     the corresponding eigensubspace so that the corresponding */
 313 /*     components of Z are zero in this new basis. */
 314
 315     *k = 0;
 316     *givptr = 0;
 317     k2 = *n + 1;
 318     i__1 = *n;
 319     for (j = 1; j <= i__1; ++j) {
 320         if (*rho * (r__1 = z__[j], dabs(r__1)) <= tol) {
 321
 322 /*           Deflate due to small z component. */
 323
 324             --k2;
 325             indxp[k2] = j;
 326             if (j == *n) {
 327                 goto L110;
 328             }
 329         } else {
 330             jlam = j;
 331             goto L80;
 332         }
 333 /* L70: */
 334     }
 335 L80:
 336     ++j;
 337     if (j > *n) {
 338         goto L100;
 339     }
 340     if (*rho * (r__1 = z__[j], dabs(r__1)) <= tol) {
 341
 342 /*        Deflate due to small z component. */
 343
 344         --k2;
 345         indxp[k2] = j;
 346     } else {
 347
 348 /*        Check if eigenvalues are close enough to allow deflation. */
 349
 350         s = z__[jlam];
 351         c__ = z__[j];
 352
 353 /*        Find sqrt(a**2+b**2) without overflow or */
 354 /*        destructive underflow. */
 355
 356         tau = slapy2_(&c__, &s);
 357         t = d__[j] - d__[jlam];
 358         c__ /= tau;
 359         s = -s / tau;
 360         if ((r__1 = t * c__ * s, dabs(r__1)) <= tol) {
 361
 362 /*           Deflation is possible. */
 363
 364             z__[j] = tau;
 365             z__[jlam] = 0.f;
 366
 367 /*           Record the appropriate Givens rotation */
 368
 369             ++(*givptr);
 370             givcol[(*givptr << 1) + 1] = indxq[indx[jlam]];
 371             givcol[(*givptr << 1) + 2] = indxq[indx[j]];
 372             givnum[(*givptr << 1) + 1] = c__;
 373             givnum[(*givptr << 1) + 2] = s;
 374             if (*icompq == 1) {
 375                 srot_(qsiz, &q[indxq[indx[jlam]] * q_dim1 + 1], &c__1, &q[
 376                         indxq[indx[j]] * q_dim1 + 1], &c__1, &c__, &s);
 377             }
 378             t = d__[jlam] * c__ * c__ + d__[j] * s * s;
 379             d__[j] = d__[jlam] * s * s + d__[j] * c__ * c__;
 380             d__[jlam] = t;
 381             --k2;
 382             i__ = 1;
 383 L90:
 384             if (k2 + i__ <= *n) {
 385                 if (d__[jlam] < d__[indxp[k2 + i__]]) {
 386                     indxp[k2 + i__ - 1] = indxp[k2 + i__];
 387                     indxp[k2 + i__] = jlam;
 388                     ++i__;
 389                     goto L90;
 390                 } else {
 391                     indxp[k2 + i__ - 1] = jlam;
 392                 }
 393             } else {
 394                 indxp[k2 + i__ - 1] = jlam;
 395             }
 396             jlam = j;
 397         } else {
 398             ++(*k);
 399             w[*k] = z__[jlam];
 400             dlamda[*k] = d__[jlam];
 401             indxp[*k] = jlam;
 402             jlam = j;
 403         }
 404     }
 405     goto L80;
 406 L100:
 407
 408 /*     Record the last eigenvalue. */
 409
 410     ++(*k);
 411     w[*k] = z__[jlam];
 412     dlamda[*k] = d__[jlam];
 413     indxp[*k] = jlam;
 414
 415 L110:
 416
 417 /*     Sort the eigenvalues and corresponding eigenvectors into DLAMDA */
 418 /*     and Q2 respectively.  The eigenvalues/vectors which were not */
 419 /*     deflated go into the first K slots of DLAMDA and Q2 respectively, */
 420 /*     while those which were deflated go into the last N - K slots. */
 421
 422     if (*icompq == 0) {
 423         i__1 = *n;
 424         for (j = 1; j <= i__1; ++j) {
 425             jp = indxp[j];
 426             dlamda[j] = d__[jp];
 427             perm[j] = indxq[indx[jp]];
 428 /* L120: */
 429         }
 430     } else {
 431         i__1 = *n;
 432         for (j = 1; j <= i__1; ++j) {
 433             jp = indxp[j];
 434             dlamda[j] = d__[jp];
 435             perm[j] = indxq[indx[jp]];
 436             scopy_(qsiz, &q[perm[j] * q_dim1 + 1], &c__1, &q2[j * q2_dim1 + 1]
 437 , &c__1);
 438 /* L130: */
 439         }
 440     }
 441
 442 /*     The deflated eigenvalues and their corresponding vectors go back */
 443 /*     into the last N - K slots of D and Q respectively. */
 444
 445     if (*k < *n) {
 446         if (*icompq == 0) {
 447             i__1 = *n - *k;
 448             scopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
 449         } else {
 450             i__1 = *n - *k;
 451             scopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
 452             i__1 = *n - *k;
 453             slacpy_("A", qsiz, &i__1, &q2[(*k + 1) * q2_dim1 + 1], ldq2, &q[(*
 454                     k + 1) * q_dim1 + 1], ldq);
 455         }
 456     }
 457
 458     return 0;
 459
 460 /*     End of SLAED8 */
 461
 462 } /* slaed8_ */