LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
claed0.f
Go to the documentation of this file.
1*> \brief \b CLAED0 used by CSTEDC. 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*> Download CLAED0 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claed0.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claed0.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claed0.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
20* IWORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER INFO, LDQ, LDQS, N, QSIZ
24* ..
25* .. Array Arguments ..
26* INTEGER IWORK( * )
27* REAL D( * ), E( * ), RWORK( * )
28* COMPLEX Q( LDQ, * ), QSTORE( LDQS, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> Using the divide and conquer method, CLAED0 computes all eigenvalues
38*> of a symmetric tridiagonal matrix which is one diagonal block of
39*> those from reducing a dense or band Hermitian matrix and
40*> corresponding eigenvectors of the dense or band matrix.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] QSIZ
47*> \verbatim
48*> QSIZ is INTEGER
49*> The dimension of the unitary matrix used to reduce
50*> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*> N is INTEGER
56*> The dimension of the symmetric tridiagonal matrix. N >= 0.
57*> \endverbatim
58*>
59*> \param[in,out] D
60*> \verbatim
61*> D is REAL array, dimension (N)
62*> On entry, the diagonal elements of the tridiagonal matrix.
63*> On exit, the eigenvalues in ascending order.
64*> \endverbatim
65*>
66*> \param[in,out] E
67*> \verbatim
68*> E is REAL array, dimension (N-1)
69*> On entry, the off-diagonal elements of the tridiagonal matrix.
70*> On exit, E has been destroyed.
71*> \endverbatim
72*>
73*> \param[in,out] Q
74*> \verbatim
75*> Q is COMPLEX array, dimension (LDQ,N)
76*> On entry, Q must contain an QSIZ x N matrix whose columns
77*> unitarily orthonormal. It is a part of the unitary matrix
78*> that reduces the full dense Hermitian matrix to a
79*> (reducible) symmetric tridiagonal matrix.
80*> \endverbatim
81*>
82*> \param[in] LDQ
83*> \verbatim
84*> LDQ is INTEGER
85*> The leading dimension of the array Q. LDQ >= max(1,N).
86*> \endverbatim
87*>
88*> \param[out] IWORK
89*> \verbatim
90*> IWORK is INTEGER array,
91*> the dimension of IWORK must be at least
92*> 6 + 6*N + 5*N*lg N
93*> ( lg( N ) = smallest integer k
94*> such that 2^k >= N )
95*> \endverbatim
96*>
97*> \param[out] RWORK
98*> \verbatim
99*> RWORK is REAL array,
100*> dimension (1 + 3*N + 2*N*lg N + 3*N**2)
101*> ( lg( N ) = smallest integer k
102*> such that 2^k >= N )
103*> \endverbatim
104*>
105*> \param[out] QSTORE
106*> \verbatim
107*> QSTORE is COMPLEX array, dimension (LDQS, N)
108*> Used to store parts of
109*> the eigenvector matrix when the updating matrix multiplies
110*> take place.
111*> \endverbatim
112*>
113*> \param[in] LDQS
114*> \verbatim
115*> LDQS is INTEGER
116*> The leading dimension of the array QSTORE.
117*> LDQS >= max(1,N).
118*> \endverbatim
119*>
120*> \param[out] INFO
121*> \verbatim
122*> INFO is INTEGER
123*> = 0: successful exit.
124*> < 0: if INFO = -i, the i-th argument had an illegal value.
125*> > 0: The algorithm failed to compute an eigenvalue while
126*> working on the submatrix lying in rows and columns
127*> INFO/(N+1) through mod(INFO,N+1).
128*> \endverbatim
129*
130* Authors:
131* ========
132*
133*> \author Univ. of Tennessee
134*> \author Univ. of California Berkeley
135*> \author Univ. of Colorado Denver
136*> \author NAG Ltd.
137*
138*> \ingroup laed0
139*
140* =====================================================================
141 SUBROUTINE claed0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
142 $ IWORK, INFO )
143*
144* -- LAPACK computational routine --
145* -- LAPACK is a software package provided by Univ. of Tennessee, --
146* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
147*
148* .. Scalar Arguments ..
149 INTEGER INFO, LDQ, LDQS, N, QSIZ
150* ..
151* .. Array Arguments ..
152 INTEGER IWORK( * )
153 REAL D( * ), E( * ), RWORK( * )
154 COMPLEX Q( LDQ, * ), QSTORE( LDQS, * )
155* ..
156*
157* =====================================================================
158*
159* Warning: N could be as big as QSIZ!
160*
161* .. Parameters ..
162 REAL TWO
163 parameter( two = 2.e+0 )
164* ..
165* .. Local Scalars ..
166 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
167 $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
168 $ j, k, lgn, ll, matsiz, msd2, smlsiz, smm1,
169 $ spm1, spm2, submat, subpbs, tlvls
170 REAL TEMP
171* ..
172* .. External Subroutines ..
173 EXTERNAL ccopy, clacrm, claed7, scopy, ssteqr,
174 $ xerbla
175* ..
176* .. External Functions ..
177 INTEGER ILAENV
178 EXTERNAL ilaenv
179* ..
180* .. Intrinsic Functions ..
181 INTRINSIC abs, int, log, max, real
182* ..
183* .. Executable Statements ..
184*
185* Test the input parameters.
186*
187 info = 0
188*
189* IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
190* INFO = -1
191* ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
192* $ THEN
193 IF( qsiz.LT.max( 0, n ) ) THEN
194 info = -1
195 ELSE IF( n.LT.0 ) THEN
196 info = -2
197 ELSE IF( ldq.LT.max( 1, n ) ) THEN
198 info = -6
199 ELSE IF( ldqs.LT.max( 1, n ) ) THEN
200 info = -8
201 END IF
202 IF( info.NE.0 ) THEN
203 CALL xerbla( 'CLAED0', -info )
204 RETURN
205 END IF
206*
207* Quick return if possible
208*
209 IF( n.EQ.0 )
210 $ RETURN
211*
212 smlsiz = ilaenv( 9, 'CLAED0', ' ', 0, 0, 0, 0 )
213*
214* Determine the size and placement of the submatrices, and save in
215* the leading elements of IWORK.
216*
217 iwork( 1 ) = n
218 subpbs = 1
219 tlvls = 0
220 10 CONTINUE
221 IF( iwork( subpbs ).GT.smlsiz ) THEN
222 DO 20 j = subpbs, 1, -1
223 iwork( 2*j ) = ( iwork( j )+1 ) / 2
224 iwork( 2*j-1 ) = iwork( j ) / 2
225 20 CONTINUE
226 tlvls = tlvls + 1
227 subpbs = 2*subpbs
228 GO TO 10
229 END IF
230 DO 30 j = 2, subpbs
231 iwork( j ) = iwork( j ) + iwork( j-1 )
232 30 CONTINUE
233*
234* Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
235* using rank-1 modifications (cuts).
236*
237 spm1 = subpbs - 1
238 DO 40 i = 1, spm1
239 submat = iwork( i ) + 1
240 smm1 = submat - 1
241 d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
242 d( submat ) = d( submat ) - abs( e( smm1 ) )
243 40 CONTINUE
244*
245 indxq = 4*n + 3
246*
247* Set up workspaces for eigenvalues only/accumulate new vectors
248* routine
249*
250 temp = log( real( n ) ) / log( two )
251 lgn = int( temp )
252 IF( 2**lgn.LT.n )
253 $ lgn = lgn + 1
254 IF( 2**lgn.LT.n )
255 $ lgn = lgn + 1
256 iprmpt = indxq + n + 1
257 iperm = iprmpt + n*lgn
258 iqptr = iperm + n*lgn
259 igivpt = iqptr + n + 2
260 igivcl = igivpt + n*lgn
261*
262 igivnm = 1
263 iq = igivnm + 2*n*lgn
264 iwrem = iq + n**2 + 1
265* Initialize pointers
266 DO 50 i = 0, subpbs
267 iwork( iprmpt+i ) = 1
268 iwork( igivpt+i ) = 1
269 50 CONTINUE
270 iwork( iqptr ) = 1
271*
272* Solve each submatrix eigenproblem at the bottom of the divide and
273* conquer tree.
274*
275 curr = 0
276 DO 70 i = 0, spm1
277 IF( i.EQ.0 ) THEN
278 submat = 1
279 matsiz = iwork( 1 )
280 ELSE
281 submat = iwork( i ) + 1
282 matsiz = iwork( i+1 ) - iwork( i )
283 END IF
284 ll = iq - 1 + iwork( iqptr+curr )
285 CALL ssteqr( 'I', matsiz, d( submat ), e( submat ),
286 $ rwork( ll ), matsiz, rwork, info )
287 CALL clacrm( qsiz, matsiz, q( 1, submat ), ldq, rwork( ll ),
288 $ matsiz, qstore( 1, submat ), ldqs,
289 $ rwork( iwrem ) )
290 iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
291 curr = curr + 1
292 IF( info.GT.0 ) THEN
293 info = submat*( n+1 ) + submat + matsiz - 1
294 RETURN
295 END IF
296 k = 1
297 DO 60 j = submat, iwork( i+1 )
298 iwork( indxq+j ) = k
299 k = k + 1
300 60 CONTINUE
301 70 CONTINUE
302*
303* Successively merge eigensystems of adjacent submatrices
304* into eigensystem for the corresponding larger matrix.
305*
306* while ( SUBPBS > 1 )
307*
308 curlvl = 1
309 80 CONTINUE
310 IF( subpbs.GT.1 ) THEN
311 spm2 = subpbs - 2
312 DO 90 i = 0, spm2, 2
313 IF( i.EQ.0 ) THEN
314 submat = 1
315 matsiz = iwork( 2 )
316 msd2 = iwork( 1 )
317 curprb = 0
318 ELSE
319 submat = iwork( i ) + 1
320 matsiz = iwork( i+2 ) - iwork( i )
321 msd2 = matsiz / 2
322 curprb = curprb + 1
323 END IF
324*
325* Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
326* into an eigensystem of size MATSIZ. CLAED7 handles the case
327* when the eigenvectors of a full or band Hermitian matrix (which
328* was reduced to tridiagonal form) are desired.
329*
330* I am free to use Q as a valuable working space until Loop 150.
331*
332 CALL claed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,
333 $ d( submat ), qstore( 1, submat ), ldqs,
334 $ e( submat+msd2-1 ), iwork( indxq+submat ),
335 $ rwork( iq ), iwork( iqptr ), iwork( iprmpt ),
336 $ iwork( iperm ), iwork( igivpt ),
337 $ iwork( igivcl ), rwork( igivnm ),
338 $ q( 1, submat ), rwork( iwrem ),
339 $ iwork( subpbs+1 ), info )
340 IF( info.GT.0 ) THEN
341 info = submat*( n+1 ) + submat + matsiz - 1
342 RETURN
343 END IF
344 iwork( i / 2+1 ) = iwork( i+2 )
345 90 CONTINUE
346 subpbs = subpbs / 2
347 curlvl = curlvl + 1
348 GO TO 80
349 END IF
350*
351* end while
352*
353* Re-merge the eigenvalues/vectors which were deflated at the final
354* merge step.
355*
356 DO 100 i = 1, n
357 j = iwork( indxq+i )
358 rwork( i ) = d( j )
359 CALL ccopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
360 100 CONTINUE
361 CALL scopy( n, rwork, 1, d, 1 )
362*
363 RETURN
364*
365* End of CLAED0
366*
367 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
clacrm
subroutine clacrm(m, n, a, lda, b, ldb, c, ldc, rwork)
CLACRM multiplies a complex matrix by a square real matrix.
Definition clacrm.f:112
claed0
subroutine claed0(qsiz, n, d, e, q, ldq, qstore, ldqs, rwork, iwork, info)
CLAED0 used by CSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition claed0.f:143
claed7
subroutine claed7(n, cutpnt, qsiz, tlvls, curlvl, curpbm, d, q, ldq, rho, indxq, qstore, qptr, prmptr, perm, givptr, givcol, givnum, work, rwork, iwork, info)
CLAED7 used by CSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition claed7.f:248
ssteqr
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:129