 LAPACK  3.10.1 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).```
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 186 of file dstedc.f.

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