 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * -- April 2011 -- 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 * DSYGS2 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 * Specifies whether the upper or lower triangular part of the 00039 * symmetric 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) 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 00081 DOUBLE PRECISION AKK, BKK, CT 00082 * .. 00083 * .. External Subroutines .. 00084 EXTERNAL DAXPY, DSCAL, DSYR2, DTRMV, DTRSV, XERBLA 00085 * .. 00086 * .. Intrinsic Functions .. 00087 INTRINSIC MAX 00088 * .. 00089 * .. External Functions .. 00090 LOGICAL LSAME 00091 EXTERNAL LSAME 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( 'DSYGS2', -INFO ) 00112 RETURN 00113 END IF 00114 * 00115 IF( ITYPE.EQ.1 ) THEN 00116 IF( UPPER ) THEN 00117 * 00118 * Compute inv(U**T)*A*inv(U) 00119 * 00120 DO 10 K = 1, N 00121 * 00122 * Update the upper triangle of A(k:n,k:n) 00123 * 00124 AKK = A( K, K ) 00125 BKK = B( K, K ) 00126 AKK = AKK / BKK**2 00127 A( K, K ) = AKK 00128 IF( K.LT.N ) THEN 00129 CALL DSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA ) 00130 CT = -HALF*AKK 00131 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 00132 $ LDA ) 00133 CALL DSYR2( UPLO, N-K, -ONE, A( K, K+1 ), LDA, 00134 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA ) 00135 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 00136 $ LDA ) 00137 CALL DTRSV( UPLO, 'Transpose', 'Non-unit', N-K, 00138 $ B( K+1, K+1 ), LDB, A( K, K+1 ), LDA ) 00139 END IF 00140 10 CONTINUE 00141 ELSE 00142 * 00143 * Compute inv(L)*A*inv(L**T) 00144 * 00145 DO 20 K = 1, N 00146 * 00147 * Update the lower triangle of A(k:n,k:n) 00148 * 00149 AKK = A( K, K ) 00150 BKK = B( K, K ) 00151 AKK = AKK / BKK**2 00152 A( K, K ) = AKK 00153 IF( K.LT.N ) THEN 00154 CALL DSCAL( N-K, ONE / BKK, A( K+1, K ), 1 ) 00155 CT = -HALF*AKK 00156 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 00157 CALL DSYR2( UPLO, N-K, -ONE, A( K+1, K ), 1, 00158 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA ) 00159 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 00160 CALL DTRSV( UPLO, 'No transpose', 'Non-unit', N-K, 00161 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 ) 00162 END IF 00163 20 CONTINUE 00164 END IF 00165 ELSE 00166 IF( UPPER ) THEN 00167 * 00168 * Compute U*A*U**T 00169 * 00170 DO 30 K = 1, N 00171 * 00172 * Update the upper triangle of A(1:k,1:k) 00173 * 00174 AKK = A( K, K ) 00175 BKK = B( K, K ) 00176 CALL DTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B, 00177 $ LDB, A( 1, K ), 1 ) 00178 CT = HALF*AKK 00179 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 00180 CALL DSYR2( UPLO, K-1, ONE, A( 1, K ), 1, B( 1, K ), 1, 00181 $ A, LDA ) 00182 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 00183 CALL DSCAL( K-1, BKK, A( 1, K ), 1 ) 00184 A( K, K ) = AKK*BKK**2 00185 30 CONTINUE 00186 ELSE 00187 * 00188 * Compute L**T *A*L 00189 * 00190 DO 40 K = 1, N 00191 * 00192 * Update the lower triangle of A(1:k,1:k) 00193 * 00194 AKK = A( K, K ) 00195 BKK = B( K, K ) 00196 CALL DTRMV( UPLO, 'Transpose', 'Non-unit', K-1, B, LDB, 00197 $ A( K, 1 ), LDA ) 00198 CT = HALF*AKK 00199 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 00200 CALL DSYR2( UPLO, K-1, ONE, A( K, 1 ), LDA, B( K, 1 ), 00201 $ LDB, A, LDA ) 00202 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 00203 CALL DSCAL( K-1, BKK, A( K, 1 ), LDA ) 00204 A( K, K ) = AKK*BKK**2 00205 40 CONTINUE 00206 END IF 00207 END IF 00208 RETURN 00209 * 00210 * End of DSYGS2 00211 * 00212 END
1.7.3 
