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

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