src/LAPACK/zhpgst.f

   1       SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          UPLO
   9       INTEGER            INFO, ITYPE, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       COMPLEX*16         AP( * ), BP( * )
  13 *     ..
  14 *
  15 *  Purpose
  16 *  =======
  17 *
  18 *  ZHPGST reduces a complex Hermitian-definite generalized
  19 *  eigenproblem to standard form, using packed storage.
  20 *
  21 *  If ITYPE = 1, the problem is A*x = lambda*B*x,
  22 *  and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
  23 *
  24 *  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
  25 *  B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
  26 *
  27 *  B must have been previously factorized as U**H*U or L*L**H by ZPPTRF.
  28 *
  29 *  Arguments
  30 *  =========
  31 *
  32 *  ITYPE   (input) INTEGER
  33 *          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
  34 *          = 2 or 3: compute U*A*U**H or L**H*A*L.
  35 *
  36 *  UPLO    (input) CHARACTER*1
  37 *          = 'U':  Upper triangle of A is stored and B is factored as
  38 *                  U**H*U;
  39 *          = 'L':  Lower triangle of A is stored and B is factored as
  40 *                  L*L**H.
  41 *
  42 *  N       (input) INTEGER
  43 *          The order of the matrices A and B.  N >= 0.
  44 *
  45 *  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
  46 *          On entry, the upper or lower triangle of the Hermitian matrix
  47 *          A, packed columnwise in a linear array.  The j-th column of A
  48 *          is stored in the array AP as follows:
  49 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  50 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  51 *
  52 *          On exit, if INFO = 0, the transformed matrix, stored in the
  53 *          same format as A.
  54 *
  55 *  BP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
  56 *          The triangular factor from the Cholesky factorization of B,
  57 *          stored in the same format as A, as returned by ZPPTRF.
  58 *
  59 *  INFO    (output) INTEGER
  60 *          = 0:  successful exit
  61 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  62 *
  63 *  =====================================================================
  64 *
  65 *     .. Parameters ..
  66       DOUBLE PRECISION   ONE, HALF
  67       PARAMETER          ( ONE = 1.0D+0, HALF = 0.5D+0 )
  68       COMPLEX*16         CONE
  69       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
  70 *     ..
  71 *     .. Local Scalars ..
  72       LOGICAL            UPPER
  73       INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
  74       DOUBLE PRECISION   AJJ, AKK, BJJ, BKK
  75       COMPLEX*16         CT
  76 *     ..
  77 *     .. External Subroutines ..
  78       EXTERNAL           XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV,
  79      $                   ZTPSV
  80 *     ..
  81 *     .. Intrinsic Functions ..
  82       INTRINSIC          DBLE
  83 *     ..
  84 *     .. External Functions ..
  85       LOGICAL            LSAME
  86       COMPLEX*16         ZDOTC
  87       EXTERNAL           LSAME, ZDOTC
  88 *     ..
  89 *     .. Executable Statements ..
  90 *
  91 *     Test the input parameters.
  92 *
  93       INFO = 0
  94       UPPER = LSAME( UPLO, 'U' )
  95       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
  96          INFO = -1
  97       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  98          INFO = -2
  99       ELSE IF( N.LT.0 ) THEN
 100          INFO = -3
 101       END IF
 102       IF( INFO.NE.0 ) THEN
 103          CALL XERBLA( 'ZHPGST', -INFO )
 104          RETURN
 105       END IF
 106 *
 107       IF( ITYPE.EQ.1 ) THEN
 108          IF( UPPER ) THEN
 109 *
 110 *           Compute inv(U')*A*inv(U)
 111 *
 112 *           J1 and JJ are the indices of A(1,j) and A(j,j)
 113 *
 114             JJ = 0
 115             DO 10 J = 1, N
 116                J1 = JJ + 1
 117                JJ = JJ + J
 118 *
 119 *              Compute the j-th column of the upper triangle of A
 120 *
 121                AP( JJ ) = DBLE( AP( JJ ) )
 122                BJJ = BP( JJ )
 123                CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J,
 124      $                     BP, AP( J1 ), 1 )
 125                CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE,
 126      $                     AP( J1 ), 1 )
 127                CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
 128                AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ),
 129      $                    1 ) ) / BJJ
 130    10       CONTINUE
 131          ELSE
 132 *
 133 *           Compute inv(L)*A*inv(L')
 134 *
 135 *           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
 136 *
 137             KK = 1
 138             DO 20 K = 1, N
 139                K1K1 = KK + N - K + 1
 140 *
 141 *              Update the lower triangle of A(k:n,k:n)
 142 *
 143                AKK = AP( KK )
 144                BKK = BP( KK )
 145                AKK = AKK / BKK**2
 146                AP( KK ) = AKK
 147                IF( K.LT.N ) THEN
 148                   CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
 149                   CT = -HALF*AKK
 150                   CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 151                   CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1,
 152      $                        BP( KK+1 ), 1, AP( K1K1 ) )
 153                   CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 154                   CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
 155      $                        BP( K1K1 ), AP( KK+1 ), 1 )
 156                END IF
 157                KK = K1K1
 158    20       CONTINUE
 159          END IF
 160       ELSE
 161          IF( UPPER ) THEN
 162 *
 163 *           Compute U*A*U'
 164 *
 165 *           K1 and KK are the indices of A(1,k) and A(k,k)
 166 *
 167             KK = 0
 168             DO 30 K = 1, N
 169                K1 = KK + 1
 170                KK = KK + K
 171 *
 172 *              Update the upper triangle of A(1:k,1:k)
 173 *
 174                AKK = AP( KK )
 175                BKK = BP( KK )
 176                CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
 177      $                     AP( K1 ), 1 )
 178                CT = HALF*AKK
 179                CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 180                CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1,
 181      $                     AP )
 182                CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 183                CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 )
 184                AP( KK ) = AKK*BKK**2
 185    30       CONTINUE
 186          ELSE
 187 *
 188 *           Compute L'*A*L
 189 *
 190 *           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
 191 *
 192             JJ = 1
 193             DO 40 J = 1, N
 194                J1J1 = JJ + N - J + 1
 195 *
 196 *              Compute the j-th column of the lower triangle of A
 197 *
 198                AJJ = AP( JJ )
 199                BJJ = BP( JJ )
 200                AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1,
 201      $                    BP( JJ+1 ), 1 )
 202                CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
 203                CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1,
 204      $                     CONE, AP( JJ+1 ), 1 )
 205                CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit',
 206      $                     N-J+1, BP( JJ ), AP( JJ ), 1 )
 207                JJ = J1J1
 208    40       CONTINUE
 209          END IF
 210       END IF
 211       RETURN
 212 *
 213 *     End of ZHPGST
 214 *
 215       END