LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 (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
November 2015
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 191 of file dstedc.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: