SRC/dsygst.f

   1 *> \brief \b DSYGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSYGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, ITYPE, LDA, LDB, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DSYGST reduces a real symmetric-definite generalized eigenproblem
  38 *> to standard form.
  39 *>
  40 *> If ITYPE = 1, the problem is A*x = lambda*B*x,
  41 *> and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
  42 *>
  43 *> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
  44 *> B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
  45 *>
  46 *> B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] ITYPE
  53 *> \verbatim
  54 *>          ITYPE is INTEGER
  55 *>          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
  56 *>          = 2 or 3: compute U*A*U**T or L**T*A*L.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] UPLO
  60 *> \verbatim
  61 *>          UPLO is CHARACTER*1
  62 *>          = 'U':  Upper triangle of A is stored and B is factored as
  63 *>                  U**T*U;
  64 *>          = 'L':  Lower triangle of A is stored and B is factored as
  65 *>                  L*L**T.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] N
  69 *> \verbatim
  70 *>          N is INTEGER
  71 *>          The order of the matrices A and B.  N >= 0.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] A
  75 *> \verbatim
  76 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  77 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  78 *>          N-by-N upper triangular part of A contains the upper
  79 *>          triangular part of the matrix A, and the strictly lower
  80 *>          triangular part of A is not referenced.  If UPLO = 'L', the
  81 *>          leading N-by-N lower triangular part of A contains the lower
  82 *>          triangular part of the matrix A, and the strictly upper
  83 *>          triangular part of A is not referenced.
  84 *>
  85 *>          On exit, if INFO = 0, the transformed matrix, stored in the
  86 *>          same format as A.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDA
  90 *> \verbatim
  91 *>          LDA is INTEGER
  92 *>          The leading dimension of the array A.  LDA >= max(1,N).
  93 *> \endverbatim
  94 *>
  95 *> \param[in] B
  96 *> \verbatim
  97 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
  98 *>          The triangular factor from the Cholesky factorization of B,
  99 *>          as returned by DPOTRF.
 100 *> \endverbatim
 101 *>
 102 *> \param[in] LDB
 103 *> \verbatim
 104 *>          LDB is INTEGER
 105 *>          The leading dimension of the array B.  LDB >= max(1,N).
 106 *> \endverbatim
 107 *>
 108 *> \param[out] INFO
 109 *> \verbatim
 110 *>          INFO is INTEGER
 111 *>          = 0:  successful exit
 112 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 113 *> \endverbatim
 114 *
 115 *  Authors:
 116 *  ========
 117 *
 118 *> \author Univ. of Tennessee
 119 *> \author Univ. of California Berkeley
 120 *> \author Univ. of Colorado Denver
 121 *> \author NAG Ltd.
 122 *
 123 *> \date November 2011
 124 *
 125 *> \ingroup doubleSYcomputational
 126 *
 127 *  =====================================================================
 128       SUBROUTINE DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 129 *
 130 *  -- LAPACK computational routine (version 3.4.0) --
 131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 133 *     November 2011
 134 *
 135 *     .. Scalar Arguments ..
 136       CHARACTER          UPLO
 137       INTEGER            INFO, ITYPE, LDA, LDB, N
 138 *     ..
 139 *     .. Array Arguments ..
 140       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
 141 *     ..
 142 *
 143 *  =====================================================================
 144 *
 145 *     .. Parameters ..
 146       DOUBLE PRECISION   ONE, HALF
 147       PARAMETER          ( ONE = 1.0D0, HALF = 0.5D0 )
 148 *     ..
 149 *     .. Local Scalars ..
 150       LOGICAL            UPPER
 151       INTEGER            K, KB, NB
 152 *     ..
 153 *     .. External Subroutines ..
 154       EXTERNAL           DSYGS2, DSYMM, DSYR2K, DTRMM, DTRSM, XERBLA
 155 *     ..
 156 *     .. Intrinsic Functions ..
 157       INTRINSIC          MAX, MIN
 158 *     ..
 159 *     .. External Functions ..
 160       LOGICAL            LSAME
 161       INTEGER            ILAENV
 162       EXTERNAL           LSAME, ILAENV
 163 *     ..
 164 *     .. Executable Statements ..
 165 *
 166 *     Test the input parameters.
 167 *
 168       INFO = 0
 169       UPPER = LSAME( UPLO, 'U' )
 170       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 171          INFO = -1
 172       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 173          INFO = -2
 174       ELSE IF( N.LT.0 ) THEN
 175          INFO = -3
 176       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 177          INFO = -5
 178       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 179          INFO = -7
 180       END IF
 181       IF( INFO.NE.0 ) THEN
 182          CALL XERBLA( 'DSYGST', -INFO )
 183          RETURN
 184       END IF
 185 *
 186 *     Quick return if possible
 187 *
 188       IF( N.EQ.0 )
 189      $   RETURN
 190 *
 191 *     Determine the block size for this environment.
 192 *
 193       NB = ILAENV( 1, 'DSYGST', UPLO, N, -1, -1, -1 )
 194 *
 195       IF( NB.LE.1 .OR. NB.GE.N ) THEN
 196 *
 197 *        Use unblocked code
 198 *
 199          CALL DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 200       ELSE
 201 *
 202 *        Use blocked code
 203 *
 204          IF( ITYPE.EQ.1 ) THEN
 205             IF( UPPER ) THEN
 206 *
 207 *              Compute inv(U**T)*A*inv(U)
 208 *
 209                DO 10 K = 1, N, NB
 210                   KB = MIN( N-K+1, NB )
 211 *
 212 *                 Update the upper triangle of A(k:n,k:n)
 213 *
 214                   CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 215      $                         B( K, K ), LDB, INFO )
 216                   IF( K+KB.LE.N ) THEN
 217                      CALL DTRSM( 'Left', UPLO, 'Transpose', 'Non-unit',
 218      $                           KB, N-K-KB+1, ONE, B( K, K ), LDB,
 219      $                           A( K, K+KB ), LDA )
 220                      CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
 221      $                           A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
 222      $                           A( K, K+KB ), LDA )
 223                      CALL DSYR2K( UPLO, 'Transpose', N-K-KB+1, KB, -ONE,
 224      $                            A( K, K+KB ), LDA, B( K, K+KB ), LDB,
 225      $                            ONE, A( K+KB, K+KB ), LDA )
 226                      CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
 227      $                           A( K, K ), LDA, B( K, K+KB ), LDB, ONE,
 228      $                           A( K, K+KB ), LDA )
 229                      CALL DTRSM( 'Right', UPLO, 'No transpose',
 230      $                           'Non-unit', KB, N-K-KB+1, ONE,
 231      $                           B( K+KB, K+KB ), LDB, A( K, K+KB ),
 232      $                           LDA )
 233                   END IF
 234    10          CONTINUE
 235             ELSE
 236 *
 237 *              Compute inv(L)*A*inv(L**T)
 238 *
 239                DO 20 K = 1, N, NB
 240                   KB = MIN( N-K+1, NB )
 241 *
 242 *                 Update the lower triangle of A(k:n,k:n)
 243 *
 244                   CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 245      $                         B( K, K ), LDB, INFO )
 246                   IF( K+KB.LE.N ) THEN
 247                      CALL DTRSM( 'Right', UPLO, 'Transpose', 'Non-unit',
 248      $                           N-K-KB+1, KB, ONE, B( K, K ), LDB,
 249      $                           A( K+KB, K ), LDA )
 250                      CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
 251      $                           A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
 252      $                           A( K+KB, K ), LDA )
 253                      CALL DSYR2K( UPLO, 'No transpose', N-K-KB+1, KB,
 254      $                            -ONE, A( K+KB, K ), LDA, B( K+KB, K ),
 255      $                            LDB, ONE, A( K+KB, K+KB ), LDA )
 256                      CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
 257      $                           A( K, K ), LDA, B( K+KB, K ), LDB, ONE,
 258      $                           A( K+KB, K ), LDA )
 259                      CALL DTRSM( 'Left', UPLO, 'No transpose',
 260      $                           'Non-unit', N-K-KB+1, KB, ONE,
 261      $                           B( K+KB, K+KB ), LDB, A( K+KB, K ),
 262      $                           LDA )
 263                   END IF
 264    20          CONTINUE
 265             END IF
 266          ELSE
 267             IF( UPPER ) THEN
 268 *
 269 *              Compute U*A*U**T
 270 *
 271                DO 30 K = 1, N, NB
 272                   KB = MIN( N-K+1, NB )
 273 *
 274 *                 Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
 275 *
 276                   CALL DTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
 277      $                        K-1, KB, ONE, B, LDB, A( 1, K ), LDA )
 278                   CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
 279      $                        LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
 280                   CALL DSYR2K( UPLO, 'No transpose', K-1, KB, ONE,
 281      $                         A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
 282      $                         LDA )
 283                   CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
 284      $                        LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA )
 285                   CALL DTRMM( 'Right', UPLO, 'Transpose', 'Non-unit',
 286      $                        K-1, KB, ONE, B( K, K ), LDB, A( 1, K ),
 287      $                        LDA )
 288                   CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 289      $                         B( K, K ), LDB, INFO )
 290    30          CONTINUE
 291             ELSE
 292 *
 293 *              Compute L**T*A*L
 294 *
 295                DO 40 K = 1, N, NB
 296                   KB = MIN( N-K+1, NB )
 297 *
 298 *                 Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
 299 *
 300                   CALL DTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
 301      $                        KB, K-1, ONE, B, LDB, A( K, 1 ), LDA )
 302                   CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
 303      $                        LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
 304                   CALL DSYR2K( UPLO, 'Transpose', K-1, KB, ONE,
 305      $                         A( K, 1 ), LDA, B( K, 1 ), LDB, ONE, A,
 306      $                         LDA )
 307                   CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
 308      $                        LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA )
 309                   CALL DTRMM( 'Left', UPLO, 'Transpose', 'Non-unit', KB,
 310      $                        K-1, ONE, B( K, K ), LDB, A( K, 1 ), LDA )
 311                   CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 312      $                         B( K, K ), LDB, INFO )
 313    40          CONTINUE
 314             END IF
 315          END IF
 316       END IF
 317       RETURN
 318 *
 319 *     End of DSYGST
 320 *
 321       END