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