![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSBGST 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSBGVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbgvd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbgvd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbgvd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, 00022 * Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBZ, UPLO 00026 * INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), W( * ), 00031 * $ WORK( * ), Z( LDZ, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> DSBGVD computes all the eigenvalues, and optionally, the eigenvectors 00041 *> of a real generalized symmetric-definite banded eigenproblem, of the 00042 *> form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric and 00043 *> banded, and B is also positive definite. If eigenvectors are 00044 *> desired, it uses a divide and conquer algorithm. 00045 *> 00046 *> The divide and conquer algorithm makes very mild assumptions about 00047 *> floating point arithmetic. It will work on machines with a guard 00048 *> digit in add/subtract, or on those binary machines without guard 00049 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00050 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines 00051 *> without guard digits, but we know of none. 00052 *> \endverbatim 00053 * 00054 * Arguments: 00055 * ========== 00056 * 00057 *> \param[in] JOBZ 00058 *> \verbatim 00059 *> JOBZ is CHARACTER*1 00060 *> = 'N': Compute eigenvalues only; 00061 *> = 'V': Compute eigenvalues and eigenvectors. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] UPLO 00065 *> \verbatim 00066 *> UPLO is CHARACTER*1 00067 *> = 'U': Upper triangles of A and B are stored; 00068 *> = 'L': Lower triangles of A and B are stored. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] N 00072 *> \verbatim 00073 *> N is INTEGER 00074 *> The order of the matrices A and B. N >= 0. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] KA 00078 *> \verbatim 00079 *> KA is INTEGER 00080 *> The number of superdiagonals of the matrix A if UPLO = 'U', 00081 *> or the number of subdiagonals if UPLO = 'L'. KA >= 0. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] KB 00085 *> \verbatim 00086 *> KB is INTEGER 00087 *> The number of superdiagonals of the matrix B if UPLO = 'U', 00088 *> or the number of subdiagonals if UPLO = 'L'. KB >= 0. 00089 *> \endverbatim 00090 *> 00091 *> \param[in,out] AB 00092 *> \verbatim 00093 *> AB is DOUBLE PRECISION array, dimension (LDAB, N) 00094 *> On entry, the upper or lower triangle of the symmetric band 00095 *> matrix A, stored in the first ka+1 rows of the array. The 00096 *> j-th column of A is stored in the j-th column of the array AB 00097 *> as follows: 00098 *> if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 00099 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 00100 *> 00101 *> On exit, the contents of AB are destroyed. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] LDAB 00105 *> \verbatim 00106 *> LDAB is INTEGER 00107 *> The leading dimension of the array AB. LDAB >= KA+1. 00108 *> \endverbatim 00109 *> 00110 *> \param[in,out] BB 00111 *> \verbatim 00112 *> BB is DOUBLE PRECISION array, dimension (LDBB, N) 00113 *> On entry, the upper or lower triangle of the symmetric band 00114 *> matrix B, stored in the first kb+1 rows of the array. The 00115 *> j-th column of B is stored in the j-th column of the array BB 00116 *> as follows: 00117 *> if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; 00118 *> if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb). 00119 *> 00120 *> On exit, the factor S from the split Cholesky factorization 00121 *> B = S**T*S, as returned by DPBSTF. 00122 *> \endverbatim 00123 *> 00124 *> \param[in] LDBB 00125 *> \verbatim 00126 *> LDBB is INTEGER 00127 *> The leading dimension of the array BB. LDBB >= KB+1. 00128 *> \endverbatim 00129 *> 00130 *> \param[out] W 00131 *> \verbatim 00132 *> W is DOUBLE PRECISION array, dimension (N) 00133 *> If INFO = 0, the eigenvalues in ascending order. 00134 *> \endverbatim 00135 *> 00136 *> \param[out] Z 00137 *> \verbatim 00138 *> Z is DOUBLE PRECISION array, dimension (LDZ, N) 00139 *> If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 00140 *> eigenvectors, with the i-th column of Z holding the 00141 *> eigenvector associated with W(i). The eigenvectors are 00142 *> normalized so Z**T*B*Z = I. 00143 *> If JOBZ = 'N', then Z is not referenced. 00144 *> \endverbatim 00145 *> 00146 *> \param[in] LDZ 00147 *> \verbatim 00148 *> LDZ is INTEGER 00149 *> The leading dimension of the array Z. LDZ >= 1, and if 00150 *> JOBZ = 'V', LDZ >= max(1,N). 00151 *> \endverbatim 00152 *> 00153 *> \param[out] WORK 00154 *> \verbatim 00155 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00156 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00157 *> \endverbatim 00158 *> 00159 *> \param[in] LWORK 00160 *> \verbatim 00161 *> LWORK is INTEGER 00162 *> The dimension of the array WORK. 00163 *> If N <= 1, LWORK >= 1. 00164 *> If JOBZ = 'N' and N > 1, LWORK >= 3*N. 00165 *> If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2. 00166 *> 00167 *> If LWORK = -1, then a workspace query is assumed; the routine 00168 *> only calculates the optimal sizes of the WORK and IWORK 00169 *> arrays, returns these values as the first entries of the WORK 00170 *> and IWORK arrays, and no error message related to LWORK or 00171 *> LIWORK is issued by XERBLA. 00172 *> \endverbatim 00173 *> 00174 *> \param[out] IWORK 00175 *> \verbatim 00176 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00177 *> On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK. 00178 *> \endverbatim 00179 *> 00180 *> \param[in] LIWORK 00181 *> \verbatim 00182 *> LIWORK is INTEGER 00183 *> The dimension of the array IWORK. 00184 *> If JOBZ = 'N' or N <= 1, LIWORK >= 1. 00185 *> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 00186 *> 00187 *> If LIWORK = -1, then a workspace query is assumed; the 00188 *> routine only calculates the optimal sizes of the WORK and 00189 *> IWORK arrays, returns these values as the first entries of 00190 *> the WORK and IWORK arrays, and no error message related to 00191 *> LWORK or LIWORK is issued by XERBLA. 00192 *> \endverbatim 00193 *> 00194 *> \param[out] INFO 00195 *> \verbatim 00196 *> INFO is INTEGER 00197 *> = 0: successful exit 00198 *> < 0: if INFO = -i, the i-th argument had an illegal value 00199 *> > 0: if INFO = i, and i is: 00200 *> <= N: the algorithm failed to converge: 00201 *> i off-diagonal elements of an intermediate 00202 *> tridiagonal form did not converge to zero; 00203 *> > N: if INFO = N + i, for 1 <= i <= N, then DPBSTF 00204 *> returned INFO = i: B is not positive definite. 00205 *> The factorization of B could not be completed and 00206 *> no eigenvalues or eigenvectors were computed. 00207 *> \endverbatim 00208 * 00209 * Authors: 00210 * ======== 00211 * 00212 *> \author Univ. of Tennessee 00213 *> \author Univ. of California Berkeley 00214 *> \author Univ. of Colorado Denver 00215 *> \author NAG Ltd. 00216 * 00217 *> \date November 2011 00218 * 00219 *> \ingroup doubleOTHEReigen 00220 * 00221 *> \par Contributors: 00222 * ================== 00223 *> 00224 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 00225 * 00226 * ===================================================================== 00227 SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, 00228 $ Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO ) 00229 * 00230 * -- LAPACK driver routine (version 3.4.0) -- 00231 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00233 * November 2011 00234 * 00235 * .. Scalar Arguments .. 00236 CHARACTER JOBZ, UPLO 00237 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N 00238 * .. 00239 * .. Array Arguments .. 00240 INTEGER IWORK( * ) 00241 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), W( * ), 00242 $ WORK( * ), Z( LDZ, * ) 00243 * .. 00244 * 00245 * ===================================================================== 00246 * 00247 * .. Parameters .. 00248 DOUBLE PRECISION ONE, ZERO 00249 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00250 * .. 00251 * .. Local Scalars .. 00252 LOGICAL LQUERY, UPPER, WANTZ 00253 CHARACTER VECT 00254 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2, 00255 $ LWMIN 00256 * .. 00257 * .. External Functions .. 00258 LOGICAL LSAME 00259 EXTERNAL LSAME 00260 * .. 00261 * .. External Subroutines .. 00262 EXTERNAL DGEMM, DLACPY, DPBSTF, DSBGST, DSBTRD, DSTEDC, 00263 $ DSTERF, XERBLA 00264 * .. 00265 * .. Executable Statements .. 00266 * 00267 * Test the input parameters. 00268 * 00269 WANTZ = LSAME( JOBZ, 'V' ) 00270 UPPER = LSAME( UPLO, 'U' ) 00271 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00272 * 00273 INFO = 0 00274 IF( N.LE.1 ) THEN 00275 LIWMIN = 1 00276 LWMIN = 1 00277 ELSE IF( WANTZ ) THEN 00278 LIWMIN = 3 + 5*N 00279 LWMIN = 1 + 5*N + 2*N**2 00280 ELSE 00281 LIWMIN = 1 00282 LWMIN = 2*N 00283 END IF 00284 * 00285 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00286 INFO = -1 00287 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 00288 INFO = -2 00289 ELSE IF( N.LT.0 ) THEN 00290 INFO = -3 00291 ELSE IF( KA.LT.0 ) THEN 00292 INFO = -4 00293 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 00294 INFO = -5 00295 ELSE IF( LDAB.LT.KA+1 ) THEN 00296 INFO = -7 00297 ELSE IF( LDBB.LT.KB+1 ) THEN 00298 INFO = -9 00299 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00300 INFO = -12 00301 END IF 00302 * 00303 IF( INFO.EQ.0 ) THEN 00304 WORK( 1 ) = LWMIN 00305 IWORK( 1 ) = LIWMIN 00306 * 00307 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00308 INFO = -14 00309 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00310 INFO = -16 00311 END IF 00312 END IF 00313 * 00314 IF( INFO.NE.0 ) THEN 00315 CALL XERBLA( 'DSBGVD', -INFO ) 00316 RETURN 00317 ELSE IF( LQUERY ) THEN 00318 RETURN 00319 END IF 00320 * 00321 * Quick return if possible 00322 * 00323 IF( N.EQ.0 ) 00324 $ RETURN 00325 * 00326 * Form a split Cholesky factorization of B. 00327 * 00328 CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO ) 00329 IF( INFO.NE.0 ) THEN 00330 INFO = N + INFO 00331 RETURN 00332 END IF 00333 * 00334 * Transform problem to standard eigenvalue problem. 00335 * 00336 INDE = 1 00337 INDWRK = INDE + N 00338 INDWK2 = INDWRK + N*N 00339 LLWRK2 = LWORK - INDWK2 + 1 00340 CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ, 00341 $ WORK( INDWRK ), IINFO ) 00342 * 00343 * Reduce to tridiagonal form. 00344 * 00345 IF( WANTZ ) THEN 00346 VECT = 'U' 00347 ELSE 00348 VECT = 'N' 00349 END IF 00350 CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ, 00351 $ WORK( INDWRK ), IINFO ) 00352 * 00353 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC. 00354 * 00355 IF( .NOT.WANTZ ) THEN 00356 CALL DSTERF( N, W, WORK( INDE ), INFO ) 00357 ELSE 00358 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, 00359 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO ) 00360 CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N, 00361 $ ZERO, WORK( INDWK2 ), N ) 00362 CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 00363 END IF 00364 * 00365 WORK( 1 ) = LWMIN 00366 IWORK( 1 ) = LIWMIN 00367 * 00368 RETURN 00369 * 00370 * End of DSBGVD 00371 * 00372 END