00001 /* chpgvd.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 integer c__1 = 1; 00019 00020 /* Subroutine */ int chpgvd_(integer *itype, char *jobz, char *uplo, integer * 00021 n, complex *ap, complex *bp, real *w, complex *z__, integer *ldz, 00022 complex *work, integer *lwork, real *rwork, integer *lrwork, integer * 00023 iwork, integer *liwork, integer *info) 00024 { 00025 /* System generated locals */ 00026 integer z_dim1, z_offset, i__1; 00027 real r__1, r__2; 00028 00029 /* Local variables */ 00030 integer j, neig; 00031 extern logical lsame_(char *, char *); 00032 integer lwmin; 00033 char trans[1]; 00034 extern /* Subroutine */ int ctpmv_(char *, char *, char *, integer *, 00035 complex *, complex *, integer *); 00036 logical upper; 00037 extern /* Subroutine */ int ctpsv_(char *, char *, char *, integer *, 00038 complex *, complex *, integer *); 00039 logical wantz; 00040 extern /* Subroutine */ int chpevd_(char *, char *, integer *, complex *, 00041 real *, complex *, integer *, complex *, integer *, real *, 00042 integer *, integer *, integer *, integer *), 00043 xerbla_(char *, integer *), chpgst_(integer *, char *, 00044 integer *, complex *, complex *, integer *), cpptrf_(char 00045 *, integer *, complex *, integer *); 00046 integer liwmin, lrwmin; 00047 logical lquery; 00048 00049 00050 /* -- LAPACK driver routine (version 3.2) -- */ 00051 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00052 /* November 2006 */ 00053 00054 /* .. Scalar Arguments .. */ 00055 /* .. */ 00056 /* .. Array Arguments .. */ 00057 /* .. */ 00058 00059 /* Purpose */ 00060 /* ======= */ 00061 00062 /* CHPGVD computes all the eigenvalues and, optionally, the eigenvectors */ 00063 /* of a complex generalized Hermitian-definite eigenproblem, of the form */ 00064 /* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and */ 00065 /* B are assumed to be Hermitian, stored in packed format, and B is also */ 00066 /* positive definite. */ 00067 /* If eigenvectors are desired, it uses a divide and conquer algorithm. */ 00068 00069 /* The divide and conquer algorithm makes very mild assumptions about */ 00070 /* floating point arithmetic. It will work on machines with a guard */ 00071 /* digit in add/subtract, or on those binary machines without guard */ 00072 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */ 00073 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */ 00074 /* without guard digits, but we know of none. */ 00075 00076 /* Arguments */ 00077 /* ========= */ 00078 00079 /* ITYPE (input) INTEGER */ 00080 /* Specifies the problem type to be solved: */ 00081 /* = 1: A*x = (lambda)*B*x */ 00082 /* = 2: A*B*x = (lambda)*x */ 00083 /* = 3: B*A*x = (lambda)*x */ 00084 00085 /* JOBZ (input) CHARACTER*1 */ 00086 /* = 'N': Compute eigenvalues only; */ 00087 /* = 'V': Compute eigenvalues and eigenvectors. */ 00088 00089 /* UPLO (input) CHARACTER*1 */ 00090 /* = 'U': Upper triangles of A and B are stored; */ 00091 /* = 'L': Lower triangles of A and B are stored. */ 00092 00093 /* N (input) INTEGER */ 00094 /* The order of the matrices A and B. N >= 0. */ 00095 00096 /* AP (input/output) COMPLEX array, dimension (N*(N+1)/2) */ 00097 /* On entry, the upper or lower triangle of the Hermitian matrix */ 00098 /* A, packed columnwise in a linear array. The j-th column of A */ 00099 /* is stored in the array AP as follows: */ 00100 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00101 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */ 00102 00103 /* On exit, the contents of AP are destroyed. */ 00104 00105 /* BP (input/output) COMPLEX array, dimension (N*(N+1)/2) */ 00106 /* On entry, the upper or lower triangle of the Hermitian matrix */ 00107 /* B, packed columnwise in a linear array. The j-th column of B */ 00108 /* is stored in the array BP as follows: */ 00109 /* if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j; */ 00110 /* if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n. */ 00111 00112 /* On exit, the triangular factor U or L from the Cholesky */ 00113 /* factorization B = U**H*U or B = L*L**H, in the same storage */ 00114 /* format as B. */ 00115 00116 /* W (output) REAL array, dimension (N) */ 00117 /* If INFO = 0, the eigenvalues in ascending order. */ 00118 00119 /* Z (output) COMPLEX array, dimension (LDZ, N) */ 00120 /* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of */ 00121 /* eigenvectors. The eigenvectors are normalized as follows: */ 00122 /* if ITYPE = 1 or 2, Z**H*B*Z = I; */ 00123 /* if ITYPE = 3, Z**H*inv(B)*Z = I. */ 00124 /* If JOBZ = 'N', then Z is not referenced. */ 00125 00126 /* LDZ (input) INTEGER */ 00127 /* The leading dimension of the array Z. LDZ >= 1, and if */ 00128 /* JOBZ = 'V', LDZ >= max(1,N). */ 00129 00130 /* WORK (workspace) COMPLEX array, dimension (MAX(1,LWORK)) */ 00131 /* On exit, if INFO = 0, WORK(1) returns the required LWORK. */ 00132 00133 /* LWORK (input) INTEGER */ 00134 /* The dimension of array WORK. */ 00135 /* If N <= 1, LWORK >= 1. */ 00136 /* If JOBZ = 'N' and N > 1, LWORK >= N. */ 00137 /* If JOBZ = 'V' and N > 1, LWORK >= 2*N. */ 00138 00139 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00140 /* only calculates the required sizes of the WORK, RWORK and */ 00141 /* IWORK arrays, returns these values as the first entries of */ 00142 /* the WORK, RWORK and IWORK arrays, and no error message */ 00143 /* related to LWORK or LRWORK or LIWORK is issued by XERBLA. */ 00144 00145 /* RWORK (workspace) REAL array, dimension (MAX(1,LRWORK)) */ 00146 /* On exit, if INFO = 0, RWORK(1) returns the required LRWORK. */ 00147 00148 /* LRWORK (input) INTEGER */ 00149 /* The dimension of array RWORK. */ 00150 /* If N <= 1, LRWORK >= 1. */ 00151 /* If JOBZ = 'N' and N > 1, LRWORK >= N. */ 00152 /* If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2. */ 00153 00154 /* If LRWORK = -1, then a workspace query is assumed; the */ 00155 /* routine only calculates the required sizes of the WORK, RWORK */ 00156 /* and IWORK arrays, returns these values as the first entries */ 00157 /* of the WORK, RWORK and IWORK arrays, and no error message */ 00158 /* related to LWORK or LRWORK or LIWORK is issued by XERBLA. */ 00159 00160 /* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */ 00161 /* On exit, if INFO = 0, IWORK(1) returns the required LIWORK. */ 00162 00163 /* LIWORK (input) INTEGER */ 00164 /* The dimension of array IWORK. */ 00165 /* If JOBZ = 'N' or N <= 1, LIWORK >= 1. */ 00166 /* If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. */ 00167 00168 /* If LIWORK = -1, then a workspace query is assumed; the */ 00169 /* routine only calculates the required sizes of the WORK, RWORK */ 00170 /* and IWORK arrays, returns these values as the first entries */ 00171 /* of the WORK, RWORK and IWORK arrays, and no error message */ 00172 /* related to LWORK or LRWORK or LIWORK is issued by XERBLA. */ 00173 00174 /* INFO (output) INTEGER */ 00175 /* = 0: successful exit */ 00176 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00177 /* > 0: CPPTRF or CHPEVD returned an error code: */ 00178 /* <= N: if INFO = i, CHPEVD failed to converge; */ 00179 /* i off-diagonal elements of an intermediate */ 00180 /* tridiagonal form did not convergeto zero; */ 00181 /* > N: if INFO = N + i, for 1 <= i <= n, then the leading */ 00182 /* minor of order i of B is not positive definite. */ 00183 /* The factorization of B could not be completed and */ 00184 /* no eigenvalues or eigenvectors were computed. */ 00185 00186 /* Further Details */ 00187 /* =============== */ 00188 00189 /* Based on contributions by */ 00190 /* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA */ 00191 00192 /* ===================================================================== */ 00193 00194 /* .. Local Scalars .. */ 00195 /* .. */ 00196 /* .. External Functions .. */ 00197 /* .. */ 00198 /* .. External Subroutines .. */ 00199 /* .. */ 00200 /* .. Intrinsic Functions .. */ 00201 /* .. */ 00202 /* .. Executable Statements .. */ 00203 00204 /* Test the input parameters. */ 00205 00206 /* Parameter adjustments */ 00207 --ap; 00208 --bp; 00209 --w; 00210 z_dim1 = *ldz; 00211 z_offset = 1 + z_dim1; 00212 z__ -= z_offset; 00213 --work; 00214 --rwork; 00215 --iwork; 00216 00217 /* Function Body */ 00218 wantz = lsame_(jobz, "V"); 00219 upper = lsame_(uplo, "U"); 00220 lquery = *lwork == -1 || *lrwork == -1 || *liwork == -1; 00221 00222 *info = 0; 00223 if (*itype < 1 || *itype > 3) { 00224 *info = -1; 00225 } else if (! (wantz || lsame_(jobz, "N"))) { 00226 *info = -2; 00227 } else if (! (upper || lsame_(uplo, "L"))) { 00228 *info = -3; 00229 } else if (*n < 0) { 00230 *info = -4; 00231 } else if (*ldz < 1 || wantz && *ldz < *n) { 00232 *info = -9; 00233 } 00234 00235 if (*info == 0) { 00236 if (*n <= 1) { 00237 lwmin = 1; 00238 liwmin = 1; 00239 lrwmin = 1; 00240 } else { 00241 if (wantz) { 00242 lwmin = *n << 1; 00243 /* Computing 2nd power */ 00244 i__1 = *n; 00245 lrwmin = *n * 5 + 1 + (i__1 * i__1 << 1); 00246 liwmin = *n * 5 + 3; 00247 } else { 00248 lwmin = *n; 00249 lrwmin = *n; 00250 liwmin = 1; 00251 } 00252 } 00253 work[1].r = (real) lwmin, work[1].i = 0.f; 00254 rwork[1] = (real) lrwmin; 00255 iwork[1] = liwmin; 00256 00257 if (*lwork < lwmin && ! lquery) { 00258 *info = -11; 00259 } else if (*lrwork < lrwmin && ! lquery) { 00260 *info = -13; 00261 } else if (*liwork < liwmin && ! lquery) { 00262 *info = -15; 00263 } 00264 } 00265 00266 if (*info != 0) { 00267 i__1 = -(*info); 00268 xerbla_("CHPGVD", &i__1); 00269 return 0; 00270 } else if (lquery) { 00271 return 0; 00272 } 00273 00274 /* Quick return if possible */ 00275 00276 if (*n == 0) { 00277 return 0; 00278 } 00279 00280 /* Form a Cholesky factorization of B. */ 00281 00282 cpptrf_(uplo, n, &bp[1], info); 00283 if (*info != 0) { 00284 *info = *n + *info; 00285 return 0; 00286 } 00287 00288 /* Transform problem to standard eigenvalue problem and solve. */ 00289 00290 chpgst_(itype, uplo, n, &ap[1], &bp[1], info); 00291 chpevd_(jobz, uplo, n, &ap[1], &w[1], &z__[z_offset], ldz, &work[1], 00292 lwork, &rwork[1], lrwork, &iwork[1], liwork, info); 00293 /* Computing MAX */ 00294 r__1 = (real) lwmin, r__2 = work[1].r; 00295 lwmin = dmax(r__1,r__2); 00296 /* Computing MAX */ 00297 r__1 = (real) lrwmin; 00298 lrwmin = dmax(r__1,rwork[1]); 00299 /* Computing MAX */ 00300 r__1 = (real) liwmin, r__2 = (real) iwork[1]; 00301 liwmin = dmax(r__1,r__2); 00302 00303 if (wantz) { 00304 00305 /* Backtransform eigenvectors to the original problem. */ 00306 00307 neig = *n; 00308 if (*info > 0) { 00309 neig = *info - 1; 00310 } 00311 if (*itype == 1 || *itype == 2) { 00312 00313 /* For A*x=(lambda)*B*x and A*B*x=(lambda)*x; */ 00314 /* backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */ 00315 00316 if (upper) { 00317 *(unsigned char *)trans = 'N'; 00318 } else { 00319 *(unsigned char *)trans = 'C'; 00320 } 00321 00322 i__1 = neig; 00323 for (j = 1; j <= i__1; ++j) { 00324 ctpsv_(uplo, trans, "Non-unit", n, &bp[1], &z__[j * z_dim1 + 00325 1], &c__1); 00326 /* L10: */ 00327 } 00328 00329 } else if (*itype == 3) { 00330 00331 /* For B*A*x=(lambda)*x; */ 00332 /* backtransform eigenvectors: x = L*y or U'*y */ 00333 00334 if (upper) { 00335 *(unsigned char *)trans = 'C'; 00336 } else { 00337 *(unsigned char *)trans = 'N'; 00338 } 00339 00340 i__1 = neig; 00341 for (j = 1; j <= i__1; ++j) { 00342 ctpmv_(uplo, trans, "Non-unit", n, &bp[1], &z__[j * z_dim1 + 00343 1], &c__1); 00344 /* L20: */ 00345 } 00346 } 00347 } 00348 00349 work[1].r = (real) lwmin, work[1].i = 0.f; 00350 rwork[1] = (real) lrwmin; 00351 iwork[1] = liwmin; 00352 return 0; 00353 00354 /* End of CHPGVD */ 00355 00356 } /* chpgvd_ */