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

◆ dstedc()

subroutine dstedc ( character compz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

DSTEDC

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

Purpose:
!>
!> DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
!> symmetric tridiagonal matrix using the divide and conquer method.
!> The eigenvectors of a full or band real symmetric matrix can also be
!> found if DSYTRD or DSPTRD or DSBTRD 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 dense symmetric
!>                  matrix also.  On entry, Z contains the orthogonal
!>                  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 DOUBLE PRECISION array, dimension (LDZ,N)
!>          On entry, if COMPZ = 'V', then Z contains the orthogonal
!>          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 symmetric 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 DOUBLE PRECISION 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 N <= 1 then LWORK must be at least 1.
!>          If COMPZ = 'V' and N > 1 then LWORK 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 then LWORK must be at least
!>                         ( 1 + 4*N + 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 LWORK need
!>          only be max(1,2*(N-1)).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK 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 then LIWORK must be at least 1.
!>          If COMPZ = 'V' and N > 1 then LIWORK must be at least
!>                         ( 6 + 6*N + 5*N*lg N ).
!>          If COMPZ = 'I' and N > 1 then 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 size of the IWORK array,
!>          returns this value as the first entry of the IWORK array, and
!>          no error message related to 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
Modified by Francoise Tisseur, University of Tennessee

Definition at line 178 of file dstedc.f.

180*
181* -- LAPACK computational routine --
182* -- LAPACK is a software package provided by Univ. of Tennessee, --
183* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184*
185* .. Scalar Arguments ..
186 CHARACTER COMPZ
187 INTEGER INFO, LDZ, LIWORK, LWORK, N
188* ..
189* .. Array Arguments ..
190 INTEGER IWORK( * )
191 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
192* ..
193*
194* =====================================================================
195*
196* .. Parameters ..
197 DOUBLE PRECISION ZERO, ONE, TWO
198 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
199* ..
200* .. Local Scalars ..
201 LOGICAL LQUERY
202 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
203 $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
204 DOUBLE PRECISION EPS, ORGNRM, P, TINY
205* ..
206* .. External Functions ..
207 LOGICAL LSAME
208 INTEGER ILAENV
209 DOUBLE PRECISION DLAMCH, DLANST
210 EXTERNAL lsame, ilaenv, dlamch, dlanst
211* ..
212* .. External Subroutines ..
213 EXTERNAL dgemm, dlacpy, dlaed0, dlascl, dlaset,
214 $ dlasrt,
215 $ dsteqr, dsterf, dswap, xerbla
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC abs, dble, int, log, max, mod, sqrt
219* ..
220* .. Executable Statements ..
221*
222* Test the input parameters.
223*
224 info = 0
225 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
226*
227 IF( lsame( compz, 'N' ) ) THEN
228 icompz = 0
229 ELSE IF( lsame( compz, 'V' ) ) THEN
230 icompz = 1
231 ELSE IF( lsame( compz, 'I' ) ) THEN
232 icompz = 2
233 ELSE
234 icompz = -1
235 END IF
236 IF( icompz.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( ( ldz.LT.1 ) .OR.
241 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
242 info = -6
243 END IF
244*
245 IF( info.EQ.0 ) THEN
246*
247* Compute the workspace requirements
248*
249 smlsiz = ilaenv( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
250 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
251 liwmin = 1
252 lwmin = 1
253 ELSE IF( n.LE.smlsiz ) THEN
254 liwmin = 1
255 lwmin = 2*( n - 1 )
256 ELSE
257 lgn = int( log( dble( n ) )/log( two ) )
258 IF( 2**lgn.LT.n )
259 $ lgn = lgn + 1
260 IF( 2**lgn.LT.n )
261 $ lgn = lgn + 1
262 IF( icompz.EQ.1 ) THEN
263 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
264 liwmin = 6 + 6*n + 5*n*lgn
265 ELSE IF( icompz.EQ.2 ) THEN
266 lwmin = 1 + 4*n + n**2
267 liwmin = 3 + 5*n
268 END IF
269 END IF
270 work( 1 ) = lwmin
271 iwork( 1 ) = liwmin
272*
273 IF( lwork.LT.lwmin .AND. .NOT. lquery ) THEN
274 info = -8
275 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery ) THEN
276 info = -10
277 END IF
278 END IF
279*
280 IF( info.NE.0 ) THEN
281 CALL xerbla( 'DSTEDC', -info )
282 RETURN
283 ELSE IF (lquery) THEN
284 RETURN
285 END IF
286*
287* Quick return if possible
288*
289 IF( n.EQ.0 )
290 $ RETURN
291 IF( n.EQ.1 ) THEN
292 IF( icompz.NE.0 )
293 $ z( 1, 1 ) = one
294 RETURN
295 END IF
296*
297* If the following conditional clause is removed, then the routine
298* will use the Divide and Conquer routine to compute only the
299* eigenvalues, which requires (3N + 3N**2) real workspace and
300* (2 + 5N + 2N lg(N)) integer workspace.
301* Since on many architectures DSTERF is much faster than any other
302* algorithm for finding eigenvalues only, it is used here
303* as the default. If the conditional clause is removed, then
304* information on the size of workspace needs to be changed.
305*
306* If COMPZ = 'N', use DSTERF to compute the eigenvalues.
307*
308 IF( icompz.EQ.0 ) THEN
309 CALL dsterf( n, d, e, info )
310 GO TO 50
311 END IF
312*
313* If N is smaller than the minimum divide size (SMLSIZ+1), then
314* solve the problem with another solver.
315*
316 IF( n.LE.smlsiz ) THEN
317*
318 CALL dsteqr( compz, n, d, e, z, ldz, work, info )
319*
320 ELSE
321*
322* If COMPZ = 'V', the Z matrix must be stored elsewhere for later
323* use.
324*
325 IF( icompz.EQ.1 ) THEN
326 storez = 1 + n*n
327 ELSE
328 storez = 1
329 END IF
330*
331 IF( icompz.EQ.2 ) THEN
332 CALL dlaset( 'Full', n, n, zero, one, z, ldz )
333 END IF
334*
335* Scale.
336*
337 orgnrm = dlanst( 'M', n, d, e )
338 IF( orgnrm.EQ.zero )
339 $ GO TO 50
340*
341 eps = dlamch( 'Epsilon' )
342*
343 start = 1
344*
345* while ( START <= N )
346*
347 10 CONTINUE
348 IF( start.LE.n ) THEN
349*
350* Let FINISH be the position of the next subdiagonal entry
351* such that E( FINISH ) <= TINY or FINISH = N if no such
352* subdiagonal exists. The matrix identified by the elements
353* between START and FINISH constitutes an independent
354* sub-problem.
355*
356 finish = start
357 20 CONTINUE
358 IF( finish.LT.n ) THEN
359 tiny = eps*sqrt( abs( d( finish ) ) )*
360 $ sqrt( abs( d( finish+1 ) ) )
361 IF( abs( e( finish ) ).GT.tiny ) THEN
362 finish = finish + 1
363 GO TO 20
364 END IF
365 END IF
366*
367* (Sub) Problem determined. Compute its size and solve it.
368*
369 m = finish - start + 1
370 IF( m.EQ.1 ) THEN
371 start = finish + 1
372 GO TO 10
373 END IF
374 IF( m.GT.smlsiz ) THEN
375*
376* Scale.
377*
378 orgnrm = dlanst( 'M', m, d( start ), e( start ) )
379 CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ),
380 $ m,
381 $ info )
382 CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1,
383 $ e( start ),
384 $ m-1, info )
385*
386 IF( icompz.EQ.1 ) THEN
387 strtrw = 1
388 ELSE
389 strtrw = start
390 END IF
391 CALL dlaed0( icompz, n, m, d( start ), e( start ),
392 $ z( strtrw, start ), ldz, work( 1 ), n,
393 $ work( storez ), iwork, info )
394 IF( info.NE.0 ) THEN
395 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
396 $ mod( info, ( m+1 ) ) + start - 1
397 GO TO 50
398 END IF
399*
400* Scale back.
401*
402 CALL dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ),
403 $ m,
404 $ info )
405*
406 ELSE
407 IF( icompz.EQ.1 ) THEN
408*
409* Since QR won't update a Z matrix which is larger than
410* the length of D, we must solve the sub-problem in a
411* workspace and then multiply back into Z.
412*
413 CALL dsteqr( 'I', m, d( start ), e( start ), work,
414 $ m,
415 $ work( m*m+1 ), info )
416 CALL dlacpy( 'A', n, m, z( 1, start ), ldz,
417 $ work( storez ), n )
418 CALL dgemm( 'N', 'N', n, m, m, one,
419 $ work( storez ), n, work, m, zero,
420 $ z( 1, start ), ldz )
421 ELSE IF( icompz.EQ.2 ) THEN
422 CALL dsteqr( 'I', m, d( start ), e( start ),
423 $ z( start, start ), ldz, work, info )
424 ELSE
425 CALL dsterf( m, d( start ), e( start ), info )
426 END IF
427 IF( info.NE.0 ) THEN
428 info = start*( n+1 ) + finish
429 GO TO 50
430 END IF
431 END IF
432*
433 start = finish + 1
434 GO TO 10
435 END IF
436*
437* endwhile
438*
439 IF( icompz.EQ.0 ) THEN
440*
441* Use Quick Sort
442*
443 CALL dlasrt( 'I', n, d, info )
444*
445 ELSE
446*
447* Use Selection Sort to minimize swaps of eigenvectors
448*
449 DO 40 ii = 2, n
450 i = ii - 1
451 k = i
452 p = d( i )
453 DO 30 j = ii, n
454 IF( d( j ).LT.p ) THEN
455 k = j
456 p = d( j )
457 END IF
458 30 CONTINUE
459 IF( k.NE.i ) THEN
460 d( k ) = d( i )
461 d( i ) = p
462 CALL dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
463 END IF
464 40 CONTINUE
465 END IF
466 END IF
467*
468 50 CONTINUE
469 work( 1 ) = lwmin
470 iwork( 1 ) = liwmin
471*
472 RETURN
473*
474* End of DSTEDC
475*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlaed0
subroutine dlaed0(icompq, qsiz, n, d, e, q, ldq, qstore, ldqs, work, iwork, info)
DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition dlaed0.f:170
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
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:98
dlascl
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:142
dlaset
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:108
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:86
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:129
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:84
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: