![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CHEGST 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CHEGST + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chegst.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chegst.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chegst.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, ITYPE, LDA, LDB, N 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX A( LDA, * ), B( LDB, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> CHEGST reduces a complex Hermitian-definite generalized 00038 *> eigenproblem to standard form. 00039 *> 00040 *> If ITYPE = 1, the problem is A*x = lambda*B*x, 00041 *> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H) 00042 *> 00043 *> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or 00044 *> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L. 00045 *> 00046 *> B must have been previously factorized as U**H*U or L*L**H by CPOTRF. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] ITYPE 00053 *> \verbatim 00054 *> ITYPE is INTEGER 00055 *> = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H); 00056 *> = 2 or 3: compute U*A*U**H or L**H*A*L. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] UPLO 00060 *> \verbatim 00061 *> UPLO is CHARACTER*1 00062 *> = 'U': Upper triangle of A is stored and B is factored as 00063 *> U**H*U; 00064 *> = 'L': Lower triangle of A is stored and B is factored as 00065 *> L*L**H. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] N 00069 *> \verbatim 00070 *> N is INTEGER 00071 *> The order of the matrices A and B. N >= 0. 00072 *> \endverbatim 00073 *> 00074 *> \param[in,out] A 00075 *> \verbatim 00076 *> A is COMPLEX array, dimension (LDA,N) 00077 *> On entry, the Hermitian matrix A. If UPLO = 'U', the leading 00078 *> N-by-N upper triangular part of A contains the upper 00079 *> triangular part of the matrix A, and the strictly lower 00080 *> triangular part of A is not referenced. If UPLO = 'L', the 00081 *> leading N-by-N lower triangular part of A contains the lower 00082 *> triangular part of the matrix A, and the strictly upper 00083 *> triangular part of A is not referenced. 00084 *> 00085 *> On exit, if INFO = 0, the transformed matrix, stored in the 00086 *> same format as A. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDA 00090 *> \verbatim 00091 *> LDA is INTEGER 00092 *> The leading dimension of the array A. LDA >= max(1,N). 00093 *> \endverbatim 00094 *> 00095 *> \param[in] B 00096 *> \verbatim 00097 *> B is COMPLEX array, dimension (LDB,N) 00098 *> The triangular factor from the Cholesky factorization of B, 00099 *> as returned by CPOTRF. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] LDB 00103 *> \verbatim 00104 *> LDB is INTEGER 00105 *> The leading dimension of the array B. LDB >= max(1,N). 00106 *> \endverbatim 00107 *> 00108 *> \param[out] INFO 00109 *> \verbatim 00110 *> INFO is INTEGER 00111 *> = 0: successful exit 00112 *> < 0: if INFO = -i, the i-th argument had an illegal value 00113 *> \endverbatim 00114 * 00115 * Authors: 00116 * ======== 00117 * 00118 *> \author Univ. of Tennessee 00119 *> \author Univ. of California Berkeley 00120 *> \author Univ. of Colorado Denver 00121 *> \author NAG Ltd. 00122 * 00123 *> \date November 2011 00124 * 00125 *> \ingroup complexHEcomputational 00126 * 00127 * ===================================================================== 00128 SUBROUTINE CHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00129 * 00130 * -- LAPACK computational routine (version 3.4.0) -- 00131 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00133 * November 2011 00134 * 00135 * .. Scalar Arguments .. 00136 CHARACTER UPLO 00137 INTEGER INFO, ITYPE, LDA, LDB, N 00138 * .. 00139 * .. Array Arguments .. 00140 COMPLEX A( LDA, * ), B( LDB, * ) 00141 * .. 00142 * 00143 * ===================================================================== 00144 * 00145 * .. Parameters .. 00146 REAL ONE 00147 PARAMETER ( ONE = 1.0E+0 ) 00148 COMPLEX CONE, HALF 00149 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 00150 $ HALF = ( 0.5E+0, 0.0E+0 ) ) 00151 * .. 00152 * .. Local Scalars .. 00153 LOGICAL UPPER 00154 INTEGER K, KB, NB 00155 * .. 00156 * .. External Subroutines .. 00157 EXTERNAL CHEGS2, CHEMM, CHER2K, CTRMM, CTRSM, XERBLA 00158 * .. 00159 * .. Intrinsic Functions .. 00160 INTRINSIC MAX, MIN 00161 * .. 00162 * .. External Functions .. 00163 LOGICAL LSAME 00164 INTEGER ILAENV 00165 EXTERNAL LSAME, ILAENV 00166 * .. 00167 * .. Executable Statements .. 00168 * 00169 * Test the input parameters. 00170 * 00171 INFO = 0 00172 UPPER = LSAME( UPLO, 'U' ) 00173 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00174 INFO = -1 00175 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00176 INFO = -2 00177 ELSE IF( N.LT.0 ) THEN 00178 INFO = -3 00179 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00180 INFO = -5 00181 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00182 INFO = -7 00183 END IF 00184 IF( INFO.NE.0 ) THEN 00185 CALL XERBLA( 'CHEGST', -INFO ) 00186 RETURN 00187 END IF 00188 * 00189 * Quick return if possible 00190 * 00191 IF( N.EQ.0 ) 00192 $ RETURN 00193 * 00194 * Determine the block size for this environment. 00195 * 00196 NB = ILAENV( 1, 'CHEGST', UPLO, N, -1, -1, -1 ) 00197 * 00198 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00199 * 00200 * Use unblocked code 00201 * 00202 CALL CHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00203 ELSE 00204 * 00205 * Use blocked code 00206 * 00207 IF( ITYPE.EQ.1 ) THEN 00208 IF( UPPER ) THEN 00209 * 00210 * Compute inv(U**H)*A*inv(U) 00211 * 00212 DO 10 K = 1, N, NB 00213 KB = MIN( N-K+1, NB ) 00214 * 00215 * Update the upper triangle of A(k:n,k:n) 00216 * 00217 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00218 $ B( K, K ), LDB, INFO ) 00219 IF( K+KB.LE.N ) THEN 00220 CALL CTRSM( 'Left', UPLO, 'Conjugate transpose', 00221 $ 'Non-unit', KB, N-K-KB+1, CONE, 00222 $ B( K, K ), LDB, A( K, K+KB ), LDA ) 00223 CALL CHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, 00224 $ A( K, K ), LDA, B( K, K+KB ), LDB, 00225 $ CONE, A( K, K+KB ), LDA ) 00226 CALL CHER2K( UPLO, 'Conjugate transpose', N-K-KB+1, 00227 $ KB, -CONE, A( K, K+KB ), LDA, 00228 $ B( K, K+KB ), LDB, ONE, 00229 $ A( K+KB, K+KB ), LDA ) 00230 CALL CHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, 00231 $ A( K, K ), LDA, B( K, K+KB ), LDB, 00232 $ CONE, A( K, K+KB ), LDA ) 00233 CALL CTRSM( 'Right', UPLO, 'No transpose', 00234 $ 'Non-unit', KB, N-K-KB+1, CONE, 00235 $ B( K+KB, K+KB ), LDB, A( K, K+KB ), 00236 $ LDA ) 00237 END IF 00238 10 CONTINUE 00239 ELSE 00240 * 00241 * Compute inv(L)*A*inv(L**H) 00242 * 00243 DO 20 K = 1, N, NB 00244 KB = MIN( N-K+1, NB ) 00245 * 00246 * Update the lower triangle of A(k:n,k:n) 00247 * 00248 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00249 $ B( K, K ), LDB, INFO ) 00250 IF( K+KB.LE.N ) THEN 00251 CALL CTRSM( 'Right', UPLO, 'Conjugate transpose', 00252 $ 'Non-unit', N-K-KB+1, KB, CONE, 00253 $ B( K, K ), LDB, A( K+KB, K ), LDA ) 00254 CALL CHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, 00255 $ A( K, K ), LDA, B( K+KB, K ), LDB, 00256 $ CONE, A( K+KB, K ), LDA ) 00257 CALL CHER2K( UPLO, 'No transpose', N-K-KB+1, KB, 00258 $ -CONE, A( K+KB, K ), LDA, 00259 $ B( K+KB, K ), LDB, ONE, 00260 $ A( K+KB, K+KB ), LDA ) 00261 CALL CHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, 00262 $ A( K, K ), LDA, B( K+KB, K ), LDB, 00263 $ CONE, A( K+KB, K ), LDA ) 00264 CALL CTRSM( 'Left', UPLO, 'No transpose', 00265 $ 'Non-unit', N-K-KB+1, KB, CONE, 00266 $ B( K+KB, K+KB ), LDB, A( K+KB, K ), 00267 $ LDA ) 00268 END IF 00269 20 CONTINUE 00270 END IF 00271 ELSE 00272 IF( UPPER ) THEN 00273 * 00274 * Compute U*A*U**H 00275 * 00276 DO 30 K = 1, N, NB 00277 KB = MIN( N-K+1, NB ) 00278 * 00279 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1) 00280 * 00281 CALL CTRMM( 'Left', UPLO, 'No transpose', 'Non-unit', 00282 $ K-1, KB, CONE, B, LDB, A( 1, K ), LDA ) 00283 CALL CHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ), 00284 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ), 00285 $ LDA ) 00286 CALL CHER2K( UPLO, 'No transpose', K-1, KB, CONE, 00287 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A, 00288 $ LDA ) 00289 CALL CHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ), 00290 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ), 00291 $ LDA ) 00292 CALL CTRMM( 'Right', UPLO, 'Conjugate transpose', 00293 $ 'Non-unit', K-1, KB, CONE, B( K, K ), LDB, 00294 $ A( 1, K ), LDA ) 00295 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00296 $ B( K, K ), LDB, INFO ) 00297 30 CONTINUE 00298 ELSE 00299 * 00300 * Compute L**H*A*L 00301 * 00302 DO 40 K = 1, N, NB 00303 KB = MIN( N-K+1, NB ) 00304 * 00305 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1) 00306 * 00307 CALL CTRMM( 'Right', UPLO, 'No transpose', 'Non-unit', 00308 $ KB, K-1, CONE, B, LDB, A( K, 1 ), LDA ) 00309 CALL CHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ), 00310 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ), 00311 $ LDA ) 00312 CALL CHER2K( UPLO, 'Conjugate transpose', K-1, KB, 00313 $ CONE, A( K, 1 ), LDA, B( K, 1 ), LDB, 00314 $ ONE, A, LDA ) 00315 CALL CHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ), 00316 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ), 00317 $ LDA ) 00318 CALL CTRMM( 'Left', UPLO, 'Conjugate transpose', 00319 $ 'Non-unit', KB, K-1, CONE, B( K, K ), LDB, 00320 $ A( K, 1 ), LDA ) 00321 CALL CHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00322 $ B( K, K ), LDB, INFO ) 00323 40 CONTINUE 00324 END IF 00325 END IF 00326 END IF 00327 RETURN 00328 * 00329 * End of CHEGST 00330 * 00331 END