1 SUBROUTINE DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 2* 3* -- LAPACK routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* February 29, 1992 7* 8* .. Scalar Arguments .. 9 CHARACTER UPLO 10 INTEGER INFO, ITYPE, LDA, LDB, N 11* .. 12* .. Array Arguments .. 13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* DSYGS2 reduces a real symmetric-definite generalized eigenproblem 20* to standard form. 21* 22* If ITYPE = 1, the problem is A*x = lambda*B*x, 23* and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L') 24* 25* If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or 26* B*A*x = lambda*x, and A is overwritten by U*A*U` or L'*A*L. 27* 28* B must have been previously factorized as U'*U or L*L' by DPOTRF. 29* 30* Arguments 31* ========= 32* 33* ITYPE (input) INTEGER 34* = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L'); 35* = 2 or 3: compute U*A*U' or L'*A*L. 36* 37* UPLO (input) CHARACTER 38* Specifies whether the upper or lower triangular part of the 39* symmetric matrix A is stored, and how B has been factorized. 40* = 'U': Upper triangular 41* = 'L': Lower triangular 42* 43* N (input) INTEGER 44* The order of the matrices A and B. N >= 0. 45* 46* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 47* On entry, the symmetric matrix A. If UPLO = 'U', the leading 48* n by n upper triangular part of A contains the upper 49* triangular part of the matrix A, and the strictly lower 50* triangular part of A is not referenced. If UPLO = 'L', the 51* leading n by n lower triangular part of A contains the lower 52* triangular part of the matrix A, and the strictly upper 53* triangular part of A is not referenced. 54* 55* On exit, if INFO = 0, the transformed matrix, stored in the 56* same format as A. 57* 58* LDA (input) INTEGER 59* The leading dimension of the array A. LDA >= max(1,N). 60* 61* B (input) DOUBLE PRECISION array, dimension (LDB,N) 62* The triangular factor from the Cholesky factorization of B, 63* as returned by DPOTRF. 64* 65* LDB (input) INTEGER 66* The leading dimension of the array B. LDB >= max(1,N). 67* 68* INFO (output) INTEGER 69* = 0: successful exit. 70* < 0: if INFO = -i, the i-th argument had an illegal value. 71* 72* ===================================================================== 73* 74* .. Parameters .. 75 DOUBLE PRECISION ONE, HALF 76 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 ) 77* .. 78* .. Local Scalars .. 79 LOGICAL UPPER 80 INTEGER K 81 DOUBLE PRECISION AKK, BKK, CT 82* .. 83* .. External Subroutines .. 84 EXTERNAL DAXPY, DSCAL, DSYR2, DTRMV, DTRSV, XERBLA 85* .. 86* .. Intrinsic Functions .. 87 INTRINSIC MAX 88* .. 89* .. External Functions .. 90 LOGICAL LSAME 91 EXTERNAL LSAME 92* .. 93* .. Executable Statements .. 94* 95* Test the input parameters. 96* 97 INFO = 0 98 UPPER = LSAME( UPLO, 'U' ) 99 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 100 INFO = -1 101 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 102 INFO = -2 103 ELSE IF( N.LT.0 ) THEN 104 INFO = -3 105 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 106 INFO = -5 107 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 108 INFO = -7 109 END IF 110 IF( INFO.NE.0 ) THEN 111 CALL XERBLA( 'DSYGS2', -INFO ) 112 RETURN 113 END IF 114* 115 IF( ITYPE.EQ.1 ) THEN 116 IF( UPPER ) THEN 117* 118* Compute inv(U')*A*inv(U) 119* 120 DO 10 K = 1, N 121* 122* Update the upper triangle of A(k:n,k:n) 123* 124 AKK = A( K, K ) 125 BKK = B( K, K ) 126 AKK = AKK / BKK**2 127 A( K, K ) = AKK 128 IF( K.LT.N ) THEN 129 CALL DSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA ) 130 CT = -HALF*AKK 131 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 132 $ LDA ) 133 CALL DSYR2( UPLO, N-K, -ONE, A( K, K+1 ), LDA, 134 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA ) 135 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 136 $ LDA ) 137 CALL DTRSV( UPLO, 'Transpose', 'Non-unit', N-K, 138 $ B( K+1, K+1 ), LDB, A( K, K+1 ), LDA ) 139 END IF 140 10 CONTINUE 141 ELSE 142* 143* Compute inv(L)*A*inv(L') 144* 145 DO 20 K = 1, N 146* 147* Update the lower triangle of A(k:n,k:n) 148* 149 AKK = A( K, K ) 150 BKK = B( K, K ) 151 AKK = AKK / BKK**2 152 A( K, K ) = AKK 153 IF( K.LT.N ) THEN 154 CALL DSCAL( N-K, ONE / BKK, A( K+1, K ), 1 ) 155 CT = -HALF*AKK 156 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 157 CALL DSYR2( UPLO, N-K, -ONE, A( K+1, K ), 1, 158 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA ) 159 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 160 CALL DTRSV( UPLO, 'No transpose', 'Non-unit', N-K, 161 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 ) 162 END IF 163 20 CONTINUE 164 END IF 165 ELSE 166 IF( UPPER ) THEN 167* 168* Compute U*A*U' 169* 170 DO 30 K = 1, N 171* 172* Update the upper triangle of A(1:k,1:k) 173* 174 AKK = A( K, K ) 175 BKK = B( K, K ) 176 CALL DTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B, 177 $ LDB, A( 1, K ), 1 ) 178 CT = HALF*AKK 179 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 180 CALL DSYR2( UPLO, K-1, ONE, A( 1, K ), 1, B( 1, K ), 1, 181 $ A, LDA ) 182 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 183 CALL DSCAL( K-1, BKK, A( 1, K ), 1 ) 184 A( K, K ) = AKK*BKK**2 185 30 CONTINUE 186 ELSE 187* 188* Compute L'*A*L 189* 190 DO 40 K = 1, N 191* 192* Update the lower triangle of A(1:k,1:k) 193* 194 AKK = A( K, K ) 195 BKK = B( K, K ) 196 CALL DTRMV( UPLO, 'Transpose', 'Non-unit', K-1, B, LDB, 197 $ A( K, 1 ), LDA ) 198 CT = HALF*AKK 199 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 200 CALL DSYR2( UPLO, K-1, ONE, A( K, 1 ), LDA, B( K, 1 ), 201 $ LDB, A, LDA ) 202 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 203 CALL DSCAL( K-1, BKK, A( K, 1 ), LDA ) 204 A( K, K ) = AKK*BKK**2 205 40 CONTINUE 206 END IF 207 END IF 208 RETURN 209* 210* End of DSYGS2 211* 212 END 213