--- rpl/lapack/lapack/dstedc.f 2010/12/21 13:53:38 1.7 +++ rpl/lapack/lapack/dstedc.f 2018/05/29 07:18:07 1.18 @@ -1,10 +1,197 @@ +*> \brief \b DSTEDC +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DSTEDC + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, +* LIWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER COMPZ +* INTEGER INFO, LDZ, LIWORK, LWORK, N +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> 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. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] COMPZ +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The dimension of the symmetric tridiagonal matrix. N >= 0. +*> \endverbatim +*> +*> \param[in,out] D +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in,out] E +*> \verbatim +*> E is DOUBLE PRECISION array, dimension (N-1) +*> On entry, the subdiagonal elements of the tridiagonal matrix. +*> On exit, E has been destroyed. +*> \endverbatim +*> +*> \param[in,out] Z +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in] LDZ +*> \verbatim +*> LDZ is INTEGER +*> The leading dimension of the array Z. LDZ >= 1. +*> If eigenvectors are desired, then LDZ >= max(1,N). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) +*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. +*> \endverbatim +*> +*> \param[in] LIWORK +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> 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). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date June 2017 +* +*> \ingroup auxOTHERcomputational +* +*> \par Contributors: +* ================== +*> +*> Jeff Rutter, Computer Science Division, University of California +*> at Berkeley, USA \n +*> Modified by Francoise Tisseur, University of Tennessee +*> +* ===================================================================== SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, $ LIWORK, INFO ) * -* -- LAPACK driver routine (version 3.2) -- +* -- LAPACK computational routine (version 3.7.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* June 2017 * * .. Scalar Arguments .. CHARACTER COMPZ @@ -15,113 +202,6 @@ DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) * .. * -* 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. -* -* Arguments -* ========= -* -* COMPZ (input) 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. -* -* N (input) INTEGER -* The dimension of the symmetric tridiagonal matrix. N >= 0. -* -* D (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the diagonal elements of the tridiagonal matrix. -* On exit, if INFO = 0, the eigenvalues in ascending order. -* -* E (input/output) DOUBLE PRECISION array, dimension (N-1) -* On entry, the subdiagonal elements of the tridiagonal matrix. -* On exit, E has been destroyed. -* -* Z (input/output) 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. -* -* LDZ (input) INTEGER -* The leading dimension of the array Z. LDZ >= 1. -* If eigenvectors are desired, then LDZ >= max(1,N). -* -* WORK (workspace/output) DOUBLE PRECISION array, -* dimension (LWORK) -* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. -* -* LWORK (input) 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 + 3*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. -* -* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) -* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. -* -* LIWORK (input) 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. -* -* INFO (output) 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). -* -* Further Details -* =============== -* -* Based on contributions by -* Jeff Rutter, Computer Science Division, University of California -* at Berkeley, USA -* Modified by Francoise Tisseur, University of Tennessee. -* * ===================================================================== * * .. Parameters .. @@ -190,7 +270,7 @@ IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IF( ICOMPZ.EQ.1 ) THEN - LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2 + LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2 LIWMIN = 6 + 6*N + 5*N*LGN ELSE IF( ICOMPZ.EQ.2 ) THEN LWMIN = 1 + 4*N + N**2 @@ -362,38 +442,32 @@ * * endwhile * -* If the problem split any number of times, then the eigenvalues -* will not be properly ordered. Here we permute the eigenvalues -* (and the associated eigenvectors) into ascending order. -* - IF( M.NE.N ) THEN - IF( ICOMPZ.EQ.0 ) THEN + IF( ICOMPZ.EQ.0 ) THEN * -* Use Quick Sort +* Use Quick Sort * - CALL DLASRT( 'I', N, D, INFO ) + CALL DLASRT( 'I', N, D, INFO ) * - ELSE + ELSE * -* Use Selection Sort to minimize swaps of eigenvectors +* Use Selection Sort to minimize swaps of eigenvectors * - DO 40 II = 2, N - I = II - 1 - K = I - P = D( I ) - DO 30 J = II, N - IF( D( J ).LT.P ) THEN - K = J - P = D( J ) - END IF - 30 CONTINUE - IF( K.NE.I ) THEN - D( K ) = D( I ) - D( I ) = P - CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) - END IF - 40 CONTINUE - END IF + DO 40 II = 2, N + I = II - 1 + K = I + P = D( I ) + DO 30 J = II, N + IF( D( J ).LT.P ) THEN + K = J + P = D( J ) + END IF + 30 CONTINUE + IF( K.NE.I ) THEN + D( K ) = D( I ) + D( I ) = P + CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) + END IF + 40 CONTINUE END IF END IF *