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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6
   7 /* Subroutine */ int slasd7_(integer *icompq, integer *nl, integer *nr,
   8         integer *sqre, integer *k, real *d__, real *z__, real *zw, real *vf,
   9         real *vfw, real *vl, real *vlw, real *alpha, real *beta, real *dsigma,
  10          integer *idx, integer *idxp, integer *idxq, integer *perm, integer *
  11         givptr, integer *givcol, integer *ldgcol, real *givnum, integer *
  12         ldgnum, real *c__, real *s, integer *info)
  13 {
  14     /* System generated locals */
  15     integer givcol_dim1, givcol_offset, givnum_dim1, givnum_offset, i__1;
  16     real r__1, r__2;
  17
  18     /* Local variables */
  19     integer i__, j, m, n, k2;
  20     real z1;
  21     integer jp;
  22     real eps, tau, tol;
  23     integer nlp1, nlp2, idxi, idxj;
  24     extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
  25             integer *, real *, real *);
  26     integer idxjp, jprev;
  27     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  28             integer *);
  29     extern doublereal slapy2_(real *, real *), slamch_(char *);
  30     extern /* Subroutine */ int xerbla_(char *, integer *), slamrg_(
  31             integer *, integer *, real *, integer *, integer *, integer *);
  32     real hlftol;
  33
  34
  35 /*  -- LAPACK auxiliary 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 /*  SLASD7 merges the two sets of singular values together into a single */
  48 /*  sorted set. Then it tries to deflate the size of the problem. There */
  49 /*  are two ways in which deflation can occur:  when two or more singular */
  50 /*  values are close together or if there is a tiny entry in the Z */
  51 /*  vector. For each such occurrence the order of the related */
  52 /*  secular equation problem is reduced by one. */
  53
  54 /*  SLASD7 is called from SLASD6. */
  55
  56 /*  Arguments */
  57 /*  ========= */
  58
  59 /*  ICOMPQ  (input) INTEGER */
  60 /*          Specifies whether singular vectors are to be computed */
  61 /*          in compact form, as follows: */
  62 /*          = 0: Compute singular values only. */
  63 /*          = 1: Compute singular vectors of upper */
  64 /*               bidiagonal matrix in compact form. */
  65
  66 /*  NL     (input) INTEGER */
  67 /*         The row dimension of the upper block. NL >= 1. */
  68
  69 /*  NR     (input) INTEGER */
  70 /*         The row dimension of the lower block. NR >= 1. */
  71
  72 /*  SQRE   (input) INTEGER */
  73 /*         = 0: the lower block is an NR-by-NR square matrix. */
  74 /*         = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
  75
  76 /*         The bidiagonal matrix has */
  77 /*         N = NL + NR + 1 rows and */
  78 /*         M = N + SQRE >= N columns. */
  79
  80 /*  K      (output) INTEGER */
  81 /*         Contains the dimension of the non-deflated matrix, this is */
  82 /*         the order of the related secular equation. 1 <= K <=N. */
  83
  84 /*  D      (input/output) REAL array, dimension ( N ) */
  85 /*         On entry D contains the singular values of the two submatrices */
  86 /*         to be combined. On exit D contains the trailing (N-K) updated */
  87 /*         singular values (those which were deflated) sorted into */
  88 /*         increasing order. */
  89
  90 /*  Z      (output) REAL array, dimension ( M ) */
  91 /*         On exit Z contains the updating row vector in the secular */
  92 /*         equation. */
  93
  94 /*  ZW     (workspace) REAL array, dimension ( M ) */
  95 /*         Workspace for Z. */
  96
  97 /*  VF     (input/output) REAL array, dimension ( M ) */
  98 /*         On entry, VF(1:NL+1) contains the first components of all */
  99 /*         right singular vectors of the upper block; and VF(NL+2:M) */
 100 /*         contains the first components of all right singular vectors */
 101 /*         of the lower block. On exit, VF contains the first components */
 102 /*         of all right singular vectors of the bidiagonal matrix. */
 103
 104 /*  VFW    (workspace) REAL array, dimension ( M ) */
 105 /*         Workspace for VF. */
 106
 107 /*  VL     (input/output) REAL array, dimension ( M ) */
 108 /*         On entry, VL(1:NL+1) contains the  last components of all */
 109 /*         right singular vectors of the upper block; and VL(NL+2:M) */
 110 /*         contains the last components of all right singular vectors */
 111 /*         of the lower block. On exit, VL contains the last components */
 112 /*         of all right singular vectors of the bidiagonal matrix. */
 113
 114 /*  VLW    (workspace) REAL array, dimension ( M ) */
 115 /*         Workspace for VL. */
 116
 117 /*  ALPHA  (input) REAL */
 118 /*         Contains the diagonal element associated with the added row. */
 119
 120 /*  BETA   (input) REAL */
 121 /*         Contains the off-diagonal element associated with the added */
 122 /*         row. */
 123
 124 /*  DSIGMA (output) REAL array, dimension ( N ) */
 125 /*         Contains a copy of the diagonal elements (K-1 singular values */
 126 /*         and one zero) in the secular equation. */
 127
 128 /*  IDX    (workspace) INTEGER array, dimension ( N ) */
 129 /*         This will contain the permutation used to sort the contents of */
 130 /*         D into ascending order. */
 131
 132 /*  IDXP   (workspace) INTEGER array, dimension ( N ) */
 133 /*         This will contain the permutation used to place deflated */
 134 /*         values of D at the end of the array. On output IDXP(2:K) */
 135 /*         points to the nondeflated D-values and IDXP(K+1:N) */
 136 /*         points to the deflated singular values. */
 137
 138 /*  IDXQ   (input) INTEGER array, dimension ( N ) */
 139 /*         This contains the permutation which separately sorts the two */
 140 /*         sub-problems in D into ascending order.  Note that entries in */
 141 /*         the first half of this permutation must first be moved one */
 142 /*         position backward; and entries in the second half */
 143 /*         must first have NL+1 added to their values. */
 144
 145 /*  PERM   (output) INTEGER array, dimension ( N ) */
 146 /*         The permutations (from deflation and sorting) to be applied */
 147 /*         to each singular block. Not referenced if ICOMPQ = 0. */
 148
 149 /*  GIVPTR (output) INTEGER */
 150 /*         The number of Givens rotations which took place in this */
 151 /*         subproblem. Not referenced if ICOMPQ = 0. */
 152
 153 /*  GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) */
 154 /*         Each pair of numbers indicates a pair of columns to take place */
 155 /*         in a Givens rotation. Not referenced if ICOMPQ = 0. */
 156
 157 /*  LDGCOL (input) INTEGER */
 158 /*         The leading dimension of GIVCOL, must be at least N. */
 159
 160 /*  GIVNUM (output) REAL array, dimension ( LDGNUM, 2 ) */
 161 /*         Each number indicates the C or S value to be used in the */
 162 /*         corresponding Givens rotation. Not referenced if ICOMPQ = 0. */
 163
 164 /*  LDGNUM (input) INTEGER */
 165 /*         The leading dimension of GIVNUM, must be at least N. */
 166
 167 /*  C      (output) REAL */
 168 /*         C contains garbage if SQRE =0 and the C-value of a Givens */
 169 /*         rotation related to the right null space if SQRE = 1. */
 170
 171 /*  S      (output) REAL */
 172 /*         S contains garbage if SQRE =0 and the S-value of a Givens */
 173 /*         rotation related to the right null space if SQRE = 1. */
 174
 175 /*  INFO   (output) INTEGER */
 176 /*         = 0:  successful exit. */
 177 /*         < 0:  if INFO = -i, the i-th argument had an illegal value. */
 178
 179 /*  Further Details */
 180 /*  =============== */
 181
 182 /*  Based on contributions by */
 183 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 184 /*     California at Berkeley, USA */
 185
 186 /*  ===================================================================== */
 187
 188 /*     .. Parameters .. */
 189 /*     .. */
 190 /*     .. Local Scalars .. */
 191
 192 /*     .. */
 193 /*     .. External Subroutines .. */
 194 /*     .. */
 195 /*     .. External Functions .. */
 196 /*     .. */
 197 /*     .. Intrinsic Functions .. */
 198 /*     .. */
 199 /*     .. Executable Statements .. */
 200
 201 /*     Test the input parameters. */
 202
 203     /* Parameter adjustments */
 204     --d__;
 205     --z__;
 206     --zw;
 207     --vf;
 208     --vfw;
 209     --vl;
 210     --vlw;
 211     --dsigma;
 212     --idx;
 213     --idxp;
 214     --idxq;
 215     --perm;
 216     givcol_dim1 = *ldgcol;
 217     givcol_offset = 1 + givcol_dim1;
 218     givcol -= givcol_offset;
 219     givnum_dim1 = *ldgnum;
 220     givnum_offset = 1 + givnum_dim1;
 221     givnum -= givnum_offset;
 222
 223     /* Function Body */
 224     *info = 0;
 225     n = *nl + *nr + 1;
 226     m = n + *sqre;
 227
 228     if (*icompq < 0 || *icompq > 1) {
 229         *info = -1;
 230     } else if (*nl < 1) {
 231         *info = -2;
 232     } else if (*nr < 1) {
 233         *info = -3;
 234     } else if (*sqre < 0 || *sqre > 1) {
 235         *info = -4;
 236     } else if (*ldgcol < n) {
 237         *info = -22;
 238     } else if (*ldgnum < n) {
 239         *info = -24;
 240     }
 241     if (*info != 0) {
 242         i__1 = -(*info);
 243         xerbla_("SLASD7", &i__1);
 244         return 0;
 245     }
 246
 247     nlp1 = *nl + 1;
 248     nlp2 = *nl + 2;
 249     if (*icompq == 1) {
 250         *givptr = 0;
 251     }
 252
 253 /*     Generate the first part of the vector Z and move the singular */
 254 /*     values in the first part of D one position backward. */
 255
 256     z1 = *alpha * vl[nlp1];
 257     vl[nlp1] = 0.f;
 258     tau = vf[nlp1];
 259     for (i__ = *nl; i__ >= 1; --i__) {
 260         z__[i__ + 1] = *alpha * vl[i__];
 261         vl[i__] = 0.f;
 262         vf[i__ + 1] = vf[i__];
 263         d__[i__ + 1] = d__[i__];
 264         idxq[i__ + 1] = idxq[i__] + 1;
 265 /* L10: */
 266     }
 267     vf[1] = tau;
 268
 269 /*     Generate the second part of the vector Z. */
 270
 271     i__1 = m;
 272     for (i__ = nlp2; i__ <= i__1; ++i__) {
 273         z__[i__] = *beta * vf[i__];
 274         vf[i__] = 0.f;
 275 /* L20: */
 276     }
 277
 278 /*     Sort the singular values into increasing order */
 279
 280     i__1 = n;
 281     for (i__ = nlp2; i__ <= i__1; ++i__) {
 282         idxq[i__] += nlp1;
 283 /* L30: */
 284     }
 285
 286 /*     DSIGMA, IDXC, IDXC, and ZW are used as storage space. */
 287
 288     i__1 = n;
 289     for (i__ = 2; i__ <= i__1; ++i__) {
 290         dsigma[i__] = d__[idxq[i__]];
 291         zw[i__] = z__[idxq[i__]];
 292         vfw[i__] = vf[idxq[i__]];
 293         vlw[i__] = vl[idxq[i__]];
 294 /* L40: */
 295     }
 296
 297     slamrg_(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
 298
 299     i__1 = n;
 300     for (i__ = 2; i__ <= i__1; ++i__) {
 301         idxi = idx[i__] + 1;
 302         d__[i__] = dsigma[idxi];
 303         z__[i__] = zw[idxi];
 304         vf[i__] = vfw[idxi];
 305         vl[i__] = vlw[idxi];
 306 /* L50: */
 307     }
 308
 309 /*     Calculate the allowable deflation tolerence */
 310
 311     eps = slamch_("Epsilon");
 312 /* Computing MAX */
 313     r__1 = dabs(*alpha), r__2 = dabs(*beta);
 314     tol = dmax(r__1,r__2);
 315 /* Computing MAX */
 316     r__2 = (r__1 = d__[n], dabs(r__1));
 317     tol = eps * 64.f * dmax(r__2,tol);
 318
 319 /*     There are 2 kinds of deflation -- first a value in the z-vector */
 320 /*     is small, second two (or more) singular values are very close */
 321 /*     together (their difference is small). */
 322
 323 /*     If the value in the z-vector is small, we simply permute the */
 324 /*     array so that the corresponding singular value is moved to the */
 325 /*     end. */
 326
 327 /*     If two values in the D-vector are close, we perform a two-sided */
 328 /*     rotation designed to make one of the corresponding z-vector */
 329 /*     entries zero, and then permute the array so that the deflated */
 330 /*     singular value is moved to the end. */
 331
 332 /*     If there are multiple singular values then the problem deflates. */
 333 /*     Here the number of equal singular values are found.  As each equal */
 334 /*     singular value is found, an elementary reflector is computed to */
 335 /*     rotate the corresponding singular subspace so that the */
 336 /*     corresponding components of Z are zero in this new basis. */
 337
 338     *k = 1;
 339     k2 = n + 1;
 340     i__1 = n;
 341     for (j = 2; j <= i__1; ++j) {
 342         if ((r__1 = z__[j], dabs(r__1)) <= tol) {
 343
 344 /*           Deflate due to small z component. */
 345
 346             --k2;
 347             idxp[k2] = j;
 348             if (j == n) {
 349                 goto L100;
 350             }
 351         } else {
 352             jprev = j;
 353             goto L70;
 354         }
 355 /* L60: */
 356     }
 357 L70:
 358     j = jprev;
 359 L80:
 360     ++j;
 361     if (j > n) {
 362         goto L90;
 363     }
 364     if ((r__1 = z__[j], dabs(r__1)) <= tol) {
 365
 366 /*        Deflate due to small z component. */
 367
 368         --k2;
 369         idxp[k2] = j;
 370     } else {
 371
 372 /*        Check if singular values are close enough to allow deflation. */
 373
 374         if ((r__1 = d__[j] - d__[jprev], dabs(r__1)) <= tol) {
 375
 376 /*           Deflation is possible. */
 377
 378             *s = z__[jprev];
 379             *c__ = z__[j];
 380
 381 /*           Find sqrt(a**2+b**2) without overflow or */
 382 /*           destructive underflow. */
 383
 384             tau = slapy2_(c__, s);
 385             z__[j] = tau;
 386             z__[jprev] = 0.f;
 387             *c__ /= tau;
 388             *s = -(*s) / tau;
 389
 390 /*           Record the appropriate Givens rotation */
 391
 392             if (*icompq == 1) {
 393                 ++(*givptr);
 394                 idxjp = idxq[idx[jprev] + 1];
 395                 idxj = idxq[idx[j] + 1];
 396                 if (idxjp <= nlp1) {
 397                     --idxjp;
 398                 }
 399                 if (idxj <= nlp1) {
 400                     --idxj;
 401                 }
 402                 givcol[*givptr + (givcol_dim1 << 1)] = idxjp;
 403                 givcol[*givptr + givcol_dim1] = idxj;
 404                 givnum[*givptr + (givnum_dim1 << 1)] = *c__;
 405                 givnum[*givptr + givnum_dim1] = *s;
 406             }
 407             srot_(&c__1, &vf[jprev], &c__1, &vf[j], &c__1, c__, s);
 408             srot_(&c__1, &vl[jprev], &c__1, &vl[j], &c__1, c__, s);
 409             --k2;
 410             idxp[k2] = jprev;
 411             jprev = j;
 412         } else {
 413             ++(*k);
 414             zw[*k] = z__[jprev];
 415             dsigma[*k] = d__[jprev];
 416             idxp[*k] = jprev;
 417             jprev = j;
 418         }
 419     }
 420     goto L80;
 421 L90:
 422
 423 /*     Record the last singular value. */
 424
 425     ++(*k);
 426     zw[*k] = z__[jprev];
 427     dsigma[*k] = d__[jprev];
 428     idxp[*k] = jprev;
 429
 430 L100:
 431
 432 /*     Sort the singular values into DSIGMA. The singular values which */
 433 /*     were not deflated go into the first K slots of DSIGMA, except */
 434 /*     that DSIGMA(1) is treated separately. */
 435
 436     i__1 = n;
 437     for (j = 2; j <= i__1; ++j) {
 438         jp = idxp[j];
 439         dsigma[j] = d__[jp];
 440         vfw[j] = vf[jp];
 441         vlw[j] = vl[jp];
 442 /* L110: */
 443     }
 444     if (*icompq == 1) {
 445         i__1 = n;
 446         for (j = 2; j <= i__1; ++j) {
 447             jp = idxp[j];
 448             perm[j] = idxq[idx[jp] + 1];
 449             if (perm[j] <= nlp1) {
 450                 --perm[j];
 451             }
 452 /* L120: */
 453         }
 454     }
 455
 456 /*     The deflated singular values go back into the last N - K slots of */
 457 /*     D. */
 458
 459     i__1 = n - *k;
 460     scopy_(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
 461
 462 /*     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and */
 463 /*     VL(M). */
 464
 465     dsigma[1] = 0.f;
 466     hlftol = tol / 2.f;
 467     if (dabs(dsigma[2]) <= hlftol) {
 468         dsigma[2] = hlftol;
 469     }
 470     if (m > n) {
 471         z__[1] = slapy2_(&z1, &z__[m]);
 472         if (z__[1] <= tol) {
 473             *c__ = 1.f;
 474             *s = 0.f;
 475             z__[1] = tol;
 476         } else {
 477             *c__ = z1 / z__[1];
 478             *s = -z__[m] / z__[1];
 479         }
 480         srot_(&c__1, &vf[m], &c__1, &vf[1], &c__1, c__, s);
 481         srot_(&c__1, &vl[m], &c__1, &vl[1], &c__1, c__, s);
 482     } else {
 483         if (dabs(z1) <= tol) {
 484             z__[1] = tol;
 485         } else {
 486             z__[1] = z1;
 487         }
 488     }
 489
 490 /*     Restore Z, VF, and VL. */
 491
 492     i__1 = *k - 1;
 493     scopy_(&i__1, &zw[2], &c__1, &z__[2], &c__1);
 494     i__1 = n - 1;
 495     scopy_(&i__1, &vfw[2], &c__1, &vf[2], &c__1);
 496     i__1 = n - 1;
 497     scopy_(&i__1, &vlw[2], &c__1, &vl[2], &c__1);
 498
 499     return 0;
 500
 501 /*     End of SLASD7 */
 502
 503 } /* slasd7_ */