|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHPGST( 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 AP( * ), BP( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CHPGST 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 CPPTRF. 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 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 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 CPPTRF. 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 REAL ONE, HALF 00068 PARAMETER ( ONE = 1.0E+0, HALF = 0.5E+0 ) 00069 COMPLEX CONE 00070 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00071 * .. 00072 * .. Local Scalars .. 00073 LOGICAL UPPER 00074 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK 00075 REAL AJJ, AKK, BJJ, BKK 00076 COMPLEX CT 00077 * .. 00078 * .. External Subroutines .. 00079 EXTERNAL CAXPY, CHPMV, CHPR2, CSSCAL, CTPMV, CTPSV, 00080 $ XERBLA 00081 * .. 00082 * .. Intrinsic Functions .. 00083 INTRINSIC REAL 00084 * .. 00085 * .. External Functions .. 00086 LOGICAL LSAME 00087 COMPLEX CDOTC 00088 EXTERNAL LSAME, CDOTC 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( 'CHPGST', -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 ) = REAL( AP( JJ ) ) 00123 BJJ = BP( JJ ) 00124 CALL CTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J, 00125 $ BP, AP( J1 ), 1 ) 00126 CALL CHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE, 00127 $ AP( J1 ), 1 ) 00128 CALL CSSCAL( J-1, ONE / BJJ, AP( J1 ), 1 ) 00129 AP( JJ ) = ( AP( JJ )-CDOTC( 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 CSSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 ) 00150 CT = -HALF*AKK 00151 CALL CAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 00152 CALL CHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1, 00153 $ BP( KK+1 ), 1, AP( K1K1 ) ) 00154 CALL CAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 00155 CALL CTPSV( 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 CTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP, 00178 $ AP( K1 ), 1 ) 00179 CT = HALF*AKK 00180 CALL CAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 00181 CALL CHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1, 00182 $ AP ) 00183 CALL CAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 00184 CALL CSSCAL( 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 + CDOTC( N-J, AP( JJ+1 ), 1, 00202 $ BP( JJ+1 ), 1 ) 00203 CALL CSSCAL( N-J, BJJ, AP( JJ+1 ), 1 ) 00204 CALL CHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1, 00205 $ CONE, AP( JJ+1 ), 1 ) 00206 CALL CTPMV( 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 CHPGST 00215 * 00216 END
1.7.3 
