![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dspgst.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dspgst.f CVS log | Up to [local] / rpl / lapack / lapack |
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