LAPACK 3.12.0
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 180 of file dstedc.f.

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