LAPACK 3.3.1 Linear Algebra PACKage

# chbgvd.f

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