1 SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
2 $ IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.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 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDQ, LDQS, N, QSIZ
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
15 COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * Using the divide and conquer method, ZLAED0 computes all eigenvalues
22 * of a symmetric tridiagonal matrix which is one diagonal block of
23 * those from reducing a dense or band Hermitian matrix and
24 * corresponding eigenvectors of the dense or band matrix.
25 *
26 * Arguments
27 * =========
28 *
29 * QSIZ (input) INTEGER
30 * The dimension of the unitary matrix used to reduce
31 * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
32 *
33 * N (input) INTEGER
34 * The dimension of the symmetric tridiagonal matrix. N >= 0.
35 *
36 * D (input/output) DOUBLE PRECISION array, dimension (N)
37 * On entry, the diagonal elements of the tridiagonal matrix.
38 * On exit, the eigenvalues in ascending order.
39 *
40 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
41 * On entry, the off-diagonal elements of the tridiagonal matrix.
42 * On exit, E has been destroyed.
43 *
44 * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
45 * On entry, Q must contain an QSIZ x N matrix whose columns
46 * unitarily orthonormal. It is a part of the unitary matrix
47 * that reduces the full dense Hermitian matrix to a
48 * (reducible) symmetric tridiagonal matrix.
49 *
50 * LDQ (input) INTEGER
51 * The leading dimension of the array Q. LDQ >= max(1,N).
52 *
53 * IWORK (workspace) INTEGER array,
54 * the dimension of IWORK must be at least
55 * 6 + 6*N + 5*N*lg N
56 * ( lg( N ) = smallest integer k
57 * such that 2^k >= N )
58 *
59 * RWORK (workspace) DOUBLE PRECISION array,
60 * dimension (1 + 3*N + 2*N*lg N + 3*N**2)
61 * ( lg( N ) = smallest integer k
62 * such that 2^k >= N )
63 *
64 * QSTORE (workspace) COMPLEX*16 array, dimension (LDQS, N)
65 * Used to store parts of
66 * the eigenvector matrix when the updating matrix multiplies
67 * take place.
68 *
69 * LDQS (input) INTEGER
70 * The leading dimension of the array QSTORE.
71 * LDQS >= max(1,N).
72 *
73 * INFO (output) INTEGER
74 * = 0: successful exit.
75 * < 0: if INFO = -i, the i-th argument had an illegal value.
76 * > 0: The algorithm failed to compute an eigenvalue while
77 * working on the submatrix lying in rows and columns
78 * INFO/(N+1) through mod(INFO,N+1).
79 *
80 * =====================================================================
81 *
82 * Warning: N could be as big as QSIZ!
83 *
84 * .. Parameters ..
85 DOUBLE PRECISION TWO
86 PARAMETER ( TWO = 2.D+0 )
87 * ..
88 * .. Local Scalars ..
89 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
90 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
91 $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
92 $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
93 DOUBLE PRECISION TEMP
94 * ..
95 * .. External Subroutines ..
96 EXTERNAL DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
97 * ..
98 * .. External Functions ..
99 INTEGER ILAENV
100 EXTERNAL ILAENV
101 * ..
102 * .. Intrinsic Functions ..
103 INTRINSIC ABS, DBLE, INT, LOG, MAX
104 * ..
105 * .. Executable Statements ..
106 *
107 * Test the input parameters.
108 *
109 INFO = 0
110 *
111 * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
112 * INFO = -1
113 * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
114 * $ THEN
115 IF( QSIZ.LT.MAX( 0, N ) ) THEN
116 INFO = -1
117 ELSE IF( N.LT.0 ) THEN
118 INFO = -2
119 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
120 INFO = -6
121 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
122 INFO = -8
123 END IF
124 IF( INFO.NE.0 ) THEN
125 CALL XERBLA( 'ZLAED0', -INFO )
126 RETURN
127 END IF
128 *
129 * Quick return if possible
130 *
131 IF( N.EQ.0 )
132 $ RETURN
133 *
134 SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
135 *
136 * Determine the size and placement of the submatrices, and save in
137 * the leading elements of IWORK.
138 *
139 IWORK( 1 ) = N
140 SUBPBS = 1
141 TLVLS = 0
142 10 CONTINUE
143 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
144 DO 20 J = SUBPBS, 1, -1
145 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
146 IWORK( 2*J-1 ) = IWORK( J ) / 2
147 20 CONTINUE
148 TLVLS = TLVLS + 1
149 SUBPBS = 2*SUBPBS
150 GO TO 10
151 END IF
152 DO 30 J = 2, SUBPBS
153 IWORK( J ) = IWORK( J ) + IWORK( J-1 )
154 30 CONTINUE
155 *
156 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
157 * using rank-1 modifications (cuts).
158 *
159 SPM1 = SUBPBS - 1
160 DO 40 I = 1, SPM1
161 SUBMAT = IWORK( I ) + 1
162 SMM1 = SUBMAT - 1
163 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
164 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
165 40 CONTINUE
166 *
167 INDXQ = 4*N + 3
168 *
169 * Set up workspaces for eigenvalues only/accumulate new vectors
170 * routine
171 *
172 TEMP = LOG( DBLE( N ) ) / LOG( TWO )
173 LGN = INT( TEMP )
174 IF( 2**LGN.LT.N )
175 $ LGN = LGN + 1
176 IF( 2**LGN.LT.N )
177 $ LGN = LGN + 1
178 IPRMPT = INDXQ + N + 1
179 IPERM = IPRMPT + N*LGN
180 IQPTR = IPERM + N*LGN
181 IGIVPT = IQPTR + N + 2
182 IGIVCL = IGIVPT + N*LGN
183 *
184 IGIVNM = 1
185 IQ = IGIVNM + 2*N*LGN
186 IWREM = IQ + N**2 + 1
187 * Initialize pointers
188 DO 50 I = 0, SUBPBS
189 IWORK( IPRMPT+I ) = 1
190 IWORK( IGIVPT+I ) = 1
191 50 CONTINUE
192 IWORK( IQPTR ) = 1
193 *
194 * Solve each submatrix eigenproblem at the bottom of the divide and
195 * conquer tree.
196 *
197 CURR = 0
198 DO 70 I = 0, SPM1
199 IF( I.EQ.0 ) THEN
200 SUBMAT = 1
201 MATSIZ = IWORK( 1 )
202 ELSE
203 SUBMAT = IWORK( I ) + 1
204 MATSIZ = IWORK( I+1 ) - IWORK( I )
205 END IF
206 LL = IQ - 1 + IWORK( IQPTR+CURR )
207 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
208 $ RWORK( LL ), MATSIZ, RWORK, INFO )
209 CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
210 $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
211 $ RWORK( IWREM ) )
212 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
213 CURR = CURR + 1
214 IF( INFO.GT.0 ) THEN
215 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
216 RETURN
217 END IF
218 K = 1
219 DO 60 J = SUBMAT, IWORK( I+1 )
220 IWORK( INDXQ+J ) = K
221 K = K + 1
222 60 CONTINUE
223 70 CONTINUE
224 *
225 * Successively merge eigensystems of adjacent submatrices
226 * into eigensystem for the corresponding larger matrix.
227 *
228 * while ( SUBPBS > 1 )
229 *
230 CURLVL = 1
231 80 CONTINUE
232 IF( SUBPBS.GT.1 ) THEN
233 SPM2 = SUBPBS - 2
234 DO 90 I = 0, SPM2, 2
235 IF( I.EQ.0 ) THEN
236 SUBMAT = 1
237 MATSIZ = IWORK( 2 )
238 MSD2 = IWORK( 1 )
239 CURPRB = 0
240 ELSE
241 SUBMAT = IWORK( I ) + 1
242 MATSIZ = IWORK( I+2 ) - IWORK( I )
243 MSD2 = MATSIZ / 2
244 CURPRB = CURPRB + 1
245 END IF
246 *
247 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
248 * into an eigensystem of size MATSIZ. ZLAED7 handles the case
249 * when the eigenvectors of a full or band Hermitian matrix (which
250 * was reduced to tridiagonal form) are desired.
251 *
252 * I am free to use Q as a valuable working space until Loop 150.
253 *
254 CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
255 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
256 $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
257 $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
258 $ IWORK( IPERM ), IWORK( IGIVPT ),
259 $ IWORK( IGIVCL ), RWORK( IGIVNM ),
260 $ Q( 1, SUBMAT ), RWORK( IWREM ),
261 $ IWORK( SUBPBS+1 ), INFO )
262 IF( INFO.GT.0 ) THEN
263 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
264 RETURN
265 END IF
266 IWORK( I / 2+1 ) = IWORK( I+2 )
267 90 CONTINUE
268 SUBPBS = SUBPBS / 2
269 CURLVL = CURLVL + 1
270 GO TO 80
271 END IF
272 *
273 * end while
274 *
275 * Re-merge the eigenvalues/vectors which were deflated at the final
276 * merge step.
277 *
278 DO 100 I = 1, N
279 J = IWORK( INDXQ+I )
280 RWORK( I ) = D( J )
281 CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
282 100 CONTINUE
283 CALL DCOPY( N, RWORK, 1, D, 1 )
284 *
285 RETURN
286 *
287 * End of ZLAED0
288 *
289 END
2 $ IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.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 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDQ, LDQS, N, QSIZ
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
15 COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * Using the divide and conquer method, ZLAED0 computes all eigenvalues
22 * of a symmetric tridiagonal matrix which is one diagonal block of
23 * those from reducing a dense or band Hermitian matrix and
24 * corresponding eigenvectors of the dense or band matrix.
25 *
26 * Arguments
27 * =========
28 *
29 * QSIZ (input) INTEGER
30 * The dimension of the unitary matrix used to reduce
31 * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
32 *
33 * N (input) INTEGER
34 * The dimension of the symmetric tridiagonal matrix. N >= 0.
35 *
36 * D (input/output) DOUBLE PRECISION array, dimension (N)
37 * On entry, the diagonal elements of the tridiagonal matrix.
38 * On exit, the eigenvalues in ascending order.
39 *
40 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
41 * On entry, the off-diagonal elements of the tridiagonal matrix.
42 * On exit, E has been destroyed.
43 *
44 * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
45 * On entry, Q must contain an QSIZ x N matrix whose columns
46 * unitarily orthonormal. It is a part of the unitary matrix
47 * that reduces the full dense Hermitian matrix to a
48 * (reducible) symmetric tridiagonal matrix.
49 *
50 * LDQ (input) INTEGER
51 * The leading dimension of the array Q. LDQ >= max(1,N).
52 *
53 * IWORK (workspace) INTEGER array,
54 * the dimension of IWORK must be at least
55 * 6 + 6*N + 5*N*lg N
56 * ( lg( N ) = smallest integer k
57 * such that 2^k >= N )
58 *
59 * RWORK (workspace) DOUBLE PRECISION array,
60 * dimension (1 + 3*N + 2*N*lg N + 3*N**2)
61 * ( lg( N ) = smallest integer k
62 * such that 2^k >= N )
63 *
64 * QSTORE (workspace) COMPLEX*16 array, dimension (LDQS, N)
65 * Used to store parts of
66 * the eigenvector matrix when the updating matrix multiplies
67 * take place.
68 *
69 * LDQS (input) INTEGER
70 * The leading dimension of the array QSTORE.
71 * LDQS >= max(1,N).
72 *
73 * INFO (output) INTEGER
74 * = 0: successful exit.
75 * < 0: if INFO = -i, the i-th argument had an illegal value.
76 * > 0: The algorithm failed to compute an eigenvalue while
77 * working on the submatrix lying in rows and columns
78 * INFO/(N+1) through mod(INFO,N+1).
79 *
80 * =====================================================================
81 *
82 * Warning: N could be as big as QSIZ!
83 *
84 * .. Parameters ..
85 DOUBLE PRECISION TWO
86 PARAMETER ( TWO = 2.D+0 )
87 * ..
88 * .. Local Scalars ..
89 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
90 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
91 $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
92 $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
93 DOUBLE PRECISION TEMP
94 * ..
95 * .. External Subroutines ..
96 EXTERNAL DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
97 * ..
98 * .. External Functions ..
99 INTEGER ILAENV
100 EXTERNAL ILAENV
101 * ..
102 * .. Intrinsic Functions ..
103 INTRINSIC ABS, DBLE, INT, LOG, MAX
104 * ..
105 * .. Executable Statements ..
106 *
107 * Test the input parameters.
108 *
109 INFO = 0
110 *
111 * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
112 * INFO = -1
113 * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
114 * $ THEN
115 IF( QSIZ.LT.MAX( 0, N ) ) THEN
116 INFO = -1
117 ELSE IF( N.LT.0 ) THEN
118 INFO = -2
119 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
120 INFO = -6
121 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
122 INFO = -8
123 END IF
124 IF( INFO.NE.0 ) THEN
125 CALL XERBLA( 'ZLAED0', -INFO )
126 RETURN
127 END IF
128 *
129 * Quick return if possible
130 *
131 IF( N.EQ.0 )
132 $ RETURN
133 *
134 SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
135 *
136 * Determine the size and placement of the submatrices, and save in
137 * the leading elements of IWORK.
138 *
139 IWORK( 1 ) = N
140 SUBPBS = 1
141 TLVLS = 0
142 10 CONTINUE
143 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
144 DO 20 J = SUBPBS, 1, -1
145 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
146 IWORK( 2*J-1 ) = IWORK( J ) / 2
147 20 CONTINUE
148 TLVLS = TLVLS + 1
149 SUBPBS = 2*SUBPBS
150 GO TO 10
151 END IF
152 DO 30 J = 2, SUBPBS
153 IWORK( J ) = IWORK( J ) + IWORK( J-1 )
154 30 CONTINUE
155 *
156 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
157 * using rank-1 modifications (cuts).
158 *
159 SPM1 = SUBPBS - 1
160 DO 40 I = 1, SPM1
161 SUBMAT = IWORK( I ) + 1
162 SMM1 = SUBMAT - 1
163 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
164 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
165 40 CONTINUE
166 *
167 INDXQ = 4*N + 3
168 *
169 * Set up workspaces for eigenvalues only/accumulate new vectors
170 * routine
171 *
172 TEMP = LOG( DBLE( N ) ) / LOG( TWO )
173 LGN = INT( TEMP )
174 IF( 2**LGN.LT.N )
175 $ LGN = LGN + 1
176 IF( 2**LGN.LT.N )
177 $ LGN = LGN + 1
178 IPRMPT = INDXQ + N + 1
179 IPERM = IPRMPT + N*LGN
180 IQPTR = IPERM + N*LGN
181 IGIVPT = IQPTR + N + 2
182 IGIVCL = IGIVPT + N*LGN
183 *
184 IGIVNM = 1
185 IQ = IGIVNM + 2*N*LGN
186 IWREM = IQ + N**2 + 1
187 * Initialize pointers
188 DO 50 I = 0, SUBPBS
189 IWORK( IPRMPT+I ) = 1
190 IWORK( IGIVPT+I ) = 1
191 50 CONTINUE
192 IWORK( IQPTR ) = 1
193 *
194 * Solve each submatrix eigenproblem at the bottom of the divide and
195 * conquer tree.
196 *
197 CURR = 0
198 DO 70 I = 0, SPM1
199 IF( I.EQ.0 ) THEN
200 SUBMAT = 1
201 MATSIZ = IWORK( 1 )
202 ELSE
203 SUBMAT = IWORK( I ) + 1
204 MATSIZ = IWORK( I+1 ) - IWORK( I )
205 END IF
206 LL = IQ - 1 + IWORK( IQPTR+CURR )
207 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
208 $ RWORK( LL ), MATSIZ, RWORK, INFO )
209 CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
210 $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
211 $ RWORK( IWREM ) )
212 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
213 CURR = CURR + 1
214 IF( INFO.GT.0 ) THEN
215 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
216 RETURN
217 END IF
218 K = 1
219 DO 60 J = SUBMAT, IWORK( I+1 )
220 IWORK( INDXQ+J ) = K
221 K = K + 1
222 60 CONTINUE
223 70 CONTINUE
224 *
225 * Successively merge eigensystems of adjacent submatrices
226 * into eigensystem for the corresponding larger matrix.
227 *
228 * while ( SUBPBS > 1 )
229 *
230 CURLVL = 1
231 80 CONTINUE
232 IF( SUBPBS.GT.1 ) THEN
233 SPM2 = SUBPBS - 2
234 DO 90 I = 0, SPM2, 2
235 IF( I.EQ.0 ) THEN
236 SUBMAT = 1
237 MATSIZ = IWORK( 2 )
238 MSD2 = IWORK( 1 )
239 CURPRB = 0
240 ELSE
241 SUBMAT = IWORK( I ) + 1
242 MATSIZ = IWORK( I+2 ) - IWORK( I )
243 MSD2 = MATSIZ / 2
244 CURPRB = CURPRB + 1
245 END IF
246 *
247 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
248 * into an eigensystem of size MATSIZ. ZLAED7 handles the case
249 * when the eigenvectors of a full or band Hermitian matrix (which
250 * was reduced to tridiagonal form) are desired.
251 *
252 * I am free to use Q as a valuable working space until Loop 150.
253 *
254 CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
255 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
256 $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
257 $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
258 $ IWORK( IPERM ), IWORK( IGIVPT ),
259 $ IWORK( IGIVCL ), RWORK( IGIVNM ),
260 $ Q( 1, SUBMAT ), RWORK( IWREM ),
261 $ IWORK( SUBPBS+1 ), INFO )
262 IF( INFO.GT.0 ) THEN
263 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
264 RETURN
265 END IF
266 IWORK( I / 2+1 ) = IWORK( I+2 )
267 90 CONTINUE
268 SUBPBS = SUBPBS / 2
269 CURLVL = CURLVL + 1
270 GO TO 80
271 END IF
272 *
273 * end while
274 *
275 * Re-merge the eigenvalues/vectors which were deflated at the final
276 * merge step.
277 *
278 DO 100 I = 1, N
279 J = IWORK( INDXQ+I )
280 RWORK( I ) = D( J )
281 CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
282 100 CONTINUE
283 CALL DCOPY( N, RWORK, 1, D, 1 )
284 *
285 RETURN
286 *
287 * End of ZLAED0
288 *
289 END