LAPACK 3.11.0
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.

 This code makes very mild assumptions about floating point
 arithmetic. It will work on machines with a guard digit in
 add/subtract, or on those binary machines without guard digits
 which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.  See DLAED3 for details.
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 210 of file zstedc.f.

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