LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ssygvd.f
Go to the documentation of this file.
00001 *> \brief \b SSYGST
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SSYGVD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssygvd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssygvd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssygvd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
00022 *                          LWORK, IWORK, LIWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBZ, UPLO
00026 *       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * )
00030 *       REAL               A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> SSYGVD computes all the eigenvalues, and optionally, the eigenvectors
00040 *> of a real generalized symmetric-definite eigenproblem, of the form
00041 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
00042 *> B are assumed to be symmetric and B is also positive definite.
00043 *> If eigenvectors are desired, it uses a divide and conquer algorithm.
00044 *>
00045 *> The divide and conquer algorithm makes very mild assumptions about
00046 *> floating point arithmetic. It will work on machines with a guard
00047 *> digit in add/subtract, or on those binary machines without guard
00048 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00049 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
00050 *> without guard digits, but we know of none.
00051 *> \endverbatim
00052 *
00053 *  Arguments:
00054 *  ==========
00055 *
00056 *> \param[in] ITYPE
00057 *> \verbatim
00058 *>          ITYPE is INTEGER
00059 *>          Specifies the problem type to be solved:
00060 *>          = 1:  A*x = (lambda)*B*x
00061 *>          = 2:  A*B*x = (lambda)*x
00062 *>          = 3:  B*A*x = (lambda)*x
00063 *> \endverbatim
00064 *>
00065 *> \param[in] JOBZ
00066 *> \verbatim
00067 *>          JOBZ is CHARACTER*1
00068 *>          = 'N':  Compute eigenvalues only;
00069 *>          = 'V':  Compute eigenvalues and eigenvectors.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] UPLO
00073 *> \verbatim
00074 *>          UPLO is CHARACTER*1
00075 *>          = 'U':  Upper triangles of A and B are stored;
00076 *>          = 'L':  Lower triangles of A and B are stored.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] N
00080 *> \verbatim
00081 *>          N is INTEGER
00082 *>          The order of the matrices A and B.  N >= 0.
00083 *> \endverbatim
00084 *>
00085 *> \param[in,out] A
00086 *> \verbatim
00087 *>          A is REAL array, dimension (LDA, N)
00088 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the
00089 *>          leading N-by-N upper triangular part of A contains the
00090 *>          upper triangular part of the matrix A.  If UPLO = 'L',
00091 *>          the leading N-by-N lower triangular part of A contains
00092 *>          the lower triangular part of the matrix A.
00093 *>
00094 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
00095 *>          matrix Z of eigenvectors.  The eigenvectors are normalized
00096 *>          as follows:
00097 *>          if ITYPE = 1 or 2, Z**T*B*Z = I;
00098 *>          if ITYPE = 3, Z**T*inv(B)*Z = I.
00099 *>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
00100 *>          or the lower triangle (if UPLO='L') of A, including the
00101 *>          diagonal, is destroyed.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] LDA
00105 *> \verbatim
00106 *>          LDA is INTEGER
00107 *>          The leading dimension of the array A.  LDA >= max(1,N).
00108 *> \endverbatim
00109 *>
00110 *> \param[in,out] B
00111 *> \verbatim
00112 *>          B is REAL array, dimension (LDB, N)
00113 *>          On entry, the symmetric matrix B.  If UPLO = 'U', the
00114 *>          leading N-by-N upper triangular part of B contains the
00115 *>          upper triangular part of the matrix B.  If UPLO = 'L',
00116 *>          the leading N-by-N lower triangular part of B contains
00117 *>          the lower triangular part of the matrix B.
00118 *>
00119 *>          On exit, if INFO <= N, the part of B containing the matrix is
00120 *>          overwritten by the triangular factor U or L from the Cholesky
00121 *>          factorization B = U**T*U or B = L*L**T.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDB
00125 *> \verbatim
00126 *>          LDB is INTEGER
00127 *>          The leading dimension of the array B.  LDB >= max(1,N).
00128 *> \endverbatim
00129 *>
00130 *> \param[out] W
00131 *> \verbatim
00132 *>          W is REAL array, dimension (N)
00133 *>          If INFO = 0, the eigenvalues in ascending order.
00134 *> \endverbatim
00135 *>
00136 *> \param[out] WORK
00137 *> \verbatim
00138 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00139 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00140 *> \endverbatim
00141 *>
00142 *> \param[in] LWORK
00143 *> \verbatim
00144 *>          LWORK is INTEGER
00145 *>          The dimension of the array WORK.
00146 *>          If N <= 1,               LWORK >= 1.
00147 *>          If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
00148 *>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
00149 *>
00150 *>          If LWORK = -1, then a workspace query is assumed; the routine
00151 *>          only calculates the optimal sizes of the WORK and IWORK
00152 *>          arrays, returns these values as the first entries of the WORK
00153 *>          and IWORK arrays, and no error message related to LWORK or
00154 *>          LIWORK is issued by XERBLA.
00155 *> \endverbatim
00156 *>
00157 *> \param[out] IWORK
00158 *> \verbatim
00159 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
00160 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00161 *> \endverbatim
00162 *>
00163 *> \param[in] LIWORK
00164 *> \verbatim
00165 *>          LIWORK is INTEGER
00166 *>          The dimension of the array IWORK.
00167 *>          If N <= 1,                LIWORK >= 1.
00168 *>          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
00169 *>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
00170 *>
00171 *>          If LIWORK = -1, then a workspace query is assumed; the
00172 *>          routine only calculates the optimal sizes of the WORK and
00173 *>          IWORK arrays, returns these values as the first entries of
00174 *>          the WORK and IWORK arrays, and no error message related to
00175 *>          LWORK or LIWORK is issued by XERBLA.
00176 *> \endverbatim
00177 *>
00178 *> \param[out] INFO
00179 *> \verbatim
00180 *>          INFO is INTEGER
00181 *>          = 0:  successful exit
00182 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00183 *>          > 0:  SPOTRF or SSYEVD returned an error code:
00184 *>             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
00185 *>                    failed to converge; i off-diagonal elements of an
00186 *>                    intermediate tridiagonal form did not converge to
00187 *>                    zero;
00188 *>                    if INFO = i and JOBZ = 'V', then the algorithm
00189 *>                    failed to compute an eigenvalue while working on
00190 *>                    the submatrix lying in rows and columns INFO/(N+1)
00191 *>                    through mod(INFO,N+1);
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 *> \endverbatim
00197 *
00198 *  Authors:
00199 *  ========
00200 *
00201 *> \author Univ. of Tennessee 
00202 *> \author Univ. of California Berkeley 
00203 *> \author Univ. of Colorado Denver 
00204 *> \author NAG Ltd. 
00205 *
00206 *> \date November 2011
00207 *
00208 *> \ingroup realSYeigen
00209 *
00210 *> \par Further Details:
00211 *  =====================
00212 *>
00213 *> \verbatim
00214 *>
00215 *>  Modified so that no backsubstitution is performed if SSYEVD fails to
00216 *>  converge (NEIG in old code could be greater than N causing out of
00217 *>  bounds reference to A - reported by Ralf Meyer).  Also corrected the
00218 *>  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
00219 *> \endverbatim
00220 *
00221 *> \par Contributors:
00222 *  ==================
00223 *>
00224 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
00225 *>
00226 *  =====================================================================
00227       SUBROUTINE SSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
00228      $                   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, ITYPE, LDA, LDB, LIWORK, LWORK, N
00238 *     ..
00239 *     .. Array Arguments ..
00240       INTEGER            IWORK( * )
00241       REAL               A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
00242 *     ..
00243 *
00244 *  =====================================================================
00245 *
00246 *     .. Parameters ..
00247       REAL               ONE
00248       PARAMETER          ( ONE = 1.0E+0 )
00249 *     ..
00250 *     .. Local Scalars ..
00251       LOGICAL            LQUERY, UPPER, WANTZ
00252       CHARACTER          TRANS
00253       INTEGER            LIOPT, LIWMIN, LOPT, LWMIN
00254 *     ..
00255 *     .. External Functions ..
00256       LOGICAL            LSAME
00257       EXTERNAL           LSAME
00258 *     ..
00259 *     .. External Subroutines ..
00260       EXTERNAL           SPOTRF, SSYEVD, SSYGST, STRMM, STRSM, XERBLA
00261 *     ..
00262 *     .. Intrinsic Functions ..
00263       INTRINSIC          MAX, REAL
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 + 6*N + 2*N**2
00280       ELSE
00281          LIWMIN = 1
00282          LWMIN = 2*N + 1
00283       END IF
00284       LOPT = LWMIN
00285       LIOPT = LIWMIN
00286       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00287          INFO = -1
00288       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00289          INFO = -2
00290       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
00291          INFO = -3
00292       ELSE IF( N.LT.0 ) THEN
00293          INFO = -4
00294       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00295          INFO = -6
00296       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00297          INFO = -8
00298       END IF
00299 *
00300       IF( INFO.EQ.0 ) THEN
00301          WORK( 1 ) = LOPT
00302          IWORK( 1 ) = LIOPT
00303 *
00304          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00305             INFO = -11
00306          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00307             INFO = -13
00308          END IF
00309       END IF
00310 *
00311       IF( INFO.NE.0 ) THEN
00312          CALL XERBLA( 'SSYGVD', -INFO )
00313          RETURN
00314       ELSE IF( LQUERY ) THEN
00315          RETURN
00316       END IF
00317 *
00318 *     Quick return if possible
00319 *
00320       IF( N.EQ.0 )
00321      $   RETURN
00322 *
00323 *     Form a Cholesky factorization of B.
00324 *
00325       CALL SPOTRF( UPLO, N, B, LDB, INFO )
00326       IF( INFO.NE.0 ) THEN
00327          INFO = N + INFO
00328          RETURN
00329       END IF
00330 *
00331 *     Transform problem to standard eigenvalue problem and solve.
00332 *
00333       CALL SSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
00334       CALL SSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
00335      $             INFO )
00336       LOPT = MAX( REAL( LOPT ), REAL( WORK( 1 ) ) )
00337       LIOPT = MAX( REAL( LIOPT ), REAL( IWORK( 1 ) ) )
00338 *
00339       IF( WANTZ .AND. INFO.EQ.0 ) THEN
00340 *
00341 *        Backtransform eigenvectors to the original problem.
00342 *
00343          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
00344 *
00345 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
00346 *           backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
00347 *
00348             IF( UPPER ) THEN
00349                TRANS = 'N'
00350             ELSE
00351                TRANS = 'T'
00352             END IF
00353 *
00354             CALL STRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
00355      $                  B, LDB, A, LDA )
00356 *
00357          ELSE IF( ITYPE.EQ.3 ) THEN
00358 *
00359 *           For B*A*x=(lambda)*x;
00360 *           backtransform eigenvectors: x = L*y or U**T*y
00361 *
00362             IF( UPPER ) THEN
00363                TRANS = 'T'
00364             ELSE
00365                TRANS = 'N'
00366             END IF
00367 *
00368             CALL STRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
00369      $                  B, LDB, A, LDA )
00370          END IF
00371       END IF
00372 *
00373       WORK( 1 ) = LOPT
00374       IWORK( 1 ) = LIOPT
00375 *
00376       RETURN
00377 *
00378 *     End of SSYGVD
00379 *
00380       END
 All Files Functions