 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CHPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, RWORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, UPLO 00011 INTEGER INFO, LDZ, N 00012 * .. 00013 * .. Array Arguments .. 00014 REAL RWORK( * ), W( * ) 00015 COMPLEX AP( * ), WORK( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CHPEV computes all the eigenvalues and, optionally, eigenvectors of a 00022 * complex Hermitian matrix in packed storage. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * JOBZ (input) CHARACTER*1 00028 * = 'N': Compute eigenvalues only; 00029 * = 'V': Compute eigenvalues and eigenvectors. 00030 * 00031 * UPLO (input) CHARACTER*1 00032 * = 'U': Upper triangle of A is stored; 00033 * = 'L': Lower triangle of A is stored. 00034 * 00035 * N (input) INTEGER 00036 * The order of the matrix A. N >= 0. 00037 * 00038 * AP (input/output) COMPLEX array, dimension (N*(N+1)/2) 00039 * On entry, the upper or lower triangle of the Hermitian matrix 00040 * A, packed columnwise in a linear array. The j-th column of A 00041 * is stored in the array AP as follows: 00042 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00043 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00044 * 00045 * On exit, AP is overwritten by values generated during the 00046 * reduction to tridiagonal form. If UPLO = 'U', the diagonal 00047 * and first superdiagonal of the tridiagonal matrix T overwrite 00048 * the corresponding elements of A, and if UPLO = 'L', the 00049 * diagonal and first subdiagonal of T overwrite the 00050 * corresponding elements of A. 00051 * 00052 * W (output) REAL array, dimension (N) 00053 * If INFO = 0, the eigenvalues in ascending order. 00054 * 00055 * Z (output) COMPLEX array, dimension (LDZ, N) 00056 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 00057 * eigenvectors of the matrix A, with the i-th column of Z 00058 * holding the eigenvector associated with W(i). 00059 * If JOBZ = 'N', then Z is not referenced. 00060 * 00061 * LDZ (input) INTEGER 00062 * The leading dimension of the array Z. LDZ >= 1, and if 00063 * JOBZ = 'V', LDZ >= max(1,N). 00064 * 00065 * WORK (workspace) COMPLEX array, dimension (max(1, 2*N-1)) 00066 * 00067 * RWORK (workspace) REAL array, dimension (max(1, 3*N-2)) 00068 * 00069 * INFO (output) INTEGER 00070 * = 0: successful exit. 00071 * < 0: if INFO = -i, the i-th argument had an illegal value. 00072 * > 0: if INFO = i, the algorithm failed to converge; i 00073 * off-diagonal elements of an intermediate tridiagonal 00074 * form did not converge to zero. 00075 * 00076 * ===================================================================== 00077 * 00078 * .. Parameters .. 00079 REAL ZERO, ONE 00080 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00081 * .. 00082 * .. Local Scalars .. 00083 LOGICAL WANTZ 00084 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK, 00085 $ ISCALE 00086 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00087 $ SMLNUM 00088 * .. 00089 * .. External Functions .. 00090 LOGICAL LSAME 00091 REAL CLANHP, SLAMCH 00092 EXTERNAL LSAME, CLANHP, SLAMCH 00093 * .. 00094 * .. External Subroutines .. 00095 EXTERNAL CHPTRD, CSSCAL, CSTEQR, CUPGTR, SSCAL, SSTERF, 00096 $ XERBLA 00097 * .. 00098 * .. Intrinsic Functions .. 00099 INTRINSIC SQRT 00100 * .. 00101 * .. Executable Statements .. 00102 * 00103 * Test the input parameters. 00104 * 00105 WANTZ = LSAME( JOBZ, 'V' ) 00106 * 00107 INFO = 0 00108 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00109 INFO = -1 00110 ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) ) 00111 $ THEN 00112 INFO = -2 00113 ELSE IF( N.LT.0 ) THEN 00114 INFO = -3 00115 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00116 INFO = -7 00117 END IF 00118 * 00119 IF( INFO.NE.0 ) THEN 00120 CALL XERBLA( 'CHPEV ', -INFO ) 00121 RETURN 00122 END IF 00123 * 00124 * Quick return if possible 00125 * 00126 IF( N.EQ.0 ) 00127 $ RETURN 00128 * 00129 IF( N.EQ.1 ) THEN 00130 W( 1 ) = AP( 1 ) 00131 RWORK( 1 ) = 1 00132 IF( WANTZ ) 00133 $ Z( 1, 1 ) = ONE 00134 RETURN 00135 END IF 00136 * 00137 * Get machine constants. 00138 * 00139 SAFMIN = SLAMCH( 'Safe minimum' ) 00140 EPS = SLAMCH( 'Precision' ) 00141 SMLNUM = SAFMIN / EPS 00142 BIGNUM = ONE / SMLNUM 00143 RMIN = SQRT( SMLNUM ) 00144 RMAX = SQRT( BIGNUM ) 00145 * 00146 * Scale matrix to allowable range, if necessary. 00147 * 00148 ANRM = CLANHP( 'M', UPLO, N, AP, RWORK ) 00149 ISCALE = 0 00150 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00151 ISCALE = 1 00152 SIGMA = RMIN / ANRM 00153 ELSE IF( ANRM.GT.RMAX ) THEN 00154 ISCALE = 1 00155 SIGMA = RMAX / ANRM 00156 END IF 00157 IF( ISCALE.EQ.1 ) THEN 00158 CALL CSSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 00159 END IF 00160 * 00161 * Call CHPTRD to reduce Hermitian packed matrix to tridiagonal form. 00162 * 00163 INDE = 1 00164 INDTAU = 1 00165 CALL CHPTRD( UPLO, N, AP, W, RWORK( INDE ), WORK( INDTAU ), 00166 $ IINFO ) 00167 * 00168 * For eigenvalues only, call SSTERF. For eigenvectors, first call 00169 * CUPGTR to generate the orthogonal matrix, then call CSTEQR. 00170 * 00171 IF( .NOT.WANTZ ) THEN 00172 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 00173 ELSE 00174 INDWRK = INDTAU + N 00175 CALL CUPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ, 00176 $ WORK( INDWRK ), IINFO ) 00177 INDRWK = INDE + N 00178 CALL CSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ, 00179 $ RWORK( INDRWK ), INFO ) 00180 END IF 00181 * 00182 * If matrix was scaled, then rescale eigenvalues appropriately. 00183 * 00184 IF( ISCALE.EQ.1 ) THEN 00185 IF( INFO.EQ.0 ) THEN 00186 IMAX = N 00187 ELSE 00188 IMAX = INFO - 1 00189 END IF 00190 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00191 END IF 00192 * 00193 RETURN 00194 * 00195 * End of CHPEV 00196 * 00197 END
1.7.3 
