00001 /* chbgvd.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static complex c_b1 = {1.f,0.f}; 00019 static complex c_b2 = {0.f,0.f}; 00020 00021 /* Subroutine */ int chbgvd_(char *jobz, char *uplo, integer *n, integer *ka, 00022 integer *kb, complex *ab, integer *ldab, complex *bb, integer *ldbb, 00023 real *w, complex *z__, integer *ldz, complex *work, integer *lwork, 00024 real *rwork, integer *lrwork, integer *iwork, integer *liwork, 00025 integer *info) 00026 { 00027 /* System generated locals */ 00028 integer ab_dim1, ab_offset, bb_dim1, bb_offset, z_dim1, z_offset, i__1; 00029 00030 /* Local variables */ 00031 integer inde; 00032 char vect[1]; 00033 integer llwk2; 00034 extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, 00035 integer *, complex *, complex *, integer *, complex *, integer *, 00036 complex *, complex *, integer *); 00037 extern logical lsame_(char *, char *); 00038 integer iinfo, lwmin; 00039 logical upper; 00040 integer llrwk; 00041 logical wantz; 00042 integer indwk2; 00043 extern /* Subroutine */ int cstedc_(char *, integer *, real *, real *, 00044 complex *, integer *, complex *, integer *, real *, integer *, 00045 integer *, integer *, integer *), chbtrd_(char *, char *, 00046 integer *, integer *, complex *, integer *, real *, real *, 00047 complex *, integer *, complex *, integer *), 00048 chbgst_(char *, char *, integer *, integer *, integer *, complex * 00049 , integer *, complex *, integer *, complex *, integer *, complex * 00050 , real *, integer *), clacpy_(char *, integer *, 00051 integer *, complex *, integer *, complex *, integer *), 00052 xerbla_(char *, integer *), cpbstf_(char *, integer *, 00053 integer *, complex *, integer *, integer *); 00054 integer indwrk, liwmin; 00055 extern /* Subroutine */ int ssterf_(integer *, real *, real *, integer *); 00056 integer lrwmin; 00057 logical lquery; 00058 00059 00060 /* -- LAPACK driver routine (version 3.2) -- */ 00061 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00062 /* November 2006 */ 00063 00064 /* .. Scalar Arguments .. */ 00065 /* .. */ 00066 /* .. Array Arguments .. */ 00067 /* .. */ 00068 00069 /* Purpose */ 00070 /* ======= */ 00071 00072 /* CHBGVD computes all the eigenvalues, and optionally, the eigenvectors */ 00073 /* of a complex generalized Hermitian-definite banded eigenproblem, of */ 00074 /* the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian */ 00075 /* and banded, and B is also positive definite. If eigenvectors are */ 00076 /* desired, it uses a divide and conquer algorithm. */ 00077 00078 /* The divide and conquer algorithm makes very mild assumptions about */ 00079 /* floating point arithmetic. It will work on machines with a guard */ 00080 /* digit in add/subtract, or on those binary machines without guard */ 00081 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */ 00082 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */ 00083 /* without guard digits, but we know of none. */ 00084 00085 /* Arguments */ 00086 /* ========= */ 00087 00088 /* JOBZ (input) CHARACTER*1 */ 00089 /* = 'N': Compute eigenvalues only; */ 00090 /* = 'V': Compute eigenvalues and eigenvectors. */ 00091 00092 /* UPLO (input) CHARACTER*1 */ 00093 /* = 'U': Upper triangles of A and B are stored; */ 00094 /* = 'L': Lower triangles of A and B are stored. */ 00095 00096 /* N (input) INTEGER */ 00097 /* The order of the matrices A and B. N >= 0. */ 00098 00099 /* KA (input) INTEGER */ 00100 /* The number of superdiagonals of the matrix A if UPLO = 'U', */ 00101 /* or the number of subdiagonals if UPLO = 'L'. KA >= 0. */ 00102 00103 /* KB (input) INTEGER */ 00104 /* The number of superdiagonals of the matrix B if UPLO = 'U', */ 00105 /* or the number of subdiagonals if UPLO = 'L'. KB >= 0. */ 00106 00107 /* AB (input/output) COMPLEX array, dimension (LDAB, N) */ 00108 /* On entry, the upper or lower triangle of the Hermitian band */ 00109 /* matrix A, stored in the first ka+1 rows of the array. The */ 00110 /* j-th column of A is stored in the j-th column of the array AB */ 00111 /* as follows: */ 00112 /* if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; */ 00113 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). */ 00114 00115 /* On exit, the contents of AB are destroyed. */ 00116 00117 /* LDAB (input) INTEGER */ 00118 /* The leading dimension of the array AB. LDAB >= KA+1. */ 00119 00120 /* BB (input/output) COMPLEX array, dimension (LDBB, N) */ 00121 /* On entry, the upper or lower triangle of the Hermitian band */ 00122 /* matrix B, stored in the first kb+1 rows of the array. The */ 00123 /* j-th column of B is stored in the j-th column of the array BB */ 00124 /* as follows: */ 00125 /* if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; */ 00126 /* if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb). */ 00127 00128 /* On exit, the factor S from the split Cholesky factorization */ 00129 /* B = S**H*S, as returned by CPBSTF. */ 00130 00131 /* LDBB (input) INTEGER */ 00132 /* The leading dimension of the array BB. LDBB >= KB+1. */ 00133 00134 /* W (output) REAL array, dimension (N) */ 00135 /* If INFO = 0, the eigenvalues in ascending order. */ 00136 00137 /* Z (output) COMPLEX array, dimension (LDZ, N) */ 00138 /* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of */ 00139 /* eigenvectors, with the i-th column of Z holding the */ 00140 /* eigenvector associated with W(i). The eigenvectors are */ 00141 /* normalized so that Z**H*B*Z = I. */ 00142 /* If JOBZ = 'N', then Z is not referenced. */ 00143 00144 /* LDZ (input) INTEGER */ 00145 /* The leading dimension of the array Z. LDZ >= 1, and if */ 00146 /* JOBZ = 'V', LDZ >= N. */ 00147 00148 /* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */ 00149 /* On exit, if INFO=0, WORK(1) returns the optimal LWORK. */ 00150 00151 /* LWORK (input) INTEGER */ 00152 /* The dimension of the array WORK. */ 00153 /* If N <= 1, LWORK >= 1. */ 00154 /* If JOBZ = 'N' and N > 1, LWORK >= N. */ 00155 /* If JOBZ = 'V' and N > 1, LWORK >= 2*N**2. */ 00156 00157 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00158 /* only calculates the optimal sizes of the WORK, RWORK and */ 00159 /* IWORK arrays, returns these values as the first entries of */ 00160 /* the WORK, RWORK and IWORK arrays, and no error message */ 00161 /* related to LWORK or LRWORK or LIWORK is issued by XERBLA. */ 00162 00163 /* RWORK (workspace/output) REAL array, dimension (MAX(1,LRWORK)) */ 00164 /* On exit, if INFO=0, RWORK(1) returns the optimal LRWORK. */ 00165 00166 /* LRWORK (input) INTEGER */ 00167 /* The dimension of array RWORK. */ 00168 /* If N <= 1, LRWORK >= 1. */ 00169 /* If JOBZ = 'N' and N > 1, LRWORK >= N. */ 00170 /* If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2. */ 00171 00172 /* If LRWORK = -1, then a workspace query is assumed; the */ 00173 /* routine only calculates the optimal sizes of the WORK, RWORK */ 00174 /* and IWORK arrays, returns these values as the first entries */ 00175 /* of the WORK, RWORK and IWORK arrays, and no error message */ 00176 /* related to LWORK or LRWORK or LIWORK is issued by XERBLA. */ 00177 00178 /* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */ 00179 /* On exit, if INFO=0, IWORK(1) returns the optimal LIWORK. */ 00180 00181 /* LIWORK (input) INTEGER */ 00182 /* The dimension of array IWORK. */ 00183 /* If JOBZ = 'N' or N <= 1, LIWORK >= 1. */ 00184 /* If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. */ 00185 00186 /* If LIWORK = -1, then a workspace query is assumed; the */ 00187 /* routine only calculates the optimal sizes of the WORK, RWORK */ 00188 /* and IWORK arrays, returns these values as the first entries */ 00189 /* of the WORK, RWORK and IWORK arrays, and no error message */ 00190 /* related to LWORK or LRWORK or LIWORK is issued by XERBLA. */ 00191 00192 /* INFO (output) INTEGER */ 00193 /* = 0: successful exit */ 00194 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00195 /* > 0: if INFO = i, and i is: */ 00196 /* <= N: the algorithm failed to converge: */ 00197 /* i off-diagonal elements of an intermediate */ 00198 /* tridiagonal form did not converge to zero; */ 00199 /* > N: if INFO = N + i, for 1 <= i <= N, then CPBSTF */ 00200 /* returned INFO = i: B is not positive definite. */ 00201 /* The factorization of B could not be completed and */ 00202 /* no eigenvalues or eigenvectors were computed. */ 00203 00204 /* Further Details */ 00205 /* =============== */ 00206 00207 /* Based on contributions by */ 00208 /* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA */ 00209 00210 /* ===================================================================== */ 00211 00212 /* .. Parameters .. */ 00213 /* .. */ 00214 /* .. Local Scalars .. */ 00215 /* .. */ 00216 /* .. External Functions .. */ 00217 /* .. */ 00218 /* .. External Subroutines .. */ 00219 /* .. */ 00220 /* .. Executable Statements .. */ 00221 00222 /* Test the input parameters. */ 00223 00224 /* Parameter adjustments */ 00225 ab_dim1 = *ldab; 00226 ab_offset = 1 + ab_dim1; 00227 ab -= ab_offset; 00228 bb_dim1 = *ldbb; 00229 bb_offset = 1 + bb_dim1; 00230 bb -= bb_offset; 00231 --w; 00232 z_dim1 = *ldz; 00233 z_offset = 1 + z_dim1; 00234 z__ -= z_offset; 00235 --work; 00236 --rwork; 00237 --iwork; 00238 00239 /* Function Body */ 00240 wantz = lsame_(jobz, "V"); 00241 upper = lsame_(uplo, "U"); 00242 lquery = *lwork == -1 || *lrwork == -1 || *liwork == -1; 00243 00244 *info = 0; 00245 if (*n <= 1) { 00246 lwmin = 1; 00247 lrwmin = 1; 00248 liwmin = 1; 00249 } else if (wantz) { 00250 /* Computing 2nd power */ 00251 i__1 = *n; 00252 lwmin = i__1 * i__1 << 1; 00253 /* Computing 2nd power */ 00254 i__1 = *n; 00255 lrwmin = *n * 5 + 1 + (i__1 * i__1 << 1); 00256 liwmin = *n * 5 + 3; 00257 } else { 00258 lwmin = *n; 00259 lrwmin = *n; 00260 liwmin = 1; 00261 } 00262 if (! (wantz || lsame_(jobz, "N"))) { 00263 *info = -1; 00264 } else if (! (upper || lsame_(uplo, "L"))) { 00265 *info = -2; 00266 } else if (*n < 0) { 00267 *info = -3; 00268 } else if (*ka < 0) { 00269 *info = -4; 00270 } else if (*kb < 0 || *kb > *ka) { 00271 *info = -5; 00272 } else if (*ldab < *ka + 1) { 00273 *info = -7; 00274 } else if (*ldbb < *kb + 1) { 00275 *info = -9; 00276 } else if (*ldz < 1 || wantz && *ldz < *n) { 00277 *info = -12; 00278 } 00279 00280 if (*info == 0) { 00281 work[1].r = (real) lwmin, work[1].i = 0.f; 00282 rwork[1] = (real) lrwmin; 00283 iwork[1] = liwmin; 00284 00285 if (*lwork < lwmin && ! lquery) { 00286 *info = -14; 00287 } else if (*lrwork < lrwmin && ! lquery) { 00288 *info = -16; 00289 } else if (*liwork < liwmin && ! lquery) { 00290 *info = -18; 00291 } 00292 } 00293 00294 if (*info != 0) { 00295 i__1 = -(*info); 00296 xerbla_("CHBGVD", &i__1); 00297 return 0; 00298 } else if (lquery) { 00299 return 0; 00300 } 00301 00302 /* Quick return if possible */ 00303 00304 if (*n == 0) { 00305 return 0; 00306 } 00307 00308 /* Form a split Cholesky factorization of B. */ 00309 00310 cpbstf_(uplo, n, kb, &bb[bb_offset], ldbb, info); 00311 if (*info != 0) { 00312 *info = *n + *info; 00313 return 0; 00314 } 00315 00316 /* Transform problem to standard eigenvalue problem. */ 00317 00318 inde = 1; 00319 indwrk = inde + *n; 00320 indwk2 = *n * *n + 1; 00321 llwk2 = *lwork - indwk2 + 2; 00322 llrwk = *lrwork - indwrk + 2; 00323 chbgst_(jobz, uplo, n, ka, kb, &ab[ab_offset], ldab, &bb[bb_offset], ldbb, 00324 &z__[z_offset], ldz, &work[1], &rwork[indwrk], &iinfo); 00325 00326 /* Reduce Hermitian band matrix to tridiagonal form. */ 00327 00328 if (wantz) { 00329 *(unsigned char *)vect = 'U'; 00330 } else { 00331 *(unsigned char *)vect = 'N'; 00332 } 00333 chbtrd_(vect, uplo, n, ka, &ab[ab_offset], ldab, &w[1], &rwork[inde], & 00334 z__[z_offset], ldz, &work[1], &iinfo); 00335 00336 /* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC. */ 00337 00338 if (! wantz) { 00339 ssterf_(n, &w[1], &rwork[inde], info); 00340 } else { 00341 cstedc_("I", n, &w[1], &rwork[inde], &work[1], n, &work[indwk2], & 00342 llwk2, &rwork[indwrk], &llrwk, &iwork[1], liwork, info); 00343 cgemm_("N", "N", n, n, n, &c_b1, &z__[z_offset], ldz, &work[1], n, & 00344 c_b2, &work[indwk2], n); 00345 clacpy_("A", n, n, &work[indwk2], n, &z__[z_offset], ldz); 00346 } 00347 00348 work[1].r = (real) lwmin, work[1].i = 0.f; 00349 rwork[1] = (real) lrwmin; 00350 iwork[1] = liwmin; 00351 return 0; 00352 00353 /* End of CHBGVD */ 00354 00355 } /* chbgvd_ */