LAPACK: SRC/zhpgst.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, ITYPE, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX*16         AP( * ), BP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZHPGST reduces a complex Hermitian-definite generalized
00020 *  eigenproblem to standard form, using packed storage.
00021 *
00022 *  If ITYPE = 1, the problem is A*x = lambda*B*x,
00023 *  and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
00024 *
00025 *  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
00026 *  B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
00027 *
00028 *  B must have been previously factorized as U**H*U or L*L**H by ZPPTRF.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  ITYPE   (input) INTEGER
00034 *          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
00035 *          = 2 or 3: compute U*A*U**H or L**H*A*L.
00036 *
00037 *  UPLO    (input) CHARACTER*1
00038 *          = 'U':  Upper triangle of A is stored and B is factored as
00039 *                  U**H*U;
00040 *          = 'L':  Lower triangle of A is stored and B is factored as
00041 *                  L*L**H.
00042 *
00043 *  N       (input) INTEGER
00044 *          The order of the matrices A and B.  N >= 0.
00045 *
00046 *  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
00047 *          On entry, the upper or lower triangle of the Hermitian matrix
00048 *          A, packed columnwise in a linear array.  The j-th column of A
00049 *          is stored in the array AP as follows:
00050 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00051 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
00052 *
00053 *          On exit, if INFO = 0, the transformed matrix, stored in the
00054 *          same format as A.
00055 *
00056 *  BP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
00057 *          The triangular factor from the Cholesky factorization of B,
00058 *          stored in the same format as A, as returned by ZPPTRF.
00059 *
00060 *  INFO    (output) INTEGER
00061 *          = 0:  successful exit
00062 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00063 *
00064 *  =====================================================================
00065 *
00066 *     .. Parameters ..
00067       DOUBLE PRECISION   ONE, HALF
00068       PARAMETER          ( ONE = 1.0D+0, HALF = 0.5D+0 )
00069       COMPLEX*16         CONE
00070       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00071 *     ..
00072 *     .. Local Scalars ..
00073       LOGICAL            UPPER
00074       INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
00075       DOUBLE PRECISION   AJJ, AKK, BJJ, BKK
00076       COMPLEX*16         CT
00077 *     ..
00078 *     .. External Subroutines ..
00079       EXTERNAL           XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV,
00080      $                   ZTPSV
00081 *     ..
00082 *     .. Intrinsic Functions ..
00083       INTRINSIC          DBLE
00084 *     ..
00085 *     .. External Functions ..
00086       LOGICAL            LSAME
00087       COMPLEX*16         ZDOTC
00088       EXTERNAL           LSAME, ZDOTC
00089 *     ..
00090 *     .. Executable Statements ..
00091 *
00092 *     Test the input parameters.
00093 *
00094       INFO = 0
00095       UPPER = LSAME( UPLO, 'U' )
00096       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00097          INFO = -1
00098       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00099          INFO = -2
00100       ELSE IF( N.LT.0 ) THEN
00101          INFO = -3
00102       END IF
00103       IF( INFO.NE.0 ) THEN
00104          CALL XERBLA( 'ZHPGST', -INFO )
00105          RETURN
00106       END IF
00107 *
00108       IF( ITYPE.EQ.1 ) THEN
00109          IF( UPPER ) THEN
00110 *
00111 *           Compute inv(U')*A*inv(U)
00112 *
00113 *           J1 and JJ are the indices of A(1,j) and A(j,j)
00114 *
00115             JJ = 0
00116             DO 10 J = 1, N
00117                J1 = JJ + 1
00118                JJ = JJ + J
00119 *
00120 *              Compute the j-th column of the upper triangle of A
00121 *
00122                AP( JJ ) = DBLE( AP( JJ ) )
00123                BJJ = BP( JJ )
00124                CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J,
00125      $                     BP, AP( J1 ), 1 )
00126                CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE,
00127      $                     AP( J1 ), 1 )
00128                CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
00129                AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ),
00130      $                    1 ) ) / BJJ
00131    10       CONTINUE
00132          ELSE
00133 *
00134 *           Compute inv(L)*A*inv(L')
00135 *
00136 *           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
00137 *
00138             KK = 1
00139             DO 20 K = 1, N
00140                K1K1 = KK + N - K + 1
00141 *
00142 *              Update the lower triangle of A(k:n,k:n)
00143 *
00144                AKK = AP( KK )
00145                BKK = BP( KK )
00146                AKK = AKK / BKK**2
00147                AP( KK ) = AKK
00148                IF( K.LT.N ) THEN
00149                   CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
00150                   CT = -HALF*AKK
00151                   CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
00152                   CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1,
00153      $                        BP( KK+1 ), 1, AP( K1K1 ) )
00154                   CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
00155                   CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
00156      $                        BP( K1K1 ), AP( KK+1 ), 1 )
00157                END IF
00158                KK = K1K1
00159    20       CONTINUE
00160          END IF
00161       ELSE
00162          IF( UPPER ) THEN
00163 *
00164 *           Compute U*A*U'
00165 *
00166 *           K1 and KK are the indices of A(1,k) and A(k,k)
00167 *
00168             KK = 0
00169             DO 30 K = 1, N
00170                K1 = KK + 1
00171                KK = KK + K
00172 *
00173 *              Update the upper triangle of A(1:k,1:k)
00174 *
00175                AKK = AP( KK )
00176                BKK = BP( KK )
00177                CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
00178      $                     AP( K1 ), 1 )
00179                CT = HALF*AKK
00180                CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
00181                CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1,
00182      $                     AP )
00183                CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
00184                CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 )
00185                AP( KK ) = AKK*BKK**2
00186    30       CONTINUE
00187          ELSE
00188 *
00189 *           Compute L'*A*L
00190 *
00191 *           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
00192 *
00193             JJ = 1
00194             DO 40 J = 1, N
00195                J1J1 = JJ + N - J + 1
00196 *
00197 *              Compute the j-th column of the lower triangle of A
00198 *
00199                AJJ = AP( JJ )
00200                BJJ = BP( JJ )
00201                AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1,
00202      $                    BP( JJ+1 ), 1 )
00203                CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
00204                CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1,
00205      $                     CONE, AP( JJ+1 ), 1 )
00206                CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit',
00207      $                     N-J+1, BP( JJ ), AP( JJ ), 1 )
00208                JJ = J1J1
00209    40       CONTINUE
00210          END IF
00211       END IF
00212       RETURN
00213 *
00214 *     End of ZHPGST
00215 *
00216       END