3 /* Table of constant values */
5 static integer c__0 = 0;
6 static real c_b7 = 1.f;
7 static integer c__1 = 1;
8 static integer c_n1 = -1;
10 /* Subroutine */ int slasd1_(integer *nl, integer *nr, integer *sqre, real *
11 d__, real *alpha, real *beta, real *u, integer *ldu, real *vt,
12 integer *ldvt, integer *idxq, integer *iwork, real *work, integer *
15 /* System generated locals */
16 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1;
20 integer i__, k, m, n, n1, n2, iq, iz, iu2, ldq, idx, ldu2, ivt2, idxc,
22 extern /* Subroutine */ int slasd2_(integer *, integer *, integer *,
23 integer *, real *, real *, real *, real *, real *, integer *,
24 real *, integer *, real *, real *, integer *, real *, integer *,
25 integer *, integer *, integer *, integer *, integer *, integer *),
26 slasd3_(integer *, integer *, integer *, integer *, real *, real
27 *, integer *, real *, real *, integer *, real *, integer *, real *
28 , integer *, real *, integer *, integer *, integer *, real *,
31 extern /* Subroutine */ int xerbla_(char *, integer *), slascl_(
32 char *, integer *, integer *, real *, real *, integer *, integer *
33 , real *, integer *, integer *), slamrg_(integer *,
34 integer *, real *, integer *, integer *, integer *);
39 /* -- LAPACK auxiliary routine (version 3.1) -- */
40 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
43 /* .. Scalar Arguments .. */
45 /* .. Array Arguments .. */
51 /* SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, */
52 /* where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0. */
54 /* A related subroutine SLASD7 handles the case in which the singular */
55 /* values (and the singular vectors in factored form) are desired. */
57 /* SLASD1 computes the SVD as follows: */
59 /* ( D1(in) 0 0 0 ) */
60 /* B = U(in) * ( Z1' a Z2' b ) * VT(in) */
61 /* ( 0 0 D2(in) 0 ) */
63 /* = U(out) * ( D(out) 0) * VT(out) */
65 /* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M */
66 /* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros */
67 /* elsewhere; and the entry b is empty if SQRE = 0. */
69 /* The left singular vectors of the original matrix are stored in U, and */
70 /* the transpose of the right singular vectors are stored in VT, and the */
71 /* singular values are in D. The algorithm consists of three stages: */
73 /* The first stage consists of deflating the size of the problem */
74 /* when there are multiple singular values or when there are zeros in */
75 /* the Z vector. For each such occurence the dimension of the */
76 /* secular equation problem is reduced by one. This stage is */
77 /* performed by the routine SLASD2. */
79 /* The second stage consists of calculating the updated */
80 /* singular values. This is done by finding the square roots of the */
81 /* roots of the secular equation via the routine SLASD4 (as called */
82 /* by SLASD3). This routine also calculates the singular vectors of */
83 /* the current problem. */
85 /* The final stage consists of computing the updated singular vectors */
86 /* directly using the updated singular values. The singular vectors */
87 /* for the current problem are multiplied with the singular vectors */
88 /* from the overall problem. */
93 /* NL (input) INTEGER */
94 /* The row dimension of the upper block. NL >= 1. */
96 /* NR (input) INTEGER */
97 /* The row dimension of the lower block. NR >= 1. */
99 /* SQRE (input) INTEGER */
100 /* = 0: the lower block is an NR-by-NR square matrix. */
101 /* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
103 /* The bidiagonal matrix has row dimension N = NL + NR + 1, */
104 /* and column dimension M = N + SQRE. */
106 /* D (input/output) REAL array, dimension (NL+NR+1). */
108 /* On entry D(1:NL,1:NL) contains the singular values of the */
109 /* upper block; and D(NL+2:N) contains the singular values of */
110 /* the lower block. On exit D(1:N) contains the singular values */
111 /* of the modified matrix. */
113 /* ALPHA (input/output) REAL */
114 /* Contains the diagonal element associated with the added row. */
116 /* BETA (input/output) REAL */
117 /* Contains the off-diagonal element associated with the added */
120 /* U (input/output) REAL array, dimension (LDU,N) */
121 /* On entry U(1:NL, 1:NL) contains the left singular vectors of */
122 /* the upper block; U(NL+2:N, NL+2:N) contains the left singular */
123 /* vectors of the lower block. On exit U contains the left */
124 /* singular vectors of the bidiagonal matrix. */
126 /* LDU (input) INTEGER */
127 /* The leading dimension of the array U. LDU >= max( 1, N ). */
129 /* VT (input/output) REAL array, dimension (LDVT,M) */
130 /* where M = N + SQRE. */
131 /* On entry VT(1:NL+1, 1:NL+1)' contains the right singular */
132 /* vectors of the upper block; VT(NL+2:M, NL+2:M)' contains */
133 /* the right singular vectors of the lower block. On exit */
134 /* VT' contains the right singular vectors of the */
135 /* bidiagonal matrix. */
137 /* LDVT (input) INTEGER */
138 /* The leading dimension of the array VT. LDVT >= max( 1, M ). */
140 /* IDXQ (output) INTEGER array, dimension (N) */
141 /* This contains the permutation which will reintegrate the */
142 /* subproblem just solved back into sorted order, i.e. */
143 /* D( IDXQ( I = 1, N ) ) will be in ascending order. */
145 /* IWORK (workspace) INTEGER array, dimension (4*N) */
147 /* WORK (workspace) REAL array, dimension (3*M**2+2*M) */
149 /* INFO (output) INTEGER */
150 /* = 0: successful exit. */
151 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
152 /* > 0: if INFO = 1, an singular value did not converge */
154 /* Further Details */
155 /* =============== */
157 /* Based on contributions by */
158 /* Ming Gu and Huan Ren, Computer Science Division, University of */
159 /* California at Berkeley, USA */
161 /* ===================================================================== */
163 /* .. Parameters .. */
166 /* .. Local Scalars .. */
168 /* .. External Subroutines .. */
170 /* .. Intrinsic Functions .. */
172 /* .. Executable Statements .. */
174 /* Test the input parameters. */
176 /* Parameter adjustments */
179 u_offset = 1 + u_dim1;
182 vt_offset = 1 + vt_dim1;
193 } else if (*nr < 1) {
195 } else if (*sqre < 0 || *sqre > 1) {
200 xerbla_("SLASD1", &i__1);
207 /* The following values are for bookkeeping purposes only. They are */
208 /* integer pointers which indicate the portion of the workspace */
209 /* used by a particular array in SLASD2 and SLASD3. */
217 ivt2 = iu2 + ldu2 * n;
218 iq = ivt2 + ldvt2 * m;
228 r__1 = dabs(*alpha), r__2 = dabs(*beta);
229 orgnrm = dmax(r__1,r__2);
232 for (i__ = 1; i__ <= i__1; ++i__) {
233 if ((r__1 = d__[i__], dabs(r__1)) > orgnrm) {
234 orgnrm = (r__1 = d__[i__], dabs(r__1));
238 slascl_("G", &c__0, &c__0, &orgnrm, &c_b7, &n, &c__1, &d__[1], &n, info);
242 /* Deflate singular values. */
244 slasd2_(nl, nr, sqre, &k, &d__[1], &work[iz], alpha, beta, &u[u_offset],
245 ldu, &vt[vt_offset], ldvt, &work[isigma], &work[iu2], &ldu2, &
246 work[ivt2], &ldvt2, &iwork[idxp], &iwork[idx], &iwork[idxc], &
247 idxq[1], &iwork[coltyp], info);
249 /* Solve Secular Equation and update singular vectors. */
252 slasd3_(nl, nr, sqre, &k, &d__[1], &work[iq], &ldq, &work[isigma], &u[
253 u_offset], ldu, &work[iu2], &ldu2, &vt[vt_offset], ldvt, &work[
254 ivt2], &ldvt2, &iwork[idxc], &iwork[coltyp], &work[iz], info);
261 slascl_("G", &c__0, &c__0, &c_b7, &orgnrm, &n, &c__1, &d__[1], &n, info);
263 /* Prepare the IDXQ sorting permutation. */
267 slamrg_(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]);