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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c__0 = 0;
   7 static doublereal c_b8 = 1.;
   8
   9 /* Subroutine */ int dlasd8_(integer *icompq, integer *k, doublereal *d__,
  10         doublereal *z__, doublereal *vf, doublereal *vl, doublereal *difl,
  11         doublereal *difr, integer *lddifr, doublereal *dsigma, doublereal *
  12         work, integer *info)
  13 {
  14     /* System generated locals */
  15     integer difr_dim1, difr_offset, i__1, i__2;
  16     doublereal d__1, d__2;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
  20
  21     /* Local variables */
  22     integer i__, j;
  23     doublereal dj, rho;
  24     integer iwk1, iwk2, iwk3;
  25     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  26             integer *);
  27     doublereal temp;
  28     extern doublereal dnrm2_(integer *, doublereal *, integer *);
  29     integer iwk2i, iwk3i;
  30     doublereal diflj, difrj, dsigj;
  31     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  32             doublereal *, integer *);
  33     extern doublereal dlamc3_(doublereal *, doublereal *);
  34     extern /* Subroutine */ int dlasd4_(integer *, integer *, doublereal *,
  35             doublereal *, doublereal *, doublereal *, doublereal *,
  36             doublereal *, integer *), dlascl_(char *, integer *, integer *,
  37             doublereal *, doublereal *, integer *, integer *, doublereal *,
  38             integer *, integer *), dlaset_(char *, integer *, integer
  39             *, doublereal *, doublereal *, doublereal *, integer *),
  40             xerbla_(char *, integer *);
  41     doublereal dsigjp;
  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 /*  DLASD8 finds the square roots of the roots of the secular equation, */
  57 /*  as defined by the values in DSIGMA and Z. It makes the appropriate */
  58 /*  calls to DLASD4, and stores, for each  element in D, the distance */
  59 /*  to its two nearest poles (elements in DSIGMA). It also updates */
  60 /*  the arrays VF and VL, the first and last components of all the */
  61 /*  right singular vectors of the original bidiagonal matrix. */
  62
  63 /*  DLASD8 is called from DLASD6. */
  64
  65 /*  Arguments */
  66 /*  ========= */
  67
  68 /*  ICOMPQ  (input) INTEGER */
  69 /*          Specifies whether singular vectors are to be computed in */
  70 /*          factored form in the calling routine: */
  71 /*          = 0: Compute singular values only. */
  72 /*          = 1: Compute singular vectors in factored form as well. */
  73
  74 /*  K       (input) INTEGER */
  75 /*          The number of terms in the rational function to be solved */
  76 /*          by DLASD4.  K >= 1. */
  77
  78 /*  D       (output) DOUBLE PRECISION array, dimension ( K ) */
  79 /*          On output, D contains the updated singular values. */
  80
  81 /*  Z       (input) DOUBLE PRECISION array, dimension ( K ) */
  82 /*          The first K elements of this array contain the components */
  83 /*          of the deflation-adjusted updating row vector. */
  84
  85 /*  VF      (input/output) DOUBLE PRECISION array, dimension ( K ) */
  86 /*          On entry, VF contains  information passed through DBEDE8. */
  87 /*          On exit, VF contains the first K components of the first */
  88 /*          components of all right singular vectors of the bidiagonal */
  89 /*          matrix. */
  90
  91 /*  VL      (input/output) DOUBLE PRECISION array, dimension ( K ) */
  92 /*          On entry, VL contains  information passed through DBEDE8. */
  93 /*          On exit, VL contains the first K components of the last */
  94 /*          components of all right singular vectors of the bidiagonal */
  95 /*          matrix. */
  96
  97 /*  DIFL    (output) DOUBLE PRECISION array, dimension ( K ) */
  98 /*          On exit, DIFL(I) = D(I) - DSIGMA(I). */
  99
 100 /*  DIFR    (output) DOUBLE PRECISION array, */
 101 /*                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and */
 102 /*                   dimension ( K ) if ICOMPQ = 0. */
 103 /*          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not */
 104 /*          defined and will not be referenced. */
 105
 106 /*          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the */
 107 /*          normalizing factors for the right singular vector matrix. */
 108
 109 /*  LDDIFR  (input) INTEGER */
 110 /*          The leading dimension of DIFR, must be at least K. */
 111
 112 /*  DSIGMA  (input) DOUBLE PRECISION array, dimension ( K ) */
 113 /*          The first K elements of this array contain the old roots */
 114 /*          of the deflated updating problem.  These are the poles */
 115 /*          of the secular equation. */
 116
 117 /*  WORK    (workspace) DOUBLE PRECISION array, dimension at least 3 * K */
 118
 119 /*  INFO    (output) INTEGER */
 120 /*          = 0:  successful exit. */
 121 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 122 /*          > 0:  if INFO = 1, an singular value did not converge */
 123
 124 /*  Further Details */
 125 /*  =============== */
 126
 127 /*  Based on contributions by */
 128 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 129 /*     California at Berkeley, USA */
 130
 131 /*  ===================================================================== */
 132
 133 /*     .. Parameters .. */
 134 /*     .. */
 135 /*     .. Local Scalars .. */
 136 /*     .. */
 137 /*     .. External Subroutines .. */
 138 /*     .. */
 139 /*     .. External Functions .. */
 140 /*     .. */
 141 /*     .. Intrinsic Functions .. */
 142 /*     .. */
 143 /*     .. Executable Statements .. */
 144
 145 /*     Test the input parameters. */
 146
 147     /* Parameter adjustments */
 148     --d__;
 149     --z__;
 150     --vf;
 151     --vl;
 152     --difl;
 153     difr_dim1 = *lddifr;
 154     difr_offset = 1 + difr_dim1;
 155     difr -= difr_offset;
 156     --dsigma;
 157     --work;
 158
 159     /* Function Body */
 160     *info = 0;
 161
 162     if (*icompq < 0 || *icompq > 1) {
 163         *info = -1;
 164     } else if (*k < 1) {
 165         *info = -2;
 166     } else if (*lddifr < *k) {
 167         *info = -9;
 168     }
 169     if (*info != 0) {
 170         i__1 = -(*info);
 171         xerbla_("DLASD8", &i__1);
 172         return 0;
 173     }
 174
 175 /*     Quick return if possible */
 176
 177     if (*k == 1) {
 178         d__[1] = abs(z__[1]);
 179         difl[1] = d__[1];
 180         if (*icompq == 1) {
 181             difl[2] = 1.;
 182             difr[(difr_dim1 << 1) + 1] = 1.;
 183         }
 184         return 0;
 185     }
 186
 187 /*     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can */
 188 /*     be computed with high relative accuracy (barring over/underflow). */
 189 /*     This is a problem on machines without a guard digit in */
 190 /*     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). */
 191 /*     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), */
 192 /*     which on any of these machines zeros out the bottommost */
 193 /*     bit of DSIGMA(I) if it is 1; this makes the subsequent */
 194 /*     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation */
 195 /*     occurs. On binary machines with a guard digit (almost all */
 196 /*     machines) it does not change DSIGMA(I) at all. On hexadecimal */
 197 /*     and decimal machines with a guard digit, it slightly */
 198 /*     changes the bottommost bits of DSIGMA(I). It does not account */
 199 /*     for hexadecimal or decimal machines without guard digits */
 200 /*     (we know of none). We use a subroutine call to compute */
 201 /*     2*DSIGMA(I) to prevent optimizing compilers from eliminating */
 202 /*     this code. */
 203
 204     i__1 = *k;
 205     for (i__ = 1; i__ <= i__1; ++i__) {
 206         dsigma[i__] = dlamc3_(&dsigma[i__], &dsigma[i__]) - dsigma[i__];
 207 /* L10: */
 208     }
 209
 210 /*     Book keeping. */
 211
 212     iwk1 = 1;
 213     iwk2 = iwk1 + *k;
 214     iwk3 = iwk2 + *k;
 215     iwk2i = iwk2 - 1;
 216     iwk3i = iwk3 - 1;
 217
 218 /*     Normalize Z. */
 219
 220     rho = dnrm2_(k, &z__[1], &c__1);
 221     dlascl_("G", &c__0, &c__0, &rho, &c_b8, k, &c__1, &z__[1], k, info);
 222     rho *= rho;
 223
 224 /*     Initialize WORK(IWK3). */
 225
 226     dlaset_("A", k, &c__1, &c_b8, &c_b8, &work[iwk3], k);
 227
 228 /*     Compute the updated singular values, the arrays DIFL, DIFR, */
 229 /*     and the updated Z. */
 230
 231     i__1 = *k;
 232     for (j = 1; j <= i__1; ++j) {
 233         dlasd4_(k, &j, &dsigma[1], &z__[1], &work[iwk1], &rho, &d__[j], &work[
 234                 iwk2], info);
 235
 236 /*        If the root finder fails, the computation is terminated. */
 237
 238         if (*info != 0) {
 239             return 0;
 240         }
 241         work[iwk3i + j] = work[iwk3i + j] * work[j] * work[iwk2i + j];
 242         difl[j] = -work[j];
 243         difr[j + difr_dim1] = -work[j + 1];
 244         i__2 = j - 1;
 245         for (i__ = 1; i__ <= i__2; ++i__) {
 246             work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
 247                     i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
 248                     j]);
 249 /* L20: */
 250         }
 251         i__2 = *k;
 252         for (i__ = j + 1; i__ <= i__2; ++i__) {
 253             work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
 254                     i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
 255                     j]);
 256 /* L30: */
 257         }
 258 /* L40: */
 259     }
 260
 261 /*     Compute updated Z. */
 262
 263     i__1 = *k;
 264     for (i__ = 1; i__ <= i__1; ++i__) {
 265         d__2 = sqrt((d__1 = work[iwk3i + i__], abs(d__1)));
 266         z__[i__] = d_sign(&d__2, &z__[i__]);
 267 /* L50: */
 268     }
 269
 270 /*     Update VF and VL. */
 271
 272     i__1 = *k;
 273     for (j = 1; j <= i__1; ++j) {
 274         diflj = difl[j];
 275         dj = d__[j];
 276         dsigj = -dsigma[j];
 277         if (j < *k) {
 278             difrj = -difr[j + difr_dim1];
 279             dsigjp = -dsigma[j + 1];
 280         }
 281         work[j] = -z__[j] / diflj / (dsigma[j] + dj);
 282         i__2 = j - 1;
 283         for (i__ = 1; i__ <= i__2; ++i__) {
 284             work[i__] = z__[i__] / (dlamc3_(&dsigma[i__], &dsigj) - diflj) / (
 285                     dsigma[i__] + dj);
 286 /* L60: */
 287         }
 288         i__2 = *k;
 289         for (i__ = j + 1; i__ <= i__2; ++i__) {
 290             work[i__] = z__[i__] / (dlamc3_(&dsigma[i__], &dsigjp) + difrj) /
 291                     (dsigma[i__] + dj);
 292 /* L70: */
 293         }
 294         temp = dnrm2_(k, &work[1], &c__1);
 295         work[iwk2i + j] = ddot_(k, &work[1], &c__1, &vf[1], &c__1) / temp;
 296         work[iwk3i + j] = ddot_(k, &work[1], &c__1, &vl[1], &c__1) / temp;
 297         if (*icompq == 1) {
 298             difr[j + (difr_dim1 << 1)] = temp;
 299         }
 300 /* L80: */
 301     }
 302
 303     dcopy_(k, &work[iwk2], &c__1, &vf[1], &c__1);
 304     dcopy_(k, &work[iwk3], &c__1, &vl[1], &c__1);
 305
 306     return 0;
 307
 308 /*     End of DLASD8 */
 309
 310 } /* dlasd8_ */