LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaed0.f
Go to the documentation of this file.
1*> \brief \b DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAED0 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
22* WORK, IWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
26* ..
27* .. Array Arguments ..
28* INTEGER IWORK( * )
29* DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
30* $ WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DLAED0 computes all eigenvalues and corresponding eigenvectors of a
40*> symmetric tridiagonal matrix using the divide and conquer method.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] ICOMPQ
47*> \verbatim
48*> ICOMPQ is INTEGER
49*> = 0: Compute eigenvalues only.
50*> = 1: Compute eigenvectors of original dense symmetric matrix
51*> also. On entry, Q contains the orthogonal matrix used
52*> to reduce the original matrix to tridiagonal form.
53*> = 2: Compute eigenvalues and eigenvectors of tridiagonal
54*> matrix.
55*> \endverbatim
56*>
57*> \param[in] QSIZ
58*> \verbatim
59*> QSIZ is INTEGER
60*> The dimension of the orthogonal matrix used to reduce
61*> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*> N is INTEGER
67*> The dimension of the symmetric tridiagonal matrix. N >= 0.
68*> \endverbatim
69*>
70*> \param[in,out] D
71*> \verbatim
72*> D is DOUBLE PRECISION array, dimension (N)
73*> On entry, the main diagonal of the tridiagonal matrix.
74*> On exit, its eigenvalues.
75*> \endverbatim
76*>
77*> \param[in] E
78*> \verbatim
79*> E is DOUBLE PRECISION array, dimension (N-1)
80*> The off-diagonal elements of the tridiagonal matrix.
81*> On exit, E has been destroyed.
82*> \endverbatim
83*>
84*> \param[in,out] Q
85*> \verbatim
86*> Q is DOUBLE PRECISION array, dimension (LDQ, N)
87*> On entry, Q must contain an N-by-N orthogonal matrix.
88*> If ICOMPQ = 0 Q is not referenced.
89*> If ICOMPQ = 1 On entry, Q is a subset of the columns of the
90*> orthogonal matrix used to reduce the full
91*> matrix to tridiagonal form corresponding to
92*> the subset of the full matrix which is being
93*> decomposed at this time.
94*> If ICOMPQ = 2 On entry, Q will be the identity matrix.
95*> On exit, Q contains the eigenvectors of the
96*> tridiagonal matrix.
97*> \endverbatim
98*>
99*> \param[in] LDQ
100*> \verbatim
101*> LDQ is INTEGER
102*> The leading dimension of the array Q. If eigenvectors are
103*> desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
104*> \endverbatim
105*>
106*> \param[out] QSTORE
107*> \verbatim
108*> QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
109*> Referenced only when ICOMPQ = 1. Used to store parts of
110*> the eigenvector matrix when the updating matrix multiplies
111*> take place.
112*> \endverbatim
113*>
114*> \param[in] LDQS
115*> \verbatim
116*> LDQS is INTEGER
117*> The leading dimension of the array QSTORE. If ICOMPQ = 1,
118*> then LDQS >= max(1,N). In any case, LDQS >= 1.
119*> \endverbatim
120*>
121*> \param[out] WORK
122*> \verbatim
123*> WORK is DOUBLE PRECISION array,
124*> If ICOMPQ = 0 or 1, the dimension of WORK must be at least
125*> 1 + 3*N + 2*N*lg N + 3*N**2
126*> ( lg( N ) = smallest integer k
127*> such that 2^k >= N )
128*> If ICOMPQ = 2, the dimension of WORK must be at least
129*> 4*N + N**2.
130*> \endverbatim
131*>
132*> \param[out] IWORK
133*> \verbatim
134*> IWORK is INTEGER array,
135*> If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
136*> 6 + 6*N + 5*N*lg N.
137*> ( lg( N ) = smallest integer k
138*> such that 2^k >= N )
139*> If ICOMPQ = 2, the dimension of IWORK must be at least
140*> 3 + 5*N.
141*> \endverbatim
142*>
143*> \param[out] INFO
144*> \verbatim
145*> INFO is INTEGER
146*> = 0: successful exit.
147*> < 0: if INFO = -i, the i-th argument had an illegal value.
148*> > 0: The algorithm failed to compute an eigenvalue while
149*> working on the submatrix lying in rows and columns
150*> INFO/(N+1) through mod(INFO,N+1).
151*> \endverbatim
152*
153* Authors:
154* ========
155*
156*> \author Univ. of Tennessee
157*> \author Univ. of California Berkeley
158*> \author Univ. of Colorado Denver
159*> \author NAG Ltd.
160*
161*> \ingroup laed0
162*
163*> \par Contributors:
164* ==================
165*>
166*> Jeff Rutter, Computer Science Division, University of California
167*> at Berkeley, USA
168*
169* =====================================================================
170 SUBROUTINE dlaed0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
171 $ WORK, IWORK, INFO )
172*
173* -- LAPACK computational routine --
174* -- LAPACK is a software package provided by Univ. of Tennessee, --
175* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176*
177* .. Scalar Arguments ..
178 INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
179* ..
180* .. Array Arguments ..
181 INTEGER IWORK( * )
182 DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
183 $ work( * )
184* ..
185*
186* =====================================================================
187*
188* .. Parameters ..
189 DOUBLE PRECISION ZERO, ONE, TWO
190 parameter( zero = 0.d0, one = 1.d0, two = 2.d0 )
191* ..
192* .. Local Scalars ..
193 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
194 $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
195 $ j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1,
196 $ spm2, submat, subpbs, tlvls
197 DOUBLE PRECISION TEMP
198* ..
199* .. External Subroutines ..
200 EXTERNAL dcopy, dgemm, dlacpy, dlaed1, dlaed7, dsteqr,
201 $ xerbla
202* ..
203* .. External Functions ..
204 INTEGER ILAENV
205 EXTERNAL ilaenv
206* ..
207* .. Intrinsic Functions ..
208 INTRINSIC abs, dble, int, log, max
209* ..
210* .. Executable Statements ..
211*
212* Test the input parameters.
213*
214 info = 0
215*
216 IF( icompq.LT.0 .OR. icompq.GT.2 ) THEN
217 info = -1
218 ELSE IF( ( icompq.EQ.1 ) .AND. ( qsiz.LT.max( 0, n ) ) ) THEN
219 info = -2
220 ELSE IF( n.LT.0 ) THEN
221 info = -3
222 ELSE IF( ldq.LT.max( 1, n ) ) THEN
223 info = -7
224 ELSE IF( ldqs.LT.max( 1, n ) ) THEN
225 info = -9
226 END IF
227 IF( info.NE.0 ) THEN
228 CALL xerbla( 'DLAED0', -info )
229 RETURN
230 END IF
231*
232* Quick return if possible
233*
234 IF( n.EQ.0 )
235 $ RETURN
236*
237 smlsiz = ilaenv( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
238*
239* Determine the size and placement of the submatrices, and save in
240* the leading elements of IWORK.
241*
242 iwork( 1 ) = n
243 subpbs = 1
244 tlvls = 0
245 10 CONTINUE
246 IF( iwork( subpbs ).GT.smlsiz ) THEN
247 DO 20 j = subpbs, 1, -1
248 iwork( 2*j ) = ( iwork( j )+1 ) / 2
249 iwork( 2*j-1 ) = iwork( j ) / 2
250 20 CONTINUE
251 tlvls = tlvls + 1
252 subpbs = 2*subpbs
253 GO TO 10
254 END IF
255 DO 30 j = 2, subpbs
256 iwork( j ) = iwork( j ) + iwork( j-1 )
257 30 CONTINUE
258*
259* Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
260* using rank-1 modifications (cuts).
261*
262 spm1 = subpbs - 1
263 DO 40 i = 1, spm1
264 submat = iwork( i ) + 1
265 smm1 = submat - 1
266 d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
267 d( submat ) = d( submat ) - abs( e( smm1 ) )
268 40 CONTINUE
269*
270 indxq = 4*n + 3
271 IF( icompq.NE.2 ) THEN
272*
273* Set up workspaces for eigenvalues only/accumulate new vectors
274* routine
275*
276 temp = log( dble( n ) ) / log( two )
277 lgn = int( temp )
278 IF( 2**lgn.LT.n )
279 $ lgn = lgn + 1
280 IF( 2**lgn.LT.n )
281 $ lgn = lgn + 1
282 iprmpt = indxq + n + 1
283 iperm = iprmpt + n*lgn
284 iqptr = iperm + n*lgn
285 igivpt = iqptr + n + 2
286 igivcl = igivpt + n*lgn
287*
288 igivnm = 1
289 iq = igivnm + 2*n*lgn
290 iwrem = iq + n**2 + 1
291*
292* Initialize pointers
293*
294 DO 50 i = 0, subpbs
295 iwork( iprmpt+i ) = 1
296 iwork( igivpt+i ) = 1
297 50 CONTINUE
298 iwork( iqptr ) = 1
299 END IF
300*
301* Solve each submatrix eigenproblem at the bottom of the divide and
302* conquer tree.
303*
304 curr = 0
305 DO 70 i = 0, spm1
306 IF( i.EQ.0 ) THEN
307 submat = 1
308 matsiz = iwork( 1 )
309 ELSE
310 submat = iwork( i ) + 1
311 matsiz = iwork( i+1 ) - iwork( i )
312 END IF
313 IF( icompq.EQ.2 ) THEN
314 CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
315 $ q( submat, submat ), ldq, work, info )
316 IF( info.NE.0 )
317 $ GO TO 130
318 ELSE
319 CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
320 $ work( iq-1+iwork( iqptr+curr ) ), matsiz, work,
321 $ info )
322 IF( info.NE.0 )
323 $ GO TO 130
324 IF( icompq.EQ.1 ) THEN
325 CALL dgemm( 'N', 'N', qsiz, matsiz, matsiz, one,
326 $ q( 1, submat ), ldq, work( iq-1+iwork( iqptr+
327 $ curr ) ), matsiz, zero, qstore( 1, submat ),
328 $ ldqs )
329 END IF
330 iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
331 curr = curr + 1
332 END IF
333 k = 1
334 DO 60 j = submat, iwork( i+1 )
335 iwork( indxq+j ) = k
336 k = k + 1
337 60 CONTINUE
338 70 CONTINUE
339*
340* Successively merge eigensystems of adjacent submatrices
341* into eigensystem for the corresponding larger matrix.
342*
343* while ( SUBPBS > 1 )
344*
345 curlvl = 1
346 80 CONTINUE
347 IF( subpbs.GT.1 ) THEN
348 spm2 = subpbs - 2
349 DO 90 i = 0, spm2, 2
350 IF( i.EQ.0 ) THEN
351 submat = 1
352 matsiz = iwork( 2 )
353 msd2 = iwork( 1 )
354 curprb = 0
355 ELSE
356 submat = iwork( i ) + 1
357 matsiz = iwork( i+2 ) - iwork( i )
358 msd2 = matsiz / 2
359 curprb = curprb + 1
360 END IF
361*
362* Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
363* into an eigensystem of size MATSIZ.
364* DLAED1 is used only for the full eigensystem of a tridiagonal
365* matrix.
366* DLAED7 handles the cases in which eigenvalues only or eigenvalues
367* and eigenvectors of a full symmetric matrix (which was reduced to
368* tridiagonal form) are desired.
369*
370 IF( icompq.EQ.2 ) THEN
371 CALL dlaed1( matsiz, d( submat ), q( submat, submat ),
372 $ ldq, iwork( indxq+submat ),
373 $ e( submat+msd2-1 ), msd2, work,
374 $ iwork( subpbs+1 ), info )
375 ELSE
376 CALL dlaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,
377 $ d( submat ), qstore( 1, submat ), ldqs,
378 $ iwork( indxq+submat ), e( submat+msd2-1 ),
379 $ msd2, work( iq ), iwork( iqptr ),
380 $ iwork( iprmpt ), iwork( iperm ),
381 $ iwork( igivpt ), iwork( igivcl ),
382 $ work( igivnm ), work( iwrem ),
383 $ iwork( subpbs+1 ), info )
384 END IF
385 IF( info.NE.0 )
386 $ GO TO 130
387 iwork( i / 2+1 ) = iwork( i+2 )
388 90 CONTINUE
389 subpbs = subpbs / 2
390 curlvl = curlvl + 1
391 GO TO 80
392 END IF
393*
394* end while
395*
396* Re-merge the eigenvalues/vectors which were deflated at the final
397* merge step.
398*
399 IF( icompq.EQ.1 ) THEN
400 DO 100 i = 1, n
401 j = iwork( indxq+i )
402 work( i ) = d( j )
403 CALL dcopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
404 100 CONTINUE
405 CALL dcopy( n, work, 1, d, 1 )
406 ELSE IF( icompq.EQ.2 ) THEN
407 DO 110 i = 1, n
408 j = iwork( indxq+i )
409 work( i ) = d( j )
410 CALL dcopy( n, q( 1, j ), 1, work( n*i+1 ), 1 )
411 110 CONTINUE
412 CALL dcopy( n, work, 1, d, 1 )
413 CALL dlacpy( 'A', n, n, work( n+1 ), n, q, ldq )
414 ELSE
415 DO 120 i = 1, n
416 j = iwork( indxq+i )
417 work( i ) = d( j )
418 120 CONTINUE
419 CALL dcopy( n, work, 1, d, 1 )
420 END IF
421 GO TO 140
422*
423 130 CONTINUE
424 info = submat*( n+1 ) + submat + matsiz - 1
425*
426 140 CONTINUE
427 RETURN
428*
429* End of DLAED0
430*
431 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dlaed0(icompq, qsiz, n, d, e, q, ldq, qstore, ldqs, work, iwork, info)
DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition dlaed0.f:172
subroutine dlaed1(n, d, q, ldq, indxq, rho, cutpnt, work, iwork, info)
DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition dlaed1.f:163
subroutine dlaed7(icompq, n, qsiz, tlvls, curlvl, curpbm, d, q, ldq, indxq, rho, cutpnt, qstore, qptr, prmptr, perm, givptr, givcol, givnum, work, iwork, info)
DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition dlaed7.f:260
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131