LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cstedc()

subroutine cstedc ( character  COMPZ,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
complex, dimension( ldz, * )  Z,
integer  LDZ,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer  LRWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

CSTEDC

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

Purpose:
 CSTEDC 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 CHETRD or CHPTRD or CHBTRD 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 SLAED3 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 REAL 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 REAL 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 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 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 REAL 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 cstedc.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  REAL D( * ), E( * ), RWORK( * )
224  COMPLEX WORK( * ), Z( LDZ, * )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  REAL ZERO, ONE, TWO
231  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
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  REAL EPS, ORGNRM, P, TINY
238 * ..
239 * .. External Functions ..
240  LOGICAL LSAME
241  INTEGER ILAENV
242  REAL SLAMCH, SLANST
243  EXTERNAL ilaenv, lsame, slamch, slanst
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL xerbla, clacpy, clacrm, claed0, csteqr, cswap,
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, int, log, max, mod, real, 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, 'CSTEDC', ' ', 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( real( 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( 'CSTEDC', -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 SSTERF 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 SSTERF to compute the eigenvalues.
344 *
345  IF( icompz.EQ.0 ) THEN
346  CALL ssterf( 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 csteqr( compz, n, d, e, z, ldz, rwork, info )
356 *
357  ELSE
358 *
359 * If COMPZ = 'I', we simply call SSTEDC instead.
360 *
361  IF( icompz.EQ.2 ) THEN
362  CALL slaset( 'Full', n, n, zero, one, rwork, n )
363  ll = n*n + 1
364  CALL sstedc( '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 = slanst( 'M', n, d, e )
380  IF( orgnrm.EQ.zero )
381  $ GO TO 70
382 *
383  eps = slamch( '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 = slanst( 'M', m, d( start ), e( start ) )
417  CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
418  $ info )
419  CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
420  $ m-1, info )
421 *
422  CALL claed0( 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 slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
433  $ info )
434 *
435  ELSE
436  CALL ssteqr( 'I', m, d( start ), e( start ), rwork, m,
437  $ rwork( m*m+1 ), info )
438  CALL clacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
439  $ rwork( m*m+1 ) )
440  CALL clacpy( '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 cswap( 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 CSTEDC
482 *
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slanst.f:100
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.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 ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:188
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
Definition: clacrm.f:114
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
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:145
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: