LAPACK  3.8.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.
Date
December 2016
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 214 of file cstedc.f.

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