3 /* Table of constant values */
5 static integer c__0 = 0;
6 static integer c__2 = 2;
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)
12 /* System generated locals */
13 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
15 /* Builtin functions */
16 integer pow_ii(integer *, integer *);
19 integer i__, j, m, i1, ic, lf, nd, ll, nl, nr, im1, ncc, nlf, nrf, iwk,
20 lvl, ndb1, nlp1, nrp1;
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 *
34 /* -- LAPACK auxiliary routine (version 3.1) -- */
35 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
38 /* .. Scalar Arguments .. */
40 /* .. Array Arguments .. */
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. */
52 /* A related subroutine, SLASDA, computes only the singular values, */
53 /* and optionally, the singular vectors in compact form. */
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. */
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; */
67 /* D (input/output) REAL array, dimension (N) */
68 /* On entry D contains the main diagonal of the bidiagonal */
70 /* On exit D, if INFO = 0, contains its singular values. */
72 /* E (input) REAL array, dimension (M-1) */
73 /* Contains the subdiagonal entries of the bidiagonal matrix. */
74 /* On exit, E has been destroyed. */
76 /* U (output) REAL array, dimension at least (LDQ, N) */
77 /* On exit, U contains the left singular vectors. */
79 /* LDU (input) INTEGER */
80 /* On entry, leading dimension of U. */
82 /* VT (output) REAL array, dimension at least (LDVT, M) */
83 /* On exit, VT' contains the right singular vectors. */
85 /* LDVT (input) INTEGER */
86 /* On entry, leading dimension of VT. */
88 /* SMLSIZ (input) INTEGER */
89 /* On entry, maximum size of the subproblems at the */
90 /* bottom of the computation tree. */
92 /* IWORK (workspace) INTEGER array, dimension (8*N) */
94 /* WORK (workspace) REAL array, dimension (3*M**2+2*M) */
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 */
101 /* Further Details */
102 /* =============== */
104 /* Based on contributions by */
105 /* Ming Gu and Huan Ren, Computer Science Division, University of */
106 /* California at Berkeley, USA */
108 /* ===================================================================== */
110 /* .. Local Scalars .. */
112 /* .. External Subroutines .. */
114 /* .. Executable Statements .. */
116 /* Test the input parameters. */
118 /* Parameter adjustments */
122 u_offset = 1 + u_dim1;
125 vt_offset = 1 + vt_dim1;
135 } else if (*sqre < 0 || *sqre > 1) {
143 } else if (*ldvt < m) {
145 } else if (*smlsiz < 3) {
150 xerbla_("SLASD0", &i__1);
154 /* If the input matrix is too small, call SLASDQ to find the SVD. */
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);
162 /* Set up the computation tree. */
169 slasdt_(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
172 /* For the nodes on bottom level of the tree, solve */
173 /* their subproblems by SLASDQ. */
178 for (i__ = ndb1; i__ <= i__1; ++i__) {
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 */
187 ic = iwork[inode + i1];
188 nl = iwork[ndiml + i1];
190 nr = iwork[ndimr + i1];
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);
201 itemp = idxq + nlf - 2;
203 for (j = 1; j <= i__2; ++j) {
204 iwork[itemp + j] = j;
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);
221 for (j = 1; j <= i__2; ++j) {
222 iwork[itemp + j - 1] = j;
228 /* Now conquer each subproblem bottom-up. */
230 for (lvl = nlvl; lvl >= 1; --lvl) {
232 /* Find the first node LF and last node LL on the */
233 /* current level LVL. */
240 lf = pow_ii(&c__2, &i__1);
244 for (i__ = lf; i__ <= i__1; ++i__) {
246 ic = iwork[inode + im1];
247 nl = iwork[ndiml + im1];
248 nr = iwork[ndimr + im1];
250 if (*sqre == 0 && i__ == ll) {
255 idxqc = idxq + nlf - 1;
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);