1 SUBROUTINE ZHBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
2 $ RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, KD, LDAB, LDZ, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * ), W( * )
15 COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZHBEV computes all the eigenvalues and, optionally, eigenvectors of
22 * a complex Hermitian band matrix A.
23 *
24 * Arguments
25 * =========
26 *
27 * JOBZ (input) CHARACTER*1
28 * = 'N': Compute eigenvalues only;
29 * = 'V': Compute eigenvalues and eigenvectors.
30 *
31 * UPLO (input) CHARACTER*1
32 * = 'U': Upper triangle of A is stored;
33 * = 'L': Lower triangle of A is stored.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * KD (input) INTEGER
39 * The number of superdiagonals of the matrix A if UPLO = 'U',
40 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
41 *
42 * AB (input/output) COMPLEX*16 array, dimension (LDAB, N)
43 * On entry, the upper or lower triangle of the Hermitian band
44 * matrix A, stored in the first KD+1 rows of the array. The
45 * j-th column of A is stored in the j-th column of the array AB
46 * as follows:
47 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
48 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
49 *
50 * On exit, AB is overwritten by values generated during the
51 * reduction to tridiagonal form. If UPLO = 'U', the first
52 * superdiagonal and the diagonal of the tridiagonal matrix T
53 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
54 * the diagonal and first subdiagonal of T are returned in the
55 * first two rows of AB.
56 *
57 * LDAB (input) INTEGER
58 * The leading dimension of the array AB. LDAB >= KD + 1.
59 *
60 * W (output) DOUBLE PRECISION array, dimension (N)
61 * If INFO = 0, the eigenvalues in ascending order.
62 *
63 * Z (output) COMPLEX*16 array, dimension (LDZ, N)
64 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
65 * eigenvectors of the matrix A, with the i-th column of Z
66 * holding the eigenvector associated with W(i).
67 * If JOBZ = 'N', then Z is not referenced.
68 *
69 * LDZ (input) INTEGER
70 * The leading dimension of the array Z. LDZ >= 1, and if
71 * JOBZ = 'V', LDZ >= max(1,N).
72 *
73 * WORK (workspace) COMPLEX*16 array, dimension (N)
74 *
75 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
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, ONE
88 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
89 * ..
90 * .. Local Scalars ..
91 LOGICAL LOWER, WANTZ
92 INTEGER IINFO, IMAX, INDE, INDRWK, ISCALE
93 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
94 $ SMLNUM
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 DOUBLE PRECISION DLAMCH, ZLANHB
99 EXTERNAL LSAME, DLAMCH, ZLANHB
100 * ..
101 * .. External Subroutines ..
102 EXTERNAL DSCAL, DSTERF, XERBLA, ZHBTRD, ZLASCL, ZSTEQR
103 * ..
104 * .. Intrinsic Functions ..
105 INTRINSIC SQRT
106 * ..
107 * .. Executable Statements ..
108 *
109 * Test the input parameters.
110 *
111 WANTZ = LSAME( JOBZ, 'V' )
112 LOWER = LSAME( UPLO, 'L' )
113 *
114 INFO = 0
115 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
116 INFO = -1
117 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
118 INFO = -2
119 ELSE IF( N.LT.0 ) THEN
120 INFO = -3
121 ELSE IF( KD.LT.0 ) THEN
122 INFO = -4
123 ELSE IF( LDAB.LT.KD+1 ) THEN
124 INFO = -6
125 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
126 INFO = -9
127 END IF
128 *
129 IF( INFO.NE.0 ) THEN
130 CALL XERBLA( 'ZHBEV ', -INFO )
131 RETURN
132 END IF
133 *
134 * Quick return if possible
135 *
136 IF( N.EQ.0 )
137 $ RETURN
138 *
139 IF( N.EQ.1 ) THEN
140 IF( LOWER ) THEN
141 W( 1 ) = AB( 1, 1 )
142 ELSE
143 W( 1 ) = AB( KD+1, 1 )
144 END IF
145 IF( WANTZ )
146 $ Z( 1, 1 ) = ONE
147 RETURN
148 END IF
149 *
150 * Get machine constants.
151 *
152 SAFMIN = DLAMCH( 'Safe minimum' )
153 EPS = DLAMCH( 'Precision' )
154 SMLNUM = SAFMIN / EPS
155 BIGNUM = ONE / SMLNUM
156 RMIN = SQRT( SMLNUM )
157 RMAX = SQRT( BIGNUM )
158 *
159 * Scale matrix to allowable range, if necessary.
160 *
161 ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
162 ISCALE = 0
163 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
164 ISCALE = 1
165 SIGMA = RMIN / ANRM
166 ELSE IF( ANRM.GT.RMAX ) THEN
167 ISCALE = 1
168 SIGMA = RMAX / ANRM
169 END IF
170 IF( ISCALE.EQ.1 ) THEN
171 IF( LOWER ) THEN
172 CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
173 ELSE
174 CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
175 END IF
176 END IF
177 *
178 * Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
179 *
180 INDE = 1
181 CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
182 $ LDZ, WORK, IINFO )
183 *
184 * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
185 *
186 IF( .NOT.WANTZ ) THEN
187 CALL DSTERF( N, W, RWORK( INDE ), INFO )
188 ELSE
189 INDRWK = INDE + N
190 CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ,
191 $ RWORK( INDRWK ), INFO )
192 END IF
193 *
194 * If matrix was scaled, then rescale eigenvalues appropriately.
195 *
196 IF( ISCALE.EQ.1 ) THEN
197 IF( INFO.EQ.0 ) THEN
198 IMAX = N
199 ELSE
200 IMAX = INFO - 1
201 END IF
202 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
203 END IF
204 *
205 RETURN
206 *
207 * End of ZHBEV
208 *
209 END
2 $ RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, KD, LDAB, LDZ, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * ), W( * )
15 COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZHBEV computes all the eigenvalues and, optionally, eigenvectors of
22 * a complex Hermitian band matrix A.
23 *
24 * Arguments
25 * =========
26 *
27 * JOBZ (input) CHARACTER*1
28 * = 'N': Compute eigenvalues only;
29 * = 'V': Compute eigenvalues and eigenvectors.
30 *
31 * UPLO (input) CHARACTER*1
32 * = 'U': Upper triangle of A is stored;
33 * = 'L': Lower triangle of A is stored.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * KD (input) INTEGER
39 * The number of superdiagonals of the matrix A if UPLO = 'U',
40 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
41 *
42 * AB (input/output) COMPLEX*16 array, dimension (LDAB, N)
43 * On entry, the upper or lower triangle of the Hermitian band
44 * matrix A, stored in the first KD+1 rows of the array. The
45 * j-th column of A is stored in the j-th column of the array AB
46 * as follows:
47 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
48 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
49 *
50 * On exit, AB is overwritten by values generated during the
51 * reduction to tridiagonal form. If UPLO = 'U', the first
52 * superdiagonal and the diagonal of the tridiagonal matrix T
53 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
54 * the diagonal and first subdiagonal of T are returned in the
55 * first two rows of AB.
56 *
57 * LDAB (input) INTEGER
58 * The leading dimension of the array AB. LDAB >= KD + 1.
59 *
60 * W (output) DOUBLE PRECISION array, dimension (N)
61 * If INFO = 0, the eigenvalues in ascending order.
62 *
63 * Z (output) COMPLEX*16 array, dimension (LDZ, N)
64 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
65 * eigenvectors of the matrix A, with the i-th column of Z
66 * holding the eigenvector associated with W(i).
67 * If JOBZ = 'N', then Z is not referenced.
68 *
69 * LDZ (input) INTEGER
70 * The leading dimension of the array Z. LDZ >= 1, and if
71 * JOBZ = 'V', LDZ >= max(1,N).
72 *
73 * WORK (workspace) COMPLEX*16 array, dimension (N)
74 *
75 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
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, ONE
88 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
89 * ..
90 * .. Local Scalars ..
91 LOGICAL LOWER, WANTZ
92 INTEGER IINFO, IMAX, INDE, INDRWK, ISCALE
93 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
94 $ SMLNUM
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 DOUBLE PRECISION DLAMCH, ZLANHB
99 EXTERNAL LSAME, DLAMCH, ZLANHB
100 * ..
101 * .. External Subroutines ..
102 EXTERNAL DSCAL, DSTERF, XERBLA, ZHBTRD, ZLASCL, ZSTEQR
103 * ..
104 * .. Intrinsic Functions ..
105 INTRINSIC SQRT
106 * ..
107 * .. Executable Statements ..
108 *
109 * Test the input parameters.
110 *
111 WANTZ = LSAME( JOBZ, 'V' )
112 LOWER = LSAME( UPLO, 'L' )
113 *
114 INFO = 0
115 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
116 INFO = -1
117 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
118 INFO = -2
119 ELSE IF( N.LT.0 ) THEN
120 INFO = -3
121 ELSE IF( KD.LT.0 ) THEN
122 INFO = -4
123 ELSE IF( LDAB.LT.KD+1 ) THEN
124 INFO = -6
125 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
126 INFO = -9
127 END IF
128 *
129 IF( INFO.NE.0 ) THEN
130 CALL XERBLA( 'ZHBEV ', -INFO )
131 RETURN
132 END IF
133 *
134 * Quick return if possible
135 *
136 IF( N.EQ.0 )
137 $ RETURN
138 *
139 IF( N.EQ.1 ) THEN
140 IF( LOWER ) THEN
141 W( 1 ) = AB( 1, 1 )
142 ELSE
143 W( 1 ) = AB( KD+1, 1 )
144 END IF
145 IF( WANTZ )
146 $ Z( 1, 1 ) = ONE
147 RETURN
148 END IF
149 *
150 * Get machine constants.
151 *
152 SAFMIN = DLAMCH( 'Safe minimum' )
153 EPS = DLAMCH( 'Precision' )
154 SMLNUM = SAFMIN / EPS
155 BIGNUM = ONE / SMLNUM
156 RMIN = SQRT( SMLNUM )
157 RMAX = SQRT( BIGNUM )
158 *
159 * Scale matrix to allowable range, if necessary.
160 *
161 ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
162 ISCALE = 0
163 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
164 ISCALE = 1
165 SIGMA = RMIN / ANRM
166 ELSE IF( ANRM.GT.RMAX ) THEN
167 ISCALE = 1
168 SIGMA = RMAX / ANRM
169 END IF
170 IF( ISCALE.EQ.1 ) THEN
171 IF( LOWER ) THEN
172 CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
173 ELSE
174 CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
175 END IF
176 END IF
177 *
178 * Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
179 *
180 INDE = 1
181 CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
182 $ LDZ, WORK, IINFO )
183 *
184 * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
185 *
186 IF( .NOT.WANTZ ) THEN
187 CALL DSTERF( N, W, RWORK( INDE ), INFO )
188 ELSE
189 INDRWK = INDE + N
190 CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ,
191 $ RWORK( INDRWK ), INFO )
192 END IF
193 *
194 * If matrix was scaled, then rescale eigenvalues appropriately.
195 *
196 IF( ISCALE.EQ.1 ) THEN
197 IF( INFO.EQ.0 ) THEN
198 IMAX = N
199 ELSE
200 IMAX = INFO - 1
201 END IF
202 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
203 END IF
204 *
205 RETURN
206 *
207 * End of ZHBEV
208 *
209 END