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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:13