1 SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
2 $ WORK, 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 ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAED0 computes all eigenvalues and corresponding eigenvectors of a
22 * symmetric tridiagonal matrix using the divide and conquer method.
23 *
24 * Arguments
25 * =========
26 *
27 * ICOMPQ (input) INTEGER
28 * = 0: Compute eigenvalues only.
29 * = 1: Compute eigenvectors of original dense symmetric matrix
30 * also. On entry, Q contains the orthogonal matrix used
31 * to reduce the original matrix to tridiagonal form.
32 * = 2: Compute eigenvalues and eigenvectors of tridiagonal
33 * matrix.
34 *
35 * QSIZ (input) INTEGER
36 * The dimension of the orthogonal matrix used to reduce
37 * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
38 *
39 * N (input) INTEGER
40 * The dimension of the symmetric tridiagonal matrix. N >= 0.
41 *
42 * D (input/output) DOUBLE PRECISION array, dimension (N)
43 * On entry, the main diagonal of the tridiagonal matrix.
44 * On exit, its eigenvalues.
45 *
46 * E (input) DOUBLE PRECISION array, dimension (N-1)
47 * The off-diagonal elements of the tridiagonal matrix.
48 * On exit, E has been destroyed.
49 *
50 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
51 * On entry, Q must contain an N-by-N orthogonal matrix.
52 * If ICOMPQ = 0 Q is not referenced.
53 * If ICOMPQ = 1 On entry, Q is a subset of the columns of the
54 * orthogonal matrix used to reduce the full
55 * matrix to tridiagonal form corresponding to
56 * the subset of the full matrix which is being
57 * decomposed at this time.
58 * If ICOMPQ = 2 On entry, Q will be the identity matrix.
59 * On exit, Q contains the eigenvectors of the
60 * tridiagonal matrix.
61 *
62 * LDQ (input) INTEGER
63 * The leading dimension of the array Q. If eigenvectors are
64 * desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
65 *
66 * QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
67 * Referenced only when ICOMPQ = 1. Used to store parts of
68 * the eigenvector matrix when the updating matrix multiplies
69 * take place.
70 *
71 * LDQS (input) INTEGER
72 * The leading dimension of the array QSTORE. If ICOMPQ = 1,
73 * then LDQS >= max(1,N). In any case, LDQS >= 1.
74 *
75 * WORK (workspace) DOUBLE PRECISION array,
76 * If ICOMPQ = 0 or 1, the dimension of WORK must be at least
77 * 1 + 3*N + 2*N*lg N + 2*N**2
78 * ( lg( N ) = smallest integer k
79 * such that 2^k >= N )
80 * If ICOMPQ = 2, the dimension of WORK must be at least
81 * 4*N + N**2.
82 *
83 * IWORK (workspace) INTEGER array,
84 * If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
85 * 6 + 6*N + 5*N*lg N.
86 * ( lg( N ) = smallest integer k
87 * such that 2^k >= N )
88 * If ICOMPQ = 2, the dimension of IWORK must be at least
89 * 3 + 5*N.
90 *
91 * INFO (output) INTEGER
92 * = 0: successful exit.
93 * < 0: if INFO = -i, the i-th argument had an illegal value.
94 * > 0: The algorithm failed to compute an eigenvalue while
95 * working on the submatrix lying in rows and columns
96 * INFO/(N+1) through mod(INFO,N+1).
97 *
98 * Further Details
99 * ===============
100 *
101 * Based on contributions by
102 * Jeff Rutter, Computer Science Division, University of California
103 * at Berkeley, USA
104 *
105 * =====================================================================
106 *
107 * .. Parameters ..
108 DOUBLE PRECISION ZERO, ONE, TWO
109 PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
110 * ..
111 * .. Local Scalars ..
112 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
113 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
114 $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
115 $ SPM2, SUBMAT, SUBPBS, TLVLS
116 DOUBLE PRECISION TEMP
117 * ..
118 * .. External Subroutines ..
119 EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
120 $ XERBLA
121 * ..
122 * .. External Functions ..
123 INTEGER ILAENV
124 EXTERNAL ILAENV
125 * ..
126 * .. Intrinsic Functions ..
127 INTRINSIC ABS, DBLE, INT, LOG, MAX
128 * ..
129 * .. Executable Statements ..
130 *
131 * Test the input parameters.
132 *
133 INFO = 0
134 *
135 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
136 INFO = -1
137 ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
138 INFO = -2
139 ELSE IF( N.LT.0 ) THEN
140 INFO = -3
141 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
142 INFO = -7
143 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
144 INFO = -9
145 END IF
146 IF( INFO.NE.0 ) THEN
147 CALL XERBLA( 'DLAED0', -INFO )
148 RETURN
149 END IF
150 *
151 * Quick return if possible
152 *
153 IF( N.EQ.0 )
154 $ RETURN
155 *
156 SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
157 *
158 * Determine the size and placement of the submatrices, and save in
159 * the leading elements of IWORK.
160 *
161 IWORK( 1 ) = N
162 SUBPBS = 1
163 TLVLS = 0
164 10 CONTINUE
165 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
166 DO 20 J = SUBPBS, 1, -1
167 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
168 IWORK( 2*J-1 ) = IWORK( J ) / 2
169 20 CONTINUE
170 TLVLS = TLVLS + 1
171 SUBPBS = 2*SUBPBS
172 GO TO 10
173 END IF
174 DO 30 J = 2, SUBPBS
175 IWORK( J ) = IWORK( J ) + IWORK( J-1 )
176 30 CONTINUE
177 *
178 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
179 * using rank-1 modifications (cuts).
180 *
181 SPM1 = SUBPBS - 1
182 DO 40 I = 1, SPM1
183 SUBMAT = IWORK( I ) + 1
184 SMM1 = SUBMAT - 1
185 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
186 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
187 40 CONTINUE
188 *
189 INDXQ = 4*N + 3
190 IF( ICOMPQ.NE.2 ) THEN
191 *
192 * Set up workspaces for eigenvalues only/accumulate new vectors
193 * routine
194 *
195 TEMP = LOG( DBLE( N ) ) / LOG( TWO )
196 LGN = INT( TEMP )
197 IF( 2**LGN.LT.N )
198 $ LGN = LGN + 1
199 IF( 2**LGN.LT.N )
200 $ LGN = LGN + 1
201 IPRMPT = INDXQ + N + 1
202 IPERM = IPRMPT + N*LGN
203 IQPTR = IPERM + N*LGN
204 IGIVPT = IQPTR + N + 2
205 IGIVCL = IGIVPT + N*LGN
206 *
207 IGIVNM = 1
208 IQ = IGIVNM + 2*N*LGN
209 IWREM = IQ + N**2 + 1
210 *
211 * Initialize pointers
212 *
213 DO 50 I = 0, SUBPBS
214 IWORK( IPRMPT+I ) = 1
215 IWORK( IGIVPT+I ) = 1
216 50 CONTINUE
217 IWORK( IQPTR ) = 1
218 END IF
219 *
220 * Solve each submatrix eigenproblem at the bottom of the divide and
221 * conquer tree.
222 *
223 CURR = 0
224 DO 70 I = 0, SPM1
225 IF( I.EQ.0 ) THEN
226 SUBMAT = 1
227 MATSIZ = IWORK( 1 )
228 ELSE
229 SUBMAT = IWORK( I ) + 1
230 MATSIZ = IWORK( I+1 ) - IWORK( I )
231 END IF
232 IF( ICOMPQ.EQ.2 ) THEN
233 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
234 $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
235 IF( INFO.NE.0 )
236 $ GO TO 130
237 ELSE
238 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
239 $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
240 $ INFO )
241 IF( INFO.NE.0 )
242 $ GO TO 130
243 IF( ICOMPQ.EQ.1 ) THEN
244 CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
245 $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
246 $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
247 $ LDQS )
248 END IF
249 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
250 CURR = CURR + 1
251 END IF
252 K = 1
253 DO 60 J = SUBMAT, IWORK( I+1 )
254 IWORK( INDXQ+J ) = K
255 K = K + 1
256 60 CONTINUE
257 70 CONTINUE
258 *
259 * Successively merge eigensystems of adjacent submatrices
260 * into eigensystem for the corresponding larger matrix.
261 *
262 * while ( SUBPBS > 1 )
263 *
264 CURLVL = 1
265 80 CONTINUE
266 IF( SUBPBS.GT.1 ) THEN
267 SPM2 = SUBPBS - 2
268 DO 90 I = 0, SPM2, 2
269 IF( I.EQ.0 ) THEN
270 SUBMAT = 1
271 MATSIZ = IWORK( 2 )
272 MSD2 = IWORK( 1 )
273 CURPRB = 0
274 ELSE
275 SUBMAT = IWORK( I ) + 1
276 MATSIZ = IWORK( I+2 ) - IWORK( I )
277 MSD2 = MATSIZ / 2
278 CURPRB = CURPRB + 1
279 END IF
280 *
281 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
282 * into an eigensystem of size MATSIZ.
283 * DLAED1 is used only for the full eigensystem of a tridiagonal
284 * matrix.
285 * DLAED7 handles the cases in which eigenvalues only or eigenvalues
286 * and eigenvectors of a full symmetric matrix (which was reduced to
287 * tridiagonal form) are desired.
288 *
289 IF( ICOMPQ.EQ.2 ) THEN
290 CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
291 $ LDQ, IWORK( INDXQ+SUBMAT ),
292 $ E( SUBMAT+MSD2-1 ), MSD2, WORK,
293 $ IWORK( SUBPBS+1 ), INFO )
294 ELSE
295 CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
296 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
297 $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
298 $ MSD2, WORK( IQ ), IWORK( IQPTR ),
299 $ IWORK( IPRMPT ), IWORK( IPERM ),
300 $ IWORK( IGIVPT ), IWORK( IGIVCL ),
301 $ WORK( IGIVNM ), WORK( IWREM ),
302 $ IWORK( SUBPBS+1 ), INFO )
303 END IF
304 IF( INFO.NE.0 )
305 $ GO TO 130
306 IWORK( I / 2+1 ) = IWORK( I+2 )
307 90 CONTINUE
308 SUBPBS = SUBPBS / 2
309 CURLVL = CURLVL + 1
310 GO TO 80
311 END IF
312 *
313 * end while
314 *
315 * Re-merge the eigenvalues/vectors which were deflated at the final
316 * merge step.
317 *
318 IF( ICOMPQ.EQ.1 ) THEN
319 DO 100 I = 1, N
320 J = IWORK( INDXQ+I )
321 WORK( I ) = D( J )
322 CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
323 100 CONTINUE
324 CALL DCOPY( N, WORK, 1, D, 1 )
325 ELSE IF( ICOMPQ.EQ.2 ) THEN
326 DO 110 I = 1, N
327 J = IWORK( INDXQ+I )
328 WORK( I ) = D( J )
329 CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
330 110 CONTINUE
331 CALL DCOPY( N, WORK, 1, D, 1 )
332 CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
333 ELSE
334 DO 120 I = 1, N
335 J = IWORK( INDXQ+I )
336 WORK( I ) = D( J )
337 120 CONTINUE
338 CALL DCOPY( N, WORK, 1, D, 1 )
339 END IF
340 GO TO 140
341 *
342 130 CONTINUE
343 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
344 *
345 140 CONTINUE
346 RETURN
347 *
348 * End of DLAED0
349 *
350 END
2 $ WORK, 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 ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAED0 computes all eigenvalues and corresponding eigenvectors of a
22 * symmetric tridiagonal matrix using the divide and conquer method.
23 *
24 * Arguments
25 * =========
26 *
27 * ICOMPQ (input) INTEGER
28 * = 0: Compute eigenvalues only.
29 * = 1: Compute eigenvectors of original dense symmetric matrix
30 * also. On entry, Q contains the orthogonal matrix used
31 * to reduce the original matrix to tridiagonal form.
32 * = 2: Compute eigenvalues and eigenvectors of tridiagonal
33 * matrix.
34 *
35 * QSIZ (input) INTEGER
36 * The dimension of the orthogonal matrix used to reduce
37 * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
38 *
39 * N (input) INTEGER
40 * The dimension of the symmetric tridiagonal matrix. N >= 0.
41 *
42 * D (input/output) DOUBLE PRECISION array, dimension (N)
43 * On entry, the main diagonal of the tridiagonal matrix.
44 * On exit, its eigenvalues.
45 *
46 * E (input) DOUBLE PRECISION array, dimension (N-1)
47 * The off-diagonal elements of the tridiagonal matrix.
48 * On exit, E has been destroyed.
49 *
50 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
51 * On entry, Q must contain an N-by-N orthogonal matrix.
52 * If ICOMPQ = 0 Q is not referenced.
53 * If ICOMPQ = 1 On entry, Q is a subset of the columns of the
54 * orthogonal matrix used to reduce the full
55 * matrix to tridiagonal form corresponding to
56 * the subset of the full matrix which is being
57 * decomposed at this time.
58 * If ICOMPQ = 2 On entry, Q will be the identity matrix.
59 * On exit, Q contains the eigenvectors of the
60 * tridiagonal matrix.
61 *
62 * LDQ (input) INTEGER
63 * The leading dimension of the array Q. If eigenvectors are
64 * desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
65 *
66 * QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
67 * Referenced only when ICOMPQ = 1. Used to store parts of
68 * the eigenvector matrix when the updating matrix multiplies
69 * take place.
70 *
71 * LDQS (input) INTEGER
72 * The leading dimension of the array QSTORE. If ICOMPQ = 1,
73 * then LDQS >= max(1,N). In any case, LDQS >= 1.
74 *
75 * WORK (workspace) DOUBLE PRECISION array,
76 * If ICOMPQ = 0 or 1, the dimension of WORK must be at least
77 * 1 + 3*N + 2*N*lg N + 2*N**2
78 * ( lg( N ) = smallest integer k
79 * such that 2^k >= N )
80 * If ICOMPQ = 2, the dimension of WORK must be at least
81 * 4*N + N**2.
82 *
83 * IWORK (workspace) INTEGER array,
84 * If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
85 * 6 + 6*N + 5*N*lg N.
86 * ( lg( N ) = smallest integer k
87 * such that 2^k >= N )
88 * If ICOMPQ = 2, the dimension of IWORK must be at least
89 * 3 + 5*N.
90 *
91 * INFO (output) INTEGER
92 * = 0: successful exit.
93 * < 0: if INFO = -i, the i-th argument had an illegal value.
94 * > 0: The algorithm failed to compute an eigenvalue while
95 * working on the submatrix lying in rows and columns
96 * INFO/(N+1) through mod(INFO,N+1).
97 *
98 * Further Details
99 * ===============
100 *
101 * Based on contributions by
102 * Jeff Rutter, Computer Science Division, University of California
103 * at Berkeley, USA
104 *
105 * =====================================================================
106 *
107 * .. Parameters ..
108 DOUBLE PRECISION ZERO, ONE, TWO
109 PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
110 * ..
111 * .. Local Scalars ..
112 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
113 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
114 $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
115 $ SPM2, SUBMAT, SUBPBS, TLVLS
116 DOUBLE PRECISION TEMP
117 * ..
118 * .. External Subroutines ..
119 EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
120 $ XERBLA
121 * ..
122 * .. External Functions ..
123 INTEGER ILAENV
124 EXTERNAL ILAENV
125 * ..
126 * .. Intrinsic Functions ..
127 INTRINSIC ABS, DBLE, INT, LOG, MAX
128 * ..
129 * .. Executable Statements ..
130 *
131 * Test the input parameters.
132 *
133 INFO = 0
134 *
135 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
136 INFO = -1
137 ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
138 INFO = -2
139 ELSE IF( N.LT.0 ) THEN
140 INFO = -3
141 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
142 INFO = -7
143 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
144 INFO = -9
145 END IF
146 IF( INFO.NE.0 ) THEN
147 CALL XERBLA( 'DLAED0', -INFO )
148 RETURN
149 END IF
150 *
151 * Quick return if possible
152 *
153 IF( N.EQ.0 )
154 $ RETURN
155 *
156 SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
157 *
158 * Determine the size and placement of the submatrices, and save in
159 * the leading elements of IWORK.
160 *
161 IWORK( 1 ) = N
162 SUBPBS = 1
163 TLVLS = 0
164 10 CONTINUE
165 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
166 DO 20 J = SUBPBS, 1, -1
167 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
168 IWORK( 2*J-1 ) = IWORK( J ) / 2
169 20 CONTINUE
170 TLVLS = TLVLS + 1
171 SUBPBS = 2*SUBPBS
172 GO TO 10
173 END IF
174 DO 30 J = 2, SUBPBS
175 IWORK( J ) = IWORK( J ) + IWORK( J-1 )
176 30 CONTINUE
177 *
178 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
179 * using rank-1 modifications (cuts).
180 *
181 SPM1 = SUBPBS - 1
182 DO 40 I = 1, SPM1
183 SUBMAT = IWORK( I ) + 1
184 SMM1 = SUBMAT - 1
185 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
186 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
187 40 CONTINUE
188 *
189 INDXQ = 4*N + 3
190 IF( ICOMPQ.NE.2 ) THEN
191 *
192 * Set up workspaces for eigenvalues only/accumulate new vectors
193 * routine
194 *
195 TEMP = LOG( DBLE( N ) ) / LOG( TWO )
196 LGN = INT( TEMP )
197 IF( 2**LGN.LT.N )
198 $ LGN = LGN + 1
199 IF( 2**LGN.LT.N )
200 $ LGN = LGN + 1
201 IPRMPT = INDXQ + N + 1
202 IPERM = IPRMPT + N*LGN
203 IQPTR = IPERM + N*LGN
204 IGIVPT = IQPTR + N + 2
205 IGIVCL = IGIVPT + N*LGN
206 *
207 IGIVNM = 1
208 IQ = IGIVNM + 2*N*LGN
209 IWREM = IQ + N**2 + 1
210 *
211 * Initialize pointers
212 *
213 DO 50 I = 0, SUBPBS
214 IWORK( IPRMPT+I ) = 1
215 IWORK( IGIVPT+I ) = 1
216 50 CONTINUE
217 IWORK( IQPTR ) = 1
218 END IF
219 *
220 * Solve each submatrix eigenproblem at the bottom of the divide and
221 * conquer tree.
222 *
223 CURR = 0
224 DO 70 I = 0, SPM1
225 IF( I.EQ.0 ) THEN
226 SUBMAT = 1
227 MATSIZ = IWORK( 1 )
228 ELSE
229 SUBMAT = IWORK( I ) + 1
230 MATSIZ = IWORK( I+1 ) - IWORK( I )
231 END IF
232 IF( ICOMPQ.EQ.2 ) THEN
233 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
234 $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
235 IF( INFO.NE.0 )
236 $ GO TO 130
237 ELSE
238 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
239 $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
240 $ INFO )
241 IF( INFO.NE.0 )
242 $ GO TO 130
243 IF( ICOMPQ.EQ.1 ) THEN
244 CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
245 $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
246 $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
247 $ LDQS )
248 END IF
249 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
250 CURR = CURR + 1
251 END IF
252 K = 1
253 DO 60 J = SUBMAT, IWORK( I+1 )
254 IWORK( INDXQ+J ) = K
255 K = K + 1
256 60 CONTINUE
257 70 CONTINUE
258 *
259 * Successively merge eigensystems of adjacent submatrices
260 * into eigensystem for the corresponding larger matrix.
261 *
262 * while ( SUBPBS > 1 )
263 *
264 CURLVL = 1
265 80 CONTINUE
266 IF( SUBPBS.GT.1 ) THEN
267 SPM2 = SUBPBS - 2
268 DO 90 I = 0, SPM2, 2
269 IF( I.EQ.0 ) THEN
270 SUBMAT = 1
271 MATSIZ = IWORK( 2 )
272 MSD2 = IWORK( 1 )
273 CURPRB = 0
274 ELSE
275 SUBMAT = IWORK( I ) + 1
276 MATSIZ = IWORK( I+2 ) - IWORK( I )
277 MSD2 = MATSIZ / 2
278 CURPRB = CURPRB + 1
279 END IF
280 *
281 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
282 * into an eigensystem of size MATSIZ.
283 * DLAED1 is used only for the full eigensystem of a tridiagonal
284 * matrix.
285 * DLAED7 handles the cases in which eigenvalues only or eigenvalues
286 * and eigenvectors of a full symmetric matrix (which was reduced to
287 * tridiagonal form) are desired.
288 *
289 IF( ICOMPQ.EQ.2 ) THEN
290 CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
291 $ LDQ, IWORK( INDXQ+SUBMAT ),
292 $ E( SUBMAT+MSD2-1 ), MSD2, WORK,
293 $ IWORK( SUBPBS+1 ), INFO )
294 ELSE
295 CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
296 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
297 $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
298 $ MSD2, WORK( IQ ), IWORK( IQPTR ),
299 $ IWORK( IPRMPT ), IWORK( IPERM ),
300 $ IWORK( IGIVPT ), IWORK( IGIVCL ),
301 $ WORK( IGIVNM ), WORK( IWREM ),
302 $ IWORK( SUBPBS+1 ), INFO )
303 END IF
304 IF( INFO.NE.0 )
305 $ GO TO 130
306 IWORK( I / 2+1 ) = IWORK( I+2 )
307 90 CONTINUE
308 SUBPBS = SUBPBS / 2
309 CURLVL = CURLVL + 1
310 GO TO 80
311 END IF
312 *
313 * end while
314 *
315 * Re-merge the eigenvalues/vectors which were deflated at the final
316 * merge step.
317 *
318 IF( ICOMPQ.EQ.1 ) THEN
319 DO 100 I = 1, N
320 J = IWORK( INDXQ+I )
321 WORK( I ) = D( J )
322 CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
323 100 CONTINUE
324 CALL DCOPY( N, WORK, 1, D, 1 )
325 ELSE IF( ICOMPQ.EQ.2 ) THEN
326 DO 110 I = 1, N
327 J = IWORK( INDXQ+I )
328 WORK( I ) = D( J )
329 CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
330 110 CONTINUE
331 CALL DCOPY( N, WORK, 1, D, 1 )
332 CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
333 ELSE
334 DO 120 I = 1, N
335 J = IWORK( INDXQ+I )
336 WORK( I ) = D( J )
337 120 CONTINUE
338 CALL DCOPY( N, WORK, 1, D, 1 )
339 END IF
340 GO TO 140
341 *
342 130 CONTINUE
343 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
344 *
345 140 CONTINUE
346 RETURN
347 *
348 * End of DLAED0
349 *
350 END