LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.

 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 DLAED3 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 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.
Date
June 2017
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 190 of file dstedc.f.

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