 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, 00002 $ Z, LDZ, WORK, LWORK, IWORK, LIWORK, 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, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 REAL AB( LDAB, * ), BB( LDBB, * ), W( * ), 00016 $ WORK( * ), Z( LDZ, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * SSBGVD computes all the eigenvalues, and optionally, the eigenvectors 00023 * of a real generalized symmetric-definite banded eigenproblem, of the 00024 * form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric and 00025 * banded, and B is also positive definite. If eigenvectors are 00026 * desired, it uses a divide and conquer algorithm. 00027 * 00028 * The divide and conquer algorithm makes very mild assumptions about 00029 * floating point arithmetic. It will work on machines with a guard 00030 * digit in add/subtract, or on those binary machines without guard 00031 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00032 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00033 * without guard digits, but we know of none. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * JOBZ (input) CHARACTER*1 00039 * = 'N': Compute eigenvalues only; 00040 * = 'V': Compute eigenvalues and eigenvectors. 00041 * 00042 * UPLO (input) CHARACTER*1 00043 * = 'U': Upper triangles of A and B are stored; 00044 * = 'L': Lower triangles of A and B are stored. 00045 * 00046 * N (input) INTEGER 00047 * The order of the matrices A and B. N >= 0. 00048 * 00049 * KA (input) INTEGER 00050 * The number of superdiagonals of the matrix A if UPLO = 'U', 00051 * or the number of subdiagonals if UPLO = 'L'. KA >= 0. 00052 * 00053 * KB (input) INTEGER 00054 * The number of superdiagonals of the matrix B if UPLO = 'U', 00055 * or the number of subdiagonals if UPLO = 'L'. KB >= 0. 00056 * 00057 * AB (input/output) REAL array, dimension (LDAB, N) 00058 * On entry, the upper or lower triangle of the symmetric band 00059 * matrix A, stored in the first ka+1 rows of the array. The 00060 * j-th column of A is stored in the j-th column of the array AB 00061 * as follows: 00062 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 00063 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 00064 * 00065 * On exit, the contents of AB are destroyed. 00066 * 00067 * LDAB (input) INTEGER 00068 * The leading dimension of the array AB. LDAB >= KA+1. 00069 * 00070 * BB (input/output) REAL array, dimension (LDBB, N) 00071 * On entry, the upper or lower triangle of the symmetric band 00072 * matrix B, stored in the first kb+1 rows of the array. The 00073 * j-th column of B is stored in the j-th column of the array BB 00074 * as follows: 00075 * if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; 00076 * if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb). 00077 * 00078 * On exit, the factor S from the split Cholesky factorization 00079 * B = S**T*S, as returned by SPBSTF. 00080 * 00081 * LDBB (input) INTEGER 00082 * The leading dimension of the array BB. LDBB >= KB+1. 00083 * 00084 * W (output) REAL array, dimension (N) 00085 * If INFO = 0, the eigenvalues in ascending order. 00086 * 00087 * Z (output) REAL array, dimension (LDZ, N) 00088 * If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 00089 * eigenvectors, with the i-th column of Z holding the 00090 * eigenvector associated with W(i). The eigenvectors are 00091 * normalized so Z**T*B*Z = I. 00092 * If JOBZ = 'N', then Z is not referenced. 00093 * 00094 * LDZ (input) INTEGER 00095 * The leading dimension of the array Z. LDZ >= 1, and if 00096 * JOBZ = 'V', LDZ >= max(1,N). 00097 * 00098 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) 00099 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00100 * 00101 * LWORK (input) INTEGER 00102 * The dimension of the array WORK. 00103 * If N <= 1, LWORK >= 1. 00104 * If JOBZ = 'N' and N > 1, LWORK >= 3*N. 00105 * If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2. 00106 * 00107 * If LWORK = -1, then a workspace query is assumed; the routine 00108 * only calculates the optimal sizes of the WORK and IWORK 00109 * arrays, returns these values as the first entries of the WORK 00110 * and IWORK arrays, and no error message related to LWORK or 00111 * LIWORK is issued by XERBLA. 00112 * 00113 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00114 * On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK. 00115 * 00116 * LIWORK (input) INTEGER 00117 * The dimension of the array IWORK. 00118 * If JOBZ = 'N' or N <= 1, LIWORK >= 1. 00119 * If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 00120 * 00121 * If LIWORK = -1, then a workspace query is assumed; the 00122 * routine only calculates the optimal sizes of the WORK and 00123 * IWORK arrays, returns these values as the first entries of 00124 * the WORK and IWORK arrays, and no error message related to 00125 * LWORK or LIWORK is issued by XERBLA. 00126 * 00127 * INFO (output) INTEGER 00128 * = 0: successful exit 00129 * < 0: if INFO = -i, the i-th argument had an illegal value 00130 * > 0: if INFO = i, and i is: 00131 * <= N: the algorithm failed to converge: 00132 * i off-diagonal elements of an intermediate 00133 * tridiagonal form did not converge to zero; 00134 * > N: if INFO = N + i, for 1 <= i <= N, then SPBSTF 00135 * returned INFO = i: B is not positive definite. 00136 * The factorization of B could not be completed and 00137 * no eigenvalues or eigenvectors were computed. 00138 * 00139 * Further Details 00140 * =============== 00141 * 00142 * Based on contributions by 00143 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 00144 * 00145 * ===================================================================== 00146 * 00147 * .. Parameters .. 00148 REAL ONE, ZERO 00149 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00150 * .. 00151 * .. Local Scalars .. 00152 LOGICAL LQUERY, UPPER, WANTZ 00153 CHARACTER VECT 00154 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2, 00155 $ LWMIN 00156 * .. 00157 * .. External Functions .. 00158 LOGICAL LSAME 00159 EXTERNAL LSAME 00160 * .. 00161 * .. External Subroutines .. 00162 EXTERNAL SGEMM, SLACPY, SPBSTF, SSBGST, SSBTRD, SSTEDC, 00163 $ SSTERF, XERBLA 00164 * .. 00165 * .. Executable Statements .. 00166 * 00167 * Test the input parameters. 00168 * 00169 WANTZ = LSAME( JOBZ, 'V' ) 00170 UPPER = LSAME( UPLO, 'U' ) 00171 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00172 * 00173 INFO = 0 00174 IF( N.LE.1 ) THEN 00175 LIWMIN = 1 00176 LWMIN = 1 00177 ELSE IF( WANTZ ) THEN 00178 LIWMIN = 3 + 5*N 00179 LWMIN = 1 + 5*N + 2*N**2 00180 ELSE 00181 LIWMIN = 1 00182 LWMIN = 2*N 00183 END IF 00184 * 00185 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00186 INFO = -1 00187 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 00188 INFO = -2 00189 ELSE IF( N.LT.0 ) THEN 00190 INFO = -3 00191 ELSE IF( KA.LT.0 ) THEN 00192 INFO = -4 00193 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 00194 INFO = -5 00195 ELSE IF( LDAB.LT.KA+1 ) THEN 00196 INFO = -7 00197 ELSE IF( LDBB.LT.KB+1 ) THEN 00198 INFO = -9 00199 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00200 INFO = -12 00201 END IF 00202 * 00203 IF( INFO.EQ.0 ) THEN 00204 WORK( 1 ) = LWMIN 00205 IWORK( 1 ) = LIWMIN 00206 * 00207 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00208 INFO = -14 00209 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00210 INFO = -16 00211 END IF 00212 END IF 00213 * 00214 IF( INFO.NE.0 ) THEN 00215 CALL XERBLA( 'SSBGVD', -INFO ) 00216 RETURN 00217 ELSE IF( LQUERY ) THEN 00218 RETURN 00219 END IF 00220 * 00221 * Quick return if possible 00222 * 00223 IF( N.EQ.0 ) 00224 $ RETURN 00225 * 00226 * Form a split Cholesky factorization of B. 00227 * 00228 CALL SPBSTF( UPLO, N, KB, BB, LDBB, INFO ) 00229 IF( INFO.NE.0 ) THEN 00230 INFO = N + INFO 00231 RETURN 00232 END IF 00233 * 00234 * Transform problem to standard eigenvalue problem. 00235 * 00236 INDE = 1 00237 INDWRK = INDE + N 00238 INDWK2 = INDWRK + N*N 00239 LLWRK2 = LWORK - INDWK2 + 1 00240 CALL SSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ, 00241 $ WORK( INDWRK ), IINFO ) 00242 * 00243 * Reduce to tridiagonal form. 00244 * 00245 IF( WANTZ ) THEN 00246 VECT = 'U' 00247 ELSE 00248 VECT = 'N' 00249 END IF 00250 CALL SSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ, 00251 $ WORK( INDWRK ), IINFO ) 00252 * 00253 * For eigenvalues only, call SSTERF. For eigenvectors, call SSTEDC. 00254 * 00255 IF( .NOT.WANTZ ) THEN 00256 CALL SSTERF( N, W, WORK( INDE ), INFO ) 00257 ELSE 00258 CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, 00259 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO ) 00260 CALL SGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N, 00261 $ ZERO, WORK( INDWK2 ), N ) 00262 CALL SLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 00263 END IF 00264 * 00265 WORK( 1 ) = LWMIN 00266 IWORK( 1 ) = LIWMIN 00267 * 00268 RETURN 00269 * 00270 * End of SSBGVD 00271 * 00272 END
1.7.3 
