3 /* Table of constant values */
5 static integer c__0 = 0;
6 static doublereal c_b7 = 1.;
7 static integer c__1 = 1;
8 static integer c_n1 = -1;
10 /* Subroutine */ int dlasd1_(integer *nl, integer *nr, integer *sqre,
11 doublereal *d__, doublereal *alpha, doublereal *beta, doublereal *u,
12 integer *ldu, doublereal *vt, integer *ldvt, integer *idxq, integer *
13 iwork, doublereal *work, integer *info)
15 /* System generated locals */
16 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1;
17 doublereal d__1, d__2;
20 integer i__, k, m, n, n1, n2, iq, iz, iu2, ldq, idx, ldu2, ivt2, idxc,
22 extern /* Subroutine */ int dlasd2_(integer *, integer *, integer *,
23 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
24 doublereal *, integer *, doublereal *, integer *, doublereal *,
25 doublereal *, integer *, doublereal *, integer *, integer *,
26 integer *, integer *, integer *, integer *, integer *), dlasd3_(
27 integer *, integer *, integer *, integer *, doublereal *,
28 doublereal *, integer *, doublereal *, doublereal *, integer *,
29 doublereal *, integer *, doublereal *, integer *, doublereal *,
30 integer *, integer *, integer *, doublereal *, integer *),
31 dlascl_(char *, integer *, integer *, doublereal *, doublereal *,
32 integer *, integer *, doublereal *, integer *, integer *),
33 dlamrg_(integer *, integer *, doublereal *, integer *, integer *,
36 extern /* Subroutine */ int xerbla_(char *, integer *);
41 /* -- LAPACK auxiliary routine (version 3.1) -- */
42 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
45 /* .. Scalar Arguments .. */
47 /* .. Array Arguments .. */
53 /* DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, */
54 /* where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0. */
56 /* A related subroutine DLASD7 handles the case in which the singular */
57 /* values (and the singular vectors in factored form) are desired. */
59 /* DLASD1 computes the SVD as follows: */
61 /* ( D1(in) 0 0 0 ) */
62 /* B = U(in) * ( Z1' a Z2' b ) * VT(in) */
63 /* ( 0 0 D2(in) 0 ) */
65 /* = U(out) * ( D(out) 0) * VT(out) */
67 /* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M */
68 /* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros */
69 /* elsewhere; and the entry b is empty if SQRE = 0. */
71 /* The left singular vectors of the original matrix are stored in U, and */
72 /* the transpose of the right singular vectors are stored in VT, and the */
73 /* singular values are in D. The algorithm consists of three stages: */
75 /* The first stage consists of deflating the size of the problem */
76 /* when there are multiple singular values or when there are zeros in */
77 /* the Z vector. For each such occurence the dimension of the */
78 /* secular equation problem is reduced by one. This stage is */
79 /* performed by the routine DLASD2. */
81 /* The second stage consists of calculating the updated */
82 /* singular values. This is done by finding the square roots of the */
83 /* roots of the secular equation via the routine DLASD4 (as called */
84 /* by DLASD3). This routine also calculates the singular vectors of */
85 /* the current problem. */
87 /* The final stage consists of computing the updated singular vectors */
88 /* directly using the updated singular values. The singular vectors */
89 /* for the current problem are multiplied with the singular vectors */
90 /* from the overall problem. */
95 /* NL (input) INTEGER */
96 /* The row dimension of the upper block. NL >= 1. */
98 /* NR (input) INTEGER */
99 /* The row dimension of the lower block. NR >= 1. */
101 /* SQRE (input) INTEGER */
102 /* = 0: the lower block is an NR-by-NR square matrix. */
103 /* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
105 /* The bidiagonal matrix has row dimension N = NL + NR + 1, */
106 /* and column dimension M = N + SQRE. */
108 /* D (input/output) DOUBLE PRECISION array, */
109 /* dimension (N = NL+NR+1). */
110 /* On entry D(1:NL,1:NL) contains the singular values of the */
111 /* upper block; and D(NL+2:N) contains the singular values of */
112 /* the lower block. On exit D(1:N) contains the singular values */
113 /* of the modified matrix. */
115 /* ALPHA (input/output) DOUBLE PRECISION */
116 /* Contains the diagonal element associated with the added row. */
118 /* BETA (input/output) DOUBLE PRECISION */
119 /* Contains the off-diagonal element associated with the added */
122 /* U (input/output) DOUBLE PRECISION array, dimension(LDU,N) */
123 /* On entry U(1:NL, 1:NL) contains the left singular vectors of */
124 /* the upper block; U(NL+2:N, NL+2:N) contains the left singular */
125 /* vectors of the lower block. On exit U contains the left */
126 /* singular vectors of the bidiagonal matrix. */
128 /* LDU (input) INTEGER */
129 /* The leading dimension of the array U. LDU >= max( 1, N ). */
131 /* VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M) */
132 /* where M = N + SQRE. */
133 /* On entry VT(1:NL+1, 1:NL+1)' contains the right singular */
134 /* vectors of the upper block; VT(NL+2:M, NL+2:M)' contains */
135 /* the right singular vectors of the lower block. On exit */
136 /* VT' contains the right singular vectors of the */
137 /* bidiagonal matrix. */
139 /* LDVT (input) INTEGER */
140 /* The leading dimension of the array VT. LDVT >= max( 1, M ). */
142 /* IDXQ (output) INTEGER array, dimension(N) */
143 /* This contains the permutation which will reintegrate the */
144 /* subproblem just solved back into sorted order, i.e. */
145 /* D( IDXQ( I = 1, N ) ) will be in ascending order. */
147 /* IWORK (workspace) INTEGER array, dimension( 4 * N ) */
149 /* WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M ) */
151 /* INFO (output) INTEGER */
152 /* = 0: successful exit. */
153 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
154 /* > 0: if INFO = 1, an singular value did not converge */
156 /* Further Details */
157 /* =============== */
159 /* Based on contributions by */
160 /* Ming Gu and Huan Ren, Computer Science Division, University of */
161 /* California at Berkeley, USA */
163 /* ===================================================================== */
165 /* .. Parameters .. */
168 /* .. Local Scalars .. */
170 /* .. External Subroutines .. */
172 /* .. Intrinsic Functions .. */
174 /* .. Executable Statements .. */
176 /* Test the input parameters. */
178 /* Parameter adjustments */
181 u_offset = 1 + u_dim1;
184 vt_offset = 1 + vt_dim1;
195 } else if (*nr < 1) {
197 } else if (*sqre < 0 || *sqre > 1) {
202 xerbla_("DLASD1", &i__1);
209 /* The following values are for bookkeeping purposes only. They are */
210 /* integer pointers which indicate the portion of the workspace */
211 /* used by a particular array in DLASD2 and DLASD3. */
219 ivt2 = iu2 + ldu2 * n;
220 iq = ivt2 + ldvt2 * m;
230 d__1 = abs(*alpha), d__2 = abs(*beta);
231 orgnrm = max(d__1,d__2);
234 for (i__ = 1; i__ <= i__1; ++i__) {
235 if ((d__1 = d__[i__], abs(d__1)) > orgnrm) {
236 orgnrm = (d__1 = d__[i__], abs(d__1));
240 dlascl_("G", &c__0, &c__0, &orgnrm, &c_b7, &n, &c__1, &d__[1], &n, info);
244 /* Deflate singular values. */
246 dlasd2_(nl, nr, sqre, &k, &d__[1], &work[iz], alpha, beta, &u[u_offset],
247 ldu, &vt[vt_offset], ldvt, &work[isigma], &work[iu2], &ldu2, &
248 work[ivt2], &ldvt2, &iwork[idxp], &iwork[idx], &iwork[idxc], &
249 idxq[1], &iwork[coltyp], info);
251 /* Solve Secular Equation and update singular vectors. */
254 dlasd3_(nl, nr, sqre, &k, &d__[1], &work[iq], &ldq, &work[isigma], &u[
255 u_offset], ldu, &work[iu2], &ldu2, &vt[vt_offset], ldvt, &work[
256 ivt2], &ldvt2, &iwork[idxc], &iwork[coltyp], &work[iz], info);
263 dlascl_("G", &c__0, &c__0, &c_b7, &orgnrm, &n, &c__1, &d__[1], &n, info);
265 /* Prepare the IDXQ sorting permutation. */
269 dlamrg_(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]);