LAPACK 3.12.1
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 202 of file cstedc.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 REAL D( * ), E( * ), RWORK( * )
216 COMPLEX WORK( * ), Z( LDZ, * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 REAL ZERO, ONE, TWO
223 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
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 REAL EPS, ORGNRM, P, TINY
230* ..
231* .. External Functions ..
232 LOGICAL LSAME
233 INTEGER ILAENV
234 REAL SLAMCH, SLANST, SROUNDUP_LWORK
235 EXTERNAL ilaenv, lsame, slamch, slanst,
236 $ sroundup_lwork
237* ..
238* .. External Subroutines ..
239 EXTERNAL xerbla, clacpy, clacrm, claed0, csteqr,
240 $ cswap,
241 $ slascl, slaset, sstedc, ssteqr, ssterf
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 ) = real( 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 ),
412 $ m,
413 $ info )
414 CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1,
415 $ e( start ),
416 $ m-1, info )
417*
418 CALL claed0( n, m, d( start ), e( start ), z( 1,
419 $ start ),
420 $ ldz, work, n, rwork, iwork, info )
421 IF( info.GT.0 ) THEN
422 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
423 $ mod( info, ( m+1 ) ) + start - 1
424 GO TO 70
425 END IF
426*
427* Scale back.
428*
429 CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ),
430 $ m,
431 $ info )
432*
433 ELSE
434 CALL ssteqr( 'I', m, d( start ), e( start ), rwork, m,
435 $ rwork( m*m+1 ), info )
436 CALL clacrm( n, m, z( 1, start ), ldz, rwork, m, work,
437 $ n,
438 $ rwork( m*m+1 ) )
439 CALL clacpy( 'A', n, m, work, n, z( 1, start ), ldz )
440 IF( info.GT.0 ) THEN
441 info = start*( n+1 ) + finish
442 GO TO 70
443 END IF
444 END IF
445*
446 start = finish + 1
447 GO TO 30
448 END IF
449*
450* endwhile
451*
452*
453* Use Selection Sort to minimize swaps of eigenvectors
454*
455 DO 60 ii = 2, n
456 i = ii - 1
457 k = i
458 p = d( i )
459 DO 50 j = ii, n
460 IF( d( j ).LT.p ) THEN
461 k = j
462 p = d( j )
463 END IF
464 50 CONTINUE
465 IF( k.NE.i ) THEN
466 d( k ) = d( i )
467 d( i ) = p
468 CALL cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
469 END IF
470 60 CONTINUE
471 END IF
472*
473 70 CONTINUE
474 work( 1 ) = sroundup_lwork(lwmin)
475 rwork( 1 ) = real( lrwmin )
476 iwork( 1 ) = liwmin
477*
478 RETURN
479*
480* End of CSTEDC
481*
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
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
clacrm
subroutine clacrm(m, n, a, lda, b, ldb, c, ldc, rwork)
CLACRM multiplies a complex matrix by a square real matrix.
Definition clacrm.f:112
claed0
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:143
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slanst
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:98
slascl
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:142
slaset
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:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
sstedc
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:180
ssteqr
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:129
csteqr
subroutine csteqr(compz, n, d, e, z, ldz, work, info)
CSTEQR
Definition csteqr.f:130
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:84
cswap
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: