1 SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
2 $ WORK, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.2.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * June 2010
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * Using a divide and conquer approach, DLASD0 computes the singular
22 * value decomposition (SVD) of a real upper bidiagonal N-by-M
23 * matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
24 * The algorithm computes orthogonal matrices U and VT such that
25 * B = U * S * VT. The singular values S are overwritten on D.
26 *
27 * A related subroutine, DLASDA, computes only the singular values,
28 * and optionally, the singular vectors in compact form.
29 *
30 * Arguments
31 * =========
32 *
33 * N (input) INTEGER
34 * On entry, the row dimension of the upper bidiagonal matrix.
35 * This is also the dimension of the main diagonal array D.
36 *
37 * SQRE (input) INTEGER
38 * Specifies the column dimension of the bidiagonal matrix.
39 * = 0: The bidiagonal matrix has column dimension M = N;
40 * = 1: The bidiagonal matrix has column dimension M = N+1;
41 *
42 * D (input/output) DOUBLE PRECISION array, dimension (N)
43 * On entry D contains the main diagonal of the bidiagonal
44 * matrix.
45 * On exit D, if INFO = 0, contains its singular values.
46 *
47 * E (input) DOUBLE PRECISION array, dimension (M-1)
48 * Contains the subdiagonal entries of the bidiagonal matrix.
49 * On exit, E has been destroyed.
50 *
51 * U (output) DOUBLE PRECISION array, dimension at least (LDQ, N)
52 * On exit, U contains the left singular vectors.
53 *
54 * LDU (input) INTEGER
55 * On entry, leading dimension of U.
56 *
57 * VT (output) DOUBLE PRECISION array, dimension at least (LDVT, M)
58 * On exit, VT**T contains the right singular vectors.
59 *
60 * LDVT (input) INTEGER
61 * On entry, leading dimension of VT.
62 *
63 * SMLSIZ (input) INTEGER
64 * On entry, maximum size of the subproblems at the
65 * bottom of the computation tree.
66 *
67 * IWORK (workspace) INTEGER work array.
68 * Dimension must be at least (8 * N)
69 *
70 * WORK (workspace) DOUBLE PRECISION work array.
71 * Dimension must be at least (3 * M**2 + 2 * M)
72 *
73 * INFO (output) INTEGER
74 * = 0: successful exit.
75 * < 0: if INFO = -i, the i-th argument had an illegal value.
76 * > 0: if INFO = 1, a singular value did not converge
77 *
78 * Further Details
79 * ===============
80 *
81 * Based on contributions by
82 * Ming Gu and Huan Ren, Computer Science Division, University of
83 * California at Berkeley, USA
84 *
85 * =====================================================================
86 *
87 * .. Local Scalars ..
88 INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK,
89 $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR,
90 $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI
91 DOUBLE PRECISION ALPHA, BETA
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL DLASD1, DLASDQ, DLASDT, XERBLA
95 * ..
96 * .. Executable Statements ..
97 *
98 * Test the input parameters.
99 *
100 INFO = 0
101 *
102 IF( N.LT.0 ) THEN
103 INFO = -1
104 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
105 INFO = -2
106 END IF
107 *
108 M = N + SQRE
109 *
110 IF( LDU.LT.N ) THEN
111 INFO = -6
112 ELSE IF( LDVT.LT.M ) THEN
113 INFO = -8
114 ELSE IF( SMLSIZ.LT.3 ) THEN
115 INFO = -9
116 END IF
117 IF( INFO.NE.0 ) THEN
118 CALL XERBLA( 'DLASD0', -INFO )
119 RETURN
120 END IF
121 *
122 * If the input matrix is too small, call DLASDQ to find the SVD.
123 *
124 IF( N.LE.SMLSIZ ) THEN
125 CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U,
126 $ LDU, WORK, INFO )
127 RETURN
128 END IF
129 *
130 * Set up the computation tree.
131 *
132 INODE = 1
133 NDIML = INODE + N
134 NDIMR = NDIML + N
135 IDXQ = NDIMR + N
136 IWK = IDXQ + N
137 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
138 $ IWORK( NDIMR ), SMLSIZ )
139 *
140 * For the nodes on bottom level of the tree, solve
141 * their subproblems by DLASDQ.
142 *
143 NDB1 = ( ND+1 ) / 2
144 NCC = 0
145 DO 30 I = NDB1, ND
146 *
147 * IC : center row of each node
148 * NL : number of rows of left subproblem
149 * NR : number of rows of right subproblem
150 * NLF: starting row of the left subproblem
151 * NRF: starting row of the right subproblem
152 *
153 I1 = I - 1
154 IC = IWORK( INODE+I1 )
155 NL = IWORK( NDIML+I1 )
156 NLP1 = NL + 1
157 NR = IWORK( NDIMR+I1 )
158 NRP1 = NR + 1
159 NLF = IC - NL
160 NRF = IC + 1
161 SQREI = 1
162 CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ),
163 $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU,
164 $ U( NLF, NLF ), LDU, WORK, INFO )
165 IF( INFO.NE.0 ) THEN
166 RETURN
167 END IF
168 ITEMP = IDXQ + NLF - 2
169 DO 10 J = 1, NL
170 IWORK( ITEMP+J ) = J
171 10 CONTINUE
172 IF( I.EQ.ND ) THEN
173 SQREI = SQRE
174 ELSE
175 SQREI = 1
176 END IF
177 NRP1 = NR + SQREI
178 CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ),
179 $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU,
180 $ U( NRF, NRF ), LDU, WORK, INFO )
181 IF( INFO.NE.0 ) THEN
182 RETURN
183 END IF
184 ITEMP = IDXQ + IC
185 DO 20 J = 1, NR
186 IWORK( ITEMP+J-1 ) = J
187 20 CONTINUE
188 30 CONTINUE
189 *
190 * Now conquer each subproblem bottom-up.
191 *
192 DO 50 LVL = NLVL, 1, -1
193 *
194 * Find the first node LF and last node LL on the
195 * current level LVL.
196 *
197 IF( LVL.EQ.1 ) THEN
198 LF = 1
199 LL = 1
200 ELSE
201 LF = 2**( LVL-1 )
202 LL = 2*LF - 1
203 END IF
204 DO 40 I = LF, LL
205 IM1 = I - 1
206 IC = IWORK( INODE+IM1 )
207 NL = IWORK( NDIML+IM1 )
208 NR = IWORK( NDIMR+IM1 )
209 NLF = IC - NL
210 IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN
211 SQREI = SQRE
212 ELSE
213 SQREI = 1
214 END IF
215 IDXQC = IDXQ + NLF - 1
216 ALPHA = D( IC )
217 BETA = E( IC )
218 CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA,
219 $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT,
220 $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO )
221 IF( INFO.NE.0 ) THEN
222 RETURN
223 END IF
224 40 CONTINUE
225 50 CONTINUE
226 *
227 RETURN
228 *
229 * End of DLASD0
230 *
231 END
2 $ WORK, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.2.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * June 2010
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * Using a divide and conquer approach, DLASD0 computes the singular
22 * value decomposition (SVD) of a real upper bidiagonal N-by-M
23 * matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
24 * The algorithm computes orthogonal matrices U and VT such that
25 * B = U * S * VT. The singular values S are overwritten on D.
26 *
27 * A related subroutine, DLASDA, computes only the singular values,
28 * and optionally, the singular vectors in compact form.
29 *
30 * Arguments
31 * =========
32 *
33 * N (input) INTEGER
34 * On entry, the row dimension of the upper bidiagonal matrix.
35 * This is also the dimension of the main diagonal array D.
36 *
37 * SQRE (input) INTEGER
38 * Specifies the column dimension of the bidiagonal matrix.
39 * = 0: The bidiagonal matrix has column dimension M = N;
40 * = 1: The bidiagonal matrix has column dimension M = N+1;
41 *
42 * D (input/output) DOUBLE PRECISION array, dimension (N)
43 * On entry D contains the main diagonal of the bidiagonal
44 * matrix.
45 * On exit D, if INFO = 0, contains its singular values.
46 *
47 * E (input) DOUBLE PRECISION array, dimension (M-1)
48 * Contains the subdiagonal entries of the bidiagonal matrix.
49 * On exit, E has been destroyed.
50 *
51 * U (output) DOUBLE PRECISION array, dimension at least (LDQ, N)
52 * On exit, U contains the left singular vectors.
53 *
54 * LDU (input) INTEGER
55 * On entry, leading dimension of U.
56 *
57 * VT (output) DOUBLE PRECISION array, dimension at least (LDVT, M)
58 * On exit, VT**T contains the right singular vectors.
59 *
60 * LDVT (input) INTEGER
61 * On entry, leading dimension of VT.
62 *
63 * SMLSIZ (input) INTEGER
64 * On entry, maximum size of the subproblems at the
65 * bottom of the computation tree.
66 *
67 * IWORK (workspace) INTEGER work array.
68 * Dimension must be at least (8 * N)
69 *
70 * WORK (workspace) DOUBLE PRECISION work array.
71 * Dimension must be at least (3 * M**2 + 2 * M)
72 *
73 * INFO (output) INTEGER
74 * = 0: successful exit.
75 * < 0: if INFO = -i, the i-th argument had an illegal value.
76 * > 0: if INFO = 1, a singular value did not converge
77 *
78 * Further Details
79 * ===============
80 *
81 * Based on contributions by
82 * Ming Gu and Huan Ren, Computer Science Division, University of
83 * California at Berkeley, USA
84 *
85 * =====================================================================
86 *
87 * .. Local Scalars ..
88 INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK,
89 $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR,
90 $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI
91 DOUBLE PRECISION ALPHA, BETA
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL DLASD1, DLASDQ, DLASDT, XERBLA
95 * ..
96 * .. Executable Statements ..
97 *
98 * Test the input parameters.
99 *
100 INFO = 0
101 *
102 IF( N.LT.0 ) THEN
103 INFO = -1
104 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
105 INFO = -2
106 END IF
107 *
108 M = N + SQRE
109 *
110 IF( LDU.LT.N ) THEN
111 INFO = -6
112 ELSE IF( LDVT.LT.M ) THEN
113 INFO = -8
114 ELSE IF( SMLSIZ.LT.3 ) THEN
115 INFO = -9
116 END IF
117 IF( INFO.NE.0 ) THEN
118 CALL XERBLA( 'DLASD0', -INFO )
119 RETURN
120 END IF
121 *
122 * If the input matrix is too small, call DLASDQ to find the SVD.
123 *
124 IF( N.LE.SMLSIZ ) THEN
125 CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U,
126 $ LDU, WORK, INFO )
127 RETURN
128 END IF
129 *
130 * Set up the computation tree.
131 *
132 INODE = 1
133 NDIML = INODE + N
134 NDIMR = NDIML + N
135 IDXQ = NDIMR + N
136 IWK = IDXQ + N
137 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
138 $ IWORK( NDIMR ), SMLSIZ )
139 *
140 * For the nodes on bottom level of the tree, solve
141 * their subproblems by DLASDQ.
142 *
143 NDB1 = ( ND+1 ) / 2
144 NCC = 0
145 DO 30 I = NDB1, ND
146 *
147 * IC : center row of each node
148 * NL : number of rows of left subproblem
149 * NR : number of rows of right subproblem
150 * NLF: starting row of the left subproblem
151 * NRF: starting row of the right subproblem
152 *
153 I1 = I - 1
154 IC = IWORK( INODE+I1 )
155 NL = IWORK( NDIML+I1 )
156 NLP1 = NL + 1
157 NR = IWORK( NDIMR+I1 )
158 NRP1 = NR + 1
159 NLF = IC - NL
160 NRF = IC + 1
161 SQREI = 1
162 CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ),
163 $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU,
164 $ U( NLF, NLF ), LDU, WORK, INFO )
165 IF( INFO.NE.0 ) THEN
166 RETURN
167 END IF
168 ITEMP = IDXQ + NLF - 2
169 DO 10 J = 1, NL
170 IWORK( ITEMP+J ) = J
171 10 CONTINUE
172 IF( I.EQ.ND ) THEN
173 SQREI = SQRE
174 ELSE
175 SQREI = 1
176 END IF
177 NRP1 = NR + SQREI
178 CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ),
179 $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU,
180 $ U( NRF, NRF ), LDU, WORK, INFO )
181 IF( INFO.NE.0 ) THEN
182 RETURN
183 END IF
184 ITEMP = IDXQ + IC
185 DO 20 J = 1, NR
186 IWORK( ITEMP+J-1 ) = J
187 20 CONTINUE
188 30 CONTINUE
189 *
190 * Now conquer each subproblem bottom-up.
191 *
192 DO 50 LVL = NLVL, 1, -1
193 *
194 * Find the first node LF and last node LL on the
195 * current level LVL.
196 *
197 IF( LVL.EQ.1 ) THEN
198 LF = 1
199 LL = 1
200 ELSE
201 LF = 2**( LVL-1 )
202 LL = 2*LF - 1
203 END IF
204 DO 40 I = LF, LL
205 IM1 = I - 1
206 IC = IWORK( INODE+IM1 )
207 NL = IWORK( NDIML+IM1 )
208 NR = IWORK( NDIMR+IM1 )
209 NLF = IC - NL
210 IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN
211 SQREI = SQRE
212 ELSE
213 SQREI = 1
214 END IF
215 IDXQC = IDXQ + NLF - 1
216 ALPHA = D( IC )
217 BETA = E( IC )
218 CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA,
219 $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT,
220 $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO )
221 IF( INFO.NE.0 ) THEN
222 RETURN
223 END IF
224 40 CONTINUE
225 50 CONTINUE
226 *
227 RETURN
228 *
229 * End of DLASD0
230 *
231 END