|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHEGST( 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 * CHEGST 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**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 CPOTRF. 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 * 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 00076 PARAMETER ( ONE = 1.0E+0 ) 00077 COMPLEX CONE, HALF 00078 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 00079 $ HALF = ( 0.5E+0, 0.0E+0 ) ) 00080 * .. 00081 * .. Local Scalars .. 00082 LOGICAL UPPER 00083 INTEGER K, KB, NB 00084 * .. 00085 * .. External Subroutines .. 00086 EXTERNAL CHEGS2, CHEMM, CHER2K, CTRMM, CTRSM, XERBLA 00087 * .. 00088 * .. Intrinsic Functions .. 00089 INTRINSIC MAX, MIN 00090 * .. 00091 * .. External Functions .. 00092 LOGICAL LSAME 00093 INTEGER ILAENV 00094 EXTERNAL LSAME, ILAENV 00095 * .. 00096 * .. Executable Statements .. 00097 * 00098 * Test the input parameters. 00099 * 00100 INFO = 0 00101 UPPER = LSAME( UPLO, 'U' ) 00102 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00103 INFO = -1 00104 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00105 INFO = -2 00106 ELSE IF( N.LT.0 ) THEN 00107 INFO = -3 00108 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00109 INFO = -5 00110 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00111 INFO = -7 00112 END IF 00113 IF( INFO.NE.0 ) THEN 00114 CALL XERBLA( 'CHEGST', -INFO ) 00115 RETURN 00116 END IF 00117 * 00118 * Quick return if possible 00119 * 00120 IF( N.EQ.0 ) 00121 $ RETURN 00122 * 00123 * Determine the block size for this environment. 00124 * 00125 NB = ILAENV( 1, 'CHEGST', UPLO, N, -1, -1, -1 ) 00126 * 00127 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00128 * 00129 * Use unblocked code 00130 * 00131 CALL CHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00132 ELSE 00133 * 00134 * Use blocked code 00135 * 00136 IF( ITYPE.EQ.1 ) THEN 00137 IF( UPPER ) THEN 00138 * 00139 * Compute inv(U')*A*inv(U) 00140 * 00141 DO 10 K = 1, N, NB 00142 KB = MIN( N-K+1, NB ) 00143 * 00144 * Update the upper triangle of A(k:n,k:n) 00145 * 00146 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00147 $ B( K, K ), LDB, INFO ) 00148 IF( K+KB.LE.N ) THEN 00149 CALL CTRSM( 'Left', UPLO, 'Conjugate transpose', 00150 $ 'Non-unit', KB, N-K-KB+1, CONE, 00151 $ B( K, K ), LDB, A( K, K+KB ), LDA ) 00152 CALL CHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, 00153 $ A( K, K ), LDA, B( K, K+KB ), LDB, 00154 $ CONE, A( K, K+KB ), LDA ) 00155 CALL CHER2K( UPLO, 'Conjugate transpose', N-K-KB+1, 00156 $ KB, -CONE, A( K, K+KB ), LDA, 00157 $ B( K, K+KB ), LDB, ONE, 00158 $ A( K+KB, K+KB ), LDA ) 00159 CALL CHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, 00160 $ A( K, K ), LDA, B( K, K+KB ), LDB, 00161 $ CONE, A( K, K+KB ), LDA ) 00162 CALL CTRSM( 'Right', UPLO, 'No transpose', 00163 $ 'Non-unit', KB, N-K-KB+1, CONE, 00164 $ B( K+KB, K+KB ), LDB, A( K, K+KB ), 00165 $ LDA ) 00166 END IF 00167 10 CONTINUE 00168 ELSE 00169 * 00170 * Compute inv(L)*A*inv(L') 00171 * 00172 DO 20 K = 1, N, NB 00173 KB = MIN( N-K+1, NB ) 00174 * 00175 * Update the lower triangle of A(k:n,k:n) 00176 * 00177 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00178 $ B( K, K ), LDB, INFO ) 00179 IF( K+KB.LE.N ) THEN 00180 CALL CTRSM( 'Right', UPLO, 'Conjugate transpose', 00181 $ 'Non-unit', N-K-KB+1, KB, CONE, 00182 $ B( K, K ), LDB, A( K+KB, K ), LDA ) 00183 CALL CHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, 00184 $ A( K, K ), LDA, B( K+KB, K ), LDB, 00185 $ CONE, A( K+KB, K ), LDA ) 00186 CALL CHER2K( UPLO, 'No transpose', N-K-KB+1, KB, 00187 $ -CONE, A( K+KB, K ), LDA, 00188 $ B( K+KB, K ), LDB, ONE, 00189 $ A( K+KB, K+KB ), LDA ) 00190 CALL CHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, 00191 $ A( K, K ), LDA, B( K+KB, K ), LDB, 00192 $ CONE, A( K+KB, K ), LDA ) 00193 CALL CTRSM( 'Left', UPLO, 'No transpose', 00194 $ 'Non-unit', N-K-KB+1, KB, CONE, 00195 $ B( K+KB, K+KB ), LDB, A( K+KB, K ), 00196 $ LDA ) 00197 END IF 00198 20 CONTINUE 00199 END IF 00200 ELSE 00201 IF( UPPER ) THEN 00202 * 00203 * Compute U*A*U' 00204 * 00205 DO 30 K = 1, N, NB 00206 KB = MIN( N-K+1, NB ) 00207 * 00208 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1) 00209 * 00210 CALL CTRMM( 'Left', UPLO, 'No transpose', 'Non-unit', 00211 $ K-1, KB, CONE, B, LDB, A( 1, K ), LDA ) 00212 CALL CHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ), 00213 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ), 00214 $ LDA ) 00215 CALL CHER2K( UPLO, 'No transpose', K-1, KB, CONE, 00216 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A, 00217 $ LDA ) 00218 CALL CHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ), 00219 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ), 00220 $ LDA ) 00221 CALL CTRMM( 'Right', UPLO, 'Conjugate transpose', 00222 $ 'Non-unit', K-1, KB, CONE, B( K, K ), LDB, 00223 $ A( 1, K ), LDA ) 00224 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00225 $ B( K, K ), LDB, INFO ) 00226 30 CONTINUE 00227 ELSE 00228 * 00229 * Compute L'*A*L 00230 * 00231 DO 40 K = 1, N, NB 00232 KB = MIN( N-K+1, NB ) 00233 * 00234 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1) 00235 * 00236 CALL CTRMM( 'Right', UPLO, 'No transpose', 'Non-unit', 00237 $ KB, K-1, CONE, B, LDB, A( K, 1 ), LDA ) 00238 CALL CHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ), 00239 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ), 00240 $ LDA ) 00241 CALL CHER2K( UPLO, 'Conjugate transpose', K-1, KB, 00242 $ CONE, A( K, 1 ), LDA, B( K, 1 ), LDB, 00243 $ ONE, A, LDA ) 00244 CALL CHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ), 00245 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ), 00246 $ LDA ) 00247 CALL CTRMM( 'Left', UPLO, 'Conjugate transpose', 00248 $ 'Non-unit', KB, K-1, CONE, B( K, K ), LDB, 00249 $ A( K, 1 ), LDA ) 00250 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00251 $ B( K, K ), LDB, INFO ) 00252 40 CONTINUE 00253 END IF 00254 END IF 00255 END IF 00256 RETURN 00257 * 00258 * End of CHEGST 00259 * 00260 END
1.7.3 
