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

◆ 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.
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 204 of file cstedc.f.

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