3 /* Table of constant values */
5 static integer c__0 = 0;
6 static integer c__2 = 2;
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 *
13 /* System generated locals */
14 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
16 /* Builtin functions */
17 integer pow_ii(integer *, integer *);
20 integer i__, j, m, i1, ic, lf, nd, ll, nl, nr, im1, ncc, nlf, nrf, iwk,
21 lvl, ndb1, nlp1, nrp1;
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_(
37 /* -- LAPACK auxiliary routine (version 3.1) -- */
38 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
41 /* .. Scalar Arguments .. */
43 /* .. Array Arguments .. */
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. */
55 /* A related subroutine, DLASDA, computes only the singular values, */
56 /* and optionally, the singular vectors in compact form. */
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. */
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; */
70 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
71 /* On entry D contains the main diagonal of the bidiagonal */
73 /* On exit D, if INFO = 0, contains its singular values. */
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. */
79 /* U (output) DOUBLE PRECISION array, dimension at least (LDQ, N) */
80 /* On exit, U contains the left singular vectors. */
82 /* LDU (input) INTEGER */
83 /* On entry, leading dimension of U. */
85 /* VT (output) DOUBLE PRECISION array, dimension at least (LDVT, M) */
86 /* On exit, VT' contains the right singular vectors. */
88 /* LDVT (input) INTEGER */
89 /* On entry, leading dimension of VT. */
91 /* SMLSIZ (input) INTEGER */
92 /* On entry, maximum size of the subproblems at the */
93 /* bottom of the computation tree. */
95 /* IWORK (workspace) INTEGER work array. */
96 /* Dimension must be at least (8 * N) */
98 /* WORK (workspace) DOUBLE PRECISION work array. */
99 /* Dimension must be at least (3 * M**2 + 2 * M) */
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 */
106 /* Further Details */
107 /* =============== */
109 /* Based on contributions by */
110 /* Ming Gu and Huan Ren, Computer Science Division, University of */
111 /* California at Berkeley, USA */
113 /* ===================================================================== */
115 /* .. Local Scalars .. */
117 /* .. External Subroutines .. */
119 /* .. Executable Statements .. */
121 /* Test the input parameters. */
123 /* Parameter adjustments */
127 u_offset = 1 + u_dim1;
130 vt_offset = 1 + vt_dim1;
140 } else if (*sqre < 0 || *sqre > 1) {
148 } else if (*ldvt < m) {
150 } else if (*smlsiz < 3) {
155 xerbla_("DLASD0", &i__1);
159 /* If the input matrix is too small, call DLASDQ to find the SVD. */
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);
167 /* Set up the computation tree. */
174 dlasdt_(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
177 /* For the nodes on bottom level of the tree, solve */
178 /* their subproblems by DLASDQ. */
183 for (i__ = ndb1; i__ <= i__1; ++i__) {
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 */
192 ic = iwork[inode + i1];
193 nl = iwork[ndiml + i1];
195 nr = iwork[ndimr + i1];
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);
206 itemp = idxq + nlf - 2;
208 for (j = 1; j <= i__2; ++j) {
209 iwork[itemp + j] = j;
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);
226 for (j = 1; j <= i__2; ++j) {
227 iwork[itemp + j - 1] = j;
233 /* Now conquer each subproblem bottom-up. */
235 for (lvl = nlvl; lvl >= 1; --lvl) {
237 /* Find the first node LF and last node LL on the */
238 /* current level LVL. */
245 lf = pow_ii(&c__2, &i__1);
249 for (i__ = lf; i__ <= i__1; ++i__) {
251 ic = iwork[inode + im1];
252 nl = iwork[ndiml + im1];
253 nr = iwork[ndimr + im1];
255 if (*sqre == 0 && i__ == ll) {
260 idxqc = idxq + nlf - 1;
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);