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