|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, 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, LDA, LDB, N 00011 * .. 00012 * .. Array Arguments .. 00013 COMPLEX A( LDA, * ), B( LDB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CHEGS2 reduces a complex Hermitian-definite generalized 00020 * eigenproblem to standard form. 00021 * 00022 * If ITYPE = 1, the problem is A*x = lambda*B*x, 00023 * and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L') 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` or L'*A*L. 00027 * 00028 * B must have been previously factorized as U'*U or L*L' by CPOTRF. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * ITYPE (input) INTEGER 00034 * = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L'); 00035 * = 2 or 3: compute U*A*U' or L'*A*L. 00036 * 00037 * UPLO (input) CHARACTER*1 00038 * Specifies whether the upper or lower triangular part of the 00039 * Hermitian matrix A is stored, and how B has been factorized. 00040 * = 'U': Upper triangular 00041 * = 'L': Lower triangular 00042 * 00043 * N (input) INTEGER 00044 * The order of the matrices A and B. N >= 0. 00045 * 00046 * A (input/output) COMPLEX array, dimension (LDA,N) 00047 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 00048 * n by n upper triangular part of A contains the upper 00049 * triangular part of the matrix A, and the strictly lower 00050 * triangular part of A is not referenced. If UPLO = 'L', the 00051 * leading n by n lower triangular part of A contains the lower 00052 * triangular part of the matrix A, and the strictly upper 00053 * triangular part of A is not referenced. 00054 * 00055 * On exit, if INFO = 0, the transformed matrix, stored in the 00056 * same format as A. 00057 * 00058 * LDA (input) INTEGER 00059 * The leading dimension of the array A. LDA >= max(1,N). 00060 * 00061 * B (input) COMPLEX array, dimension (LDB,N) 00062 * The triangular factor from the Cholesky factorization of B, 00063 * as returned by CPOTRF. 00064 * 00065 * LDB (input) INTEGER 00066 * The leading dimension of the array B. LDB >= max(1,N). 00067 * 00068 * INFO (output) INTEGER 00069 * = 0: successful exit. 00070 * < 0: if INFO = -i, the i-th argument had an illegal value. 00071 * 00072 * ===================================================================== 00073 * 00074 * .. Parameters .. 00075 REAL ONE, HALF 00076 PARAMETER ( ONE = 1.0E+0, HALF = 0.5E+0 ) 00077 COMPLEX CONE 00078 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00079 * .. 00080 * .. Local Scalars .. 00081 LOGICAL UPPER 00082 INTEGER K 00083 REAL AKK, BKK 00084 COMPLEX CT 00085 * .. 00086 * .. External Subroutines .. 00087 EXTERNAL CAXPY, CHER2, CLACGV, CSSCAL, CTRMV, CTRSV, 00088 $ XERBLA 00089 * .. 00090 * .. Intrinsic Functions .. 00091 INTRINSIC MAX 00092 * .. 00093 * .. External Functions .. 00094 LOGICAL LSAME 00095 EXTERNAL LSAME 00096 * .. 00097 * .. Executable Statements .. 00098 * 00099 * Test the input parameters. 00100 * 00101 INFO = 0 00102 UPPER = LSAME( UPLO, 'U' ) 00103 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00104 INFO = -1 00105 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00106 INFO = -2 00107 ELSE IF( N.LT.0 ) THEN 00108 INFO = -3 00109 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00110 INFO = -5 00111 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00112 INFO = -7 00113 END IF 00114 IF( INFO.NE.0 ) THEN 00115 CALL XERBLA( 'CHEGS2', -INFO ) 00116 RETURN 00117 END IF 00118 * 00119 IF( ITYPE.EQ.1 ) THEN 00120 IF( UPPER ) THEN 00121 * 00122 * Compute inv(U')*A*inv(U) 00123 * 00124 DO 10 K = 1, N 00125 * 00126 * Update the upper triangle of A(k:n,k:n) 00127 * 00128 AKK = A( K, K ) 00129 BKK = B( K, K ) 00130 AKK = AKK / BKK**2 00131 A( K, K ) = AKK 00132 IF( K.LT.N ) THEN 00133 CALL CSSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA ) 00134 CT = -HALF*AKK 00135 CALL CLACGV( N-K, A( K, K+1 ), LDA ) 00136 CALL CLACGV( N-K, B( K, K+1 ), LDB ) 00137 CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 00138 $ LDA ) 00139 CALL CHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA, 00140 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA ) 00141 CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 00142 $ LDA ) 00143 CALL CLACGV( N-K, B( K, K+1 ), LDB ) 00144 CALL CTRSV( UPLO, 'Conjugate transpose', 'Non-unit', 00145 $ N-K, B( K+1, K+1 ), LDB, A( K, K+1 ), 00146 $ LDA ) 00147 CALL CLACGV( N-K, A( K, K+1 ), LDA ) 00148 END IF 00149 10 CONTINUE 00150 ELSE 00151 * 00152 * Compute inv(L)*A*inv(L') 00153 * 00154 DO 20 K = 1, N 00155 * 00156 * Update the lower triangle of A(k:n,k:n) 00157 * 00158 AKK = A( K, K ) 00159 BKK = B( K, K ) 00160 AKK = AKK / BKK**2 00161 A( K, K ) = AKK 00162 IF( K.LT.N ) THEN 00163 CALL CSSCAL( N-K, ONE / BKK, A( K+1, K ), 1 ) 00164 CT = -HALF*AKK 00165 CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 00166 CALL CHER2( UPLO, N-K, -CONE, A( K+1, K ), 1, 00167 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA ) 00168 CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 00169 CALL CTRSV( UPLO, 'No transpose', 'Non-unit', N-K, 00170 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 ) 00171 END IF 00172 20 CONTINUE 00173 END IF 00174 ELSE 00175 IF( UPPER ) THEN 00176 * 00177 * Compute U*A*U' 00178 * 00179 DO 30 K = 1, N 00180 * 00181 * Update the upper triangle of A(1:k,1:k) 00182 * 00183 AKK = A( K, K ) 00184 BKK = B( K, K ) 00185 CALL CTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B, 00186 $ LDB, A( 1, K ), 1 ) 00187 CT = HALF*AKK 00188 CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 00189 CALL CHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1, 00190 $ A, LDA ) 00191 CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 00192 CALL CSSCAL( K-1, BKK, A( 1, K ), 1 ) 00193 A( K, K ) = AKK*BKK**2 00194 30 CONTINUE 00195 ELSE 00196 * 00197 * Compute L'*A*L 00198 * 00199 DO 40 K = 1, N 00200 * 00201 * Update the lower triangle of A(1:k,1:k) 00202 * 00203 AKK = A( K, K ) 00204 BKK = B( K, K ) 00205 CALL CLACGV( K-1, A( K, 1 ), LDA ) 00206 CALL CTRMV( UPLO, 'Conjugate transpose', 'Non-unit', K-1, 00207 $ B, LDB, A( K, 1 ), LDA ) 00208 CT = HALF*AKK 00209 CALL CLACGV( K-1, B( K, 1 ), LDB ) 00210 CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 00211 CALL CHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ), 00212 $ LDB, A, LDA ) 00213 CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 00214 CALL CLACGV( K-1, B( K, 1 ), LDB ) 00215 CALL CSSCAL( K-1, BKK, A( K, 1 ), LDA ) 00216 CALL CLACGV( K-1, A( K, 1 ), LDA ) 00217 A( K, K ) = AKK*BKK**2 00218 40 CONTINUE 00219 END IF 00220 END IF 00221 RETURN 00222 * 00223 * End of CHEGS2 00224 * 00225 END
1.7.3 
