git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlasd3.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_b13 = 1.;
   8 static doublereal c_b26 = 0.;
   9
  10 /* Subroutine */ int dlasd3_(integer *nl, integer *nr, integer *sqre, integer
  11         *k, doublereal *d__, doublereal *q, integer *ldq, doublereal *dsigma,
  12         doublereal *u, integer *ldu, doublereal *u2, integer *ldu2,
  13         doublereal *vt, integer *ldvt, doublereal *vt2, integer *ldvt2,
  14         integer *idxc, integer *ctot, doublereal *z__, integer *info)
  15 {
  16     /* System generated locals */
  17     integer q_dim1, q_offset, u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1,
  18             vt_offset, vt2_dim1, vt2_offset, i__1, i__2;
  19     doublereal d__1, d__2;
  20
  21     /* Builtin functions */
  22     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
  23
  24     /* Local variables */
  25     integer i__, j, m, n, jc;
  26     doublereal rho;
  27     integer nlp1, nlp2, nrp1;
  28     doublereal temp;
  29     extern doublereal dnrm2_(integer *, doublereal *, integer *);
  30     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
  31             integer *, doublereal *, doublereal *, integer *, doublereal *,
  32             integer *, doublereal *, doublereal *, integer *);
  33     integer ctemp;
  34     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  35             doublereal *, integer *);
  36     integer ktemp;
  37     extern doublereal dlamc3_(doublereal *, doublereal *);
  38     extern /* Subroutine */ int dlasd4_(integer *, integer *, doublereal *,
  39             doublereal *, doublereal *, doublereal *, doublereal *,
  40             doublereal *, integer *), dlascl_(char *, integer *, integer *,
  41             doublereal *, doublereal *, integer *, integer *, doublereal *,
  42             integer *, integer *), dlacpy_(char *, integer *, integer
  43             *, doublereal *, integer *, doublereal *, integer *),
  44             xerbla_(char *, integer *);
  45
  46
  47 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  48 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  49 /*     November 2006 */
  50
  51 /*     .. Scalar Arguments .. */
  52 /*     .. */
  53 /*     .. Array Arguments .. */
  54 /*     .. */
  55
  56 /*  Purpose */
  57 /*  ======= */
  58
  59 /*  DLASD3 finds all the square roots of the roots of the secular */
  60 /*  equation, as defined by the values in D and Z.  It makes the */
  61 /*  appropriate calls to DLASD4 and then updates the singular */
  62 /*  vectors by matrix multiplication. */
  63
  64 /*  This code makes very mild assumptions about floating point */
  65 /*  arithmetic. It will work on machines with a guard digit in */
  66 /*  add/subtract, or on those binary machines without guard digits */
  67 /*  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. */
  68 /*  It could conceivably fail on hexadecimal or decimal machines */
  69 /*  without guard digits, but we know of none. */
  70
  71 /*  DLASD3 is called from DLASD1. */
  72
  73 /*  Arguments */
  74 /*  ========= */
  75
  76 /*  NL     (input) INTEGER */
  77 /*         The row dimension of the upper block.  NL >= 1. */
  78
  79 /*  NR     (input) INTEGER */
  80 /*         The row dimension of the lower block.  NR >= 1. */
  81
  82 /*  SQRE   (input) INTEGER */
  83 /*         = 0: the lower block is an NR-by-NR square matrix. */
  84 /*         = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
  85
  86 /*         The bidiagonal matrix has N = NL + NR + 1 rows and */
  87 /*         M = N + SQRE >= N columns. */
  88
  89 /*  K      (input) INTEGER */
  90 /*         The size of the secular equation, 1 =< K = < N. */
  91
  92 /*  D      (output) DOUBLE PRECISION array, dimension(K) */
  93 /*         On exit the square roots of the roots of the secular equation, */
  94 /*         in ascending order. */
  95
  96 /*  Q      (workspace) DOUBLE PRECISION array, */
  97 /*                     dimension at least (LDQ,K). */
  98
  99 /*  LDQ    (input) INTEGER */
 100 /*         The leading dimension of the array Q.  LDQ >= K. */
 101
 102 /*  DSIGMA (input) DOUBLE PRECISION array, dimension(K) */
 103 /*         The first K elements of this array contain the old roots */
 104 /*         of the deflated updating problem.  These are the poles */
 105 /*         of the secular equation. */
 106
 107 /*  U      (output) DOUBLE PRECISION array, dimension (LDU, N) */
 108 /*         The last N - K columns of this matrix contain the deflated */
 109 /*         left singular vectors. */
 110
 111 /*  LDU    (input) INTEGER */
 112 /*         The leading dimension of the array U.  LDU >= N. */
 113
 114 /*  U2     (input/output) DOUBLE PRECISION array, dimension (LDU2, N) */
 115 /*         The first K columns of this matrix contain the non-deflated */
 116 /*         left singular vectors for the split problem. */
 117
 118 /*  LDU2   (input) INTEGER */
 119 /*         The leading dimension of the array U2.  LDU2 >= N. */
 120
 121 /*  VT     (output) DOUBLE PRECISION array, dimension (LDVT, M) */
 122 /*         The last M - K columns of VT' contain the deflated */
 123 /*         right singular vectors. */
 124
 125 /*  LDVT   (input) INTEGER */
 126 /*         The leading dimension of the array VT.  LDVT >= N. */
 127
 128 /*  VT2    (input/output) DOUBLE PRECISION array, dimension (LDVT2, N) */
 129 /*         The first K columns of VT2' contain the non-deflated */
 130 /*         right singular vectors for the split problem. */
 131
 132 /*  LDVT2  (input) INTEGER */
 133 /*         The leading dimension of the array VT2.  LDVT2 >= N. */
 134
 135 /*  IDXC   (input) INTEGER array, dimension ( N ) */
 136 /*         The permutation used to arrange the columns of U (and rows of */
 137 /*         VT) into three groups:  the first group contains non-zero */
 138 /*         entries only at and above (or before) NL +1; the second */
 139 /*         contains non-zero entries only at and below (or after) NL+2; */
 140 /*         and the third is dense. The first column of U and the row of */
 141 /*         VT are treated separately, however. */
 142
 143 /*         The rows of the singular vectors found by DLASD4 */
 144 /*         must be likewise permuted before the matrix multiplies can */
 145 /*         take place. */
 146
 147 /*  CTOT   (input) INTEGER array, dimension ( 4 ) */
 148 /*         A count of the total number of the various types of columns */
 149 /*         in U (or rows in VT), as described in IDXC. The fourth column */
 150 /*         type is any column which has been deflated. */
 151
 152 /*  Z      (input) DOUBLE PRECISION array, dimension (K) */
 153 /*         The first K elements of this array contain the components */
 154 /*         of the deflation-adjusted updating row vector. */
 155
 156 /*  INFO   (output) INTEGER */
 157 /*         = 0:  successful exit. */
 158 /*         < 0:  if INFO = -i, the i-th argument had an illegal value. */
 159 /*         > 0:  if INFO = 1, an singular value did not converge */
 160
 161 /*  Further Details */
 162 /*  =============== */
 163
 164 /*  Based on contributions by */
 165 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 166 /*     California at Berkeley, USA */
 167
 168 /*  ===================================================================== */
 169
 170 /*     .. Parameters .. */
 171 /*     .. */
 172 /*     .. Local Scalars .. */
 173 /*     .. */
 174 /*     .. External Functions .. */
 175 /*     .. */
 176 /*     .. External Subroutines .. */
 177 /*     .. */
 178 /*     .. Intrinsic Functions .. */
 179 /*     .. */
 180 /*     .. Executable Statements .. */
 181
 182 /*     Test the input parameters. */
 183
 184     /* Parameter adjustments */
 185     --d__;
 186     q_dim1 = *ldq;
 187     q_offset = 1 + q_dim1;
 188     q -= q_offset;
 189     --dsigma;
 190     u_dim1 = *ldu;
 191     u_offset = 1 + u_dim1;
 192     u -= u_offset;
 193     u2_dim1 = *ldu2;
 194     u2_offset = 1 + u2_dim1;
 195     u2 -= u2_offset;
 196     vt_dim1 = *ldvt;
 197     vt_offset = 1 + vt_dim1;
 198     vt -= vt_offset;
 199     vt2_dim1 = *ldvt2;
 200     vt2_offset = 1 + vt2_dim1;
 201     vt2 -= vt2_offset;
 202     --idxc;
 203     --ctot;
 204     --z__;
 205
 206     /* Function Body */
 207     *info = 0;
 208
 209     if (*nl < 1) {
 210         *info = -1;
 211     } else if (*nr < 1) {
 212         *info = -2;
 213     } else if (*sqre != 1 && *sqre != 0) {
 214         *info = -3;
 215     }
 216
 217     n = *nl + *nr + 1;
 218     m = n + *sqre;
 219     nlp1 = *nl + 1;
 220     nlp2 = *nl + 2;
 221
 222     if (*k < 1 || *k > n) {
 223         *info = -4;
 224     } else if (*ldq < *k) {
 225         *info = -7;
 226     } else if (*ldu < n) {
 227         *info = -10;
 228     } else if (*ldu2 < n) {
 229         *info = -12;
 230     } else if (*ldvt < m) {
 231         *info = -14;
 232     } else if (*ldvt2 < m) {
 233         *info = -16;
 234     }
 235     if (*info != 0) {
 236         i__1 = -(*info);
 237         xerbla_("DLASD3", &i__1);
 238         return 0;
 239     }
 240
 241 /*     Quick return if possible */
 242
 243     if (*k == 1) {
 244         d__[1] = abs(z__[1]);
 245         dcopy_(&m, &vt2[vt2_dim1 + 1], ldvt2, &vt[vt_dim1 + 1], ldvt);
 246         if (z__[1] > 0.) {
 247             dcopy_(&n, &u2[u2_dim1 + 1], &c__1, &u[u_dim1 + 1], &c__1);
 248         } else {
 249             i__1 = n;
 250             for (i__ = 1; i__ <= i__1; ++i__) {
 251                 u[i__ + u_dim1] = -u2[i__ + u2_dim1];
 252 /* L10: */
 253             }
 254         }
 255         return 0;
 256     }
 257
 258 /*     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can */
 259 /*     be computed with high relative accuracy (barring over/underflow). */
 260 /*     This is a problem on machines without a guard digit in */
 261 /*     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). */
 262 /*     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), */
 263 /*     which on any of these machines zeros out the bottommost */
 264 /*     bit of DSIGMA(I) if it is 1; this makes the subsequent */
 265 /*     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation */
 266 /*     occurs. On binary machines with a guard digit (almost all */
 267 /*     machines) it does not change DSIGMA(I) at all. On hexadecimal */
 268 /*     and decimal machines with a guard digit, it slightly */
 269 /*     changes the bottommost bits of DSIGMA(I). It does not account */
 270 /*     for hexadecimal or decimal machines without guard digits */
 271 /*     (we know of none). We use a subroutine call to compute */
 272 /*     2*DSIGMA(I) to prevent optimizing compilers from eliminating */
 273 /*     this code. */
 274
 275     i__1 = *k;
 276     for (i__ = 1; i__ <= i__1; ++i__) {
 277         dsigma[i__] = dlamc3_(&dsigma[i__], &dsigma[i__]) - dsigma[i__];
 278 /* L20: */
 279     }
 280
 281 /*     Keep a copy of Z. */
 282
 283     dcopy_(k, &z__[1], &c__1, &q[q_offset], &c__1);
 284
 285 /*     Normalize Z. */
 286
 287     rho = dnrm2_(k, &z__[1], &c__1);
 288     dlascl_("G", &c__0, &c__0, &rho, &c_b13, k, &c__1, &z__[1], k, info);
 289     rho *= rho;
 290
 291 /*     Find the new singular values. */
 292
 293     i__1 = *k;
 294     for (j = 1; j <= i__1; ++j) {
 295         dlasd4_(k, &j, &dsigma[1], &z__[1], &u[j * u_dim1 + 1], &rho, &d__[j],
 296                  &vt[j * vt_dim1 + 1], info);
 297
 298 /*        If the zero finder fails, the computation is terminated. */
 299
 300         if (*info != 0) {
 301             return 0;
 302         }
 303 /* L30: */
 304     }
 305
 306 /*     Compute updated Z. */
 307
 308     i__1 = *k;
 309     for (i__ = 1; i__ <= i__1; ++i__) {
 310         z__[i__] = u[i__ + *k * u_dim1] * vt[i__ + *k * vt_dim1];
 311         i__2 = i__ - 1;
 312         for (j = 1; j <= i__2; ++j) {
 313             z__[i__] *= u[i__ + j * u_dim1] * vt[i__ + j * vt_dim1] / (dsigma[
 314                     i__] - dsigma[j]) / (dsigma[i__] + dsigma[j]);
 315 /* L40: */
 316         }
 317         i__2 = *k - 1;
 318         for (j = i__; j <= i__2; ++j) {
 319             z__[i__] *= u[i__ + j * u_dim1] * vt[i__ + j * vt_dim1] / (dsigma[
 320                     i__] - dsigma[j + 1]) / (dsigma[i__] + dsigma[j + 1]);
 321 /* L50: */
 322         }
 323         d__2 = sqrt((d__1 = z__[i__], abs(d__1)));
 324         z__[i__] = d_sign(&d__2, &q[i__ + q_dim1]);
 325 /* L60: */
 326     }
 327
 328 /*     Compute left singular vectors of the modified diagonal matrix, */
 329 /*     and store related information for the right singular vectors. */
 330
 331     i__1 = *k;
 332     for (i__ = 1; i__ <= i__1; ++i__) {
 333         vt[i__ * vt_dim1 + 1] = z__[1] / u[i__ * u_dim1 + 1] / vt[i__ *
 334                 vt_dim1 + 1];
 335         u[i__ * u_dim1 + 1] = -1.;
 336         i__2 = *k;
 337         for (j = 2; j <= i__2; ++j) {
 338             vt[j + i__ * vt_dim1] = z__[j] / u[j + i__ * u_dim1] / vt[j + i__
 339                     * vt_dim1];
 340             u[j + i__ * u_dim1] = dsigma[j] * vt[j + i__ * vt_dim1];
 341 /* L70: */
 342         }
 343         temp = dnrm2_(k, &u[i__ * u_dim1 + 1], &c__1);
 344         q[i__ * q_dim1 + 1] = u[i__ * u_dim1 + 1] / temp;
 345         i__2 = *k;
 346         for (j = 2; j <= i__2; ++j) {
 347             jc = idxc[j];
 348             q[j + i__ * q_dim1] = u[jc + i__ * u_dim1] / temp;
 349 /* L80: */
 350         }
 351 /* L90: */
 352     }
 353
 354 /*     Update the left singular vector matrix. */
 355
 356     if (*k == 2) {
 357         dgemm_("N", "N", &n, k, k, &c_b13, &u2[u2_offset], ldu2, &q[q_offset],
 358                  ldq, &c_b26, &u[u_offset], ldu);
 359         goto L100;
 360     }
 361     if (ctot[1] > 0) {
 362         dgemm_("N", "N", nl, k, &ctot[1], &c_b13, &u2[(u2_dim1 << 1) + 1],
 363                 ldu2, &q[q_dim1 + 2], ldq, &c_b26, &u[u_dim1 + 1], ldu);
 364         if (ctot[3] > 0) {
 365             ktemp = ctot[1] + 2 + ctot[2];
 366             dgemm_("N", "N", nl, k, &ctot[3], &c_b13, &u2[ktemp * u2_dim1 + 1]
 367 , ldu2, &q[ktemp + q_dim1], ldq, &c_b13, &u[u_dim1 + 1],
 368                     ldu);
 369         }
 370     } else if (ctot[3] > 0) {
 371         ktemp = ctot[1] + 2 + ctot[2];
 372         dgemm_("N", "N", nl, k, &ctot[3], &c_b13, &u2[ktemp * u2_dim1 + 1],
 373                 ldu2, &q[ktemp + q_dim1], ldq, &c_b26, &u[u_dim1 + 1], ldu);
 374     } else {
 375         dlacpy_("F", nl, k, &u2[u2_offset], ldu2, &u[u_offset], ldu);
 376     }
 377     dcopy_(k, &q[q_dim1 + 1], ldq, &u[nlp1 + u_dim1], ldu);
 378     ktemp = ctot[1] + 2;
 379     ctemp = ctot[2] + ctot[3];
 380     dgemm_("N", "N", nr, k, &ctemp, &c_b13, &u2[nlp2 + ktemp * u2_dim1], ldu2,
 381              &q[ktemp + q_dim1], ldq, &c_b26, &u[nlp2 + u_dim1], ldu);
 382
 383 /*     Generate the right singular vectors. */
 384
 385 L100:
 386     i__1 = *k;
 387     for (i__ = 1; i__ <= i__1; ++i__) {
 388         temp = dnrm2_(k, &vt[i__ * vt_dim1 + 1], &c__1);
 389         q[i__ + q_dim1] = vt[i__ * vt_dim1 + 1] / temp;
 390         i__2 = *k;
 391         for (j = 2; j <= i__2; ++j) {
 392             jc = idxc[j];
 393             q[i__ + j * q_dim1] = vt[jc + i__ * vt_dim1] / temp;
 394 /* L110: */
 395         }
 396 /* L120: */
 397     }
 398
 399 /*     Update the right singular vector matrix. */
 400
 401     if (*k == 2) {
 402         dgemm_("N", "N", k, &m, k, &c_b13, &q[q_offset], ldq, &vt2[vt2_offset]
 403 , ldvt2, &c_b26, &vt[vt_offset], ldvt);
 404         return 0;
 405     }
 406     ktemp = ctot[1] + 1;
 407     dgemm_("N", "N", k, &nlp1, &ktemp, &c_b13, &q[q_dim1 + 1], ldq, &vt2[
 408             vt2_dim1 + 1], ldvt2, &c_b26, &vt[vt_dim1 + 1], ldvt);
 409     ktemp = ctot[1] + 2 + ctot[2];
 410     if (ktemp <= *ldvt2) {
 411         dgemm_("N", "N", k, &nlp1, &ctot[3], &c_b13, &q[ktemp * q_dim1 + 1],
 412                 ldq, &vt2[ktemp + vt2_dim1], ldvt2, &c_b13, &vt[vt_dim1 + 1],
 413                 ldvt);
 414     }
 415
 416     ktemp = ctot[1] + 1;
 417     nrp1 = *nr + *sqre;
 418     if (ktemp > 1) {
 419         i__1 = *k;
 420         for (i__ = 1; i__ <= i__1; ++i__) {
 421             q[i__ + ktemp * q_dim1] = q[i__ + q_dim1];
 422 /* L130: */
 423         }
 424         i__1 = m;
 425         for (i__ = nlp2; i__ <= i__1; ++i__) {
 426             vt2[ktemp + i__ * vt2_dim1] = vt2[i__ * vt2_dim1 + 1];
 427 /* L140: */
 428         }
 429     }
 430     ctemp = ctot[2] + 1 + ctot[3];
 431     dgemm_("N", "N", k, &nrp1, &ctemp, &c_b13, &q[ktemp * q_dim1 + 1], ldq, &
 432             vt2[ktemp + nlp2 * vt2_dim1], ldvt2, &c_b26, &vt[nlp2 * vt_dim1 +
 433             1], ldvt);
 434
 435     return 0;
 436
 437 /*     End of DLASD3 */
 438
 439 } /* dlasd3_ */