SRC/dsygvd.f

   1 *> \brief \b DSYGVD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSYGVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygvd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygvd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygvd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
  22 *                          LWORK, IWORK, LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DSYGVD computes all the eigenvalues, and optionally, the eigenvectors
  40 *> of a real generalized symmetric-definite eigenproblem, of the form
  41 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
  42 *> B are assumed to be symmetric and B is also positive definite.
  43 *> If eigenvectors are desired, it uses a divide and conquer algorithm.
  44 *>
  45 *> The divide and conquer algorithm makes very mild assumptions about
  46 *> floating point arithmetic. It will work on machines with a guard
  47 *> digit in add/subtract, or on those binary machines without guard
  48 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
  49 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
  50 *> without guard digits, but we know of none.
  51 *> \endverbatim
  52 *
  53 *  Arguments:
  54 *  ==========
  55 *
  56 *> \param[in] ITYPE
  57 *> \verbatim
  58 *>          ITYPE is INTEGER
  59 *>          Specifies the problem type to be solved:
  60 *>          = 1:  A*x = (lambda)*B*x
  61 *>          = 2:  A*B*x = (lambda)*x
  62 *>          = 3:  B*A*x = (lambda)*x
  63 *> \endverbatim
  64 *>
  65 *> \param[in] JOBZ
  66 *> \verbatim
  67 *>          JOBZ is CHARACTER*1
  68 *>          = 'N':  Compute eigenvalues only;
  69 *>          = 'V':  Compute eigenvalues and eigenvectors.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] UPLO
  73 *> \verbatim
  74 *>          UPLO is CHARACTER*1
  75 *>          = 'U':  Upper triangles of A and B are stored;
  76 *>          = 'L':  Lower triangles of A and B are stored.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] N
  80 *> \verbatim
  81 *>          N is INTEGER
  82 *>          The order of the matrices A and B.  N >= 0.
  83 *> \endverbatim
  84 *>
  85 *> \param[in,out] A
  86 *> \verbatim
  87 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
  88 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the
  89 *>          leading N-by-N upper triangular part of A contains the
  90 *>          upper triangular part of the matrix A.  If UPLO = 'L',
  91 *>          the leading N-by-N lower triangular part of A contains
  92 *>          the lower triangular part of the matrix A.
  93 *>
  94 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  95 *>          matrix Z of eigenvectors.  The eigenvectors are normalized
  96 *>          as follows:
  97 *>          if ITYPE = 1 or 2, Z**T*B*Z = I;
  98 *>          if ITYPE = 3, Z**T*inv(B)*Z = I.
  99 *>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
 100 *>          or the lower triangle (if UPLO='L') of A, including the
 101 *>          diagonal, is destroyed.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] LDA
 105 *> \verbatim
 106 *>          LDA is INTEGER
 107 *>          The leading dimension of the array A.  LDA >= max(1,N).
 108 *> \endverbatim
 109 *>
 110 *> \param[in,out] B
 111 *> \verbatim
 112 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
 113 *>          On entry, the symmetric matrix B.  If UPLO = 'U', the
 114 *>          leading N-by-N upper triangular part of B contains the
 115 *>          upper triangular part of the matrix B.  If UPLO = 'L',
 116 *>          the leading N-by-N lower triangular part of B contains
 117 *>          the lower triangular part of the matrix B.
 118 *>
 119 *>          On exit, if INFO <= N, the part of B containing the matrix is
 120 *>          overwritten by the triangular factor U or L from the Cholesky
 121 *>          factorization B = U**T*U or B = L*L**T.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] LDB
 125 *> \verbatim
 126 *>          LDB is INTEGER
 127 *>          The leading dimension of the array B.  LDB >= max(1,N).
 128 *> \endverbatim
 129 *>
 130 *> \param[out] W
 131 *> \verbatim
 132 *>          W is DOUBLE PRECISION array, dimension (N)
 133 *>          If INFO = 0, the eigenvalues in ascending order.
 134 *> \endverbatim
 135 *>
 136 *> \param[out] WORK
 137 *> \verbatim
 138 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 139 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 140 *> \endverbatim
 141 *>
 142 *> \param[in] LWORK
 143 *> \verbatim
 144 *>          LWORK is INTEGER
 145 *>          The dimension of the array WORK.
 146 *>          If N <= 1,               LWORK >= 1.
 147 *>          If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
 148 *>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
 149 *>
 150 *>          If LWORK = -1, then a workspace query is assumed; the routine
 151 *>          only calculates the optimal sizes of the WORK and IWORK
 152 *>          arrays, returns these values as the first entries of the WORK
 153 *>          and IWORK arrays, and no error message related to LWORK or
 154 *>          LIWORK is issued by XERBLA.
 155 *> \endverbatim
 156 *>
 157 *> \param[out] IWORK
 158 *> \verbatim
 159 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 160 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 161 *> \endverbatim
 162 *>
 163 *> \param[in] LIWORK
 164 *> \verbatim
 165 *>          LIWORK is INTEGER
 166 *>          The dimension of the array IWORK.
 167 *>          If N <= 1,                LIWORK >= 1.
 168 *>          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
 169 *>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
 170 *>
 171 *>          If LIWORK = -1, then a workspace query is assumed; the
 172 *>          routine only calculates the optimal sizes of the WORK and
 173 *>          IWORK arrays, returns these values as the first entries of
 174 *>          the WORK and IWORK arrays, and no error message related to
 175 *>          LWORK or LIWORK is issued by XERBLA.
 176 *> \endverbatim
 177 *>
 178 *> \param[out] INFO
 179 *> \verbatim
 180 *>          INFO is INTEGER
 181 *>          = 0:  successful exit
 182 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 183 *>          > 0:  DPOTRF or DSYEVD returned an error code:
 184 *>             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
 185 *>                    failed to converge; i off-diagonal elements of an
 186 *>                    intermediate tridiagonal form did not converge to
 187 *>                    zero;
 188 *>                    if INFO = i and JOBZ = 'V', then the algorithm
 189 *>                    failed to compute an eigenvalue while working on
 190 *>                    the submatrix lying in rows and columns INFO/(N+1)
 191 *>                    through mod(INFO,N+1);
 192 *>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
 193 *>                    minor of order i of B is not positive definite.
 194 *>                    The factorization of B could not be completed and
 195 *>                    no eigenvalues or eigenvectors were computed.
 196 *> \endverbatim
 197 *
 198 *  Authors:
 199 *  ========
 200 *
 201 *> \author Univ. of Tennessee
 202 *> \author Univ. of California Berkeley
 203 *> \author Univ. of Colorado Denver
 204 *> \author NAG Ltd.
 205 *
 206 *> \date November 2015
 207 *
 208 *> \ingroup doubleSYeigen
 209 *
 210 *> \par Further Details:
 211 *  =====================
 212 *>
 213 *> \verbatim
 214 *>
 215 *>  Modified so that no backsubstitution is performed if DSYEVD fails to
 216 *>  converge (NEIG in old code could be greater than N causing out of
 217 *>  bounds reference to A - reported by Ralf Meyer).  Also corrected the
 218 *>  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
 219 *> \endverbatim
 220 *
 221 *> \par Contributors:
 222 *  ==================
 223 *>
 224 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 225 *>
 226 *  =====================================================================
 227       SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
 228      $                   LWORK, IWORK, LIWORK, INFO )
 229 *
 230 *  -- LAPACK driver routine (version 3.6.0) --
 231 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 232 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 233 *     November 2015
 234 *
 235 *     .. Scalar Arguments ..
 236       CHARACTER          JOBZ, UPLO
 237       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
 238 *     ..
 239 *     .. Array Arguments ..
 240       INTEGER            IWORK( * )
 241       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
 242 *     ..
 243 *
 244 *  =====================================================================
 245 *
 246 *     .. Parameters ..
 247       DOUBLE PRECISION   ONE
 248       PARAMETER          ( ONE = 1.0D+0 )
 249 *     ..
 250 *     .. Local Scalars ..
 251       LOGICAL            LQUERY, UPPER, WANTZ
 252       CHARACTER          TRANS
 253       INTEGER            LIOPT, LIWMIN, LOPT, LWMIN
 254 *     ..
 255 *     .. External Functions ..
 256       LOGICAL            LSAME
 257       EXTERNAL           LSAME
 258 *     ..
 259 *     .. External Subroutines ..
 260       EXTERNAL           DPOTRF, DSYEVD, DSYGST, DTRMM, DTRSM, XERBLA
 261 *     ..
 262 *     .. Intrinsic Functions ..
 263       INTRINSIC          DBLE, MAX
 264 *     ..
 265 *     .. Executable Statements ..
 266 *
 267 *     Test the input parameters.
 268 *
 269       WANTZ = LSAME( JOBZ, 'V' )
 270       UPPER = LSAME( UPLO, 'U' )
 271       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 272 *
 273       INFO = 0
 274       IF( N.LE.1 ) THEN
 275          LIWMIN = 1
 276          LWMIN = 1
 277       ELSE IF( WANTZ ) THEN
 278          LIWMIN = 3 + 5*N
 279          LWMIN = 1 + 6*N + 2*N**2
 280       ELSE
 281          LIWMIN = 1
 282          LWMIN = 2*N + 1
 283       END IF
 284       LOPT = LWMIN
 285       LIOPT = LIWMIN
 286       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 287          INFO = -1
 288       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 289          INFO = -2
 290       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 291          INFO = -3
 292       ELSE IF( N.LT.0 ) THEN
 293          INFO = -4
 294       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 295          INFO = -6
 296       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 297          INFO = -8
 298       END IF
 299 *
 300       IF( INFO.EQ.0 ) THEN
 301          WORK( 1 ) = LOPT
 302          IWORK( 1 ) = LIOPT
 303 *
 304          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 305             INFO = -11
 306          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 307             INFO = -13
 308          END IF
 309       END IF
 310 *
 311       IF( INFO.NE.0 ) THEN
 312          CALL XERBLA( 'DSYGVD', -INFO )
 313          RETURN
 314       ELSE IF( LQUERY ) THEN
 315          RETURN
 316       END IF
 317 *
 318 *     Quick return if possible
 319 *
 320       IF( N.EQ.0 )
 321      $   RETURN
 322 *
 323 *     Form a Cholesky factorization of B.
 324 *
 325       CALL DPOTRF( UPLO, N, B, LDB, INFO )
 326       IF( INFO.NE.0 ) THEN
 327          INFO = N + INFO
 328          RETURN
 329       END IF
 330 *
 331 *     Transform problem to standard eigenvalue problem and solve.
 332 *
 333       CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 334       CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
 335      $             INFO )
 336       LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
 337       LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
 338 *
 339       IF( WANTZ .AND. INFO.EQ.0 ) THEN
 340 *
 341 *        Backtransform eigenvectors to the original problem.
 342 *
 343          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
 344 *
 345 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
 346 *           backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
 347 *
 348             IF( UPPER ) THEN
 349                TRANS = 'N'
 350             ELSE
 351                TRANS = 'T'
 352             END IF
 353 *
 354             CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
 355      $                  B, LDB, A, LDA )
 356 *
 357          ELSE IF( ITYPE.EQ.3 ) THEN
 358 *
 359 *           For B*A*x=(lambda)*x;
 360 *           backtransform eigenvectors: x = L*y or U**T*y
 361 *
 362             IF( UPPER ) THEN
 363                TRANS = 'T'
 364             ELSE
 365                TRANS = 'N'
 366             END IF
 367 *
 368             CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
 369      $                  B, LDB, A, LDA )
 370          END IF
 371       END IF
 372 *
 373       WORK( 1 ) = LOPT
 374       IWORK( 1 ) = LIOPT
 375 *
 376       RETURN
 377 *
 378 *     End of DSYGVD
 379 *
 380       END