1 SUBROUTINE WFN1_AD_DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, 2 $ INFO ) 3 USE WFN1_AD1 4 IMPLICIT NONE 5#include "blas/double/intf_wfn1_ad_dscal.fh" 6#include "lapack/double/intf_wfn1_ad_dlansy.fh" 7#include "lapack/double/intf_wfn1_ad_dlascl.fh" 8#include "lapack/double/intf_wfn1_ad_dorgtr.fh" 9#include "lapack/double/intf_wfn1_ad_dsteqr.fh" 10#include "lapack/double/intf_wfn1_ad_dsterf.fh" 11#include "lapack/double/intf_wfn1_ad_dsytrd.fh" 12* 13* -- LAPACK driver routine (version 3.2) -- 14* -- LAPACK is a software package provided by Univ. of Tennessee, -- 15* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 16* November 2006 17* 18* .. Scalar Arguments .. 19 CHARACTER JOBZ, UPLO 20 INTEGER INFO, LDA, LWORK, N 21* .. 22* .. Array Arguments .. 23 TYPE(WFN1_AD_DBLE) :: A( LDA, * ), W( * ), WORK( * ) 24* .. 25* 26* Purpose 27* ======= 28* 29* DSYEV computes all eigenvalues and, optionally, eigenvectors of a 30* real symmetric matrix A. 31* 32* Arguments 33* ========= 34* 35* JOBZ (input) CHARACTER*1 36* = 'N': Compute eigenvalues only; 37* = 'V': Compute eigenvalues and eigenvectors. 38* 39* UPLO (input) CHARACTER*1 40* = 'U': Upper triangle of A is stored; 41* = 'L': Lower triangle of A is stored. 42* 43* N (input) INTEGER 44* The order of the matrix A. N >= 0. 45* 46* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 47* On entry, the symmetric matrix A. If UPLO = 'U', the 48* leading N-by-N upper triangular part of A contains the 49* upper triangular part of the matrix A. If UPLO = 'L', 50* the leading N-by-N lower triangular part of A contains 51* the lower triangular part of the matrix A. 52* On exit, if JOBZ = 'V', then if INFO = 0, A contains the 53* orthonormal eigenvectors of the matrix A. 54* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 55* or the upper triangle (if UPLO='U') of A, including the 56* diagonal, is destroyed. 57* 58* LDA (input) INTEGER 59* The leading dimension of the array A. LDA >= max(1,N). 60* 61* W (output) DOUBLE PRECISION array, dimension (N) 62* If INFO = 0, the eigenvalues in ascending order. 63* 64* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 65* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 66* 67* LWORK (input) INTEGER 68* The length of the array WORK. LWORK >= max(1,3*N-1). 69* For optimal efficiency, LWORK >= (NB+2)*N, 70* where NB is the blocksize for DSYTRD returned by ILAENV. 71* 72* If LWORK = -1, then a workspace query is assumed; the routine 73* only calculates the optimal size of the WORK array, returns 74* this value as the first entry of the WORK array, and no error 75* message related to LWORK is issued by XERBLA. 76* 77* INFO (output) INTEGER 78* = 0: successful exit 79* < 0: if INFO = -i, the i-th argument had an illegal value 80* > 0: if INFO = i, the algorithm failed to converge; i 81* off-diagonal elements of an intermediate tridiagonal 82* form did not converge to zero. 83* 84* ===================================================================== 85* 86* .. Parameters .. 87 DOUBLE PRECISION ZERO 88 TYPE(WFN1_AD_DBLE) :: ONE 89 PARAMETER ( ZERO = 0.0D0 ) 90* .. 91* .. Local Scalars .. 92 LOGICAL LOWER, LQUERY, WANTZ 93 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, 94 $ LLWORK, LWKOPT, NB 95 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SMLNUM 96 TYPE(WFN1_AD_DBLE) :: ANRM, SIGMA 97* .. 98* .. External Functions .. 99 LOGICAL LSAME 100 INTEGER ILAENV 101 DOUBLE PRECISION DLAMCH 102 EXTERNAL LSAME, ILAENV, DLAMCH 103* .. 104* .. External Subroutines .. 105 EXTERNAL XERBLA 106* .. 107* .. Intrinsic Functions .. 108c INTRINSIC MAX, SQRT 109* .. 110* .. Executable Statements .. 111* 112* Test the input parameters. 113* 114 ONE = 1.0d0 115 WANTZ = LSAME( JOBZ, 'V' ) 116 LOWER = LSAME( UPLO, 'L' ) 117 LQUERY = ( LWORK.EQ.-1 ) 118* 119 INFO = 0 120 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 121 INFO = -1 122 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 123 INFO = -2 124 ELSE IF( N.LT.0 ) THEN 125 INFO = -3 126 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 127 INFO = -5 128 END IF 129* 130 IF( INFO.EQ.0 ) THEN 131 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) 132 LWKOPT = MAX( 1, ( NB+2 )*N ) 133 WORK( 1 ) = LWKOPT*1.0d0 134* 135 IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY ) 136 $ INFO = -8 137 END IF 138* 139 IF( INFO.NE.0 ) THEN 140 CALL XERBLA( 'DSYEV ', -INFO ) 141 RETURN 142 ELSE IF( LQUERY ) THEN 143 RETURN 144 END IF 145* 146* Quick return if possible 147* 148 IF( N.EQ.0 ) THEN 149 RETURN 150 END IF 151* 152 IF( N.EQ.1 ) THEN 153 W( 1 ) = A( 1, 1 ) 154 WORK( 1 ) = 2.0d0 155 IF( WANTZ ) 156 $ A( 1, 1 ) = ONE 157 RETURN 158 END IF 159* 160* Get machine constants. 161* 162 SAFMIN = DLAMCH( 'Safe minimum' ) 163 EPS = DLAMCH( 'Precision' ) 164 SMLNUM = SAFMIN / EPS 165 BIGNUM = 1.0d0 / SMLNUM 166 RMIN = SQRT( SMLNUM ) 167 RMAX = SQRT( BIGNUM ) 168* 169* Scale matrix to allowable range, if necessary. 170* 171 ANRM = WFN1_AD_DLANSY( 'M', UPLO, N, A, LDA, WORK ) 172 ISCALE = 0 173 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 174 ISCALE = 1 175 SIGMA = RMIN / ANRM 176 ELSE IF( ANRM.GT.RMAX ) THEN 177 ISCALE = 1 178 SIGMA = RMAX / ANRM 179 END IF 180 IF( ISCALE.EQ.1 ) 181 $ CALL WFN1_AD_DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, 182 $ INFO ) 183* 184* Call DSYTRD to reduce symmetric matrix to tridiagonal form. 185* 186 INDE = 1 187 INDTAU = INDE + N 188 INDWRK = INDTAU + N 189 LLWORK = LWORK - INDWRK + 1 190 CALL WFN1_AD_DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), 191 $ WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO ) 192* 193* For eigenvalues only, call DSTERF. For eigenvectors, first call 194* DORGTR to generate the orthogonal matrix, then call DSTEQR. 195* 196 IF( .NOT.WANTZ ) THEN 197 CALL WFN1_AD_DSTERF( N, W, WORK( INDE ), INFO ) 198 ELSE 199 CALL WFN1_AD_DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), 200 $ WORK( INDWRK ), LLWORK, IINFO ) 201 CALL WFN1_AD_DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, 202 $ WORK( INDTAU ), INFO ) 203 END IF 204* 205* If matrix was scaled, then rescale eigenvalues appropriately. 206* 207 IF( ISCALE.EQ.1 ) THEN 208 IF( INFO.EQ.0 ) THEN 209 IMAX = N 210 ELSE 211 IMAX = INFO - 1 212 END IF 213 CALL WFN1_AD_DSCAL( IMAX, ONE / SIGMA, W, 1 ) 214 END IF 215* 216* Set WORK(1) to optimal workspace size. 217* 218 WORK( 1 ) = LWKOPT*1.0d0 219* 220 RETURN 221* 222* End of DSYEV 223* 224 END 225