Annotation of rpl/lapack/lapack/dsygs2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 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*1
! 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsygs2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack