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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__9 = 9;
   6 static integer c__0 = 0;
   7 static doublereal c_b15 = 1.;
   8 static integer c__1 = 1;
   9 static doublereal c_b29 = 0.;
  10
  11 /* Subroutine */ int dbdsdc_(char *uplo, char *compq, integer *n, doublereal *
  12         d__, doublereal *e, doublereal *u, integer *ldu, doublereal *vt,
  13         integer *ldvt, doublereal *q, integer *iq, doublereal *work, integer *
  14         iwork, integer *info)
  15 {
  16     /* System generated locals */
  17     integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  18     doublereal d__1;
  19
  20     /* Builtin functions */
  21     double d_sign(doublereal *, doublereal *), log(doublereal);
  22
  23     /* Local variables */
  24     integer i__, j, k;
  25     doublereal p, r__;
  26     integer z__, ic, ii, kk;
  27     doublereal cs;
  28     integer is, iu;
  29     doublereal sn;
  30     integer nm1;
  31     doublereal eps;
  32     integer ivt, difl, difr, ierr, perm, mlvl, sqre;
  33     extern logical lsame_(char *, char *);
  34     extern /* Subroutine */ int dlasr_(char *, char *, char *, integer *,
  35             integer *, doublereal *, doublereal *, doublereal *, integer *), dcopy_(integer *, doublereal *, integer *
  36 , doublereal *, integer *), dswap_(integer *, doublereal *,
  37             integer *, doublereal *, integer *);
  38     integer poles, iuplo, nsize, start;
  39     extern /* Subroutine */ int dlasd0_(integer *, integer *, doublereal *,
  40             doublereal *, doublereal *, integer *, doublereal *, integer *,
  41             integer *, integer *, doublereal *, integer *);
  42     extern doublereal dlamch_(char *);
  43     extern /* Subroutine */ int dlasda_(integer *, integer *, integer *,
  44             integer *, doublereal *, doublereal *, doublereal *, integer *,
  45             doublereal *, integer *, doublereal *, doublereal *, doublereal *,
  46              doublereal *, integer *, integer *, integer *, integer *,
  47             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
  48              integer *), dlascl_(char *, integer *, integer *, doublereal *,
  49             doublereal *, integer *, integer *, doublereal *, integer *,
  50             integer *), dlasdq_(char *, integer *, integer *, integer
  51             *, integer *, integer *, doublereal *, doublereal *, doublereal *,
  52              integer *, doublereal *, integer *, doublereal *, integer *,
  53             doublereal *, integer *), dlaset_(char *, integer *,
  54             integer *, doublereal *, doublereal *, doublereal *, integer *), dlartg_(doublereal *, doublereal *, doublereal *,
  55             doublereal *, doublereal *);
  56     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  57             integer *, integer *);
  58     extern /* Subroutine */ int xerbla_(char *, integer *);
  59     integer givcol;
  60     extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);
  61     integer icompq;
  62     doublereal orgnrm;
  63     integer givnum, givptr, qstart, smlsiz, wstart, smlszp;
  64
  65
  66 /*  -- LAPACK routine (version 3.1) -- */
  67 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  68 /*     November 2006 */
  69
  70 /*     .. Scalar Arguments .. */
  71 /*     .. */
  72 /*     .. Array Arguments .. */
  73 /*     .. */
  74
  75 /*  Purpose */
  76 /*  ======= */
  77
  78 /*  DBDSDC computes the singular value decomposition (SVD) of a real */
  79 /*  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT, */
  80 /*  using a divide and conquer method, where S is a diagonal matrix */
  81 /*  with non-negative diagonal elements (the singular values of B), and */
  82 /*  U and VT are orthogonal matrices of left and right singular vectors, */
  83 /*  respectively. DBDSDC can be used to compute all singular values, */
  84 /*  and optionally, singular vectors or singular vectors in compact form. */
  85
  86 /*  This code makes very mild assumptions about floating point */
  87 /*  arithmetic. It will work on machines with a guard digit in */
  88 /*  add/subtract, or on those binary machines without guard digits */
  89 /*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. */
  90 /*  It could conceivably fail on hexadecimal or decimal machines */
  91 /*  without guard digits, but we know of none.  See DLASD3 for details. */
  92
  93 /*  The code currently calls DLASDQ if singular values only are desired. */
  94 /*  However, it can be slightly modified to compute singular values */
  95 /*  using the divide and conquer method. */
  96
  97 /*  Arguments */
  98 /*  ========= */
  99
 100 /*  UPLO    (input) CHARACTER*1 */
 101 /*          = 'U':  B is upper bidiagonal. */
 102 /*          = 'L':  B is lower bidiagonal. */
 103
 104 /*  COMPQ   (input) CHARACTER*1 */
 105 /*          Specifies whether singular vectors are to be computed */
 106 /*          as follows: */
 107 /*          = 'N':  Compute singular values only; */
 108 /*          = 'P':  Compute singular values and compute singular */
 109 /*                  vectors in compact form; */
 110 /*          = 'I':  Compute singular values and singular vectors. */
 111
 112 /*  N       (input) INTEGER */
 113 /*          The order of the matrix B.  N >= 0. */
 114
 115 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
 116 /*          On entry, the n diagonal elements of the bidiagonal matrix B. */
 117 /*          On exit, if INFO=0, the singular values of B. */
 118
 119 /*  E       (input/output) DOUBLE PRECISION array, dimension (N-1) */
 120 /*          On entry, the elements of E contain the offdiagonal */
 121 /*          elements of the bidiagonal matrix whose SVD is desired. */
 122 /*          On exit, E has been destroyed. */
 123
 124 /*  U       (output) DOUBLE PRECISION array, dimension (LDU,N) */
 125 /*          If  COMPQ = 'I', then: */
 126 /*             On exit, if INFO = 0, U contains the left singular vectors */
 127 /*             of the bidiagonal matrix. */
 128 /*          For other values of COMPQ, U is not referenced. */
 129
 130 /*  LDU     (input) INTEGER */
 131 /*          The leading dimension of the array U.  LDU >= 1. */
 132 /*          If singular vectors are desired, then LDU >= max( 1, N ). */
 133
 134 /*  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N) */
 135 /*          If  COMPQ = 'I', then: */
 136 /*             On exit, if INFO = 0, VT' contains the right singular */
 137 /*             vectors of the bidiagonal matrix. */
 138 /*          For other values of COMPQ, VT is not referenced. */
 139
 140 /*  LDVT    (input) INTEGER */
 141 /*          The leading dimension of the array VT.  LDVT >= 1. */
 142 /*          If singular vectors are desired, then LDVT >= max( 1, N ). */
 143
 144 /*  Q       (output) DOUBLE PRECISION array, dimension (LDQ) */
 145 /*          If  COMPQ = 'P', then: */
 146 /*             On exit, if INFO = 0, Q and IQ contain the left */
 147 /*             and right singular vectors in a compact form, */
 148 /*             requiring O(N log N) space instead of 2*N**2. */
 149 /*             In particular, Q contains all the DOUBLE PRECISION data in */
 150 /*             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) */
 151 /*             words of memory, where SMLSIZ is returned by ILAENV and */
 152 /*             is equal to the maximum size of the subproblems at the */
 153 /*             bottom of the computation tree (usually about 25). */
 154 /*          For other values of COMPQ, Q is not referenced. */
 155
 156 /*  IQ      (output) INTEGER array, dimension (LDIQ) */
 157 /*          If  COMPQ = 'P', then: */
 158 /*             On exit, if INFO = 0, Q and IQ contain the left */
 159 /*             and right singular vectors in a compact form, */
 160 /*             requiring O(N log N) space instead of 2*N**2. */
 161 /*             In particular, IQ contains all INTEGER data in */
 162 /*             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) */
 163 /*             words of memory, where SMLSIZ is returned by ILAENV and */
 164 /*             is equal to the maximum size of the subproblems at the */
 165 /*             bottom of the computation tree (usually about 25). */
 166 /*          For other values of COMPQ, IQ is not referenced. */
 167
 168 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
 169 /*          If COMPQ = 'N' then LWORK >= (4 * N). */
 170 /*          If COMPQ = 'P' then LWORK >= (6 * N). */
 171 /*          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). */
 172
 173 /*  IWORK   (workspace) INTEGER array, dimension (8*N) */
 174
 175 /*  INFO    (output) INTEGER */
 176 /*          = 0:  successful exit. */
 177 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 178 /*          > 0:  The algorithm failed to compute an singular value. */
 179 /*                The update process of divide and conquer failed. */
 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 /*  Changed dimension statement in comment describing E from (N) to */
 190 /*  (N-1).  Sven, 17 Feb 05. */
 191 /*  ===================================================================== */
 192
 193 /*     .. Parameters .. */
 194 /*     .. */
 195 /*     .. Local Scalars .. */
 196 /*     .. */
 197 /*     .. External Functions .. */
 198 /*     .. */
 199 /*     .. External Subroutines .. */
 200 /*     .. */
 201 /*     .. Intrinsic Functions .. */
 202 /*     .. */
 203 /*     .. Executable Statements .. */
 204
 205 /*     Test the input parameters. */
 206
 207     /* Parameter adjustments */
 208     --d__;
 209     --e;
 210     u_dim1 = *ldu;
 211     u_offset = 1 + u_dim1;
 212     u -= u_offset;
 213     vt_dim1 = *ldvt;
 214     vt_offset = 1 + vt_dim1;
 215     vt -= vt_offset;
 216     --q;
 217     --iq;
 218     --work;
 219     --iwork;
 220
 221     /* Function Body */
 222     *info = 0;
 223
 224     iuplo = 0;
 225     if (lsame_(uplo, "U")) {
 226         iuplo = 1;
 227     }
 228     if (lsame_(uplo, "L")) {
 229         iuplo = 2;
 230     }
 231     if (lsame_(compq, "N")) {
 232         icompq = 0;
 233     } else if (lsame_(compq, "P")) {
 234         icompq = 1;
 235     } else if (lsame_(compq, "I")) {
 236         icompq = 2;
 237     } else {
 238         icompq = -1;
 239     }
 240     if (iuplo == 0) {
 241         *info = -1;
 242     } else if (icompq < 0) {
 243         *info = -2;
 244     } else if (*n < 0) {
 245         *info = -3;
 246     } else if (*ldu < 1 || icompq == 2 && *ldu < *n) {
 247         *info = -7;
 248     } else if (*ldvt < 1 || icompq == 2 && *ldvt < *n) {
 249         *info = -9;
 250     }
 251     if (*info != 0) {
 252         i__1 = -(*info);
 253         xerbla_("DBDSDC", &i__1);
 254         return 0;
 255     }
 256
 257 /*     Quick return if possible */
 258
 259     if (*n == 0) {
 260         return 0;
 261     }
 262     smlsiz = ilaenv_(&c__9, "DBDSDC", " ", &c__0, &c__0, &c__0, &c__0);
 263     if (*n == 1) {
 264         if (icompq == 1) {
 265             q[1] = d_sign(&c_b15, &d__[1]);
 266             q[smlsiz * *n + 1] = 1.;
 267         } else if (icompq == 2) {
 268             u[u_dim1 + 1] = d_sign(&c_b15, &d__[1]);
 269             vt[vt_dim1 + 1] = 1.;
 270         }
 271         d__[1] = abs(d__[1]);
 272         return 0;
 273     }
 274     nm1 = *n - 1;
 275
 276 /*     If matrix lower bidiagonal, rotate to be upper bidiagonal */
 277 /*     by applying Givens rotations on the left */
 278
 279     wstart = 1;
 280     qstart = 3;
 281     if (icompq == 1) {
 282         dcopy_(n, &d__[1], &c__1, &q[1], &c__1);
 283         i__1 = *n - 1;
 284         dcopy_(&i__1, &e[1], &c__1, &q[*n + 1], &c__1);
 285     }
 286     if (iuplo == 2) {
 287         qstart = 5;
 288         wstart = (*n << 1) - 1;
 289         i__1 = *n - 1;
 290         for (i__ = 1; i__ <= i__1; ++i__) {
 291             dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
 292             d__[i__] = r__;
 293             e[i__] = sn * d__[i__ + 1];
 294             d__[i__ + 1] = cs * d__[i__ + 1];
 295             if (icompq == 1) {
 296                 q[i__ + (*n << 1)] = cs;
 297                 q[i__ + *n * 3] = sn;
 298             } else if (icompq == 2) {
 299                 work[i__] = cs;
 300                 work[nm1 + i__] = -sn;
 301             }
 302 /* L10: */
 303         }
 304     }
 305
 306 /*     If ICOMPQ = 0, use DLASDQ to compute the singular values. */
 307
 308     if (icompq == 0) {
 309         dlasdq_("U", &c__0, n, &c__0, &c__0, &c__0, &d__[1], &e[1], &vt[
 310                 vt_offset], ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[
 311                 wstart], info);
 312         goto L40;
 313     }
 314
 315 /*     If N is smaller than the minimum divide size SMLSIZ, then solve */
 316 /*     the problem with another solver. */
 317
 318     if (*n <= smlsiz) {
 319         if (icompq == 2) {
 320             dlaset_("A", n, n, &c_b29, &c_b15, &u[u_offset], ldu);
 321             dlaset_("A", n, n, &c_b29, &c_b15, &vt[vt_offset], ldvt);
 322             dlasdq_("U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &vt[vt_offset]
 323 , ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[
 324                     wstart], info);
 325         } else if (icompq == 1) {
 326             iu = 1;
 327             ivt = iu + *n;
 328             dlaset_("A", n, n, &c_b29, &c_b15, &q[iu + (qstart - 1) * *n], n);
 329             dlaset_("A", n, n, &c_b29, &c_b15, &q[ivt + (qstart - 1) * *n], n);
 330             dlasdq_("U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &q[ivt + (
 331                     qstart - 1) * *n], n, &q[iu + (qstart - 1) * *n], n, &q[
 332                     iu + (qstart - 1) * *n], n, &work[wstart], info);
 333         }
 334         goto L40;
 335     }
 336
 337     if (icompq == 2) {
 338         dlaset_("A", n, n, &c_b29, &c_b15, &u[u_offset], ldu);
 339         dlaset_("A", n, n, &c_b29, &c_b15, &vt[vt_offset], ldvt);
 340     }
 341
 342 /*     Scale. */
 343
 344     orgnrm = dlanst_("M", n, &d__[1], &e[1]);
 345     if (orgnrm == 0.) {
 346         return 0;
 347     }
 348     dlascl_("G", &c__0, &c__0, &orgnrm, &c_b15, n, &c__1, &d__[1], n, &ierr);
 349     dlascl_("G", &c__0, &c__0, &orgnrm, &c_b15, &nm1, &c__1, &e[1], &nm1, &
 350             ierr);
 351
 352     eps = dlamch_("Epsilon");
 353
 354     mlvl = (integer) (log((doublereal) (*n) / (doublereal) (smlsiz + 1)) /
 355             log(2.)) + 1;
 356     smlszp = smlsiz + 1;
 357
 358     if (icompq == 1) {
 359         iu = 1;
 360         ivt = smlsiz + 1;
 361         difl = ivt + smlszp;
 362         difr = difl + mlvl;
 363         z__ = difr + (mlvl << 1);
 364         ic = z__ + mlvl;
 365         is = ic + 1;
 366         poles = is + 1;
 367         givnum = poles + (mlvl << 1);
 368
 369         k = 1;
 370         givptr = 2;
 371         perm = 3;
 372         givcol = perm + mlvl;
 373     }
 374
 375     i__1 = *n;
 376     for (i__ = 1; i__ <= i__1; ++i__) {
 377         if ((d__1 = d__[i__], abs(d__1)) < eps) {
 378             d__[i__] = d_sign(&eps, &d__[i__]);
 379         }
 380 /* L20: */
 381     }
 382
 383     start = 1;
 384     sqre = 0;
 385
 386     i__1 = nm1;
 387     for (i__ = 1; i__ <= i__1; ++i__) {
 388         if ((d__1 = e[i__], abs(d__1)) < eps || i__ == nm1) {
 389
 390 /*        Subproblem found. First determine its size and then */
 391 /*        apply divide and conquer on it. */
 392
 393             if (i__ < nm1) {
 394
 395 /*        A subproblem with E(I) small for I < NM1. */
 396
 397                 nsize = i__ - start + 1;
 398             } else if ((d__1 = e[i__], abs(d__1)) >= eps) {
 399
 400 /*        A subproblem with E(NM1) not too small but I = NM1. */
 401
 402                 nsize = *n - start + 1;
 403             } else {
 404
 405 /*        A subproblem with E(NM1) small. This implies an */
 406 /*        1-by-1 subproblem at D(N). Solve this 1-by-1 problem */
 407 /*        first. */
 408
 409                 nsize = i__ - start + 1;
 410                 if (icompq == 2) {
 411                     u[*n + *n * u_dim1] = d_sign(&c_b15, &d__[*n]);
 412                     vt[*n + *n * vt_dim1] = 1.;
 413                 } else if (icompq == 1) {
 414                     q[*n + (qstart - 1) * *n] = d_sign(&c_b15, &d__[*n]);
 415                     q[*n + (smlsiz + qstart - 1) * *n] = 1.;
 416                 }
 417                 d__[*n] = (d__1 = d__[*n], abs(d__1));
 418             }
 419             if (icompq == 2) {
 420                 dlasd0_(&nsize, &sqre, &d__[start], &e[start], &u[start +
 421                         start * u_dim1], ldu, &vt[start + start * vt_dim1],
 422                         ldvt, &smlsiz, &iwork[1], &work[wstart], info);
 423             } else {
 424                 dlasda_(&icompq, &smlsiz, &nsize, &sqre, &d__[start], &e[
 425                         start], &q[start + (iu + qstart - 2) * *n], n, &q[
 426                         start + (ivt + qstart - 2) * *n], &iq[start + k * *n],
 427                          &q[start + (difl + qstart - 2) * *n], &q[start + (
 428                         difr + qstart - 2) * *n], &q[start + (z__ + qstart -
 429                         2) * *n], &q[start + (poles + qstart - 2) * *n], &iq[
 430                         start + givptr * *n], &iq[start + givcol * *n], n, &
 431                         iq[start + perm * *n], &q[start + (givnum + qstart -
 432                         2) * *n], &q[start + (ic + qstart - 2) * *n], &q[
 433                         start + (is + qstart - 2) * *n], &work[wstart], &
 434                         iwork[1], info);
 435                 if (*info != 0) {
 436                     return 0;
 437                 }
 438             }
 439             start = i__ + 1;
 440         }
 441 /* L30: */
 442     }
 443
 444 /*     Unscale */
 445
 446     dlascl_("G", &c__0, &c__0, &c_b15, &orgnrm, n, &c__1, &d__[1], n, &ierr);
 447 L40:
 448
 449 /*     Use Selection Sort to minimize swaps of singular vectors */
 450
 451     i__1 = *n;
 452     for (ii = 2; ii <= i__1; ++ii) {
 453         i__ = ii - 1;
 454         kk = i__;
 455         p = d__[i__];
 456         i__2 = *n;
 457         for (j = ii; j <= i__2; ++j) {
 458             if (d__[j] > p) {
 459                 kk = j;
 460                 p = d__[j];
 461             }
 462 /* L50: */
 463         }
 464         if (kk != i__) {
 465             d__[kk] = d__[i__];
 466             d__[i__] = p;
 467             if (icompq == 1) {
 468                 iq[i__] = kk;
 469             } else if (icompq == 2) {
 470                 dswap_(n, &u[i__ * u_dim1 + 1], &c__1, &u[kk * u_dim1 + 1], &
 471                         c__1);
 472                 dswap_(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
 473             }
 474         } else if (icompq == 1) {
 475             iq[i__] = i__;
 476         }
 477 /* L60: */
 478     }
 479
 480 /*     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO */
 481
 482     if (icompq == 1) {
 483         if (iuplo == 1) {
 484             iq[*n] = 1;
 485         } else {
 486             iq[*n] = 0;
 487         }
 488     }
 489
 490 /*     If B is lower bidiagonal, update U by those Givens rotations */
 491 /*     which rotated B to be upper bidiagonal */
 492
 493     if (iuplo == 2 && icompq == 2) {
 494         dlasr_("L", "V", "B", n, n, &work[1], &work[*n], &u[u_offset], ldu);
 495     }
 496
 497     return 0;
 498
 499 /*     End of DBDSDC */
 500
 501 } /* dbdsdc_ */