LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zstedc()

subroutine zstedc ( character compz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
complex*16, dimension( ldz, * ) z,
integer ldz,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer lrwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

ZSTEDC

Download ZSTEDC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZSTEDC computes all eigenvalues and, optionally, eigenvectors of a
!> symmetric tridiagonal matrix using the divide and conquer method.
!> The eigenvectors of a full or band complex Hermitian matrix can also
!> be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this
!> matrix to tridiagonal form.
!>
!>
Parameters
[in]COMPZ
!>          COMPZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only.
!>          = 'I':  Compute eigenvectors of tridiagonal matrix also.
!>          = 'V':  Compute eigenvectors of original Hermitian matrix
!>                  also.  On entry, Z contains the unitary matrix used
!>                  to reduce the original matrix to tridiagonal form.
!>
[in]N
!>          N is INTEGER
!>          The dimension of the symmetric tridiagonal matrix.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          On entry, the diagonal elements of the tridiagonal matrix.
!>          On exit, if INFO = 0, the eigenvalues in ascending order.
!>
[in,out]E
!>          E is DOUBLE PRECISION array, dimension (N-1)
!>          On entry, the subdiagonal elements of the tridiagonal matrix.
!>          On exit, E has been destroyed.
!>
[in,out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ,N)
!>          On entry, if COMPZ = 'V', then Z contains the unitary
!>          matrix used in the reduction to tridiagonal form.
!>          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
!>          orthonormal eigenvectors of the original Hermitian matrix,
!>          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
!>          of the symmetric tridiagonal matrix.
!>          If  COMPZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1.
!>          If eigenvectors are desired, then LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1.
!>          If COMPZ = 'V' and N > 1, LWORK must be at least N*N.
!>          Note that for COMPZ = 'V', then if N is less than or
!>          equal to the minimum divide size, usually 25, then LWORK need
!>          only be 1.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK, RWORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK.
!>          If COMPZ = 'N' or N <= 1, LRWORK must be at least 1.
!>          If COMPZ = 'V' and N > 1, LRWORK must be at least
!>                         1 + 3*N + 2*N*lg N + 4*N**2 ,
!>                         where lg( N ) = smallest integer k such
!>                         that 2**k >= N.
!>          If COMPZ = 'I' and N > 1, LRWORK must be at least
!>                         1 + 4*N + 2*N**2 .
!>          Note that for COMPZ = 'I' or 'V', then if N is less than or
!>          equal to the minimum divide size, usually 25, then LRWORK
!>          need only be max(1,2*(N-1)).
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If COMPZ = 'N' or N <= 1, LIWORK must be at least 1.
!>          If COMPZ = 'V' or N > 1,  LIWORK must be at least
!>                                    6 + 6*N + 5*N*lg N.
!>          If COMPZ = 'I' or N > 1,  LIWORK must be at least
!>                                    3 + 5*N .
!>          Note that for COMPZ = 'I' or 'V', then if N is less than or
!>          equal to the minimum divide size, usually 25, then LIWORK
!>          need only be 1.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK, RWORK
!>          and IWORK arrays, returns these values as the first entries
!>          of the WORK, RWORK and IWORK arrays, and no error message
!>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  The algorithm failed to compute an eigenvalue while
!>                working on the submatrix lying in rows and columns
!>                INFO/(N+1) through mod(INFO,N+1).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 202 of file zstedc.f.

204*
205* -- LAPACK computational routine --
206* -- LAPACK is a software package provided by Univ. of Tennessee, --
207* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
208*
209* .. Scalar Arguments ..
210 CHARACTER COMPZ
211 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
212* ..
213* .. Array Arguments ..
214 INTEGER IWORK( * )
215 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
216 COMPLEX*16 WORK( * ), Z( LDZ, * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 DOUBLE PRECISION ZERO, ONE, TWO
223 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
224* ..
225* .. Local Scalars ..
226 LOGICAL LQUERY
227 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL,
228 $ LRWMIN, LWMIN, M, SMLSIZ, START
229 DOUBLE PRECISION EPS, ORGNRM, P, TINY
230* ..
231* .. External Functions ..
232 LOGICAL LSAME
233 INTEGER ILAENV
234 DOUBLE PRECISION DLAMCH, DLANST
235 EXTERNAL lsame, ilaenv, dlamch, dlanst
236* ..
237* .. External Subroutines ..
238 EXTERNAL dlascl, dlaset, dstedc, dsteqr, dsterf,
239 $ xerbla,
240 $ zlacpy, zlacrm, zlaed0, zsteqr, zswap
241* ..
242* .. Intrinsic Functions ..
243 INTRINSIC abs, dble, int, log, max, mod, sqrt
244* ..
245* .. Executable Statements ..
246*
247* Test the input parameters.
248*
249 info = 0
250 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
251*
252 IF( lsame( compz, 'N' ) ) THEN
253 icompz = 0
254 ELSE IF( lsame( compz, 'V' ) ) THEN
255 icompz = 1
256 ELSE IF( lsame( compz, 'I' ) ) THEN
257 icompz = 2
258 ELSE
259 icompz = -1
260 END IF
261 IF( icompz.LT.0 ) THEN
262 info = -1
263 ELSE IF( n.LT.0 ) THEN
264 info = -2
265 ELSE IF( ( ldz.LT.1 ) .OR.
266 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
267 info = -6
268 END IF
269*
270 IF( info.EQ.0 ) THEN
271*
272* Compute the workspace requirements
273*
274 smlsiz = ilaenv( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 )
275 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
276 lwmin = 1
277 liwmin = 1
278 lrwmin = 1
279 ELSE IF( n.LE.smlsiz ) THEN
280 lwmin = 1
281 liwmin = 1
282 lrwmin = 2*( n - 1 )
283 ELSE IF( icompz.EQ.1 ) THEN
284 lgn = int( log( dble( n ) ) / log( two ) )
285 IF( 2**lgn.LT.n )
286 $ lgn = lgn + 1
287 IF( 2**lgn.LT.n )
288 $ lgn = lgn + 1
289 lwmin = n*n
290 lrwmin = 1 + 3*n + 2*n*lgn + 4*n**2
291 liwmin = 6 + 6*n + 5*n*lgn
292 ELSE IF( icompz.EQ.2 ) THEN
293 lwmin = 1
294 lrwmin = 1 + 4*n + 2*n**2
295 liwmin = 3 + 5*n
296 END IF
297 work( 1 ) = lwmin
298 rwork( 1 ) = real( lrwmin )
299 iwork( 1 ) = liwmin
300*
301 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
302 info = -8
303 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
304 info = -10
305 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
306 info = -12
307 END IF
308 END IF
309*
310 IF( info.NE.0 ) THEN
311 CALL xerbla( 'ZSTEDC', -info )
312 RETURN
313 ELSE IF( lquery ) THEN
314 RETURN
315 END IF
316*
317* Quick return if possible
318*
319 IF( n.EQ.0 )
320 $ RETURN
321 IF( n.EQ.1 ) THEN
322 IF( icompz.NE.0 )
323 $ z( 1, 1 ) = one
324 RETURN
325 END IF
326*
327* If the following conditional clause is removed, then the routine
328* will use the Divide and Conquer routine to compute only the
329* eigenvalues, which requires (3N + 3N**2) real workspace and
330* (2 + 5N + 2N lg(N)) integer workspace.
331* Since on many architectures DSTERF is much faster than any other
332* algorithm for finding eigenvalues only, it is used here
333* as the default. If the conditional clause is removed, then
334* information on the size of workspace needs to be changed.
335*
336* If COMPZ = 'N', use DSTERF to compute the eigenvalues.
337*
338 IF( icompz.EQ.0 ) THEN
339 CALL dsterf( n, d, e, info )
340 GO TO 70
341 END IF
342*
343* If N is smaller than the minimum divide size (SMLSIZ+1), then
344* solve the problem with another solver.
345*
346 IF( n.LE.smlsiz ) THEN
347*
348 CALL zsteqr( compz, n, d, e, z, ldz, rwork, info )
349*
350 ELSE
351*
352* If COMPZ = 'I', we simply call DSTEDC instead.
353*
354 IF( icompz.EQ.2 ) THEN
355 CALL dlaset( 'Full', n, n, zero, one, rwork, n )
356 ll = n*n + 1
357 CALL dstedc( 'I', n, d, e, rwork, n,
358 $ rwork( ll ), lrwork-ll+1, iwork, liwork, info )
359 DO 20 j = 1, n
360 DO 10 i = 1, n
361 z( i, j ) = rwork( ( j-1 )*n+i )
362 10 CONTINUE
363 20 CONTINUE
364 GO TO 70
365 END IF
366*
367* From now on, only option left to be handled is COMPZ = 'V',
368* i.e. ICOMPZ = 1.
369*
370* Scale.
371*
372 orgnrm = dlanst( 'M', n, d, e )
373 IF( orgnrm.EQ.zero )
374 $ GO TO 70
375*
376 eps = dlamch( 'Epsilon' )
377*
378 start = 1
379*
380* while ( START <= N )
381*
382 30 CONTINUE
383 IF( start.LE.n ) THEN
384*
385* Let FINISH be the position of the next subdiagonal entry
386* such that E( FINISH ) <= TINY or FINISH = N if no such
387* subdiagonal exists. The matrix identified by the elements
388* between START and FINISH constitutes an independent
389* sub-problem.
390*
391 finish = start
392 40 CONTINUE
393 IF( finish.LT.n ) THEN
394 tiny = eps*sqrt( abs( d( finish ) ) )*
395 $ sqrt( abs( d( finish+1 ) ) )
396 IF( abs( e( finish ) ).GT.tiny ) THEN
397 finish = finish + 1
398 GO TO 40
399 END IF
400 END IF
401*
402* (Sub) Problem determined. Compute its size and solve it.
403*
404 m = finish - start + 1
405 IF( m.GT.smlsiz ) THEN
406*
407* Scale.
408*
409 orgnrm = dlanst( 'M', m, d( start ), e( start ) )
410 CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ),
411 $ m,
412 $ info )
413 CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1,
414 $ e( start ),
415 $ m-1, info )
416*
417 CALL zlaed0( n, m, d( start ), e( start ), z( 1,
418 $ start ),
419 $ ldz, work, n, rwork, iwork, info )
420 IF( info.GT.0 ) THEN
421 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
422 $ mod( info, ( m+1 ) ) + start - 1
423 GO TO 70
424 END IF
425*
426* Scale back.
427*
428 CALL dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ),
429 $ m,
430 $ info )
431*
432 ELSE
433 CALL dsteqr( 'I', m, d( start ), e( start ), rwork, m,
434 $ rwork( m*m+1 ), info )
435 CALL zlacrm( n, m, z( 1, start ), ldz, rwork, m, work,
436 $ n,
437 $ rwork( m*m+1 ) )
438 CALL zlacpy( 'A', n, m, work, n, z( 1, start ), ldz )
439 IF( info.GT.0 ) THEN
440 info = start*( n+1 ) + finish
441 GO TO 70
442 END IF
443 END IF
444*
445 start = finish + 1
446 GO TO 30
447 END IF
448*
449* endwhile
450*
451*
452* Use Selection Sort to minimize swaps of eigenvectors
453*
454 DO 60 ii = 2, n
455 i = ii - 1
456 k = i
457 p = d( i )
458 DO 50 j = ii, n
459 IF( d( j ).LT.p ) THEN
460 k = j
461 p = d( j )
462 END IF
463 50 CONTINUE
464 IF( k.NE.i ) THEN
465 d( k ) = d( i )
466 d( i ) = p
467 CALL zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
468 END IF
469 60 CONTINUE
470 END IF
471*
472 70 CONTINUE
473 work( 1 ) = lwmin
474 rwork( 1 ) = real( lrwmin )
475 iwork( 1 ) = liwmin
476*
477 RETURN
478*
479* End of ZSTEDC
480*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
zlacrm
subroutine zlacrm(m, n, a, lda, b, ldb, c, ldc, rwork)
ZLACRM multiplies a complex matrix by a square real matrix.
Definition zlacrm.f:112
zlaed0
subroutine zlaed0(qsiz, n, d, e, q, ldq, qstore, ldqs, rwork, iwork, info)
ZLAED0 used by ZSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition zlaed0.f:143
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
double precision function dlanst(norm, n, d, e)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlanst.f:98
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:142
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dstedc
subroutine dstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
DSTEDC
Definition dstedc.f:180
zsteqr
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:130
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:129
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:84
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: