1 SUBROUTINE DSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
2 $ LWORK, IWORK, LIWORK, 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, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSBEVD computes all the eigenvalues and, optionally, eigenvectors of
22 * a real symmetric band matrix A. If eigenvectors are desired, it uses
23 * a divide and conquer algorithm.
24 *
25 * The divide and conquer algorithm makes very mild assumptions about
26 * floating point arithmetic. It will work on machines with a guard
27 * digit in add/subtract, or on those binary machines without guard
28 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30 * without guard digits, but we know of none.
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 * KD (input) INTEGER
47 * The number of superdiagonals of the matrix A if UPLO = 'U',
48 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
49 *
50 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
51 * On entry, the upper or lower triangle of the symmetric band
52 * matrix A, stored in the first KD+1 rows of the array. The
53 * j-th column of A is stored in the j-th column of the array AB
54 * as follows:
55 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
56 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
57 *
58 * On exit, AB is overwritten by values generated during the
59 * reduction to tridiagonal form. If UPLO = 'U', the first
60 * superdiagonal and the diagonal of the tridiagonal matrix T
61 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
62 * the diagonal and first subdiagonal of T are returned in the
63 * first two rows of AB.
64 *
65 * LDAB (input) INTEGER
66 * The leading dimension of the array AB. LDAB >= KD + 1.
67 *
68 * W (output) DOUBLE PRECISION array, dimension (N)
69 * If INFO = 0, the eigenvalues in ascending order.
70 *
71 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
72 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
73 * eigenvectors of the matrix A, with the i-th column of Z
74 * holding the eigenvector associated with W(i).
75 * If JOBZ = 'N', then Z is not referenced.
76 *
77 * LDZ (input) INTEGER
78 * The leading dimension of the array Z. LDZ >= 1, and if
79 * JOBZ = 'V', LDZ >= max(1,N).
80 *
81 * WORK (workspace/output) DOUBLE PRECISION array,
82 * dimension (LWORK)
83 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84 *
85 * LWORK (input) INTEGER
86 * The dimension of the array WORK.
87 * IF N <= 1, LWORK must be at least 1.
88 * If JOBZ = 'N' and N > 2, LWORK must be at least 2*N.
89 * If JOBZ = 'V' and N > 2, LWORK must be at least
90 * ( 1 + 5*N + 2*N**2 ).
91 *
92 * If LWORK = -1, then a workspace query is assumed; the routine
93 * only calculates the optimal sizes of the WORK and IWORK
94 * arrays, returns these values as the first entries of the WORK
95 * and IWORK arrays, and no error message related to LWORK or
96 * LIWORK is issued by XERBLA.
97 *
98 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
99 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
100 *
101 * LIWORK (input) INTEGER
102 * The dimension of the array LIWORK.
103 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
104 * If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
105 *
106 * If LIWORK = -1, then a workspace query is assumed; the
107 * routine only calculates the optimal sizes of the WORK and
108 * IWORK arrays, returns these values as the first entries of
109 * the WORK and IWORK arrays, and no error message related to
110 * LWORK or LIWORK is issued by XERBLA.
111 *
112 * INFO (output) INTEGER
113 * = 0: successful exit
114 * < 0: if INFO = -i, the i-th argument had an illegal value
115 * > 0: if INFO = i, the algorithm failed to converge; i
116 * off-diagonal elements of an intermediate tridiagonal
117 * form did not converge to zero.
118 *
119 * =====================================================================
120 *
121 * .. Parameters ..
122 DOUBLE PRECISION ZERO, ONE
123 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
124 * ..
125 * .. Local Scalars ..
126 LOGICAL LOWER, LQUERY, WANTZ
127 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
128 $ LLWRK2, LWMIN
129 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
130 $ SMLNUM
131 * ..
132 * .. External Functions ..
133 LOGICAL LSAME
134 DOUBLE PRECISION DLAMCH, DLANSB
135 EXTERNAL LSAME, DLAMCH, DLANSB
136 * ..
137 * .. External Subroutines ..
138 EXTERNAL DGEMM, DLACPY, DLASCL, DSBTRD, DSCAL, DSTEDC,
139 $ DSTERF, XERBLA
140 * ..
141 * .. Intrinsic Functions ..
142 INTRINSIC SQRT
143 * ..
144 * .. Executable Statements ..
145 *
146 * Test the input parameters.
147 *
148 WANTZ = LSAME( JOBZ, 'V' )
149 LOWER = LSAME( UPLO, 'L' )
150 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
151 *
152 INFO = 0
153 IF( N.LE.1 ) THEN
154 LIWMIN = 1
155 LWMIN = 1
156 ELSE
157 IF( WANTZ ) THEN
158 LIWMIN = 3 + 5*N
159 LWMIN = 1 + 5*N + 2*N**2
160 ELSE
161 LIWMIN = 1
162 LWMIN = 2*N
163 END IF
164 END IF
165 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
166 INFO = -1
167 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
168 INFO = -2
169 ELSE IF( N.LT.0 ) THEN
170 INFO = -3
171 ELSE IF( KD.LT.0 ) THEN
172 INFO = -4
173 ELSE IF( LDAB.LT.KD+1 ) THEN
174 INFO = -6
175 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
176 INFO = -9
177 END IF
178 *
179 IF( INFO.EQ.0 ) THEN
180 WORK( 1 ) = LWMIN
181 IWORK( 1 ) = LIWMIN
182 *
183 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
184 INFO = -11
185 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
186 INFO = -13
187 END IF
188 END IF
189 *
190 IF( INFO.NE.0 ) THEN
191 CALL XERBLA( 'DSBEVD', -INFO )
192 RETURN
193 ELSE IF( LQUERY ) THEN
194 RETURN
195 END IF
196 *
197 * Quick return if possible
198 *
199 IF( N.EQ.0 )
200 $ RETURN
201 *
202 IF( N.EQ.1 ) THEN
203 W( 1 ) = AB( 1, 1 )
204 IF( WANTZ )
205 $ Z( 1, 1 ) = ONE
206 RETURN
207 END IF
208 *
209 * Get machine constants.
210 *
211 SAFMIN = DLAMCH( 'Safe minimum' )
212 EPS = DLAMCH( 'Precision' )
213 SMLNUM = SAFMIN / EPS
214 BIGNUM = ONE / SMLNUM
215 RMIN = SQRT( SMLNUM )
216 RMAX = SQRT( BIGNUM )
217 *
218 * Scale matrix to allowable range, if necessary.
219 *
220 ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
221 ISCALE = 0
222 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
223 ISCALE = 1
224 SIGMA = RMIN / ANRM
225 ELSE IF( ANRM.GT.RMAX ) THEN
226 ISCALE = 1
227 SIGMA = RMAX / ANRM
228 END IF
229 IF( ISCALE.EQ.1 ) THEN
230 IF( LOWER ) THEN
231 CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
232 ELSE
233 CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
234 END IF
235 END IF
236 *
237 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
238 *
239 INDE = 1
240 INDWRK = INDE + N
241 INDWK2 = INDWRK + N*N
242 LLWRK2 = LWORK - INDWK2 + 1
243 CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
244 $ WORK( INDWRK ), IINFO )
245 *
246 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
247 *
248 IF( .NOT.WANTZ ) THEN
249 CALL DSTERF( N, W, WORK( INDE ), INFO )
250 ELSE
251 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
252 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
253 CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
254 $ ZERO, WORK( INDWK2 ), N )
255 CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
256 END IF
257 *
258 * If matrix was scaled, then rescale eigenvalues appropriately.
259 *
260 IF( ISCALE.EQ.1 )
261 $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
262 *
263 WORK( 1 ) = LWMIN
264 IWORK( 1 ) = LIWMIN
265 RETURN
266 *
267 * End of DSBEVD
268 *
269 END
2 $ LWORK, IWORK, LIWORK, 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, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSBEVD computes all the eigenvalues and, optionally, eigenvectors of
22 * a real symmetric band matrix A. If eigenvectors are desired, it uses
23 * a divide and conquer algorithm.
24 *
25 * The divide and conquer algorithm makes very mild assumptions about
26 * floating point arithmetic. It will work on machines with a guard
27 * digit in add/subtract, or on those binary machines without guard
28 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30 * without guard digits, but we know of none.
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 * KD (input) INTEGER
47 * The number of superdiagonals of the matrix A if UPLO = 'U',
48 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
49 *
50 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
51 * On entry, the upper or lower triangle of the symmetric band
52 * matrix A, stored in the first KD+1 rows of the array. The
53 * j-th column of A is stored in the j-th column of the array AB
54 * as follows:
55 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
56 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
57 *
58 * On exit, AB is overwritten by values generated during the
59 * reduction to tridiagonal form. If UPLO = 'U', the first
60 * superdiagonal and the diagonal of the tridiagonal matrix T
61 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
62 * the diagonal and first subdiagonal of T are returned in the
63 * first two rows of AB.
64 *
65 * LDAB (input) INTEGER
66 * The leading dimension of the array AB. LDAB >= KD + 1.
67 *
68 * W (output) DOUBLE PRECISION array, dimension (N)
69 * If INFO = 0, the eigenvalues in ascending order.
70 *
71 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
72 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
73 * eigenvectors of the matrix A, with the i-th column of Z
74 * holding the eigenvector associated with W(i).
75 * If JOBZ = 'N', then Z is not referenced.
76 *
77 * LDZ (input) INTEGER
78 * The leading dimension of the array Z. LDZ >= 1, and if
79 * JOBZ = 'V', LDZ >= max(1,N).
80 *
81 * WORK (workspace/output) DOUBLE PRECISION array,
82 * dimension (LWORK)
83 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84 *
85 * LWORK (input) INTEGER
86 * The dimension of the array WORK.
87 * IF N <= 1, LWORK must be at least 1.
88 * If JOBZ = 'N' and N > 2, LWORK must be at least 2*N.
89 * If JOBZ = 'V' and N > 2, LWORK must be at least
90 * ( 1 + 5*N + 2*N**2 ).
91 *
92 * If LWORK = -1, then a workspace query is assumed; the routine
93 * only calculates the optimal sizes of the WORK and IWORK
94 * arrays, returns these values as the first entries of the WORK
95 * and IWORK arrays, and no error message related to LWORK or
96 * LIWORK is issued by XERBLA.
97 *
98 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
99 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
100 *
101 * LIWORK (input) INTEGER
102 * The dimension of the array LIWORK.
103 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
104 * If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
105 *
106 * If LIWORK = -1, then a workspace query is assumed; the
107 * routine only calculates the optimal sizes of the WORK and
108 * IWORK arrays, returns these values as the first entries of
109 * the WORK and IWORK arrays, and no error message related to
110 * LWORK or LIWORK is issued by XERBLA.
111 *
112 * INFO (output) INTEGER
113 * = 0: successful exit
114 * < 0: if INFO = -i, the i-th argument had an illegal value
115 * > 0: if INFO = i, the algorithm failed to converge; i
116 * off-diagonal elements of an intermediate tridiagonal
117 * form did not converge to zero.
118 *
119 * =====================================================================
120 *
121 * .. Parameters ..
122 DOUBLE PRECISION ZERO, ONE
123 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
124 * ..
125 * .. Local Scalars ..
126 LOGICAL LOWER, LQUERY, WANTZ
127 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
128 $ LLWRK2, LWMIN
129 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
130 $ SMLNUM
131 * ..
132 * .. External Functions ..
133 LOGICAL LSAME
134 DOUBLE PRECISION DLAMCH, DLANSB
135 EXTERNAL LSAME, DLAMCH, DLANSB
136 * ..
137 * .. External Subroutines ..
138 EXTERNAL DGEMM, DLACPY, DLASCL, DSBTRD, DSCAL, DSTEDC,
139 $ DSTERF, XERBLA
140 * ..
141 * .. Intrinsic Functions ..
142 INTRINSIC SQRT
143 * ..
144 * .. Executable Statements ..
145 *
146 * Test the input parameters.
147 *
148 WANTZ = LSAME( JOBZ, 'V' )
149 LOWER = LSAME( UPLO, 'L' )
150 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
151 *
152 INFO = 0
153 IF( N.LE.1 ) THEN
154 LIWMIN = 1
155 LWMIN = 1
156 ELSE
157 IF( WANTZ ) THEN
158 LIWMIN = 3 + 5*N
159 LWMIN = 1 + 5*N + 2*N**2
160 ELSE
161 LIWMIN = 1
162 LWMIN = 2*N
163 END IF
164 END IF
165 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
166 INFO = -1
167 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
168 INFO = -2
169 ELSE IF( N.LT.0 ) THEN
170 INFO = -3
171 ELSE IF( KD.LT.0 ) THEN
172 INFO = -4
173 ELSE IF( LDAB.LT.KD+1 ) THEN
174 INFO = -6
175 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
176 INFO = -9
177 END IF
178 *
179 IF( INFO.EQ.0 ) THEN
180 WORK( 1 ) = LWMIN
181 IWORK( 1 ) = LIWMIN
182 *
183 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
184 INFO = -11
185 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
186 INFO = -13
187 END IF
188 END IF
189 *
190 IF( INFO.NE.0 ) THEN
191 CALL XERBLA( 'DSBEVD', -INFO )
192 RETURN
193 ELSE IF( LQUERY ) THEN
194 RETURN
195 END IF
196 *
197 * Quick return if possible
198 *
199 IF( N.EQ.0 )
200 $ RETURN
201 *
202 IF( N.EQ.1 ) THEN
203 W( 1 ) = AB( 1, 1 )
204 IF( WANTZ )
205 $ Z( 1, 1 ) = ONE
206 RETURN
207 END IF
208 *
209 * Get machine constants.
210 *
211 SAFMIN = DLAMCH( 'Safe minimum' )
212 EPS = DLAMCH( 'Precision' )
213 SMLNUM = SAFMIN / EPS
214 BIGNUM = ONE / SMLNUM
215 RMIN = SQRT( SMLNUM )
216 RMAX = SQRT( BIGNUM )
217 *
218 * Scale matrix to allowable range, if necessary.
219 *
220 ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
221 ISCALE = 0
222 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
223 ISCALE = 1
224 SIGMA = RMIN / ANRM
225 ELSE IF( ANRM.GT.RMAX ) THEN
226 ISCALE = 1
227 SIGMA = RMAX / ANRM
228 END IF
229 IF( ISCALE.EQ.1 ) THEN
230 IF( LOWER ) THEN
231 CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
232 ELSE
233 CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
234 END IF
235 END IF
236 *
237 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
238 *
239 INDE = 1
240 INDWRK = INDE + N
241 INDWK2 = INDWRK + N*N
242 LLWRK2 = LWORK - INDWK2 + 1
243 CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
244 $ WORK( INDWRK ), IINFO )
245 *
246 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
247 *
248 IF( .NOT.WANTZ ) THEN
249 CALL DSTERF( N, W, WORK( INDE ), INFO )
250 ELSE
251 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
252 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
253 CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
254 $ ZERO, WORK( INDWK2 ), N )
255 CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
256 END IF
257 *
258 * If matrix was scaled, then rescale eigenvalues appropriately.
259 *
260 IF( ISCALE.EQ.1 )
261 $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
262 *
263 WORK( 1 ) = LWMIN
264 IWORK( 1 ) = LIWMIN
265 RETURN
266 *
267 * End of DSBEVD
268 *
269 END