1 SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
2 $ Z, LDZ, WORK, 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, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), W( * ),
16 $ WORK( * ), Z( LDZ, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DSBGVD computes all the eigenvalues, and optionally, the eigenvectors
23 * of a real generalized symmetric-definite banded eigenproblem, of the
24 * form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric and
25 * banded, and B is also positive definite. If eigenvectors are
26 * desired, it uses a divide and conquer algorithm.
27 *
28 * The divide and conquer algorithm makes very mild assumptions about
29 * floating point arithmetic. It will work on machines with a guard
30 * digit in add/subtract, or on those binary machines without guard
31 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
32 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
33 * without guard digits, but we know of none.
34 *
35 * Arguments
36 * =========
37 *
38 * JOBZ (input) CHARACTER*1
39 * = 'N': Compute eigenvalues only;
40 * = 'V': Compute eigenvalues and eigenvectors.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangles of A and B are stored;
44 * = 'L': Lower triangles of A and B are stored.
45 *
46 * N (input) INTEGER
47 * The order of the matrices A and B. N >= 0.
48 *
49 * KA (input) INTEGER
50 * The number of superdiagonals of the matrix A if UPLO = 'U',
51 * or the number of subdiagonals if UPLO = 'L'. KA >= 0.
52 *
53 * KB (input) INTEGER
54 * The number of superdiagonals of the matrix B if UPLO = 'U',
55 * or the number of subdiagonals if UPLO = 'L'. KB >= 0.
56 *
57 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
58 * On entry, the upper or lower triangle of the symmetric band
59 * matrix A, stored in the first ka+1 rows of the array. The
60 * j-th column of A is stored in the j-th column of the array AB
61 * as follows:
62 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
63 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
64 *
65 * On exit, the contents of AB are destroyed.
66 *
67 * LDAB (input) INTEGER
68 * The leading dimension of the array AB. LDAB >= KA+1.
69 *
70 * BB (input/output) DOUBLE PRECISION array, dimension (LDBB, N)
71 * On entry, the upper or lower triangle of the symmetric band
72 * matrix B, stored in the first kb+1 rows of the array. The
73 * j-th column of B is stored in the j-th column of the array BB
74 * as follows:
75 * if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
76 * if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb).
77 *
78 * On exit, the factor S from the split Cholesky factorization
79 * B = S**T*S, as returned by DPBSTF.
80 *
81 * LDBB (input) INTEGER
82 * The leading dimension of the array BB. LDBB >= KB+1.
83 *
84 * W (output) DOUBLE PRECISION array, dimension (N)
85 * If INFO = 0, the eigenvalues in ascending order.
86 *
87 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
88 * If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
89 * eigenvectors, with the i-th column of Z holding the
90 * eigenvector associated with W(i). The eigenvectors are
91 * normalized so Z**T*B*Z = I.
92 * If JOBZ = 'N', then Z is not referenced.
93 *
94 * LDZ (input) INTEGER
95 * The leading dimension of the array Z. LDZ >= 1, and if
96 * JOBZ = 'V', LDZ >= max(1,N).
97 *
98 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
99 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
100 *
101 * LWORK (input) INTEGER
102 * The dimension of the array WORK.
103 * If N <= 1, LWORK >= 1.
104 * If JOBZ = 'N' and N > 1, LWORK >= 3*N.
105 * If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2.
106 *
107 * If LWORK = -1, then a workspace query is assumed; the routine
108 * only calculates the optimal sizes of the WORK and IWORK
109 * arrays, returns these values as the first entries of the WORK
110 * and IWORK arrays, and no error message related to LWORK or
111 * LIWORK is issued by XERBLA.
112 *
113 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
114 * On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
115 *
116 * LIWORK (input) INTEGER
117 * The dimension of the array IWORK.
118 * If JOBZ = 'N' or N <= 1, LIWORK >= 1.
119 * If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
120 *
121 * If LIWORK = -1, then a workspace query is assumed; the
122 * routine only calculates the optimal sizes of the WORK and
123 * IWORK arrays, returns these values as the first entries of
124 * the WORK and IWORK arrays, and no error message related to
125 * LWORK or LIWORK is issued by XERBLA.
126 *
127 * INFO (output) INTEGER
128 * = 0: successful exit
129 * < 0: if INFO = -i, the i-th argument had an illegal value
130 * > 0: if INFO = i, and i is:
131 * <= N: the algorithm failed to converge:
132 * i off-diagonal elements of an intermediate
133 * tridiagonal form did not converge to zero;
134 * > N: if INFO = N + i, for 1 <= i <= N, then DPBSTF
135 * returned INFO = i: B is not positive definite.
136 * The factorization of B could not be completed and
137 * no eigenvalues or eigenvectors were computed.
138 *
139 * Further Details
140 * ===============
141 *
142 * Based on contributions by
143 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148 DOUBLE PRECISION ONE, ZERO
149 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
150 * ..
151 * .. Local Scalars ..
152 LOGICAL LQUERY, UPPER, WANTZ
153 CHARACTER VECT
154 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2,
155 $ LWMIN
156 * ..
157 * .. External Functions ..
158 LOGICAL LSAME
159 EXTERNAL LSAME
160 * ..
161 * .. External Subroutines ..
162 EXTERNAL DGEMM, DLACPY, DPBSTF, DSBGST, DSBTRD, DSTEDC,
163 $ DSTERF, XERBLA
164 * ..
165 * .. Executable Statements ..
166 *
167 * Test the input parameters.
168 *
169 WANTZ = LSAME( JOBZ, 'V' )
170 UPPER = LSAME( UPLO, 'U' )
171 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
172 *
173 INFO = 0
174 IF( N.LE.1 ) THEN
175 LIWMIN = 1
176 LWMIN = 1
177 ELSE IF( WANTZ ) THEN
178 LIWMIN = 3 + 5*N
179 LWMIN = 1 + 5*N + 2*N**2
180 ELSE
181 LIWMIN = 1
182 LWMIN = 2*N
183 END IF
184 *
185 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
186 INFO = -1
187 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
188 INFO = -2
189 ELSE IF( N.LT.0 ) THEN
190 INFO = -3
191 ELSE IF( KA.LT.0 ) THEN
192 INFO = -4
193 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
194 INFO = -5
195 ELSE IF( LDAB.LT.KA+1 ) THEN
196 INFO = -7
197 ELSE IF( LDBB.LT.KB+1 ) THEN
198 INFO = -9
199 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
200 INFO = -12
201 END IF
202 *
203 IF( INFO.EQ.0 ) THEN
204 WORK( 1 ) = LWMIN
205 IWORK( 1 ) = LIWMIN
206 *
207 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
208 INFO = -14
209 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
210 INFO = -16
211 END IF
212 END IF
213 *
214 IF( INFO.NE.0 ) THEN
215 CALL XERBLA( 'DSBGVD', -INFO )
216 RETURN
217 ELSE IF( LQUERY ) THEN
218 RETURN
219 END IF
220 *
221 * Quick return if possible
222 *
223 IF( N.EQ.0 )
224 $ RETURN
225 *
226 * Form a split Cholesky factorization of B.
227 *
228 CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO )
229 IF( INFO.NE.0 ) THEN
230 INFO = N + INFO
231 RETURN
232 END IF
233 *
234 * Transform problem to standard eigenvalue problem.
235 *
236 INDE = 1
237 INDWRK = INDE + N
238 INDWK2 = INDWRK + N*N
239 LLWRK2 = LWORK - INDWK2 + 1
240 CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
241 $ WORK( INDWRK ), IINFO )
242 *
243 * Reduce to tridiagonal form.
244 *
245 IF( WANTZ ) THEN
246 VECT = 'U'
247 ELSE
248 VECT = 'N'
249 END IF
250 CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ,
251 $ WORK( INDWRK ), IINFO )
252 *
253 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
254 *
255 IF( .NOT.WANTZ ) THEN
256 CALL DSTERF( N, W, WORK( INDE ), INFO )
257 ELSE
258 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
259 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
260 CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
261 $ ZERO, WORK( INDWK2 ), N )
262 CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
263 END IF
264 *
265 WORK( 1 ) = LWMIN
266 IWORK( 1 ) = LIWMIN
267 *
268 RETURN
269 *
270 * End of DSBGVD
271 *
272 END
2 $ Z, LDZ, WORK, 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, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), W( * ),
16 $ WORK( * ), Z( LDZ, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DSBGVD computes all the eigenvalues, and optionally, the eigenvectors
23 * of a real generalized symmetric-definite banded eigenproblem, of the
24 * form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric and
25 * banded, and B is also positive definite. If eigenvectors are
26 * desired, it uses a divide and conquer algorithm.
27 *
28 * The divide and conquer algorithm makes very mild assumptions about
29 * floating point arithmetic. It will work on machines with a guard
30 * digit in add/subtract, or on those binary machines without guard
31 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
32 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
33 * without guard digits, but we know of none.
34 *
35 * Arguments
36 * =========
37 *
38 * JOBZ (input) CHARACTER*1
39 * = 'N': Compute eigenvalues only;
40 * = 'V': Compute eigenvalues and eigenvectors.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangles of A and B are stored;
44 * = 'L': Lower triangles of A and B are stored.
45 *
46 * N (input) INTEGER
47 * The order of the matrices A and B. N >= 0.
48 *
49 * KA (input) INTEGER
50 * The number of superdiagonals of the matrix A if UPLO = 'U',
51 * or the number of subdiagonals if UPLO = 'L'. KA >= 0.
52 *
53 * KB (input) INTEGER
54 * The number of superdiagonals of the matrix B if UPLO = 'U',
55 * or the number of subdiagonals if UPLO = 'L'. KB >= 0.
56 *
57 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
58 * On entry, the upper or lower triangle of the symmetric band
59 * matrix A, stored in the first ka+1 rows of the array. The
60 * j-th column of A is stored in the j-th column of the array AB
61 * as follows:
62 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
63 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
64 *
65 * On exit, the contents of AB are destroyed.
66 *
67 * LDAB (input) INTEGER
68 * The leading dimension of the array AB. LDAB >= KA+1.
69 *
70 * BB (input/output) DOUBLE PRECISION array, dimension (LDBB, N)
71 * On entry, the upper or lower triangle of the symmetric band
72 * matrix B, stored in the first kb+1 rows of the array. The
73 * j-th column of B is stored in the j-th column of the array BB
74 * as follows:
75 * if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
76 * if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb).
77 *
78 * On exit, the factor S from the split Cholesky factorization
79 * B = S**T*S, as returned by DPBSTF.
80 *
81 * LDBB (input) INTEGER
82 * The leading dimension of the array BB. LDBB >= KB+1.
83 *
84 * W (output) DOUBLE PRECISION array, dimension (N)
85 * If INFO = 0, the eigenvalues in ascending order.
86 *
87 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
88 * If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
89 * eigenvectors, with the i-th column of Z holding the
90 * eigenvector associated with W(i). The eigenvectors are
91 * normalized so Z**T*B*Z = I.
92 * If JOBZ = 'N', then Z is not referenced.
93 *
94 * LDZ (input) INTEGER
95 * The leading dimension of the array Z. LDZ >= 1, and if
96 * JOBZ = 'V', LDZ >= max(1,N).
97 *
98 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
99 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
100 *
101 * LWORK (input) INTEGER
102 * The dimension of the array WORK.
103 * If N <= 1, LWORK >= 1.
104 * If JOBZ = 'N' and N > 1, LWORK >= 3*N.
105 * If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2.
106 *
107 * If LWORK = -1, then a workspace query is assumed; the routine
108 * only calculates the optimal sizes of the WORK and IWORK
109 * arrays, returns these values as the first entries of the WORK
110 * and IWORK arrays, and no error message related to LWORK or
111 * LIWORK is issued by XERBLA.
112 *
113 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
114 * On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
115 *
116 * LIWORK (input) INTEGER
117 * The dimension of the array IWORK.
118 * If JOBZ = 'N' or N <= 1, LIWORK >= 1.
119 * If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
120 *
121 * If LIWORK = -1, then a workspace query is assumed; the
122 * routine only calculates the optimal sizes of the WORK and
123 * IWORK arrays, returns these values as the first entries of
124 * the WORK and IWORK arrays, and no error message related to
125 * LWORK or LIWORK is issued by XERBLA.
126 *
127 * INFO (output) INTEGER
128 * = 0: successful exit
129 * < 0: if INFO = -i, the i-th argument had an illegal value
130 * > 0: if INFO = i, and i is:
131 * <= N: the algorithm failed to converge:
132 * i off-diagonal elements of an intermediate
133 * tridiagonal form did not converge to zero;
134 * > N: if INFO = N + i, for 1 <= i <= N, then DPBSTF
135 * returned INFO = i: B is not positive definite.
136 * The factorization of B could not be completed and
137 * no eigenvalues or eigenvectors were computed.
138 *
139 * Further Details
140 * ===============
141 *
142 * Based on contributions by
143 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148 DOUBLE PRECISION ONE, ZERO
149 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
150 * ..
151 * .. Local Scalars ..
152 LOGICAL LQUERY, UPPER, WANTZ
153 CHARACTER VECT
154 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2,
155 $ LWMIN
156 * ..
157 * .. External Functions ..
158 LOGICAL LSAME
159 EXTERNAL LSAME
160 * ..
161 * .. External Subroutines ..
162 EXTERNAL DGEMM, DLACPY, DPBSTF, DSBGST, DSBTRD, DSTEDC,
163 $ DSTERF, XERBLA
164 * ..
165 * .. Executable Statements ..
166 *
167 * Test the input parameters.
168 *
169 WANTZ = LSAME( JOBZ, 'V' )
170 UPPER = LSAME( UPLO, 'U' )
171 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
172 *
173 INFO = 0
174 IF( N.LE.1 ) THEN
175 LIWMIN = 1
176 LWMIN = 1
177 ELSE IF( WANTZ ) THEN
178 LIWMIN = 3 + 5*N
179 LWMIN = 1 + 5*N + 2*N**2
180 ELSE
181 LIWMIN = 1
182 LWMIN = 2*N
183 END IF
184 *
185 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
186 INFO = -1
187 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
188 INFO = -2
189 ELSE IF( N.LT.0 ) THEN
190 INFO = -3
191 ELSE IF( KA.LT.0 ) THEN
192 INFO = -4
193 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
194 INFO = -5
195 ELSE IF( LDAB.LT.KA+1 ) THEN
196 INFO = -7
197 ELSE IF( LDBB.LT.KB+1 ) THEN
198 INFO = -9
199 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
200 INFO = -12
201 END IF
202 *
203 IF( INFO.EQ.0 ) THEN
204 WORK( 1 ) = LWMIN
205 IWORK( 1 ) = LIWMIN
206 *
207 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
208 INFO = -14
209 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
210 INFO = -16
211 END IF
212 END IF
213 *
214 IF( INFO.NE.0 ) THEN
215 CALL XERBLA( 'DSBGVD', -INFO )
216 RETURN
217 ELSE IF( LQUERY ) THEN
218 RETURN
219 END IF
220 *
221 * Quick return if possible
222 *
223 IF( N.EQ.0 )
224 $ RETURN
225 *
226 * Form a split Cholesky factorization of B.
227 *
228 CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO )
229 IF( INFO.NE.0 ) THEN
230 INFO = N + INFO
231 RETURN
232 END IF
233 *
234 * Transform problem to standard eigenvalue problem.
235 *
236 INDE = 1
237 INDWRK = INDE + N
238 INDWK2 = INDWRK + N*N
239 LLWRK2 = LWORK - INDWK2 + 1
240 CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
241 $ WORK( INDWRK ), IINFO )
242 *
243 * Reduce to tridiagonal form.
244 *
245 IF( WANTZ ) THEN
246 VECT = 'U'
247 ELSE
248 VECT = 'N'
249 END IF
250 CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ,
251 $ WORK( INDWRK ), IINFO )
252 *
253 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
254 *
255 IF( .NOT.WANTZ ) THEN
256 CALL DSTERF( N, W, WORK( INDE ), INFO )
257 ELSE
258 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
259 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
260 CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
261 $ ZERO, WORK( INDWK2 ), N )
262 CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
263 END IF
264 *
265 WORK( 1 ) = LWMIN
266 IWORK( 1 ) = LIWMIN
267 *
268 RETURN
269 *
270 * End of DSBGVD
271 *
272 END