LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsbgvd.f
Go to the documentation of this file.
00001 *> \brief \b DSBGST
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DSBGVD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbgvd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbgvd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbgvd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
00022 *                          Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBZ, UPLO
00026 *       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * )
00030 *       DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), W( * ),
00031 *      $                   WORK( * ), Z( LDZ, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DSBGVD computes all the eigenvalues, and optionally, the eigenvectors
00041 *> of a real generalized symmetric-definite banded eigenproblem, of the
00042 *> form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric and
00043 *> banded, and B is also positive definite.  If eigenvectors are
00044 *> desired, it uses a divide and conquer algorithm.
00045 *>
00046 *> The divide and conquer algorithm makes very mild assumptions about
00047 *> floating point arithmetic. It will work on machines with a guard
00048 *> digit in add/subtract, or on those binary machines without guard
00049 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00050 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
00051 *> without guard digits, but we know of none.
00052 *> \endverbatim
00053 *
00054 *  Arguments:
00055 *  ==========
00056 *
00057 *> \param[in] JOBZ
00058 *> \verbatim
00059 *>          JOBZ is CHARACTER*1
00060 *>          = 'N':  Compute eigenvalues only;
00061 *>          = 'V':  Compute eigenvalues and eigenvectors.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] UPLO
00065 *> \verbatim
00066 *>          UPLO is CHARACTER*1
00067 *>          = 'U':  Upper triangles of A and B are stored;
00068 *>          = 'L':  Lower triangles of A and B are stored.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] N
00072 *> \verbatim
00073 *>          N is INTEGER
00074 *>          The order of the matrices A and B.  N >= 0.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] KA
00078 *> \verbatim
00079 *>          KA is INTEGER
00080 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
00081 *>          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] KB
00085 *> \verbatim
00086 *>          KB is INTEGER
00087 *>          The number of superdiagonals of the matrix B if UPLO = 'U',
00088 *>          or the number of subdiagonals if UPLO = 'L'.  KB >= 0.
00089 *> \endverbatim
00090 *>
00091 *> \param[in,out] AB
00092 *> \verbatim
00093 *>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
00094 *>          On entry, the upper or lower triangle of the symmetric band
00095 *>          matrix A, stored in the first ka+1 rows of the array.  The
00096 *>          j-th column of A is stored in the j-th column of the array AB
00097 *>          as follows:
00098 *>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
00099 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
00100 *>
00101 *>          On exit, the contents of AB are destroyed.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] LDAB
00105 *> \verbatim
00106 *>          LDAB is INTEGER
00107 *>          The leading dimension of the array AB.  LDAB >= KA+1.
00108 *> \endverbatim
00109 *>
00110 *> \param[in,out] BB
00111 *> \verbatim
00112 *>          BB is DOUBLE PRECISION array, dimension (LDBB, N)
00113 *>          On entry, the upper or lower triangle of the symmetric band
00114 *>          matrix B, stored in the first kb+1 rows of the array.  The
00115 *>          j-th column of B is stored in the j-th column of the array BB
00116 *>          as follows:
00117 *>          if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
00118 *>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
00119 *>
00120 *>          On exit, the factor S from the split Cholesky factorization
00121 *>          B = S**T*S, as returned by DPBSTF.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDBB
00125 *> \verbatim
00126 *>          LDBB is INTEGER
00127 *>          The leading dimension of the array BB.  LDBB >= KB+1.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] W
00131 *> \verbatim
00132 *>          W is DOUBLE PRECISION array, dimension (N)
00133 *>          If INFO = 0, the eigenvalues in ascending order.
00134 *> \endverbatim
00135 *>
00136 *> \param[out] Z
00137 *> \verbatim
00138 *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
00139 *>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
00140 *>          eigenvectors, with the i-th column of Z holding the
00141 *>          eigenvector associated with W(i).  The eigenvectors are
00142 *>          normalized so Z**T*B*Z = I.
00143 *>          If JOBZ = 'N', then Z is not referenced.
00144 *> \endverbatim
00145 *>
00146 *> \param[in] LDZ
00147 *> \verbatim
00148 *>          LDZ is INTEGER
00149 *>          The leading dimension of the array Z.  LDZ >= 1, and if
00150 *>          JOBZ = 'V', LDZ >= max(1,N).
00151 *> \endverbatim
00152 *>
00153 *> \param[out] WORK
00154 *> \verbatim
00155 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00156 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] LWORK
00160 *> \verbatim
00161 *>          LWORK is INTEGER
00162 *>          The dimension of the array WORK.
00163 *>          If N <= 1,               LWORK >= 1.
00164 *>          If JOBZ = 'N' and N > 1, LWORK >= 3*N.
00165 *>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2.
00166 *>
00167 *>          If LWORK = -1, then a workspace query is assumed; the routine
00168 *>          only calculates the optimal sizes of the WORK and IWORK
00169 *>          arrays, returns these values as the first entries of the WORK
00170 *>          and IWORK arrays, and no error message related to LWORK or
00171 *>          LIWORK is issued by XERBLA.
00172 *> \endverbatim
00173 *>
00174 *> \param[out] IWORK
00175 *> \verbatim
00176 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
00177 *>          On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
00178 *> \endverbatim
00179 *>
00180 *> \param[in] LIWORK
00181 *> \verbatim
00182 *>          LIWORK is INTEGER
00183 *>          The dimension of the array IWORK.
00184 *>          If JOBZ  = 'N' or N <= 1, LIWORK >= 1.
00185 *>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
00186 *>
00187 *>          If LIWORK = -1, then a workspace query is assumed; the
00188 *>          routine only calculates the optimal sizes of the WORK and
00189 *>          IWORK arrays, returns these values as the first entries of
00190 *>          the WORK and IWORK arrays, and no error message related to
00191 *>          LWORK or LIWORK is issued by XERBLA.
00192 *> \endverbatim
00193 *>
00194 *> \param[out] INFO
00195 *> \verbatim
00196 *>          INFO is INTEGER
00197 *>          = 0:  successful exit
00198 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00199 *>          > 0:  if INFO = i, and i is:
00200 *>             <= N:  the algorithm failed to converge:
00201 *>                    i off-diagonal elements of an intermediate
00202 *>                    tridiagonal form did not converge to zero;
00203 *>             > N:   if INFO = N + i, for 1 <= i <= N, then DPBSTF
00204 *>                    returned INFO = i: B is not positive definite.
00205 *>                    The factorization of B could not be completed and
00206 *>                    no eigenvalues or eigenvectors were computed.
00207 *> \endverbatim
00208 *
00209 *  Authors:
00210 *  ========
00211 *
00212 *> \author Univ. of Tennessee 
00213 *> \author Univ. of California Berkeley 
00214 *> \author Univ. of Colorado Denver 
00215 *> \author NAG Ltd. 
00216 *
00217 *> \date November 2011
00218 *
00219 *> \ingroup doubleOTHEReigen
00220 *
00221 *> \par Contributors:
00222 *  ==================
00223 *>
00224 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
00225 *
00226 *  =====================================================================
00227       SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
00228      $                   Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO )
00229 *
00230 *  -- LAPACK driver routine (version 3.4.0) --
00231 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00232 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00233 *     November 2011
00234 *
00235 *     .. Scalar Arguments ..
00236       CHARACTER          JOBZ, UPLO
00237       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
00238 *     ..
00239 *     .. Array Arguments ..
00240       INTEGER            IWORK( * )
00241       DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), W( * ),
00242      $                   WORK( * ), Z( LDZ, * )
00243 *     ..
00244 *
00245 *  =====================================================================
00246 *
00247 *     .. Parameters ..
00248       DOUBLE PRECISION   ONE, ZERO
00249       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00250 *     ..
00251 *     .. Local Scalars ..
00252       LOGICAL            LQUERY, UPPER, WANTZ
00253       CHARACTER          VECT
00254       INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2,
00255      $                   LWMIN
00256 *     ..
00257 *     .. External Functions ..
00258       LOGICAL            LSAME
00259       EXTERNAL           LSAME
00260 *     ..
00261 *     .. External Subroutines ..
00262       EXTERNAL           DGEMM, DLACPY, DPBSTF, DSBGST, DSBTRD, DSTEDC,
00263      $                   DSTERF, XERBLA
00264 *     ..
00265 *     .. Executable Statements ..
00266 *
00267 *     Test the input parameters.
00268 *
00269       WANTZ = LSAME( JOBZ, 'V' )
00270       UPPER = LSAME( UPLO, 'U' )
00271       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00272 *
00273       INFO = 0
00274       IF( N.LE.1 ) THEN
00275          LIWMIN = 1
00276          LWMIN = 1
00277       ELSE IF( WANTZ ) THEN
00278          LIWMIN = 3 + 5*N
00279          LWMIN = 1 + 5*N + 2*N**2
00280       ELSE
00281          LIWMIN = 1
00282          LWMIN = 2*N
00283       END IF
00284 *
00285       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00286          INFO = -1
00287       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
00288          INFO = -2
00289       ELSE IF( N.LT.0 ) THEN
00290          INFO = -3
00291       ELSE IF( KA.LT.0 ) THEN
00292          INFO = -4
00293       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
00294          INFO = -5
00295       ELSE IF( LDAB.LT.KA+1 ) THEN
00296          INFO = -7
00297       ELSE IF( LDBB.LT.KB+1 ) THEN
00298          INFO = -9
00299       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
00300          INFO = -12
00301       END IF
00302 *
00303       IF( INFO.EQ.0 ) THEN
00304          WORK( 1 ) = LWMIN
00305          IWORK( 1 ) = LIWMIN
00306 *
00307          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00308             INFO = -14
00309          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00310             INFO = -16
00311          END IF
00312       END IF
00313 *
00314       IF( INFO.NE.0 ) THEN
00315          CALL XERBLA( 'DSBGVD', -INFO )
00316          RETURN
00317       ELSE IF( LQUERY ) THEN
00318          RETURN
00319       END IF
00320 *
00321 *     Quick return if possible
00322 *
00323       IF( N.EQ.0 )
00324      $   RETURN
00325 *
00326 *     Form a split Cholesky factorization of B.
00327 *
00328       CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO )
00329       IF( INFO.NE.0 ) THEN
00330          INFO = N + INFO
00331          RETURN
00332       END IF
00333 *
00334 *     Transform problem to standard eigenvalue problem.
00335 *
00336       INDE = 1
00337       INDWRK = INDE + N
00338       INDWK2 = INDWRK + N*N
00339       LLWRK2 = LWORK - INDWK2 + 1
00340       CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
00341      $             WORK( INDWRK ), IINFO )
00342 *
00343 *     Reduce to tridiagonal form.
00344 *
00345       IF( WANTZ ) THEN
00346          VECT = 'U'
00347       ELSE
00348          VECT = 'N'
00349       END IF
00350       CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ,
00351      $             WORK( INDWRK ), IINFO )
00352 *
00353 *     For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
00354 *
00355       IF( .NOT.WANTZ ) THEN
00356          CALL DSTERF( N, W, WORK( INDE ), INFO )
00357       ELSE
00358          CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
00359      $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
00360          CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
00361      $               ZERO, WORK( INDWK2 ), N )
00362          CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
00363       END IF
00364 *
00365       WORK( 1 ) = LWMIN
00366       IWORK( 1 ) = LIWMIN
00367 *
00368       RETURN
00369 *
00370 *     End of DSBGVD
00371 *
00372       END
 All Files Functions