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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static doublereal c_b30 = 0.;
   7
   8 /* Subroutine */ int dlasd2_(integer *nl, integer *nr, integer *sqre, integer
   9         *k, doublereal *d__, doublereal *z__, doublereal *alpha, doublereal *
  10         beta, doublereal *u, integer *ldu, doublereal *vt, integer *ldvt,
  11         doublereal *dsigma, doublereal *u2, integer *ldu2, doublereal *vt2,
  12         integer *ldvt2, integer *idxp, integer *idx, integer *idxc, integer *
  13         idxq, integer *coltyp, integer *info)
  14 {
  15     /* System generated locals */
  16     integer u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1, vt_offset,
  17             vt2_dim1, vt2_offset, i__1;
  18     doublereal d__1, d__2;
  19
  20     /* Local variables */
  21     doublereal c__;
  22     integer i__, j, m, n;
  23     doublereal s;
  24     integer k2;
  25     doublereal z1;
  26     integer ct, jp;
  27     doublereal eps, tau, tol;
  28     integer psm[4], nlp1, nlp2, idxi, idxj;
  29     extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
  30             doublereal *, integer *, doublereal *, doublereal *);
  31     integer ctot[4], idxjp;
  32     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  33             doublereal *, integer *);
  34     integer jprev;
  35     extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);
  36     extern /* Subroutine */ int dlamrg_(integer *, integer *, doublereal *,
  37             integer *, integer *, integer *), dlacpy_(char *, integer *,
  38             integer *, doublereal *, integer *, doublereal *, integer *), dlaset_(char *, integer *, integer *, doublereal *,
  39             doublereal *, doublereal *, integer *), xerbla_(char *,
  40             integer *);
  41     doublereal hlftol;
  42
  43
  44 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  45 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  46 /*     November 2006 */
  47
  48 /*     .. Scalar Arguments .. */
  49 /*     .. */
  50 /*     .. Array Arguments .. */
  51 /*     .. */
  52
  53 /*  Purpose */
  54 /*  ======= */
  55
  56 /*  DLASD2 merges the two sets of singular values together into a single */
  57 /*  sorted set.  Then it tries to deflate the size of the problem. */
  58 /*  There are two ways in which deflation can occur:  when two or more */
  59 /*  singular values are close together or if there is a tiny entry in the */
  60 /*  Z vector.  For each such occurrence the order of the related secular */
  61 /*  equation problem is reduced by one. */
  62
  63 /*  DLASD2 is called from DLASD1. */
  64
  65 /*  Arguments */
  66 /*  ========= */
  67
  68 /*  NL     (input) INTEGER */
  69 /*         The row dimension of the upper block.  NL >= 1. */
  70
  71 /*  NR     (input) INTEGER */
  72 /*         The row dimension of the lower block.  NR >= 1. */
  73
  74 /*  SQRE   (input) INTEGER */
  75 /*         = 0: the lower block is an NR-by-NR square matrix. */
  76 /*         = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
  77
  78 /*         The bidiagonal matrix has N = NL + NR + 1 rows and */
  79 /*         M = N + SQRE >= N columns. */
  80
  81 /*  K      (output) INTEGER */
  82 /*         Contains the dimension of the non-deflated matrix, */
  83 /*         This is the order of the related secular equation. 1 <= K <=N. */
  84
  85 /*  D      (input/output) DOUBLE PRECISION array, dimension(N) */
  86 /*         On entry D contains the singular values of the two submatrices */
  87 /*         to be combined.  On exit D contains the trailing (N-K) updated */
  88 /*         singular values (those which were deflated) sorted into */
  89 /*         increasing order. */
  90
  91 /*  Z      (output) DOUBLE PRECISION array, dimension(N) */
  92 /*         On exit Z contains the updating row vector in the secular */
  93 /*         equation. */
  94
  95 /*  ALPHA  (input) DOUBLE PRECISION */
  96 /*         Contains the diagonal element associated with the added row. */
  97
  98 /*  BETA   (input) DOUBLE PRECISION */
  99 /*         Contains the off-diagonal element associated with the added */
 100 /*         row. */
 101
 102 /*  U      (input/output) DOUBLE PRECISION array, dimension(LDU,N) */
 103 /*         On entry U contains the left singular vectors of two */
 104 /*         submatrices in the two square blocks with corners at (1,1), */
 105 /*         (NL, NL), and (NL+2, NL+2), (N,N). */
 106 /*         On exit U contains the trailing (N-K) updated left singular */
 107 /*         vectors (those which were deflated) in its last N-K columns. */
 108
 109 /*  LDU    (input) INTEGER */
 110 /*         The leading dimension of the array U.  LDU >= N. */
 111
 112 /*  VT     (input/output) DOUBLE PRECISION array, dimension(LDVT,M) */
 113 /*         On entry VT' contains the right singular vectors of two */
 114 /*         submatrices in the two square blocks with corners at (1,1), */
 115 /*         (NL+1, NL+1), and (NL+2, NL+2), (M,M). */
 116 /*         On exit VT' contains the trailing (N-K) updated right singular */
 117 /*         vectors (those which were deflated) in its last N-K columns. */
 118 /*         In case SQRE =1, the last row of VT spans the right null */
 119 /*         space. */
 120
 121 /*  LDVT   (input) INTEGER */
 122 /*         The leading dimension of the array VT.  LDVT >= M. */
 123
 124 /*  DSIGMA (output) DOUBLE PRECISION 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 /*  U2     (output) DOUBLE PRECISION array, dimension(LDU2,N) */
 129 /*         Contains a copy of the first K-1 left singular vectors which */
 130 /*         will be used by DLASD3 in a matrix multiply (DGEMM) to solve */
 131 /*         for the new left singular vectors. U2 is arranged into four */
 132 /*         blocks. The first block contains a column with 1 at NL+1 and */
 133 /*         zero everywhere else; the second block contains non-zero */
 134 /*         entries only at and above NL; the third contains non-zero */
 135 /*         entries only below NL+1; and the fourth is dense. */
 136
 137 /*  LDU2   (input) INTEGER */
 138 /*         The leading dimension of the array U2.  LDU2 >= N. */
 139
 140 /*  VT2    (output) DOUBLE PRECISION array, dimension(LDVT2,N) */
 141 /*         VT2' contains a copy of the first K right singular vectors */
 142 /*         which will be used by DLASD3 in a matrix multiply (DGEMM) to */
 143 /*         solve for the new right singular vectors. VT2 is arranged into */
 144 /*         three blocks. The first block contains a row that corresponds */
 145 /*         to the special 0 diagonal element in SIGMA; the second block */
 146 /*         contains non-zeros only at and before NL +1; the third block */
 147 /*         contains non-zeros only at and after  NL +2. */
 148
 149 /*  LDVT2  (input) INTEGER */
 150 /*         The leading dimension of the array VT2.  LDVT2 >= M. */
 151
 152 /*  IDXP   (workspace) INTEGER array dimension(N) */
 153 /*         This will contain the permutation used to place deflated */
 154 /*         values of D at the end of the array. On output IDXP(2:K) */
 155 /*         points to the nondeflated D-values and IDXP(K+1:N) */
 156 /*         points to the deflated singular values. */
 157
 158 /*  IDX    (workspace) INTEGER array dimension(N) */
 159 /*         This will contain the permutation used to sort the contents of */
 160 /*         D into ascending order. */
 161
 162 /*  IDXC   (output) INTEGER array dimension(N) */
 163 /*         This will contain the permutation used to arrange the columns */
 164 /*         of the deflated U matrix into three groups:  the first group */
 165 /*         contains non-zero entries only at and above NL, the second */
 166 /*         contains non-zero entries only below NL+2, and the third is */
 167 /*         dense. */
 168
 169 /*  IDXQ   (input/output) INTEGER array dimension(N) */
 170 /*         This contains the permutation which separately sorts the two */
 171 /*         sub-problems in D into ascending order.  Note that entries in */
 172 /*         the first hlaf of this permutation must first be moved one */
 173 /*         position backward; and entries in the second half */
 174 /*         must first have NL+1 added to their values. */
 175
 176 /*  COLTYP (workspace/output) INTEGER array dimension(N) */
 177 /*         As workspace, this will contain a label which will indicate */
 178 /*         which of the following types a column in the U2 matrix or a */
 179 /*         row in the VT2 matrix is: */
 180 /*         1 : non-zero in the upper half only */
 181 /*         2 : non-zero in the lower half only */
 182 /*         3 : dense */
 183 /*         4 : deflated */
 184
 185 /*         On exit, it is an array of dimension 4, with COLTYP(I) being */
 186 /*         the dimension of the I-th type columns. */
 187
 188 /*  INFO   (output) INTEGER */
 189 /*          = 0:  successful exit. */
 190 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 191
 192 /*  Further Details */
 193 /*  =============== */
 194
 195 /*  Based on contributions by */
 196 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 197 /*     California at Berkeley, USA */
 198
 199 /*  ===================================================================== */
 200
 201 /*     .. Parameters .. */
 202 /*     .. */
 203 /*     .. Local Arrays .. */
 204 /*     .. */
 205 /*     .. Local Scalars .. */
 206 /*     .. */
 207 /*     .. External Functions .. */
 208 /*     .. */
 209 /*     .. External Subroutines .. */
 210 /*     .. */
 211 /*     .. Intrinsic Functions .. */
 212 /*     .. */
 213 /*     .. Executable Statements .. */
 214
 215 /*     Test the input parameters. */
 216
 217     /* Parameter adjustments */
 218     --d__;
 219     --z__;
 220     u_dim1 = *ldu;
 221     u_offset = 1 + u_dim1;
 222     u -= u_offset;
 223     vt_dim1 = *ldvt;
 224     vt_offset = 1 + vt_dim1;
 225     vt -= vt_offset;
 226     --dsigma;
 227     u2_dim1 = *ldu2;
 228     u2_offset = 1 + u2_dim1;
 229     u2 -= u2_offset;
 230     vt2_dim1 = *ldvt2;
 231     vt2_offset = 1 + vt2_dim1;
 232     vt2 -= vt2_offset;
 233     --idxp;
 234     --idx;
 235     --idxc;
 236     --idxq;
 237     --coltyp;
 238
 239     /* Function Body */
 240     *info = 0;
 241
 242     if (*nl < 1) {
 243         *info = -1;
 244     } else if (*nr < 1) {
 245         *info = -2;
 246     } else if (*sqre != 1 && *sqre != 0) {
 247         *info = -3;
 248     }
 249
 250     n = *nl + *nr + 1;
 251     m = n + *sqre;
 252
 253     if (*ldu < n) {
 254         *info = -10;
 255     } else if (*ldvt < m) {
 256         *info = -12;
 257     } else if (*ldu2 < n) {
 258         *info = -15;
 259     } else if (*ldvt2 < m) {
 260         *info = -17;
 261     }
 262     if (*info != 0) {
 263         i__1 = -(*info);
 264         xerbla_("DLASD2", &i__1);
 265         return 0;
 266     }
 267
 268     nlp1 = *nl + 1;
 269     nlp2 = *nl + 2;
 270
 271 /*     Generate the first part of the vector Z; and move the singular */
 272 /*     values in the first part of D one position backward. */
 273
 274     z1 = *alpha * vt[nlp1 + nlp1 * vt_dim1];
 275     z__[1] = z1;
 276     for (i__ = *nl; i__ >= 1; --i__) {
 277         z__[i__ + 1] = *alpha * vt[i__ + nlp1 * vt_dim1];
 278         d__[i__ + 1] = d__[i__];
 279         idxq[i__ + 1] = idxq[i__] + 1;
 280 /* L10: */
 281     }
 282
 283 /*     Generate the second part of the vector Z. */
 284
 285     i__1 = m;
 286     for (i__ = nlp2; i__ <= i__1; ++i__) {
 287         z__[i__] = *beta * vt[i__ + nlp2 * vt_dim1];
 288 /* L20: */
 289     }
 290
 291 /*     Initialize some reference arrays. */
 292
 293     i__1 = nlp1;
 294     for (i__ = 2; i__ <= i__1; ++i__) {
 295         coltyp[i__] = 1;
 296 /* L30: */
 297     }
 298     i__1 = n;
 299     for (i__ = nlp2; i__ <= i__1; ++i__) {
 300         coltyp[i__] = 2;
 301 /* L40: */
 302     }
 303
 304 /*     Sort the singular values into increasing order */
 305
 306     i__1 = n;
 307     for (i__ = nlp2; i__ <= i__1; ++i__) {
 308         idxq[i__] += nlp1;
 309 /* L50: */
 310     }
 311
 312 /*     DSIGMA, IDXC, IDXC, and the first column of U2 */
 313 /*     are used as storage space. */
 314
 315     i__1 = n;
 316     for (i__ = 2; i__ <= i__1; ++i__) {
 317         dsigma[i__] = d__[idxq[i__]];
 318         u2[i__ + u2_dim1] = z__[idxq[i__]];
 319         idxc[i__] = coltyp[idxq[i__]];
 320 /* L60: */
 321     }
 322
 323     dlamrg_(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
 324
 325     i__1 = n;
 326     for (i__ = 2; i__ <= i__1; ++i__) {
 327         idxi = idx[i__] + 1;
 328         d__[i__] = dsigma[idxi];
 329         z__[i__] = u2[idxi + u2_dim1];
 330         coltyp[i__] = idxc[idxi];
 331 /* L70: */
 332     }
 333
 334 /*     Calculate the allowable deflation tolerance */
 335
 336     eps = dlamch_("Epsilon");
 337 /* Computing MAX */
 338     d__1 = abs(*alpha), d__2 = abs(*beta);
 339     tol = max(d__1,d__2);
 340 /* Computing MAX */
 341     d__2 = (d__1 = d__[n], abs(d__1));
 342     tol = eps * 8. * max(d__2,tol);
 343
 344 /*     There are 2 kinds of deflation -- first a value in the z-vector */
 345 /*     is small, second two (or more) singular values are very close */
 346 /*     together (their difference is small). */
 347
 348 /*     If the value in the z-vector is small, we simply permute the */
 349 /*     array so that the corresponding singular value is moved to the */
 350 /*     end. */
 351
 352 /*     If two values in the D-vector are close, we perform a two-sided */
 353 /*     rotation designed to make one of the corresponding z-vector */
 354 /*     entries zero, and then permute the array so that the deflated */
 355 /*     singular value is moved to the end. */
 356
 357 /*     If there are multiple singular values then the problem deflates. */
 358 /*     Here the number of equal singular values are found.  As each equal */
 359 /*     singular value is found, an elementary reflector is computed to */
 360 /*     rotate the corresponding singular subspace so that the */
 361 /*     corresponding components of Z are zero in this new basis. */
 362
 363     *k = 1;
 364     k2 = n + 1;
 365     i__1 = n;
 366     for (j = 2; j <= i__1; ++j) {
 367         if ((d__1 = z__[j], abs(d__1)) <= tol) {
 368
 369 /*           Deflate due to small z component. */
 370
 371             --k2;
 372             idxp[k2] = j;
 373             coltyp[j] = 4;
 374             if (j == n) {
 375                 goto L120;
 376             }
 377         } else {
 378             jprev = j;
 379             goto L90;
 380         }
 381 /* L80: */
 382     }
 383 L90:
 384     j = jprev;
 385 L100:
 386     ++j;
 387     if (j > n) {
 388         goto L110;
 389     }
 390     if ((d__1 = z__[j], abs(d__1)) <= tol) {
 391
 392 /*        Deflate due to small z component. */
 393
 394         --k2;
 395         idxp[k2] = j;
 396         coltyp[j] = 4;
 397     } else {
 398
 399 /*        Check if singular values are close enough to allow deflation. */
 400
 401         if ((d__1 = d__[j] - d__[jprev], abs(d__1)) <= tol) {
 402
 403 /*           Deflation is possible. */
 404
 405             s = z__[jprev];
 406             c__ = z__[j];
 407
 408 /*           Find sqrt(a**2+b**2) without overflow or */
 409 /*           destructive underflow. */
 410
 411             tau = dlapy2_(&c__, &s);
 412             c__ /= tau;
 413             s = -s / tau;
 414             z__[j] = tau;
 415             z__[jprev] = 0.;
 416
 417 /*           Apply back the Givens rotation to the left and right */
 418 /*           singular vector matrices. */
 419
 420             idxjp = idxq[idx[jprev] + 1];
 421             idxj = idxq[idx[j] + 1];
 422             if (idxjp <= nlp1) {
 423                 --idxjp;
 424             }
 425             if (idxj <= nlp1) {
 426                 --idxj;
 427             }
 428             drot_(&n, &u[idxjp * u_dim1 + 1], &c__1, &u[idxj * u_dim1 + 1], &
 429                     c__1, &c__, &s);
 430             drot_(&m, &vt[idxjp + vt_dim1], ldvt, &vt[idxj + vt_dim1], ldvt, &
 431                     c__, &s);
 432             if (coltyp[j] != coltyp[jprev]) {
 433                 coltyp[j] = 3;
 434             }
 435             coltyp[jprev] = 4;
 436             --k2;
 437             idxp[k2] = jprev;
 438             jprev = j;
 439         } else {
 440             ++(*k);
 441             u2[*k + u2_dim1] = z__[jprev];
 442             dsigma[*k] = d__[jprev];
 443             idxp[*k] = jprev;
 444             jprev = j;
 445         }
 446     }
 447     goto L100;
 448 L110:
 449
 450 /*     Record the last singular value. */
 451
 452     ++(*k);
 453     u2[*k + u2_dim1] = z__[jprev];
 454     dsigma[*k] = d__[jprev];
 455     idxp[*k] = jprev;
 456
 457 L120:
 458
 459 /*     Count up the total number of the various types of columns, then */
 460 /*     form a permutation which positions the four column types into */
 461 /*     four groups of uniform structure (although one or more of these */
 462 /*     groups may be empty). */
 463
 464     for (j = 1; j <= 4; ++j) {
 465         ctot[j - 1] = 0;
 466 /* L130: */
 467     }
 468     i__1 = n;
 469     for (j = 2; j <= i__1; ++j) {
 470         ct = coltyp[j];
 471         ++ctot[ct - 1];
 472 /* L140: */
 473     }
 474
 475 /*     PSM(*) = Position in SubMatrix (of types 1 through 4) */
 476
 477     psm[0] = 2;
 478     psm[1] = ctot[0] + 2;
 479     psm[2] = psm[1] + ctot[1];
 480     psm[3] = psm[2] + ctot[2];
 481
 482 /*     Fill out the IDXC array so that the permutation which it induces */
 483 /*     will place all type-1 columns first, all type-2 columns next, */
 484 /*     then all type-3's, and finally all type-4's, starting from the */
 485 /*     second column. This applies similarly to the rows of VT. */
 486
 487     i__1 = n;
 488     for (j = 2; j <= i__1; ++j) {
 489         jp = idxp[j];
 490         ct = coltyp[jp];
 491         idxc[psm[ct - 1]] = j;
 492         ++psm[ct - 1];
 493 /* L150: */
 494     }
 495
 496 /*     Sort the singular values and corresponding singular vectors into */
 497 /*     DSIGMA, U2, and VT2 respectively.  The singular values/vectors */
 498 /*     which were not deflated go into the first K slots of DSIGMA, U2, */
 499 /*     and VT2 respectively, while those which were deflated go into the */
 500 /*     last N - K slots, except that the first column/row will be treated */
 501 /*     separately. */
 502
 503     i__1 = n;
 504     for (j = 2; j <= i__1; ++j) {
 505         jp = idxp[j];
 506         dsigma[j] = d__[jp];
 507         idxj = idxq[idx[idxp[idxc[j]]] + 1];
 508         if (idxj <= nlp1) {
 509             --idxj;
 510         }
 511         dcopy_(&n, &u[idxj * u_dim1 + 1], &c__1, &u2[j * u2_dim1 + 1], &c__1);
 512         dcopy_(&m, &vt[idxj + vt_dim1], ldvt, &vt2[j + vt2_dim1], ldvt2);
 513 /* L160: */
 514     }
 515
 516 /*     Determine DSIGMA(1), DSIGMA(2) and Z(1) */
 517
 518     dsigma[1] = 0.;
 519     hlftol = tol / 2.;
 520     if (abs(dsigma[2]) <= hlftol) {
 521         dsigma[2] = hlftol;
 522     }
 523     if (m > n) {
 524         z__[1] = dlapy2_(&z1, &z__[m]);
 525         if (z__[1] <= tol) {
 526             c__ = 1.;
 527             s = 0.;
 528             z__[1] = tol;
 529         } else {
 530             c__ = z1 / z__[1];
 531             s = z__[m] / z__[1];
 532         }
 533     } else {
 534         if (abs(z1) <= tol) {
 535             z__[1] = tol;
 536         } else {
 537             z__[1] = z1;
 538         }
 539     }
 540
 541 /*     Move the rest of the updating row to Z. */
 542
 543     i__1 = *k - 1;
 544     dcopy_(&i__1, &u2[u2_dim1 + 2], &c__1, &z__[2], &c__1);
 545
 546 /*     Determine the first column of U2, the first row of VT2 and the */
 547 /*     last row of VT. */
 548
 549     dlaset_("A", &n, &c__1, &c_b30, &c_b30, &u2[u2_offset], ldu2);
 550     u2[nlp1 + u2_dim1] = 1.;
 551     if (m > n) {
 552         i__1 = nlp1;
 553         for (i__ = 1; i__ <= i__1; ++i__) {
 554             vt[m + i__ * vt_dim1] = -s * vt[nlp1 + i__ * vt_dim1];
 555             vt2[i__ * vt2_dim1 + 1] = c__ * vt[nlp1 + i__ * vt_dim1];
 556 /* L170: */
 557         }
 558         i__1 = m;
 559         for (i__ = nlp2; i__ <= i__1; ++i__) {
 560             vt2[i__ * vt2_dim1 + 1] = s * vt[m + i__ * vt_dim1];
 561             vt[m + i__ * vt_dim1] = c__ * vt[m + i__ * vt_dim1];
 562 /* L180: */
 563         }
 564     } else {
 565         dcopy_(&m, &vt[nlp1 + vt_dim1], ldvt, &vt2[vt2_dim1 + 1], ldvt2);
 566     }
 567     if (m > n) {
 568         dcopy_(&m, &vt[m + vt_dim1], ldvt, &vt2[m + vt2_dim1], ldvt2);
 569     }
 570
 571 /*     The deflated singular values and their corresponding vectors go */
 572 /*     into the back of D, U, and V respectively. */
 573
 574     if (n > *k) {
 575         i__1 = n - *k;
 576         dcopy_(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
 577         i__1 = n - *k;
 578         dlacpy_("A", &n, &i__1, &u2[(*k + 1) * u2_dim1 + 1], ldu2, &u[(*k + 1)
 579                  * u_dim1 + 1], ldu);
 580         i__1 = n - *k;
 581         dlacpy_("A", &i__1, &m, &vt2[*k + 1 + vt2_dim1], ldvt2, &vt[*k + 1 +
 582                 vt_dim1], ldvt);
 583     }
 584
 585 /*     Copy CTOT into COLTYP for referencing in DLASD3. */
 586
 587     for (j = 1; j <= 4; ++j) {
 588         coltyp[j] = ctot[j - 1];
 589 /* L190: */
 590     }
 591
 592     return 0;
 593
 594 /*     End of DLASD2 */
 595
 596 } /* dlasd2_ */