SRC/dsygv.f

   1 *> \brief \b DSYGV
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSYGV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
  22 *                         LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, ITYPE, LDA, LDB, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DSYGV computes all the eigenvalues, and optionally, the eigenvectors
  39 *> of a real generalized symmetric-definite eigenproblem, of the form
  40 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
  41 *> Here A and B are assumed to be symmetric and B is also
  42 *> positive definite.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] ITYPE
  49 *> \verbatim
  50 *>          ITYPE is INTEGER
  51 *>          Specifies the problem type to be solved:
  52 *>          = 1:  A*x = (lambda)*B*x
  53 *>          = 2:  A*B*x = (lambda)*x
  54 *>          = 3:  B*A*x = (lambda)*x
  55 *> \endverbatim
  56 *>
  57 *> \param[in] JOBZ
  58 *> \verbatim
  59 *>          JOBZ is CHARACTER*1
  60 *>          = 'N':  Compute eigenvalues only;
  61 *>          = 'V':  Compute eigenvalues and eigenvectors.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] UPLO
  65 *> \verbatim
  66 *>          UPLO is CHARACTER*1
  67 *>          = 'U':  Upper triangles of A and B are stored;
  68 *>          = 'L':  Lower triangles of A and B are stored.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] N
  72 *> \verbatim
  73 *>          N is INTEGER
  74 *>          The order of the matrices A and B.  N >= 0.
  75 *> \endverbatim
  76 *>
  77 *> \param[in,out] A
  78 *> \verbatim
  79 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
  80 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the
  81 *>          leading N-by-N upper triangular part of A contains the
  82 *>          upper triangular part of the matrix A.  If UPLO = 'L',
  83 *>          the leading N-by-N lower triangular part of A contains
  84 *>          the lower triangular part of the matrix A.
  85 *>
  86 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  87 *>          matrix Z of eigenvectors.  The eigenvectors are normalized
  88 *>          as follows:
  89 *>          if ITYPE = 1 or 2, Z**T*B*Z = I;
  90 *>          if ITYPE = 3, Z**T*inv(B)*Z = I.
  91 *>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
  92 *>          or the lower triangle (if UPLO='L') of A, including the
  93 *>          diagonal, is destroyed.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] LDA
  97 *> \verbatim
  98 *>          LDA is INTEGER
  99 *>          The leading dimension of the array A.  LDA >= max(1,N).
 100 *> \endverbatim
 101 *>
 102 *> \param[in,out] B
 103 *> \verbatim
 104 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
 105 *>          On entry, the symmetric positive definite matrix B.
 106 *>          If UPLO = 'U', the leading N-by-N upper triangular part of B
 107 *>          contains the upper triangular part of the matrix B.
 108 *>          If UPLO = 'L', the leading N-by-N lower triangular part of B
 109 *>          contains the lower triangular part of the matrix B.
 110 *>
 111 *>          On exit, if INFO <= N, the part of B containing the matrix is
 112 *>          overwritten by the triangular factor U or L from the Cholesky
 113 *>          factorization B = U**T*U or B = L*L**T.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] LDB
 117 *> \verbatim
 118 *>          LDB is INTEGER
 119 *>          The leading dimension of the array B.  LDB >= max(1,N).
 120 *> \endverbatim
 121 *>
 122 *> \param[out] W
 123 *> \verbatim
 124 *>          W is DOUBLE PRECISION array, dimension (N)
 125 *>          If INFO = 0, the eigenvalues in ascending order.
 126 *> \endverbatim
 127 *>
 128 *> \param[out] WORK
 129 *> \verbatim
 130 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 131 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] LWORK
 135 *> \verbatim
 136 *>          LWORK is INTEGER
 137 *>          The length of the array WORK.  LWORK >= max(1,3*N-1).
 138 *>          For optimal efficiency, LWORK >= (NB+2)*N,
 139 *>          where NB is the blocksize for DSYTRD returned by ILAENV.
 140 *>
 141 *>          If LWORK = -1, then a workspace query is assumed; the routine
 142 *>          only calculates the optimal size of the WORK array, returns
 143 *>          this value as the first entry of the WORK array, and no error
 144 *>          message related to LWORK is issued by XERBLA.
 145 *> \endverbatim
 146 *>
 147 *> \param[out] INFO
 148 *> \verbatim
 149 *>          INFO is INTEGER
 150 *>          = 0:  successful exit
 151 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 152 *>          > 0:  DPOTRF or DSYEV returned an error code:
 153 *>             <= N:  if INFO = i, DSYEV failed to converge;
 154 *>                    i off-diagonal elements of an intermediate
 155 *>                    tridiagonal form did not converge to zero;
 156 *>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
 157 *>                    minor of order i of B is not positive definite.
 158 *>                    The factorization of B could not be completed and
 159 *>                    no eigenvalues or eigenvectors were computed.
 160 *> \endverbatim
 161 *
 162 *  Authors:
 163 *  ========
 164 *
 165 *> \author Univ. of Tennessee
 166 *> \author Univ. of California Berkeley
 167 *> \author Univ. of Colorado Denver
 168 *> \author NAG Ltd.
 169 *
 170 *> \date November 2015
 171 *
 172 *> \ingroup doubleSYeigen
 173 *
 174 *  =====================================================================
 175       SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
 176      $                  LWORK, INFO )
 177 *
 178 *  -- LAPACK driver routine (version 3.6.0) --
 179 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 180 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 181 *     November 2015
 182 *
 183 *     .. Scalar Arguments ..
 184       CHARACTER          JOBZ, UPLO
 185       INTEGER            INFO, ITYPE, LDA, LDB, LWORK, N
 186 *     ..
 187 *     .. Array Arguments ..
 188       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
 189 *     ..
 190 *
 191 *  =====================================================================
 192 *
 193 *     .. Parameters ..
 194       DOUBLE PRECISION   ONE
 195       PARAMETER          ( ONE = 1.0D+0 )
 196 *     ..
 197 *     .. Local Scalars ..
 198       LOGICAL            LQUERY, UPPER, WANTZ
 199       CHARACTER          TRANS
 200       INTEGER            LWKMIN, LWKOPT, NB, NEIG
 201 *     ..
 202 *     .. External Functions ..
 203       LOGICAL            LSAME
 204       INTEGER            ILAENV
 205       EXTERNAL           LSAME, ILAENV
 206 *     ..
 207 *     .. External Subroutines ..
 208       EXTERNAL           DPOTRF, DSYEV, DSYGST, DTRMM, DTRSM, XERBLA
 209 *     ..
 210 *     .. Intrinsic Functions ..
 211       INTRINSIC          MAX
 212 *     ..
 213 *     .. Executable Statements ..
 214 *
 215 *     Test the input parameters.
 216 *
 217       WANTZ = LSAME( JOBZ, 'V' )
 218       UPPER = LSAME( UPLO, 'U' )
 219       LQUERY = ( LWORK.EQ.-1 )
 220 *
 221       INFO = 0
 222       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 223          INFO = -1
 224       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 225          INFO = -2
 226       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 227          INFO = -3
 228       ELSE IF( N.LT.0 ) THEN
 229          INFO = -4
 230       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 231          INFO = -6
 232       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 233          INFO = -8
 234       END IF
 235 *
 236       IF( INFO.EQ.0 ) THEN
 237          LWKMIN = MAX( 1, 3*N - 1 )
 238          NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
 239          LWKOPT = MAX( LWKMIN, ( NB + 2 )*N )
 240          WORK( 1 ) = LWKOPT
 241 *
 242          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
 243             INFO = -11
 244          END IF
 245       END IF
 246 *
 247       IF( INFO.NE.0 ) THEN
 248          CALL XERBLA( 'DSYGV ', -INFO )
 249          RETURN
 250       ELSE IF( LQUERY ) THEN
 251          RETURN
 252       END IF
 253 *
 254 *     Quick return if possible
 255 *
 256       IF( N.EQ.0 )
 257      $   RETURN
 258 *
 259 *     Form a Cholesky factorization of B.
 260 *
 261       CALL DPOTRF( UPLO, N, B, LDB, INFO )
 262       IF( INFO.NE.0 ) THEN
 263          INFO = N + INFO
 264          RETURN
 265       END IF
 266 *
 267 *     Transform problem to standard eigenvalue problem and solve.
 268 *
 269       CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 270       CALL DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
 271 *
 272       IF( WANTZ ) THEN
 273 *
 274 *        Backtransform eigenvectors to the original problem.
 275 *
 276          NEIG = N
 277          IF( INFO.GT.0 )
 278      $      NEIG = INFO - 1
 279          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
 280 *
 281 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
 282 *           backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
 283 *
 284             IF( UPPER ) THEN
 285                TRANS = 'N'
 286             ELSE
 287                TRANS = 'T'
 288             END IF
 289 *
 290             CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
 291      $                  B, LDB, A, LDA )
 292 *
 293          ELSE IF( ITYPE.EQ.3 ) THEN
 294 *
 295 *           For B*A*x=(lambda)*x;
 296 *           backtransform eigenvectors: x = L*y or U**T*y
 297 *
 298             IF( UPPER ) THEN
 299                TRANS = 'T'
 300             ELSE
 301                TRANS = 'N'
 302             END IF
 303 *
 304             CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
 305      $                  B, LDB, A, LDA )
 306          END IF
 307       END IF
 308 *
 309       WORK( 1 ) = LWKOPT
 310       RETURN
 311 *
 312 *     End of DSYGV
 313 *
 314       END