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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__0 = 0;
   6 static integer c__2 = 2;
   7
   8 /* Subroutine */ int dlasd0_(integer *n, integer *sqre, doublereal *d__,
   9         doublereal *e, doublereal *u, integer *ldu, doublereal *vt, integer *
  10         ldvt, integer *smlsiz, integer *iwork, doublereal *work, integer *
  11         info)
  12 {
  13     /* System generated locals */
  14     integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  15
  16     /* Builtin functions */
  17     integer pow_ii(integer *, integer *);
  18
  19     /* Local variables */
  20     integer i__, j, m, i1, ic, lf, nd, ll, nl, nr, im1, ncc, nlf, nrf, iwk,
  21             lvl, ndb1, nlp1, nrp1;
  22     doublereal beta;
  23     integer idxq, nlvl;
  24     doublereal alpha;
  25     integer inode, ndiml, idxqc, ndimr, itemp, sqrei;
  26     extern /* Subroutine */ int dlasd1_(integer *, integer *, integer *,
  27             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
  28              doublereal *, integer *, integer *, integer *, doublereal *,
  29             integer *), dlasdq_(char *, integer *, integer *, integer *,
  30             integer *, integer *, doublereal *, doublereal *, doublereal *,
  31             integer *, doublereal *, integer *, doublereal *, integer *,
  32             doublereal *, integer *), dlasdt_(integer *, integer *,
  33             integer *, integer *, integer *, integer *, integer *), xerbla_(
  34             char *, integer *);
  35
  36
  37 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  38 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  39 /*     November 2006 */
  40
  41 /*     .. Scalar Arguments .. */
  42 /*     .. */
  43 /*     .. Array Arguments .. */
  44 /*     .. */
  45
  46 /*  Purpose */
  47 /*  ======= */
  48
  49 /*  Using a divide and conquer approach, DLASD0 computes the singular */
  50 /*  value decomposition (SVD) of a real upper bidiagonal N-by-M */
  51 /*  matrix B with diagonal D and offdiagonal E, where M = N + SQRE. */
  52 /*  The algorithm computes orthogonal matrices U and VT such that */
  53 /*  B = U * S * VT. The singular values S are overwritten on D. */
  54
  55 /*  A related subroutine, DLASDA, computes only the singular values, */
  56 /*  and optionally, the singular vectors in compact form. */
  57
  58 /*  Arguments */
  59 /*  ========= */
  60
  61 /*  N      (input) INTEGER */
  62 /*         On entry, the row dimension of the upper bidiagonal matrix. */
  63 /*         This is also the dimension of the main diagonal array D. */
  64
  65 /*  SQRE   (input) INTEGER */
  66 /*         Specifies the column dimension of the bidiagonal matrix. */
  67 /*         = 0: The bidiagonal matrix has column dimension M = N; */
  68 /*         = 1: The bidiagonal matrix has column dimension M = N+1; */
  69
  70 /*  D      (input/output) DOUBLE PRECISION array, dimension (N) */
  71 /*         On entry D contains the main diagonal of the bidiagonal */
  72 /*         matrix. */
  73 /*         On exit D, if INFO = 0, contains its singular values. */
  74
  75 /*  E      (input) DOUBLE PRECISION array, dimension (M-1) */
  76 /*         Contains the subdiagonal entries of the bidiagonal matrix. */
  77 /*         On exit, E has been destroyed. */
  78
  79 /*  U      (output) DOUBLE PRECISION array, dimension at least (LDQ, N) */
  80 /*         On exit, U contains the left singular vectors. */
  81
  82 /*  LDU    (input) INTEGER */
  83 /*         On entry, leading dimension of U. */
  84
  85 /*  VT     (output) DOUBLE PRECISION array, dimension at least (LDVT, M) */
  86 /*         On exit, VT' contains the right singular vectors. */
  87
  88 /*  LDVT   (input) INTEGER */
  89 /*         On entry, leading dimension of VT. */
  90
  91 /*  SMLSIZ (input) INTEGER */
  92 /*         On entry, maximum size of the subproblems at the */
  93 /*         bottom of the computation tree. */
  94
  95 /*  IWORK  (workspace) INTEGER work array. */
  96 /*         Dimension must be at least (8 * N) */
  97
  98 /*  WORK   (workspace) DOUBLE PRECISION work array. */
  99 /*         Dimension must be at least (3 * M**2 + 2 * M) */
 100
 101 /*  INFO   (output) INTEGER */
 102 /*          = 0:  successful exit. */
 103 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 104 /*          > 0:  if INFO = 1, an singular value did not converge */
 105
 106 /*  Further Details */
 107 /*  =============== */
 108
 109 /*  Based on contributions by */
 110 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 111 /*     California at Berkeley, USA */
 112
 113 /*  ===================================================================== */
 114
 115 /*     .. Local Scalars .. */
 116 /*     .. */
 117 /*     .. External Subroutines .. */
 118 /*     .. */
 119 /*     .. Executable Statements .. */
 120
 121 /*     Test the input parameters. */
 122
 123     /* Parameter adjustments */
 124     --d__;
 125     --e;
 126     u_dim1 = *ldu;
 127     u_offset = 1 + u_dim1;
 128     u -= u_offset;
 129     vt_dim1 = *ldvt;
 130     vt_offset = 1 + vt_dim1;
 131     vt -= vt_offset;
 132     --iwork;
 133     --work;
 134
 135     /* Function Body */
 136     *info = 0;
 137
 138     if (*n < 0) {
 139         *info = -1;
 140     } else if (*sqre < 0 || *sqre > 1) {
 141         *info = -2;
 142     }
 143
 144     m = *n + *sqre;
 145
 146     if (*ldu < *n) {
 147         *info = -6;
 148     } else if (*ldvt < m) {
 149         *info = -8;
 150     } else if (*smlsiz < 3) {
 151         *info = -9;
 152     }
 153     if (*info != 0) {
 154         i__1 = -(*info);
 155         xerbla_("DLASD0", &i__1);
 156         return 0;
 157     }
 158
 159 /*     If the input matrix is too small, call DLASDQ to find the SVD. */
 160
 161     if (*n <= *smlsiz) {
 162         dlasdq_("U", sqre, n, &m, n, &c__0, &d__[1], &e[1], &vt[vt_offset],
 163                 ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[1], info);
 164         return 0;
 165     }
 166
 167 /*     Set up the computation tree. */
 168
 169     inode = 1;
 170     ndiml = inode + *n;
 171     ndimr = ndiml + *n;
 172     idxq = ndimr + *n;
 173     iwk = idxq + *n;
 174     dlasdt_(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
 175             smlsiz);
 176
 177 /*     For the nodes on bottom level of the tree, solve */
 178 /*     their subproblems by DLASDQ. */
 179
 180     ndb1 = (nd + 1) / 2;
 181     ncc = 0;
 182     i__1 = nd;
 183     for (i__ = ndb1; i__ <= i__1; ++i__) {
 184
 185 /*     IC : center row of each node */
 186 /*     NL : number of rows of left  subproblem */
 187 /*     NR : number of rows of right subproblem */
 188 /*     NLF: starting row of the left   subproblem */
 189 /*     NRF: starting row of the right  subproblem */
 190
 191         i1 = i__ - 1;
 192         ic = iwork[inode + i1];
 193         nl = iwork[ndiml + i1];
 194         nlp1 = nl + 1;
 195         nr = iwork[ndimr + i1];
 196         nrp1 = nr + 1;
 197         nlf = ic - nl;
 198         nrf = ic + 1;
 199         sqrei = 1;
 200         dlasdq_("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &vt[
 201                 nlf + nlf * vt_dim1], ldvt, &u[nlf + nlf * u_dim1], ldu, &u[
 202                 nlf + nlf * u_dim1], ldu, &work[1], info);
 203         if (*info != 0) {
 204             return 0;
 205         }
 206         itemp = idxq + nlf - 2;
 207         i__2 = nl;
 208         for (j = 1; j <= i__2; ++j) {
 209             iwork[itemp + j] = j;
 210 /* L10: */
 211         }
 212         if (i__ == nd) {
 213             sqrei = *sqre;
 214         } else {
 215             sqrei = 1;
 216         }
 217         nrp1 = nr + sqrei;
 218         dlasdq_("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &vt[
 219                 nrf + nrf * vt_dim1], ldvt, &u[nrf + nrf * u_dim1], ldu, &u[
 220                 nrf + nrf * u_dim1], ldu, &work[1], info);
 221         if (*info != 0) {
 222             return 0;
 223         }
 224         itemp = idxq + ic;
 225         i__2 = nr;
 226         for (j = 1; j <= i__2; ++j) {
 227             iwork[itemp + j - 1] = j;
 228 /* L20: */
 229         }
 230 /* L30: */
 231     }
 232
 233 /*     Now conquer each subproblem bottom-up. */
 234
 235     for (lvl = nlvl; lvl >= 1; --lvl) {
 236
 237 /*        Find the first node LF and last node LL on the */
 238 /*        current level LVL. */
 239
 240         if (lvl == 1) {
 241             lf = 1;
 242             ll = 1;
 243         } else {
 244             i__1 = lvl - 1;
 245             lf = pow_ii(&c__2, &i__1);
 246             ll = (lf << 1) - 1;
 247         }
 248         i__1 = ll;
 249         for (i__ = lf; i__ <= i__1; ++i__) {
 250             im1 = i__ - 1;
 251             ic = iwork[inode + im1];
 252             nl = iwork[ndiml + im1];
 253             nr = iwork[ndimr + im1];
 254             nlf = ic - nl;
 255             if (*sqre == 0 && i__ == ll) {
 256                 sqrei = *sqre;
 257             } else {
 258                 sqrei = 1;
 259             }
 260             idxqc = idxq + nlf - 1;
 261             alpha = d__[ic];
 262             beta = e[ic];
 263             dlasd1_(&nl, &nr, &sqrei, &d__[nlf], &alpha, &beta, &u[nlf + nlf *
 264                      u_dim1], ldu, &vt[nlf + nlf * vt_dim1], ldvt, &iwork[
 265                     idxqc], &iwork[iwk], &work[1], info);
 266             if (*info != 0) {
 267                 return 0;
 268             }
 269 /* L40: */
 270         }
 271 /* L50: */
 272     }
 273
 274     return 0;
 275
 276 /*     End of DLASD0 */
 277
 278 } /* dlasd0_ */