SRC/zhegst.f

   1 *> \brief \b ZHEGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHEGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHEGST( 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 *       COMPLEX*16         A( LDA, * ), B( LDB, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZHEGST reduces a complex Hermitian-definite generalized
  38 *> eigenproblem to standard form.
  39 *>
  40 *> If ITYPE = 1, the problem is A*x = lambda*B*x,
  41 *> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
  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**H or L**H*A*L.
  45 *>
  46 *> B must have been previously factorized as U**H*U or L*L**H by ZPOTRF.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] ITYPE
  53 *> \verbatim
  54 *>          ITYPE is INTEGER
  55 *>          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
  56 *>          = 2 or 3: compute U*A*U**H or L**H*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**H*U;
  64 *>          = 'L':  Lower triangle of A is stored and B is factored as
  65 *>                  L*L**H.
  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 COMPLEX*16 array, dimension (LDA,N)
  77 *>          On entry, the Hermitian 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,out] B
  96 *> \verbatim
  97 *>          B is COMPLEX*16 array, dimension (LDB,N)
  98 *>          The triangular factor from the Cholesky factorization of B,
  99 *>          as returned by ZPOTRF.
 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 September 2012
 124 *
 125 *> \ingroup complex16HEcomputational
 126 *
 127 *  =====================================================================
 128       SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 129 *
 130 *  -- LAPACK computational routine (version 3.4.2) --
 131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 133 *     September 2012
 134 *
 135 *     .. Scalar Arguments ..
 136       CHARACTER          UPLO
 137       INTEGER            INFO, ITYPE, LDA, LDB, N
 138 *     ..
 139 *     .. Array Arguments ..
 140       COMPLEX*16         A( LDA, * ), B( LDB, * )
 141 *     ..
 142 *
 143 *  =====================================================================
 144 *
 145 *     .. Parameters ..
 146       DOUBLE PRECISION   ONE
 147       PARAMETER          ( ONE = 1.0D+0 )
 148       COMPLEX*16         CONE, HALF
 149       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
 150      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
 151 *     ..
 152 *     .. Local Scalars ..
 153       LOGICAL            UPPER
 154       INTEGER            K, KB, NB
 155 *     ..
 156 *     .. External Subroutines ..
 157       EXTERNAL           XERBLA, ZHEGS2, ZHEMM, ZHER2K, ZTRMM, ZTRSM
 158 *     ..
 159 *     .. Intrinsic Functions ..
 160       INTRINSIC          MAX, MIN
 161 *     ..
 162 *     .. External Functions ..
 163       LOGICAL            LSAME
 164       INTEGER            ILAENV
 165       EXTERNAL           LSAME, ILAENV
 166 *     ..
 167 *     .. Executable Statements ..
 168 *
 169 *     Test the input parameters.
 170 *
 171       INFO = 0
 172       UPPER = LSAME( UPLO, 'U' )
 173       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 174          INFO = -1
 175       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 176          INFO = -2
 177       ELSE IF( N.LT.0 ) THEN
 178          INFO = -3
 179       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 180          INFO = -5
 181       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 182          INFO = -7
 183       END IF
 184       IF( INFO.NE.0 ) THEN
 185          CALL XERBLA( 'ZHEGST', -INFO )
 186          RETURN
 187       END IF
 188 *
 189 *     Quick return if possible
 190 *
 191       IF( N.EQ.0 )
 192      $   RETURN
 193 *
 194 *     Determine the block size for this environment.
 195 *
 196       NB = ILAENV( 1, 'ZHEGST', UPLO, N, -1, -1, -1 )
 197 *
 198       IF( NB.LE.1 .OR. NB.GE.N ) THEN
 199 *
 200 *        Use unblocked code
 201 *
 202          CALL ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 203       ELSE
 204 *
 205 *        Use blocked code
 206 *
 207          IF( ITYPE.EQ.1 ) THEN
 208             IF( UPPER ) THEN
 209 *
 210 *              Compute inv(U**H)*A*inv(U)
 211 *
 212                DO 10 K = 1, N, NB
 213                   KB = MIN( N-K+1, NB )
 214 *
 215 *                 Update the upper triangle of A(k:n,k:n)
 216 *
 217                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 218      $                         B( K, K ), LDB, INFO )
 219                   IF( K+KB.LE.N ) THEN
 220                      CALL ZTRSM( 'Left', UPLO, 'Conjugate transpose',
 221      $                           'Non-unit', KB, N-K-KB+1, CONE,
 222      $                           B( K, K ), LDB, A( K, K+KB ), LDA )
 223                      CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
 224      $                           A( K, K ), LDA, B( K, K+KB ), LDB,
 225      $                           CONE, A( K, K+KB ), LDA )
 226                      CALL ZHER2K( UPLO, 'Conjugate transpose', N-K-KB+1,
 227      $                            KB, -CONE, A( K, K+KB ), LDA,
 228      $                            B( K, K+KB ), LDB, ONE,
 229      $                            A( K+KB, K+KB ), LDA )
 230                      CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
 231      $                           A( K, K ), LDA, B( K, K+KB ), LDB,
 232      $                           CONE, A( K, K+KB ), LDA )
 233                      CALL ZTRSM( 'Right', UPLO, 'No transpose',
 234      $                           'Non-unit', KB, N-K-KB+1, CONE,
 235      $                           B( K+KB, K+KB ), LDB, A( K, K+KB ),
 236      $                           LDA )
 237                   END IF
 238    10          CONTINUE
 239             ELSE
 240 *
 241 *              Compute inv(L)*A*inv(L**H)
 242 *
 243                DO 20 K = 1, N, NB
 244                   KB = MIN( N-K+1, NB )
 245 *
 246 *                 Update the lower triangle of A(k:n,k:n)
 247 *
 248                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 249      $                         B( K, K ), LDB, INFO )
 250                   IF( K+KB.LE.N ) THEN
 251                      CALL ZTRSM( 'Right', UPLO, 'Conjugate transpose',
 252      $                           'Non-unit', N-K-KB+1, KB, CONE,
 253      $                           B( K, K ), LDB, A( K+KB, K ), LDA )
 254                      CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
 255      $                           A( K, K ), LDA, B( K+KB, K ), LDB,
 256      $                           CONE, A( K+KB, K ), LDA )
 257                      CALL ZHER2K( UPLO, 'No transpose', N-K-KB+1, KB,
 258      $                            -CONE, A( K+KB, K ), LDA,
 259      $                            B( K+KB, K ), LDB, ONE,
 260      $                            A( K+KB, K+KB ), LDA )
 261                      CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
 262      $                           A( K, K ), LDA, B( K+KB, K ), LDB,
 263      $                           CONE, A( K+KB, K ), LDA )
 264                      CALL ZTRSM( 'Left', UPLO, 'No transpose',
 265      $                           'Non-unit', N-K-KB+1, KB, CONE,
 266      $                           B( K+KB, K+KB ), LDB, A( K+KB, K ),
 267      $                           LDA )
 268                   END IF
 269    20          CONTINUE
 270             END IF
 271          ELSE
 272             IF( UPPER ) THEN
 273 *
 274 *              Compute U*A*U**H
 275 *
 276                DO 30 K = 1, N, NB
 277                   KB = MIN( N-K+1, NB )
 278 *
 279 *                 Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
 280 *
 281                   CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
 282      $                        K-1, KB, CONE, B, LDB, A( 1, K ), LDA )
 283                   CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
 284      $                        LDA, B( 1, K ), LDB, CONE, A( 1, K ),
 285      $                        LDA )
 286                   CALL ZHER2K( UPLO, 'No transpose', K-1, KB, CONE,
 287      $                         A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
 288      $                         LDA )
 289                   CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
 290      $                        LDA, B( 1, K ), LDB, CONE, A( 1, K ),
 291      $                        LDA )
 292                   CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose',
 293      $                        'Non-unit', K-1, KB, CONE, B( K, K ), LDB,
 294      $                        A( 1, K ), LDA )
 295                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 296      $                         B( K, K ), LDB, INFO )
 297    30          CONTINUE
 298             ELSE
 299 *
 300 *              Compute L**H*A*L
 301 *
 302                DO 40 K = 1, N, NB
 303                   KB = MIN( N-K+1, NB )
 304 *
 305 *                 Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
 306 *
 307                   CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
 308      $                        KB, K-1, CONE, B, LDB, A( K, 1 ), LDA )
 309                   CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
 310      $                        LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
 311      $                        LDA )
 312                   CALL ZHER2K( UPLO, 'Conjugate transpose', K-1, KB,
 313      $                         CONE, A( K, 1 ), LDA, B( K, 1 ), LDB,
 314      $                         ONE, A, LDA )
 315                   CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
 316      $                        LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
 317      $                        LDA )
 318                   CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose',
 319      $                        'Non-unit', KB, K-1, CONE, B( K, K ), LDB,
 320      $                        A( K, 1 ), LDA )
 321                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
 322      $                         B( K, K ), LDB, INFO )
 323    40          CONTINUE
 324             END IF
 325          END IF
 326       END IF
 327       RETURN
 328 *
 329 *     End of ZHEGST
 330 *
 331       END