1 SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
2 $ IDXQ, IWORK, WORK, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDU, LDVT, NL, NR, SQRE
11 DOUBLE PRECISION ALPHA, BETA
12 * ..
13 * .. Array Arguments ..
14 INTEGER IDXQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
22 * where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
23 *
24 * A related subroutine DLASD7 handles the case in which the singular
25 * values (and the singular vectors in factored form) are desired.
26 *
27 * DLASD1 computes the SVD as follows:
28 *
29 * ( D1(in) 0 0 0 )
30 * B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
31 * ( 0 0 D2(in) 0 )
32 *
33 * = U(out) * ( D(out) 0) * VT(out)
34 *
35 * where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
36 * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
37 * elsewhere; and the entry b is empty if SQRE = 0.
38 *
39 * The left singular vectors of the original matrix are stored in U, and
40 * the transpose of the right singular vectors are stored in VT, and the
41 * singular values are in D. The algorithm consists of three stages:
42 *
43 * The first stage consists of deflating the size of the problem
44 * when there are multiple singular values or when there are zeros in
45 * the Z vector. For each such occurence the dimension of the
46 * secular equation problem is reduced by one. This stage is
47 * performed by the routine DLASD2.
48 *
49 * The second stage consists of calculating the updated
50 * singular values. This is done by finding the square roots of the
51 * roots of the secular equation via the routine DLASD4 (as called
52 * by DLASD3). This routine also calculates the singular vectors of
53 * the current problem.
54 *
55 * The final stage consists of computing the updated singular vectors
56 * directly using the updated singular values. The singular vectors
57 * for the current problem are multiplied with the singular vectors
58 * from the overall problem.
59 *
60 * Arguments
61 * =========
62 *
63 * NL (input) INTEGER
64 * The row dimension of the upper block. NL >= 1.
65 *
66 * NR (input) INTEGER
67 * The row dimension of the lower block. NR >= 1.
68 *
69 * SQRE (input) INTEGER
70 * = 0: the lower block is an NR-by-NR square matrix.
71 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
72 *
73 * The bidiagonal matrix has row dimension N = NL + NR + 1,
74 * and column dimension M = N + SQRE.
75 *
76 * D (input/output) DOUBLE PRECISION array,
77 * dimension (N = NL+NR+1).
78 * On entry D(1:NL,1:NL) contains the singular values of the
79 * upper block; and D(NL+2:N) contains the singular values of
80 * the lower block. On exit D(1:N) contains the singular values
81 * of the modified matrix.
82 *
83 * ALPHA (input/output) DOUBLE PRECISION
84 * Contains the diagonal element associated with the added row.
85 *
86 * BETA (input/output) DOUBLE PRECISION
87 * Contains the off-diagonal element associated with the added
88 * row.
89 *
90 * U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
91 * On entry U(1:NL, 1:NL) contains the left singular vectors of
92 * the upper block; U(NL+2:N, NL+2:N) contains the left singular
93 * vectors of the lower block. On exit U contains the left
94 * singular vectors of the bidiagonal matrix.
95 *
96 * LDU (input) INTEGER
97 * The leading dimension of the array U. LDU >= max( 1, N ).
98 *
99 * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
100 * where M = N + SQRE.
101 * On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
102 * vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
103 * the right singular vectors of the lower block. On exit
104 * VT**T contains the right singular vectors of the
105 * bidiagonal matrix.
106 *
107 * LDVT (input) INTEGER
108 * The leading dimension of the array VT. LDVT >= max( 1, M ).
109 *
110 * IDXQ (output) INTEGER array, dimension(N)
111 * This contains the permutation which will reintegrate the
112 * subproblem just solved back into sorted order, i.e.
113 * D( IDXQ( I = 1, N ) ) will be in ascending order.
114 *
115 * IWORK (workspace) INTEGER array, dimension( 4 * N )
116 *
117 * WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
118 *
119 * INFO (output) INTEGER
120 * = 0: successful exit.
121 * < 0: if INFO = -i, the i-th argument had an illegal value.
122 * > 0: if INFO = 1, a singular value did not converge
123 *
124 * Further Details
125 * ===============
126 *
127 * Based on contributions by
128 * Ming Gu and Huan Ren, Computer Science Division, University of
129 * California at Berkeley, USA
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134 *
135 DOUBLE PRECISION ONE, ZERO
136 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
137 * ..
138 * .. Local Scalars ..
139 INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
140 $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
141 DOUBLE PRECISION ORGNRM
142 * ..
143 * .. External Subroutines ..
144 EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
145 * ..
146 * .. Intrinsic Functions ..
147 INTRINSIC ABS, MAX
148 * ..
149 * .. Executable Statements ..
150 *
151 * Test the input parameters.
152 *
153 INFO = 0
154 *
155 IF( NL.LT.1 ) THEN
156 INFO = -1
157 ELSE IF( NR.LT.1 ) THEN
158 INFO = -2
159 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
160 INFO = -3
161 END IF
162 IF( INFO.NE.0 ) THEN
163 CALL XERBLA( 'DLASD1', -INFO )
164 RETURN
165 END IF
166 *
167 N = NL + NR + 1
168 M = N + SQRE
169 *
170 * The following values are for bookkeeping purposes only. They are
171 * integer pointers which indicate the portion of the workspace
172 * used by a particular array in DLASD2 and DLASD3.
173 *
174 LDU2 = N
175 LDVT2 = M
176 *
177 IZ = 1
178 ISIGMA = IZ + M
179 IU2 = ISIGMA + N
180 IVT2 = IU2 + LDU2*N
181 IQ = IVT2 + LDVT2*M
182 *
183 IDX = 1
184 IDXC = IDX + N
185 COLTYP = IDXC + N
186 IDXP = COLTYP + N
187 *
188 * Scale.
189 *
190 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
191 D( NL+1 ) = ZERO
192 DO 10 I = 1, N
193 IF( ABS( D( I ) ).GT.ORGNRM ) THEN
194 ORGNRM = ABS( D( I ) )
195 END IF
196 10 CONTINUE
197 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
198 ALPHA = ALPHA / ORGNRM
199 BETA = BETA / ORGNRM
200 *
201 * Deflate singular values.
202 *
203 CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
204 $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
205 $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
206 $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
207 *
208 * Solve Secular Equation and update singular vectors.
209 *
210 LDQ = K
211 CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
212 $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
213 $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
214 $ INFO )
215 IF( INFO.NE.0 ) THEN
216 RETURN
217 END IF
218 *
219 * Unscale.
220 *
221 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
222 *
223 * Prepare the IDXQ sorting permutation.
224 *
225 N1 = K
226 N2 = N - K
227 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
228 *
229 RETURN
230 *
231 * End of DLASD1
232 *
233 END
2 $ IDXQ, IWORK, WORK, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDU, LDVT, NL, NR, SQRE
11 DOUBLE PRECISION ALPHA, BETA
12 * ..
13 * .. Array Arguments ..
14 INTEGER IDXQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
22 * where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
23 *
24 * A related subroutine DLASD7 handles the case in which the singular
25 * values (and the singular vectors in factored form) are desired.
26 *
27 * DLASD1 computes the SVD as follows:
28 *
29 * ( D1(in) 0 0 0 )
30 * B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
31 * ( 0 0 D2(in) 0 )
32 *
33 * = U(out) * ( D(out) 0) * VT(out)
34 *
35 * where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
36 * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
37 * elsewhere; and the entry b is empty if SQRE = 0.
38 *
39 * The left singular vectors of the original matrix are stored in U, and
40 * the transpose of the right singular vectors are stored in VT, and the
41 * singular values are in D. The algorithm consists of three stages:
42 *
43 * The first stage consists of deflating the size of the problem
44 * when there are multiple singular values or when there are zeros in
45 * the Z vector. For each such occurence the dimension of the
46 * secular equation problem is reduced by one. This stage is
47 * performed by the routine DLASD2.
48 *
49 * The second stage consists of calculating the updated
50 * singular values. This is done by finding the square roots of the
51 * roots of the secular equation via the routine DLASD4 (as called
52 * by DLASD3). This routine also calculates the singular vectors of
53 * the current problem.
54 *
55 * The final stage consists of computing the updated singular vectors
56 * directly using the updated singular values. The singular vectors
57 * for the current problem are multiplied with the singular vectors
58 * from the overall problem.
59 *
60 * Arguments
61 * =========
62 *
63 * NL (input) INTEGER
64 * The row dimension of the upper block. NL >= 1.
65 *
66 * NR (input) INTEGER
67 * The row dimension of the lower block. NR >= 1.
68 *
69 * SQRE (input) INTEGER
70 * = 0: the lower block is an NR-by-NR square matrix.
71 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
72 *
73 * The bidiagonal matrix has row dimension N = NL + NR + 1,
74 * and column dimension M = N + SQRE.
75 *
76 * D (input/output) DOUBLE PRECISION array,
77 * dimension (N = NL+NR+1).
78 * On entry D(1:NL,1:NL) contains the singular values of the
79 * upper block; and D(NL+2:N) contains the singular values of
80 * the lower block. On exit D(1:N) contains the singular values
81 * of the modified matrix.
82 *
83 * ALPHA (input/output) DOUBLE PRECISION
84 * Contains the diagonal element associated with the added row.
85 *
86 * BETA (input/output) DOUBLE PRECISION
87 * Contains the off-diagonal element associated with the added
88 * row.
89 *
90 * U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
91 * On entry U(1:NL, 1:NL) contains the left singular vectors of
92 * the upper block; U(NL+2:N, NL+2:N) contains the left singular
93 * vectors of the lower block. On exit U contains the left
94 * singular vectors of the bidiagonal matrix.
95 *
96 * LDU (input) INTEGER
97 * The leading dimension of the array U. LDU >= max( 1, N ).
98 *
99 * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
100 * where M = N + SQRE.
101 * On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
102 * vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
103 * the right singular vectors of the lower block. On exit
104 * VT**T contains the right singular vectors of the
105 * bidiagonal matrix.
106 *
107 * LDVT (input) INTEGER
108 * The leading dimension of the array VT. LDVT >= max( 1, M ).
109 *
110 * IDXQ (output) INTEGER array, dimension(N)
111 * This contains the permutation which will reintegrate the
112 * subproblem just solved back into sorted order, i.e.
113 * D( IDXQ( I = 1, N ) ) will be in ascending order.
114 *
115 * IWORK (workspace) INTEGER array, dimension( 4 * N )
116 *
117 * WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
118 *
119 * INFO (output) INTEGER
120 * = 0: successful exit.
121 * < 0: if INFO = -i, the i-th argument had an illegal value.
122 * > 0: if INFO = 1, a singular value did not converge
123 *
124 * Further Details
125 * ===============
126 *
127 * Based on contributions by
128 * Ming Gu and Huan Ren, Computer Science Division, University of
129 * California at Berkeley, USA
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134 *
135 DOUBLE PRECISION ONE, ZERO
136 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
137 * ..
138 * .. Local Scalars ..
139 INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
140 $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
141 DOUBLE PRECISION ORGNRM
142 * ..
143 * .. External Subroutines ..
144 EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
145 * ..
146 * .. Intrinsic Functions ..
147 INTRINSIC ABS, MAX
148 * ..
149 * .. Executable Statements ..
150 *
151 * Test the input parameters.
152 *
153 INFO = 0
154 *
155 IF( NL.LT.1 ) THEN
156 INFO = -1
157 ELSE IF( NR.LT.1 ) THEN
158 INFO = -2
159 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
160 INFO = -3
161 END IF
162 IF( INFO.NE.0 ) THEN
163 CALL XERBLA( 'DLASD1', -INFO )
164 RETURN
165 END IF
166 *
167 N = NL + NR + 1
168 M = N + SQRE
169 *
170 * The following values are for bookkeeping purposes only. They are
171 * integer pointers which indicate the portion of the workspace
172 * used by a particular array in DLASD2 and DLASD3.
173 *
174 LDU2 = N
175 LDVT2 = M
176 *
177 IZ = 1
178 ISIGMA = IZ + M
179 IU2 = ISIGMA + N
180 IVT2 = IU2 + LDU2*N
181 IQ = IVT2 + LDVT2*M
182 *
183 IDX = 1
184 IDXC = IDX + N
185 COLTYP = IDXC + N
186 IDXP = COLTYP + N
187 *
188 * Scale.
189 *
190 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
191 D( NL+1 ) = ZERO
192 DO 10 I = 1, N
193 IF( ABS( D( I ) ).GT.ORGNRM ) THEN
194 ORGNRM = ABS( D( I ) )
195 END IF
196 10 CONTINUE
197 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
198 ALPHA = ALPHA / ORGNRM
199 BETA = BETA / ORGNRM
200 *
201 * Deflate singular values.
202 *
203 CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
204 $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
205 $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
206 $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
207 *
208 * Solve Secular Equation and update singular vectors.
209 *
210 LDQ = K
211 CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
212 $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
213 $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
214 $ INFO )
215 IF( INFO.NE.0 ) THEN
216 RETURN
217 END IF
218 *
219 * Unscale.
220 *
221 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
222 *
223 * Prepare the IDXQ sorting permutation.
224 *
225 N1 = K
226 N2 = N - K
227 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
228 *
229 RETURN
230 *
231 * End of DLASD1
232 *
233 END