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