SRC/zhpgst.f

   1 *> \brief \b ZHPGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHPGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpgst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpgst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpgst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, ITYPE, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX*16         AP( * ), BP( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZHPGST reduces a complex Hermitian-definite generalized
  38 *> eigenproblem to standard form, using packed storage.
  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 ZPPTRF.
  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] AP
  75 *> \verbatim
  76 *>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
  77 *>          On entry, the upper or lower triangle of the Hermitian matrix
  78 *>          A, packed columnwise in a linear array.  The j-th column of A
  79 *>          is stored in the array AP as follows:
  80 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  81 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  82 *>
  83 *>          On exit, if INFO = 0, the transformed matrix, stored in the
  84 *>          same format as A.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] BP
  88 *> \verbatim
  89 *>          BP is COMPLEX*16 array, dimension (N*(N+1)/2)
  90 *>          The triangular factor from the Cholesky factorization of B,
  91 *>          stored in the same format as A, as returned by ZPPTRF.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] INFO
  95 *> \verbatim
  96 *>          INFO is INTEGER
  97 *>          = 0:  successful exit
  98 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  99 *> \endverbatim
 100 *
 101 *  Authors:
 102 *  ========
 103 *
 104 *> \author Univ. of Tennessee
 105 *> \author Univ. of California Berkeley
 106 *> \author Univ. of Colorado Denver
 107 *> \author NAG Ltd.
 108 *
 109 *> \date November 2011
 110 *
 111 *> \ingroup complex16OTHERcomputational
 112 *
 113 *  =====================================================================
 114       SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
 115 *
 116 *  -- LAPACK computational routine (version 3.4.0) --
 117 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 118 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 119 *     November 2011
 120 *
 121 *     .. Scalar Arguments ..
 122       CHARACTER          UPLO
 123       INTEGER            INFO, ITYPE, N
 124 *     ..
 125 *     .. Array Arguments ..
 126       COMPLEX*16         AP( * ), BP( * )
 127 *     ..
 128 *
 129 *  =====================================================================
 130 *
 131 *     .. Parameters ..
 132       DOUBLE PRECISION   ONE, HALF
 133       PARAMETER          ( ONE = 1.0D+0, HALF = 0.5D+0 )
 134       COMPLEX*16         CONE
 135       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
 136 *     ..
 137 *     .. Local Scalars ..
 138       LOGICAL            UPPER
 139       INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
 140       DOUBLE PRECISION   AJJ, AKK, BJJ, BKK
 141       COMPLEX*16         CT
 142 *     ..
 143 *     .. External Subroutines ..
 144       EXTERNAL           XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV,
 145      $                   ZTPSV
 146 *     ..
 147 *     .. Intrinsic Functions ..
 148       INTRINSIC          DBLE
 149 *     ..
 150 *     .. External Functions ..
 151       LOGICAL            LSAME
 152       COMPLEX*16         ZDOTC
 153       EXTERNAL           LSAME, ZDOTC
 154 *     ..
 155 *     .. Executable Statements ..
 156 *
 157 *     Test the input parameters.
 158 *
 159       INFO = 0
 160       UPPER = LSAME( UPLO, 'U' )
 161       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 162          INFO = -1
 163       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 164          INFO = -2
 165       ELSE IF( N.LT.0 ) THEN
 166          INFO = -3
 167       END IF
 168       IF( INFO.NE.0 ) THEN
 169          CALL XERBLA( 'ZHPGST', -INFO )
 170          RETURN
 171       END IF
 172 *
 173       IF( ITYPE.EQ.1 ) THEN
 174          IF( UPPER ) THEN
 175 *
 176 *           Compute inv(U**H)*A*inv(U)
 177 *
 178 *           J1 and JJ are the indices of A(1,j) and A(j,j)
 179 *
 180             JJ = 0
 181             DO 10 J = 1, N
 182                J1 = JJ + 1
 183                JJ = JJ + J
 184 *
 185 *              Compute the j-th column of the upper triangle of A
 186 *
 187                AP( JJ ) = DBLE( AP( JJ ) )
 188                BJJ = BP( JJ )
 189                CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J,
 190      $                     BP, AP( J1 ), 1 )
 191                CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE,
 192      $                     AP( J1 ), 1 )
 193                CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
 194                AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ),
 195      $                    1 ) ) / BJJ
 196    10       CONTINUE
 197          ELSE
 198 *
 199 *           Compute inv(L)*A*inv(L**H)
 200 *
 201 *           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
 202 *
 203             KK = 1
 204             DO 20 K = 1, N
 205                K1K1 = KK + N - K + 1
 206 *
 207 *              Update the lower triangle of A(k:n,k:n)
 208 *
 209                AKK = AP( KK )
 210                BKK = BP( KK )
 211                AKK = AKK / BKK**2
 212                AP( KK ) = AKK
 213                IF( K.LT.N ) THEN
 214                   CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
 215                   CT = -HALF*AKK
 216                   CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 217                   CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1,
 218      $                        BP( KK+1 ), 1, AP( K1K1 ) )
 219                   CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 220                   CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
 221      $                        BP( K1K1 ), AP( KK+1 ), 1 )
 222                END IF
 223                KK = K1K1
 224    20       CONTINUE
 225          END IF
 226       ELSE
 227          IF( UPPER ) THEN
 228 *
 229 *           Compute U*A*U**H
 230 *
 231 *           K1 and KK are the indices of A(1,k) and A(k,k)
 232 *
 233             KK = 0
 234             DO 30 K = 1, N
 235                K1 = KK + 1
 236                KK = KK + K
 237 *
 238 *              Update the upper triangle of A(1:k,1:k)
 239 *
 240                AKK = AP( KK )
 241                BKK = BP( KK )
 242                CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
 243      $                     AP( K1 ), 1 )
 244                CT = HALF*AKK
 245                CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 246                CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1,
 247      $                     AP )
 248                CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 249                CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 )
 250                AP( KK ) = AKK*BKK**2
 251    30       CONTINUE
 252          ELSE
 253 *
 254 *           Compute L**H *A*L
 255 *
 256 *           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
 257 *
 258             JJ = 1
 259             DO 40 J = 1, N
 260                J1J1 = JJ + N - J + 1
 261 *
 262 *              Compute the j-th column of the lower triangle of A
 263 *
 264                AJJ = AP( JJ )
 265                BJJ = BP( JJ )
 266                AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1,
 267      $                    BP( JJ+1 ), 1 )
 268                CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
 269                CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1,
 270      $                     CONE, AP( JJ+1 ), 1 )
 271                CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit',
 272      $                     N-J+1, BP( JJ ), AP( JJ ), 1 )
 273                JJ = J1J1
 274    40       CONTINUE
 275          END IF
 276       END IF
 277       RETURN
 278 *
 279 *     End of ZHPGST
 280 *
 281       END