Annotation of rpl/lapack/lapack/dspgst.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSPGST( ITYPE, UPLO, N, AP, BP, 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, N
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION AP( * ), BP( * )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * DSPGST reduces a real symmetric-definite generalized eigenproblem
! 20: * to standard form, using packed storage.
! 21: *
! 22: * If ITYPE = 1, the problem is A*x = lambda*B*x,
! 23: * and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
! 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**T or L**T*A*L.
! 27: *
! 28: * B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
! 29: *
! 30: * Arguments
! 31: * =========
! 32: *
! 33: * ITYPE (input) INTEGER
! 34: * = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
! 35: * = 2 or 3: compute U*A*U**T or L**T*A*L.
! 36: *
! 37: * UPLO (input) CHARACTER*1
! 38: * = 'U': Upper triangle of A is stored and B is factored as
! 39: * U**T*U;
! 40: * = 'L': Lower triangle of A is stored and B is factored as
! 41: * L*L**T.
! 42: *
! 43: * N (input) INTEGER
! 44: * The order of the matrices A and B. N >= 0.
! 45: *
! 46: * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
! 47: * On entry, the upper or lower triangle of the symmetric matrix
! 48: * A, packed columnwise in a linear array. The j-th column of A
! 49: * is stored in the array AP as follows:
! 50: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
! 51: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
! 52: *
! 53: * On exit, if INFO = 0, the transformed matrix, stored in the
! 54: * same format as A.
! 55: *
! 56: * BP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
! 57: * The triangular factor from the Cholesky factorization of B,
! 58: * stored in the same format as A, as returned by DPPTRF.
! 59: *
! 60: * INFO (output) INTEGER
! 61: * = 0: successful exit
! 62: * < 0: if INFO = -i, the i-th argument had an illegal value
! 63: *
! 64: * =====================================================================
! 65: *
! 66: * .. Parameters ..
! 67: DOUBLE PRECISION ONE, HALF
! 68: PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 )
! 69: * ..
! 70: * .. Local Scalars ..
! 71: LOGICAL UPPER
! 72: INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK
! 73: DOUBLE PRECISION AJJ, AKK, BJJ, BKK, CT
! 74: * ..
! 75: * .. External Subroutines ..
! 76: EXTERNAL DAXPY, DSCAL, DSPMV, DSPR2, DTPMV, DTPSV,
! 77: $ XERBLA
! 78: * ..
! 79: * .. External Functions ..
! 80: LOGICAL LSAME
! 81: DOUBLE PRECISION DDOT
! 82: EXTERNAL LSAME, DDOT
! 83: * ..
! 84: * .. Executable Statements ..
! 85: *
! 86: * Test the input parameters.
! 87: *
! 88: INFO = 0
! 89: UPPER = LSAME( UPLO, 'U' )
! 90: IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
! 91: INFO = -1
! 92: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 93: INFO = -2
! 94: ELSE IF( N.LT.0 ) THEN
! 95: INFO = -3
! 96: END IF
! 97: IF( INFO.NE.0 ) THEN
! 98: CALL XERBLA( 'DSPGST', -INFO )
! 99: RETURN
! 100: END IF
! 101: *
! 102: IF( ITYPE.EQ.1 ) THEN
! 103: IF( UPPER ) THEN
! 104: *
! 105: * Compute inv(U')*A*inv(U)
! 106: *
! 107: * J1 and JJ are the indices of A(1,j) and A(j,j)
! 108: *
! 109: JJ = 0
! 110: DO 10 J = 1, N
! 111: J1 = JJ + 1
! 112: JJ = JJ + J
! 113: *
! 114: * Compute the j-th column of the upper triangle of A
! 115: *
! 116: BJJ = BP( JJ )
! 117: CALL DTPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
! 118: $ AP( J1 ), 1 )
! 119: CALL DSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
! 120: $ AP( J1 ), 1 )
! 121: CALL DSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
! 122: AP( JJ ) = ( AP( JJ )-DDOT( J-1, AP( J1 ), 1, BP( J1 ),
! 123: $ 1 ) ) / BJJ
! 124: 10 CONTINUE
! 125: ELSE
! 126: *
! 127: * Compute inv(L)*A*inv(L')
! 128: *
! 129: * KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
! 130: *
! 131: KK = 1
! 132: DO 20 K = 1, N
! 133: K1K1 = KK + N - K + 1
! 134: *
! 135: * Update the lower triangle of A(k:n,k:n)
! 136: *
! 137: AKK = AP( KK )
! 138: BKK = BP( KK )
! 139: AKK = AKK / BKK**2
! 140: AP( KK ) = AKK
! 141: IF( K.LT.N ) THEN
! 142: CALL DSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
! 143: CT = -HALF*AKK
! 144: CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
! 145: CALL DSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
! 146: $ BP( KK+1 ), 1, AP( K1K1 ) )
! 147: CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
! 148: CALL DTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
! 149: $ BP( K1K1 ), AP( KK+1 ), 1 )
! 150: END IF
! 151: KK = K1K1
! 152: 20 CONTINUE
! 153: END IF
! 154: ELSE
! 155: IF( UPPER ) THEN
! 156: *
! 157: * Compute U*A*U'
! 158: *
! 159: * K1 and KK are the indices of A(1,k) and A(k,k)
! 160: *
! 161: KK = 0
! 162: DO 30 K = 1, N
! 163: K1 = KK + 1
! 164: KK = KK + K
! 165: *
! 166: * Update the upper triangle of A(1:k,1:k)
! 167: *
! 168: AKK = AP( KK )
! 169: BKK = BP( KK )
! 170: CALL DTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
! 171: $ AP( K1 ), 1 )
! 172: CT = HALF*AKK
! 173: CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
! 174: CALL DSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
! 175: $ AP )
! 176: CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
! 177: CALL DSCAL( K-1, BKK, AP( K1 ), 1 )
! 178: AP( KK ) = AKK*BKK**2
! 179: 30 CONTINUE
! 180: ELSE
! 181: *
! 182: * Compute L'*A*L
! 183: *
! 184: * JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
! 185: *
! 186: JJ = 1
! 187: DO 40 J = 1, N
! 188: J1J1 = JJ + N - J + 1
! 189: *
! 190: * Compute the j-th column of the lower triangle of A
! 191: *
! 192: AJJ = AP( JJ )
! 193: BJJ = BP( JJ )
! 194: AP( JJ ) = AJJ*BJJ + DDOT( N-J, AP( JJ+1 ), 1,
! 195: $ BP( JJ+1 ), 1 )
! 196: CALL DSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
! 197: CALL DSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
! 198: $ ONE, AP( JJ+1 ), 1 )
! 199: CALL DTPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
! 200: $ BP( JJ ), AP( JJ ), 1 )
! 201: JJ = J1J1
! 202: 40 CONTINUE
! 203: END IF
! 204: END IF
! 205: RETURN
! 206: *
! 207: * End of DSPGST
! 208: *
! 209: END
CVSweb interface <joel.bertrand@systella.fr>