|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSYGST( 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DSYGST reduces a real symmetric-definite generalized eigenproblem 00020 * to standard form. 00021 * 00022 * If ITYPE = 1, the problem is A*x = lambda*B*x, 00023 * and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) 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**T or L**T*A*L. 00027 * 00028 * B must have been previously factorized as U**T*U or L*L**T by DPOTRF. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * ITYPE (input) INTEGER 00034 * = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T); 00035 * = 2 or 3: compute U*A*U**T or L**T*A*L. 00036 * 00037 * UPLO (input) CHARACTER*1 00038 * = 'U': Upper triangle of A is stored and B is factored as 00039 * U**T*U; 00040 * = 'L': Lower triangle of A is stored and B is factored as 00041 * L*L**T. 00042 * 00043 * N (input) INTEGER 00044 * The order of the matrices A and B. N >= 0. 00045 * 00046 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00047 * On entry, the symmetric 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) DOUBLE PRECISION array, dimension (LDB,N) 00062 * The triangular factor from the Cholesky factorization of B, 00063 * as returned by DPOTRF. 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 DOUBLE PRECISION ONE, HALF 00076 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 ) 00077 * .. 00078 * .. Local Scalars .. 00079 LOGICAL UPPER 00080 INTEGER K, KB, NB 00081 * .. 00082 * .. External Subroutines .. 00083 EXTERNAL DSYGS2, DSYMM, DSYR2K, DTRMM, DTRSM, XERBLA 00084 * .. 00085 * .. Intrinsic Functions .. 00086 INTRINSIC MAX, MIN 00087 * .. 00088 * .. External Functions .. 00089 LOGICAL LSAME 00090 INTEGER ILAENV 00091 EXTERNAL LSAME, ILAENV 00092 * .. 00093 * .. Executable Statements .. 00094 * 00095 * Test the input parameters. 00096 * 00097 INFO = 0 00098 UPPER = LSAME( UPLO, 'U' ) 00099 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00100 INFO = -1 00101 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00102 INFO = -2 00103 ELSE IF( N.LT.0 ) THEN 00104 INFO = -3 00105 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00106 INFO = -5 00107 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00108 INFO = -7 00109 END IF 00110 IF( INFO.NE.0 ) THEN 00111 CALL XERBLA( 'DSYGST', -INFO ) 00112 RETURN 00113 END IF 00114 * 00115 * Quick return if possible 00116 * 00117 IF( N.EQ.0 ) 00118 $ RETURN 00119 * 00120 * Determine the block size for this environment. 00121 * 00122 NB = ILAENV( 1, 'DSYGST', UPLO, N, -1, -1, -1 ) 00123 * 00124 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00125 * 00126 * Use unblocked code 00127 * 00128 CALL DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00129 ELSE 00130 * 00131 * Use blocked code 00132 * 00133 IF( ITYPE.EQ.1 ) THEN 00134 IF( UPPER ) THEN 00135 * 00136 * Compute inv(U')*A*inv(U) 00137 * 00138 DO 10 K = 1, N, NB 00139 KB = MIN( N-K+1, NB ) 00140 * 00141 * Update the upper triangle of A(k:n,k:n) 00142 * 00143 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00144 $ B( K, K ), LDB, INFO ) 00145 IF( K+KB.LE.N ) THEN 00146 CALL DTRSM( 'Left', UPLO, 'Transpose', 'Non-unit', 00147 $ KB, N-K-KB+1, ONE, B( K, K ), LDB, 00148 $ A( K, K+KB ), LDA ) 00149 CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, 00150 $ A( K, K ), LDA, B( K, K+KB ), LDB, ONE, 00151 $ A( K, K+KB ), LDA ) 00152 CALL DSYR2K( UPLO, 'Transpose', N-K-KB+1, KB, -ONE, 00153 $ A( K, K+KB ), LDA, B( K, K+KB ), LDB, 00154 $ ONE, A( K+KB, K+KB ), LDA ) 00155 CALL DSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, 00156 $ A( K, K ), LDA, B( K, K+KB ), LDB, ONE, 00157 $ A( K, K+KB ), LDA ) 00158 CALL DTRSM( 'Right', UPLO, 'No transpose', 00159 $ 'Non-unit', KB, N-K-KB+1, ONE, 00160 $ B( K+KB, K+KB ), LDB, A( K, K+KB ), 00161 $ LDA ) 00162 END IF 00163 10 CONTINUE 00164 ELSE 00165 * 00166 * Compute inv(L)*A*inv(L') 00167 * 00168 DO 20 K = 1, N, NB 00169 KB = MIN( N-K+1, NB ) 00170 * 00171 * Update the lower triangle of A(k:n,k:n) 00172 * 00173 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00174 $ B( K, K ), LDB, INFO ) 00175 IF( K+KB.LE.N ) THEN 00176 CALL DTRSM( 'Right', UPLO, 'Transpose', 'Non-unit', 00177 $ N-K-KB+1, KB, ONE, B( K, K ), LDB, 00178 $ A( K+KB, K ), LDA ) 00179 CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, 00180 $ A( K, K ), LDA, B( K+KB, K ), LDB, ONE, 00181 $ A( K+KB, K ), LDA ) 00182 CALL DSYR2K( UPLO, 'No transpose', N-K-KB+1, KB, 00183 $ -ONE, A( K+KB, K ), LDA, B( K+KB, K ), 00184 $ LDB, ONE, A( K+KB, K+KB ), LDA ) 00185 CALL DSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, 00186 $ A( K, K ), LDA, B( K+KB, K ), LDB, ONE, 00187 $ A( K+KB, K ), LDA ) 00188 CALL DTRSM( 'Left', UPLO, 'No transpose', 00189 $ 'Non-unit', N-K-KB+1, KB, ONE, 00190 $ B( K+KB, K+KB ), LDB, A( K+KB, K ), 00191 $ LDA ) 00192 END IF 00193 20 CONTINUE 00194 END IF 00195 ELSE 00196 IF( UPPER ) THEN 00197 * 00198 * Compute U*A*U' 00199 * 00200 DO 30 K = 1, N, NB 00201 KB = MIN( N-K+1, NB ) 00202 * 00203 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1) 00204 * 00205 CALL DTRMM( 'Left', UPLO, 'No transpose', 'Non-unit', 00206 $ K-1, KB, ONE, B, LDB, A( 1, K ), LDA ) 00207 CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ), 00208 $ LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA ) 00209 CALL DSYR2K( UPLO, 'No transpose', K-1, KB, ONE, 00210 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A, 00211 $ LDA ) 00212 CALL DSYMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ), 00213 $ LDA, B( 1, K ), LDB, ONE, A( 1, K ), LDA ) 00214 CALL DTRMM( 'Right', UPLO, 'Transpose', 'Non-unit', 00215 $ K-1, KB, ONE, B( K, K ), LDB, A( 1, K ), 00216 $ LDA ) 00217 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00218 $ B( K, K ), LDB, INFO ) 00219 30 CONTINUE 00220 ELSE 00221 * 00222 * Compute L'*A*L 00223 * 00224 DO 40 K = 1, N, NB 00225 KB = MIN( N-K+1, NB ) 00226 * 00227 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1) 00228 * 00229 CALL DTRMM( 'Right', UPLO, 'No transpose', 'Non-unit', 00230 $ KB, K-1, ONE, B, LDB, A( K, 1 ), LDA ) 00231 CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ), 00232 $ LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA ) 00233 CALL DSYR2K( UPLO, 'Transpose', K-1, KB, ONE, 00234 $ A( K, 1 ), LDA, B( K, 1 ), LDB, ONE, A, 00235 $ LDA ) 00236 CALL DSYMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ), 00237 $ LDA, B( K, 1 ), LDB, ONE, A( K, 1 ), LDA ) 00238 CALL DTRMM( 'Left', UPLO, 'Transpose', 'Non-unit', KB, 00239 $ K-1, ONE, B( K, K ), LDB, A( K, 1 ), LDA ) 00240 CALL DSYGS2( ITYPE, UPLO, KB, A( K, K ), LDA, 00241 $ B( K, K ), LDB, INFO ) 00242 40 CONTINUE 00243 END IF 00244 END IF 00245 END IF 00246 RETURN 00247 * 00248 * End of DSYGST 00249 * 00250 END
1.7.3 
