dspgvx.c
Go to the documentation of this file.
00001 /* dspgvx.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 dspgvx_(integer *itype, char *jobz, char *range, char *
00021         uplo, integer *n, doublereal *ap, doublereal *bp, doublereal *vl, 
00022         doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer 
00023         *m, doublereal *w, doublereal *z__, integer *ldz, doublereal *work, 
00024         integer *iwork, integer *ifail, integer *info)
00025 {
00026     /* System generated locals */
00027     integer z_dim1, z_offset, i__1;
00028 
00029     /* Local variables */
00030     integer j;
00031     extern logical lsame_(char *, char *);
00032     char trans[1];
00033     logical upper;
00034     extern /* Subroutine */ int dtpmv_(char *, char *, char *, integer *, 
00035             doublereal *, doublereal *, integer *), 
00036             dtpsv_(char *, char *, char *, integer *, doublereal *, 
00037             doublereal *, integer *);
00038     logical wantz, alleig, indeig, valeig;
00039     extern /* Subroutine */ int xerbla_(char *, integer *), dpptrf_(
00040             char *, integer *, doublereal *, integer *), dspgst_(
00041             integer *, char *, integer *, doublereal *, doublereal *, integer 
00042             *), dspevx_(char *, char *, char *, integer *, doublereal 
00043             *, doublereal *, doublereal *, integer *, integer *, doublereal *, 
00044              integer *, doublereal *, doublereal *, integer *, doublereal *, 
00045             integer *, integer *, integer *);
00046 
00047 
00048 /*  -- LAPACK driver routine (version 3.2) -- */
00049 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00050 /*     November 2006 */
00051 
00052 /*     .. Scalar Arguments .. */
00053 /*     .. */
00054 /*     .. Array Arguments .. */
00055 /*     .. */
00056 
00057 /*  Purpose */
00058 /*  ======= */
00059 
00060 /*  DSPGVX computes selected eigenvalues, and optionally, eigenvectors */
00061 /*  of a real generalized symmetric-definite eigenproblem, of the form */
00062 /*  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A */
00063 /*  and B are assumed to be symmetric, stored in packed storage, and B */
00064 /*  is also positive definite.  Eigenvalues and eigenvectors can be */
00065 /*  selected by specifying either a range of values or a range of indices */
00066 /*  for the desired eigenvalues. */
00067 
00068 /*  Arguments */
00069 /*  ========= */
00070 
00071 /*  ITYPE   (input) INTEGER */
00072 /*          Specifies the problem type to be solved: */
00073 /*          = 1:  A*x = (lambda)*B*x */
00074 /*          = 2:  A*B*x = (lambda)*x */
00075 /*          = 3:  B*A*x = (lambda)*x */
00076 
00077 /*  JOBZ    (input) CHARACTER*1 */
00078 /*          = 'N':  Compute eigenvalues only; */
00079 /*          = 'V':  Compute eigenvalues and eigenvectors. */
00080 
00081 /*  RANGE   (input) CHARACTER*1 */
00082 /*          = 'A': all eigenvalues will be found. */
00083 /*          = 'V': all eigenvalues in the half-open interval (VL,VU] */
00084 /*                 will be found. */
00085 /*          = 'I': the IL-th through IU-th eigenvalues will be found. */
00086 
00087 /*  UPLO    (input) CHARACTER*1 */
00088 /*          = 'U':  Upper triangle of A and B are stored; */
00089 /*          = 'L':  Lower triangle of A and B are stored. */
00090 
00091 /*  N       (input) INTEGER */
00092 /*          The order of the matrix pencil (A,B).  N >= 0. */
00093 
00094 /*  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
00095 /*          On entry, the upper or lower triangle of the symmetric matrix */
00096 /*          A, packed columnwise in a linear array.  The j-th column of A */
00097 /*          is stored in the array AP as follows: */
00098 /*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
00099 /*          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */
00100 
00101 /*          On exit, the contents of AP are destroyed. */
00102 
00103 /*  BP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
00104 /*          On entry, the upper or lower triangle of the symmetric matrix */
00105 /*          B, packed columnwise in a linear array.  The j-th column of B */
00106 /*          is stored in the array BP as follows: */
00107 /*          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j; */
00108 /*          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n. */
00109 
00110 /*          On exit, the triangular factor U or L from the Cholesky */
00111 /*          factorization B = U**T*U or B = L*L**T, in the same storage */
00112 /*          format as B. */
00113 
00114 /*  VL      (input) DOUBLE PRECISION */
00115 /*  VU      (input) DOUBLE PRECISION */
00116 /*          If RANGE='V', the lower and upper bounds of the interval to */
00117 /*          be searched for eigenvalues. VL < VU. */
00118 /*          Not referenced if RANGE = 'A' or 'I'. */
00119 
00120 /*  IL      (input) INTEGER */
00121 /*  IU      (input) INTEGER */
00122 /*          If RANGE='I', the indices (in ascending order) of the */
00123 /*          smallest and largest eigenvalues to be returned. */
00124 /*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
00125 /*          Not referenced if RANGE = 'A' or 'V'. */
00126 
00127 /*  ABSTOL  (input) DOUBLE PRECISION */
00128 /*          The absolute error tolerance for the eigenvalues. */
00129 /*          An approximate eigenvalue is accepted as converged */
00130 /*          when it is determined to lie in an interval [a,b] */
00131 /*          of width less than or equal to */
00132 
00133 /*                  ABSTOL + EPS *   max( |a|,|b| ) , */
00134 
00135 /*          where EPS is the machine precision.  If ABSTOL is less than */
00136 /*          or equal to zero, then  EPS*|T|  will be used in its place, */
00137 /*          where |T| is the 1-norm of the tridiagonal matrix obtained */
00138 /*          by reducing A to tridiagonal form. */
00139 
00140 /*          Eigenvalues will be computed most accurately when ABSTOL is */
00141 /*          set to twice the underflow threshold 2*DLAMCH('S'), not zero. */
00142 /*          If this routine returns with INFO>0, indicating that some */
00143 /*          eigenvectors did not converge, try setting ABSTOL to */
00144 /*          2*DLAMCH('S'). */
00145 
00146 /*  M       (output) INTEGER */
00147 /*          The total number of eigenvalues found.  0 <= M <= N. */
00148 /*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
00149 
00150 /*  W       (output) DOUBLE PRECISION array, dimension (N) */
00151 /*          On normal exit, the first M elements contain the selected */
00152 /*          eigenvalues in ascending order. */
00153 
00154 /*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M)) */
00155 /*          If JOBZ = 'N', then Z is not referenced. */
00156 /*          If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
00157 /*          contain the orthonormal eigenvectors of the matrix A */
00158 /*          corresponding to the selected eigenvalues, with the i-th */
00159 /*          column of Z holding the eigenvector associated with W(i). */
00160 /*          The eigenvectors are normalized as follows: */
00161 /*          if ITYPE = 1 or 2, Z**T*B*Z = I; */
00162 /*          if ITYPE = 3, Z**T*inv(B)*Z = I. */
00163 
00164 /*          If an eigenvector fails to converge, then that column of Z */
00165 /*          contains the latest approximation to the eigenvector, and the */
00166 /*          index of the eigenvector is returned in IFAIL. */
00167 /*          Note: the user must ensure that at least max(1,M) columns are */
00168 /*          supplied in the array Z; if RANGE = 'V', the exact value of M */
00169 /*          is not known in advance and an upper bound must be used. */
00170 
00171 /*  LDZ     (input) INTEGER */
00172 /*          The leading dimension of the array Z.  LDZ >= 1, and if */
00173 /*          JOBZ = 'V', LDZ >= max(1,N). */
00174 
00175 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (8*N) */
00176 
00177 /*  IWORK   (workspace) INTEGER array, dimension (5*N) */
00178 
00179 /*  IFAIL   (output) INTEGER array, dimension (N) */
00180 /*          If JOBZ = 'V', then if INFO = 0, the first M elements of */
00181 /*          IFAIL are zero.  If INFO > 0, then IFAIL contains the */
00182 /*          indices of the eigenvectors that failed to converge. */
00183 /*          If JOBZ = 'N', then IFAIL is not referenced. */
00184 
00185 /*  INFO    (output) INTEGER */
00186 /*          = 0:  successful exit */
00187 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00188 /*          > 0:  DPPTRF or DSPEVX returned an error code: */
00189 /*             <= N:  if INFO = i, DSPEVX failed to converge; */
00190 /*                    i eigenvectors failed to converge.  Their indices */
00191 /*                    are stored in array IFAIL. */
00192 /*             > N:   if INFO = N + i, for 1 <= i <= N, then the leading */
00193 /*                    minor of order i of B is not positive definite. */
00194 /*                    The factorization of B could not be completed and */
00195 /*                    no eigenvalues or eigenvectors were computed. */
00196 
00197 /*  Further Details */
00198 /*  =============== */
00199 
00200 /*  Based on contributions by */
00201 /*     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA */
00202 
00203 /* ===================================================================== */
00204 
00205 /*     .. Local Scalars .. */
00206 /*     .. */
00207 /*     .. External Functions .. */
00208 /*     .. */
00209 /*     .. External Subroutines .. */
00210 /*     .. */
00211 /*     .. Intrinsic Functions .. */
00212 /*     .. */
00213 /*     .. Executable Statements .. */
00214 
00215 /*     Test the input parameters. */
00216 
00217     /* Parameter adjustments */
00218     --ap;
00219     --bp;
00220     --w;
00221     z_dim1 = *ldz;
00222     z_offset = 1 + z_dim1;
00223     z__ -= z_offset;
00224     --work;
00225     --iwork;
00226     --ifail;
00227 
00228     /* Function Body */
00229     upper = lsame_(uplo, "U");
00230     wantz = lsame_(jobz, "V");
00231     alleig = lsame_(range, "A");
00232     valeig = lsame_(range, "V");
00233     indeig = lsame_(range, "I");
00234 
00235     *info = 0;
00236     if (*itype < 1 || *itype > 3) {
00237         *info = -1;
00238     } else if (! (wantz || lsame_(jobz, "N"))) {
00239         *info = -2;
00240     } else if (! (alleig || valeig || indeig)) {
00241         *info = -3;
00242     } else if (! (upper || lsame_(uplo, "L"))) {
00243         *info = -4;
00244     } else if (*n < 0) {
00245         *info = -5;
00246     } else {
00247         if (valeig) {
00248             if (*n > 0 && *vu <= *vl) {
00249                 *info = -9;
00250             }
00251         } else if (indeig) {
00252             if (*il < 1) {
00253                 *info = -10;
00254             } else if (*iu < min(*n,*il) || *iu > *n) {
00255                 *info = -11;
00256             }
00257         }
00258     }
00259     if (*info == 0) {
00260         if (*ldz < 1 || wantz && *ldz < *n) {
00261             *info = -16;
00262         }
00263     }
00264 
00265     if (*info != 0) {
00266         i__1 = -(*info);
00267         xerbla_("DSPGVX", &i__1);
00268         return 0;
00269     }
00270 
00271 /*     Quick return if possible */
00272 
00273     *m = 0;
00274     if (*n == 0) {
00275         return 0;
00276     }
00277 
00278 /*     Form a Cholesky factorization of B. */
00279 
00280     dpptrf_(uplo, n, &bp[1], info);
00281     if (*info != 0) {
00282         *info = *n + *info;
00283         return 0;
00284     }
00285 
00286 /*     Transform problem to standard eigenvalue problem and solve. */
00287 
00288     dspgst_(itype, uplo, n, &ap[1], &bp[1], info);
00289     dspevx_(jobz, range, uplo, n, &ap[1], vl, vu, il, iu, abstol, m, &w[1], &
00290             z__[z_offset], ldz, &work[1], &iwork[1], &ifail[1], info);
00291 
00292     if (wantz) {
00293 
00294 /*        Backtransform eigenvectors to the original problem. */
00295 
00296         if (*info > 0) {
00297             *m = *info - 1;
00298         }
00299         if (*itype == 1 || *itype == 2) {
00300 
00301 /*           For A*x=(lambda)*B*x and A*B*x=(lambda)*x; */
00302 /*           backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */
00303 
00304             if (upper) {
00305                 *(unsigned char *)trans = 'N';
00306             } else {
00307                 *(unsigned char *)trans = 'T';
00308             }
00309 
00310             i__1 = *m;
00311             for (j = 1; j <= i__1; ++j) {
00312                 dtpsv_(uplo, trans, "Non-unit", n, &bp[1], &z__[j * z_dim1 + 
00313                         1], &c__1);
00314 /* L10: */
00315             }
00316 
00317         } else if (*itype == 3) {
00318 
00319 /*           For B*A*x=(lambda)*x; */
00320 /*           backtransform eigenvectors: x = L*y or U'*y */
00321 
00322             if (upper) {
00323                 *(unsigned char *)trans = 'T';
00324             } else {
00325                 *(unsigned char *)trans = 'N';
00326             }
00327 
00328             i__1 = *m;
00329             for (j = 1; j <= i__1; ++j) {
00330                 dtpmv_(uplo, trans, "Non-unit", n, &bp[1], &z__[j * z_dim1 + 
00331                         1], &c__1);
00332 /* L20: */
00333             }
00334         }
00335     }
00336 
00337     return 0;
00338 
00339 /*     End of DSPGVX */
00340 
00341 } /* dspgvx_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:48