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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static real c_b6 = 0.f;
   7 static integer c__0 = 0;
   8 static real c_b11 = 1.f;
   9
  10 /* Subroutine */ int slalsd_(char *uplo, integer *smlsiz, integer *n, integer
  11         *nrhs, real *d__, real *e, real *b, integer *ldb, real *rcond,
  12         integer *rank, real *work, integer *iwork, integer *info)
  13 {
  14     /* System generated locals */
  15     integer b_dim1, b_offset, i__1, i__2;
  16     real r__1;
  17
  18     /* Builtin functions */
  19     double log(doublereal), r_sign(real *, real *);
  20
  21     /* Local variables */
  22     integer c__, i__, j, k;
  23     real r__;
  24     integer s, u, z__;
  25     real cs;
  26     integer bx;
  27     real sn;
  28     integer st, vt, nm1, st1;
  29     real eps;
  30     integer iwk;
  31     real tol;
  32     integer difl, difr;
  33     real rcnd;
  34     integer perm, nsub, nlvl, sqre, bxst;
  35     extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
  36             integer *, real *, real *), sgemm_(char *, char *, integer *,
  37             integer *, integer *, real *, real *, integer *, real *, integer *
  38 , real *, real *, integer *);
  39     integer poles, sizei, nsize;
  40     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  41             integer *);
  42     integer nwork, icmpq1, icmpq2;
  43     extern doublereal slamch_(char *);
  44     extern /* Subroutine */ int slasda_(integer *, integer *, integer *,
  45             integer *, real *, real *, real *, integer *, real *, integer *,
  46             real *, real *, real *, real *, integer *, integer *, integer *,
  47             integer *, real *, real *, real *, real *, integer *, integer *),
  48             xerbla_(char *, integer *), slalsa_(integer *, integer *,
  49             integer *, integer *, real *, integer *, real *, integer *, real *
  50 , integer *, real *, integer *, real *, real *, real *, real *,
  51             integer *, integer *, integer *, integer *, real *, real *, real *
  52 , real *, integer *, integer *), slascl_(char *, integer *,
  53             integer *, real *, real *, integer *, integer *, real *, integer *
  54 , integer *);
  55     integer givcol;
  56     extern integer isamax_(integer *, real *, integer *);
  57     extern /* Subroutine */ int slasdq_(char *, integer *, integer *, integer
  58             *, integer *, integer *, real *, real *, real *, integer *, real *
  59 , integer *, real *, integer *, real *, integer *),
  60             slacpy_(char *, integer *, integer *, real *, integer *, real *,
  61             integer *), slartg_(real *, real *, real *, real *, real *
  62 ), slaset_(char *, integer *, integer *, real *, real *, real *,
  63             integer *);
  64     real orgnrm;
  65     integer givnum;
  66     extern doublereal slanst_(char *, integer *, real *, real *);
  67     extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);
  68     integer givptr, smlszp;
  69
  70
  71 /*  -- LAPACK routine (version 3.1) -- */
  72 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  73 /*     November 2006 */
  74
  75 /*     .. Scalar Arguments .. */
  76 /*     .. */
  77 /*     .. Array Arguments .. */
  78 /*     .. */
  79
  80 /*  Purpose */
  81 /*  ======= */
  82
  83 /*  SLALSD uses the singular value decomposition of A to solve the least */
  84 /*  squares problem of finding X to minimize the Euclidean norm of each */
  85 /*  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B */
  86 /*  are N-by-NRHS. The solution X overwrites B. */
  87
  88 /*  The singular values of A smaller than RCOND times the largest */
  89 /*  singular value are treated as zero in solving the least squares */
  90 /*  problem; in this case a minimum norm solution is returned. */
  91 /*  The actual singular values are returned in D in ascending order. */
  92
  93 /*  This code makes very mild assumptions about floating point */
  94 /*  arithmetic. It will work on machines with a guard digit in */
  95 /*  add/subtract, or on those binary machines without guard digits */
  96 /*  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. */
  97 /*  It could conceivably fail on hexadecimal or decimal machines */
  98 /*  without guard digits, but we know of none. */
  99
 100 /*  Arguments */
 101 /*  ========= */
 102
 103 /*  UPLO   (input) CHARACTER*1 */
 104 /*         = 'U': D and E define an upper bidiagonal matrix. */
 105 /*         = 'L': D and E define a  lower bidiagonal matrix. */
 106
 107 /*  SMLSIZ (input) INTEGER */
 108 /*         The maximum size of the subproblems at the bottom of the */
 109 /*         computation tree. */
 110
 111 /*  N      (input) INTEGER */
 112 /*         The dimension of the  bidiagonal matrix.  N >= 0. */
 113
 114 /*  NRHS   (input) INTEGER */
 115 /*         The number of columns of B. NRHS must be at least 1. */
 116
 117 /*  D      (input/output) REAL array, dimension (N) */
 118 /*         On entry D contains the main diagonal of the bidiagonal */
 119 /*         matrix. On exit, if INFO = 0, D contains its singular values. */
 120
 121 /*  E      (input/output) REAL array, dimension (N-1) */
 122 /*         Contains the super-diagonal entries of the bidiagonal matrix. */
 123 /*         On exit, E has been destroyed. */
 124
 125 /*  B      (input/output) REAL array, dimension (LDB,NRHS) */
 126 /*         On input, B contains the right hand sides of the least */
 127 /*         squares problem. On output, B contains the solution X. */
 128
 129 /*  LDB    (input) INTEGER */
 130 /*         The leading dimension of B in the calling subprogram. */
 131 /*         LDB must be at least max(1,N). */
 132
 133 /*  RCOND  (input) REAL */
 134 /*         The singular values of A less than or equal to RCOND times */
 135 /*         the largest singular value are treated as zero in solving */
 136 /*         the least squares problem. If RCOND is negative, */
 137 /*         machine precision is used instead. */
 138 /*         For example, if diag(S)*X=B were the least squares problem, */
 139 /*         where diag(S) is a diagonal matrix of singular values, the */
 140 /*         solution would be X(i) = B(i) / S(i) if S(i) is greater than */
 141 /*         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to */
 142 /*         RCOND*max(S). */
 143
 144 /*  RANK   (output) INTEGER */
 145 /*         The number of singular values of A greater than RCOND times */
 146 /*         the largest singular value. */
 147
 148 /*  WORK   (workspace) REAL array, dimension at least */
 149 /*         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2), */
 150 /*         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1). */
 151
 152 /*  IWORK  (workspace) INTEGER array, dimension at least */
 153 /*         (3*N*NLVL + 11*N) */
 154
 155 /*  INFO   (output) INTEGER */
 156 /*         = 0:  successful exit. */
 157 /*         < 0:  if INFO = -i, the i-th argument had an illegal value. */
 158 /*         > 0:  The algorithm failed to compute an singular value while */
 159 /*               working on the submatrix lying in rows and columns */
 160 /*               INFO/(N+1) through MOD(INFO,N+1). */
 161
 162 /*  Further Details */
 163 /*  =============== */
 164
 165 /*  Based on contributions by */
 166 /*     Ming Gu and Ren-Cang Li, Computer Science Division, University of */
 167 /*       California at Berkeley, USA */
 168 /*     Osni Marques, LBNL/NERSC, USA */
 169
 170 /*  ===================================================================== */
 171
 172 /*     .. Parameters .. */
 173 /*     .. */
 174 /*     .. Local Scalars .. */
 175 /*     .. */
 176 /*     .. External Functions .. */
 177 /*     .. */
 178 /*     .. External Subroutines .. */
 179 /*     .. */
 180 /*     .. Intrinsic Functions .. */
 181 /*     .. */
 182 /*     .. Executable Statements .. */
 183
 184 /*     Test the input parameters. */
 185
 186     /* Parameter adjustments */
 187     --d__;
 188     --e;
 189     b_dim1 = *ldb;
 190     b_offset = 1 + b_dim1;
 191     b -= b_offset;
 192     --work;
 193     --iwork;
 194
 195     /* Function Body */
 196     *info = 0;
 197
 198     if (*n < 0) {
 199         *info = -3;
 200     } else if (*nrhs < 1) {
 201         *info = -4;
 202     } else if (*ldb < 1 || *ldb < *n) {
 203         *info = -8;
 204     }
 205     if (*info != 0) {
 206         i__1 = -(*info);
 207         xerbla_("SLALSD", &i__1);
 208         return 0;
 209     }
 210
 211     eps = slamch_("Epsilon");
 212
 213 /*     Set up the tolerance. */
 214
 215     if (*rcond <= 0.f || *rcond >= 1.f) {
 216         rcnd = eps;
 217     } else {
 218         rcnd = *rcond;
 219     }
 220
 221     *rank = 0;
 222
 223 /*     Quick return if possible. */
 224
 225     if (*n == 0) {
 226         return 0;
 227     } else if (*n == 1) {
 228         if (d__[1] == 0.f) {
 229             slaset_("A", &c__1, nrhs, &c_b6, &c_b6, &b[b_offset], ldb);
 230         } else {
 231             *rank = 1;
 232             slascl_("G", &c__0, &c__0, &d__[1], &c_b11, &c__1, nrhs, &b[
 233                     b_offset], ldb, info);
 234             d__[1] = dabs(d__[1]);
 235         }
 236         return 0;
 237     }
 238
 239 /*     Rotate the matrix if it is lower bidiagonal. */
 240
 241     if (*(unsigned char *)uplo == 'L') {
 242         i__1 = *n - 1;
 243         for (i__ = 1; i__ <= i__1; ++i__) {
 244             slartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
 245             d__[i__] = r__;
 246             e[i__] = sn * d__[i__ + 1];
 247             d__[i__ + 1] = cs * d__[i__ + 1];
 248             if (*nrhs == 1) {
 249                 srot_(&c__1, &b[i__ + b_dim1], &c__1, &b[i__ + 1 + b_dim1], &
 250                         c__1, &cs, &sn);
 251             } else {
 252                 work[(i__ << 1) - 1] = cs;
 253                 work[i__ * 2] = sn;
 254             }
 255 /* L10: */
 256         }
 257         if (*nrhs > 1) {
 258             i__1 = *nrhs;
 259             for (i__ = 1; i__ <= i__1; ++i__) {
 260                 i__2 = *n - 1;
 261                 for (j = 1; j <= i__2; ++j) {
 262                     cs = work[(j << 1) - 1];
 263                     sn = work[j * 2];
 264                     srot_(&c__1, &b[j + i__ * b_dim1], &c__1, &b[j + 1 + i__ *
 265                              b_dim1], &c__1, &cs, &sn);
 266 /* L20: */
 267                 }
 268 /* L30: */
 269             }
 270         }
 271     }
 272
 273 /*     Scale. */
 274
 275     nm1 = *n - 1;
 276     orgnrm = slanst_("M", n, &d__[1], &e[1]);
 277     if (orgnrm == 0.f) {
 278         slaset_("A", n, nrhs, &c_b6, &c_b6, &b[b_offset], ldb);
 279         return 0;
 280     }
 281
 282     slascl_("G", &c__0, &c__0, &orgnrm, &c_b11, n, &c__1, &d__[1], n, info);
 283     slascl_("G", &c__0, &c__0, &orgnrm, &c_b11, &nm1, &c__1, &e[1], &nm1,
 284             info);
 285
 286 /*     If N is smaller than the minimum divide size SMLSIZ, then solve */
 287 /*     the problem with another solver. */
 288
 289     if (*n <= *smlsiz) {
 290         nwork = *n * *n + 1;
 291         slaset_("A", n, n, &c_b6, &c_b11, &work[1], n);
 292         slasdq_("U", &c__0, n, n, &c__0, nrhs, &d__[1], &e[1], &work[1], n, &
 293                 work[1], n, &b[b_offset], ldb, &work[nwork], info);
 294         if (*info != 0) {
 295             return 0;
 296         }
 297         tol = rcnd * (r__1 = d__[isamax_(n, &d__[1], &c__1)], dabs(r__1));
 298         i__1 = *n;
 299         for (i__ = 1; i__ <= i__1; ++i__) {
 300             if (d__[i__] <= tol) {
 301                 slaset_("A", &c__1, nrhs, &c_b6, &c_b6, &b[i__ + b_dim1], ldb);
 302             } else {
 303                 slascl_("G", &c__0, &c__0, &d__[i__], &c_b11, &c__1, nrhs, &b[
 304                         i__ + b_dim1], ldb, info);
 305                 ++(*rank);
 306             }
 307 /* L40: */
 308         }
 309         sgemm_("T", "N", n, nrhs, n, &c_b11, &work[1], n, &b[b_offset], ldb, &
 310                 c_b6, &work[nwork], n);
 311         slacpy_("A", n, nrhs, &work[nwork], n, &b[b_offset], ldb);
 312
 313 /*        Unscale. */
 314
 315         slascl_("G", &c__0, &c__0, &c_b11, &orgnrm, n, &c__1, &d__[1], n,
 316                 info);
 317         slasrt_("D", n, &d__[1], info);
 318         slascl_("G", &c__0, &c__0, &orgnrm, &c_b11, n, nrhs, &b[b_offset],
 319                 ldb, info);
 320
 321         return 0;
 322     }
 323
 324 /*     Book-keeping and setting up some constants. */
 325
 326     nlvl = (integer) (log((real) (*n) / (real) (*smlsiz + 1)) / log(2.f)) + 1;
 327
 328     smlszp = *smlsiz + 1;
 329
 330     u = 1;
 331     vt = *smlsiz * *n + 1;
 332     difl = vt + smlszp * *n;
 333     difr = difl + nlvl * *n;
 334     z__ = difr + (nlvl * *n << 1);
 335     c__ = z__ + nlvl * *n;
 336     s = c__ + *n;
 337     poles = s + *n;
 338     givnum = poles + (nlvl << 1) * *n;
 339     bx = givnum + (nlvl << 1) * *n;
 340     nwork = bx + *n * *nrhs;
 341
 342     sizei = *n + 1;
 343     k = sizei + *n;
 344     givptr = k + *n;
 345     perm = givptr + *n;
 346     givcol = perm + nlvl * *n;
 347     iwk = givcol + (nlvl * *n << 1);
 348
 349     st = 1;
 350     sqre = 0;
 351     icmpq1 = 1;
 352     icmpq2 = 0;
 353     nsub = 0;
 354
 355     i__1 = *n;
 356     for (i__ = 1; i__ <= i__1; ++i__) {
 357         if ((r__1 = d__[i__], dabs(r__1)) < eps) {
 358             d__[i__] = r_sign(&eps, &d__[i__]);
 359         }
 360 /* L50: */
 361     }
 362
 363     i__1 = nm1;
 364     for (i__ = 1; i__ <= i__1; ++i__) {
 365         if ((r__1 = e[i__], dabs(r__1)) < eps || i__ == nm1) {
 366             ++nsub;
 367             iwork[nsub] = st;
 368
 369 /*           Subproblem found. First determine its size and then */
 370 /*           apply divide and conquer on it. */
 371
 372             if (i__ < nm1) {
 373
 374 /*              A subproblem with E(I) small for I < NM1. */
 375
 376                 nsize = i__ - st + 1;
 377                 iwork[sizei + nsub - 1] = nsize;
 378             } else if ((r__1 = e[i__], dabs(r__1)) >= eps) {
 379
 380 /*              A subproblem with E(NM1) not too small but I = NM1. */
 381
 382                 nsize = *n - st + 1;
 383                 iwork[sizei + nsub - 1] = nsize;
 384             } else {
 385
 386 /*              A subproblem with E(NM1) small. This implies an */
 387 /*              1-by-1 subproblem at D(N), which is not solved */
 388 /*              explicitly. */
 389
 390                 nsize = i__ - st + 1;
 391                 iwork[sizei + nsub - 1] = nsize;
 392                 ++nsub;
 393                 iwork[nsub] = *n;
 394                 iwork[sizei + nsub - 1] = 1;
 395                 scopy_(nrhs, &b[*n + b_dim1], ldb, &work[bx + nm1], n);
 396             }
 397             st1 = st - 1;
 398             if (nsize == 1) {
 399
 400 /*              This is a 1-by-1 subproblem and is not solved */
 401 /*              explicitly. */
 402
 403                 scopy_(nrhs, &b[st + b_dim1], ldb, &work[bx + st1], n);
 404             } else if (nsize <= *smlsiz) {
 405
 406 /*              This is a small subproblem and is solved by SLASDQ. */
 407
 408                 slaset_("A", &nsize, &nsize, &c_b6, &c_b11, &work[vt + st1],
 409                         n);
 410                 slasdq_("U", &c__0, &nsize, &nsize, &c__0, nrhs, &d__[st], &e[
 411                         st], &work[vt + st1], n, &work[nwork], n, &b[st +
 412                         b_dim1], ldb, &work[nwork], info);
 413                 if (*info != 0) {
 414                     return 0;
 415                 }
 416                 slacpy_("A", &nsize, nrhs, &b[st + b_dim1], ldb, &work[bx +
 417                         st1], n);
 418             } else {
 419
 420 /*              A large problem. Solve it using divide and conquer. */
 421
 422                 slasda_(&icmpq1, smlsiz, &nsize, &sqre, &d__[st], &e[st], &
 423                         work[u + st1], n, &work[vt + st1], &iwork[k + st1], &
 424                         work[difl + st1], &work[difr + st1], &work[z__ + st1],
 425                          &work[poles + st1], &iwork[givptr + st1], &iwork[
 426                         givcol + st1], n, &iwork[perm + st1], &work[givnum +
 427                         st1], &work[c__ + st1], &work[s + st1], &work[nwork],
 428                         &iwork[iwk], info);
 429                 if (*info != 0) {
 430                     return 0;
 431                 }
 432                 bxst = bx + st1;
 433                 slalsa_(&icmpq2, smlsiz, &nsize, nrhs, &b[st + b_dim1], ldb, &
 434                         work[bxst], n, &work[u + st1], n, &work[vt + st1], &
 435                         iwork[k + st1], &work[difl + st1], &work[difr + st1],
 436                         &work[z__ + st1], &work[poles + st1], &iwork[givptr +
 437                         st1], &iwork[givcol + st1], n, &iwork[perm + st1], &
 438                         work[givnum + st1], &work[c__ + st1], &work[s + st1],
 439                         &work[nwork], &iwork[iwk], info);
 440                 if (*info != 0) {
 441                     return 0;
 442                 }
 443             }
 444             st = i__ + 1;
 445         }
 446 /* L60: */
 447     }
 448
 449 /*     Apply the singular values and treat the tiny ones as zero. */
 450
 451     tol = rcnd * (r__1 = d__[isamax_(n, &d__[1], &c__1)], dabs(r__1));
 452
 453     i__1 = *n;
 454     for (i__ = 1; i__ <= i__1; ++i__) {
 455
 456 /*        Some of the elements in D can be negative because 1-by-1 */
 457 /*        subproblems were not solved explicitly. */
 458
 459         if ((r__1 = d__[i__], dabs(r__1)) <= tol) {
 460             slaset_("A", &c__1, nrhs, &c_b6, &c_b6, &work[bx + i__ - 1], n);
 461         } else {
 462             ++(*rank);
 463             slascl_("G", &c__0, &c__0, &d__[i__], &c_b11, &c__1, nrhs, &work[
 464                     bx + i__ - 1], n, info);
 465         }
 466         d__[i__] = (r__1 = d__[i__], dabs(r__1));
 467 /* L70: */
 468     }
 469
 470 /*     Now apply back the right singular vectors. */
 471
 472     icmpq2 = 1;
 473     i__1 = nsub;
 474     for (i__ = 1; i__ <= i__1; ++i__) {
 475         st = iwork[i__];
 476         st1 = st - 1;
 477         nsize = iwork[sizei + i__ - 1];
 478         bxst = bx + st1;
 479         if (nsize == 1) {
 480             scopy_(nrhs, &work[bxst], n, &b[st + b_dim1], ldb);
 481         } else if (nsize <= *smlsiz) {
 482             sgemm_("T", "N", &nsize, nrhs, &nsize, &c_b11, &work[vt + st1], n,
 483                      &work[bxst], n, &c_b6, &b[st + b_dim1], ldb);
 484         } else {
 485             slalsa_(&icmpq2, smlsiz, &nsize, nrhs, &work[bxst], n, &b[st +
 486                     b_dim1], ldb, &work[u + st1], n, &work[vt + st1], &iwork[
 487                     k + st1], &work[difl + st1], &work[difr + st1], &work[z__
 488                     + st1], &work[poles + st1], &iwork[givptr + st1], &iwork[
 489                     givcol + st1], n, &iwork[perm + st1], &work[givnum + st1],
 490                      &work[c__ + st1], &work[s + st1], &work[nwork], &iwork[
 491                     iwk], info);
 492             if (*info != 0) {
 493                 return 0;
 494             }
 495         }
 496 /* L80: */
 497     }
 498
 499 /*     Unscale and sort the singular values. */
 500
 501     slascl_("G", &c__0, &c__0, &c_b11, &orgnrm, n, &c__1, &d__[1], n, info);
 502     slasrt_("D", n, &d__[1], info);
 503     slascl_("G", &c__0, &c__0, &orgnrm, &c_b11, n, nrhs, &b[b_offset], ldb,
 504             info);
 505
 506     return 0;
 507
 508 /*     End of SLALSD */
 509
 510 } /* slalsd_ */