1: SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
2: $ WORK, INFO )
3: *
4: * -- LAPACK PROTOTYPE routine (version 3.3.0) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2010
8: *
9: * -- Written by Julie Langou of the Univ. of TN --
10: *
11: * .. Scalar Arguments ..
12: CHARACTER UPLO
13: INTEGER INFO, LDA, LDB, N, NRHS
14: * ..
15: * .. Array Arguments ..
16: INTEGER IPIV( * )
17: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * DSYTRS2 solves a system of linear equations A*X = B with a real
24: * symmetric matrix A using the factorization A = U*D*U**T or
25: * A = L*D*L**T computed by DSYTRF and converted by DSYCONV.
26: *
27: * Arguments
28: * =========
29: *
30: * UPLO (input) CHARACTER*1
31: * Specifies whether the details of the factorization are stored
32: * as an upper or lower triangular matrix.
33: * = 'U': Upper triangular, form is A = U*D*U**T;
34: * = 'L': Lower triangular, form is A = L*D*L**T.
35: *
36: * N (input) INTEGER
37: * The order of the matrix A. N >= 0.
38: *
39: * NRHS (input) INTEGER
40: * The number of right hand sides, i.e., the number of columns
41: * of the matrix B. NRHS >= 0.
42: *
43: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
44: * The block diagonal matrix D and the multipliers used to
45: * obtain the factor U or L as computed by DSYTRF.
46: *
47: * LDA (input) INTEGER
48: * The leading dimension of the array A. LDA >= max(1,N).
49: *
50: * IPIV (input) INTEGER array, dimension (N)
51: * Details of the interchanges and the block structure of D
52: * as determined by DSYTRF.
53: *
54: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
55: * On entry, the right hand side matrix B.
56: * On exit, the solution matrix X.
57: *
58: * LDB (input) INTEGER
59: * The leading dimension of the array B. LDB >= max(1,N).
60: *
61: * WORK (workspace) REAL array, dimension (N)
62: *
63: * INFO (output) INTEGER
64: * = 0: successful exit
65: * < 0: if INFO = -i, the i-th argument had an illegal value
66: *
67: * =====================================================================
68: *
69: * .. Parameters ..
70: DOUBLE PRECISION ONE
71: PARAMETER ( ONE = 1.0D+0 )
72: * ..
73: * .. Local Scalars ..
74: LOGICAL UPPER
75: INTEGER I, IINFO, J, K, KP
76: DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
77: * ..
78: * .. External Functions ..
79: LOGICAL LSAME
80: EXTERNAL LSAME
81: * ..
82: * .. External Subroutines ..
83: EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA
84: * ..
85: * .. Intrinsic Functions ..
86: INTRINSIC MAX
87: * ..
88: * .. Executable Statements ..
89: *
90: INFO = 0
91: UPPER = LSAME( UPLO, 'U' )
92: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
93: INFO = -1
94: ELSE IF( N.LT.0 ) THEN
95: INFO = -2
96: ELSE IF( NRHS.LT.0 ) THEN
97: INFO = -3
98: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
99: INFO = -5
100: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
101: INFO = -8
102: END IF
103: IF( INFO.NE.0 ) THEN
104: CALL XERBLA( 'DSYTRS2', -INFO )
105: RETURN
106: END IF
107: *
108: * Quick return if possible
109: *
110: IF( N.EQ.0 .OR. NRHS.EQ.0 )
111: $ RETURN
112: *
113: * Convert A
114: *
115: CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
116: *
117: IF( UPPER ) THEN
118: *
119: * Solve A*X = B, where A = U*D*U'.
120: *
121: * P' * B
122: K=N
123: DO WHILE ( K .GE. 1 )
124: IF( IPIV( K ).GT.0 ) THEN
125: * 1 x 1 diagonal block
126: * Interchange rows K and IPIV(K).
127: KP = IPIV( K )
128: IF( KP.NE.K )
129: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
130: K=K-1
131: ELSE
132: * 2 x 2 diagonal block
133: * Interchange rows K-1 and -IPIV(K).
134: KP = -IPIV( K )
135: IF( KP.EQ.-IPIV( K-1 ) )
136: $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
137: K=K-2
138: END IF
139: END DO
140: *
141: * Compute (U \P' * B) -> B [ (U \P' * B) ]
142: *
143: CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
144: *
145: * Compute D \ B -> B [ D \ (U \P' * B) ]
146: *
147: I=N
148: DO WHILE ( I .GE. 1 )
149: IF( IPIV(I) .GT. 0 ) THEN
150: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
151: ELSEIF ( I .GT. 1) THEN
152: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
153: AKM1K = WORK(I)
154: AKM1 = A( I-1, I-1 ) / AKM1K
155: AK = A( I, I ) / AKM1K
156: DENOM = AKM1*AK - ONE
157: DO 15 J = 1, NRHS
158: BKM1 = B( I-1, J ) / AKM1K
159: BK = B( I, J ) / AKM1K
160: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
161: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
162: 15 CONTINUE
163: I = I - 1
164: ENDIF
165: ENDIF
166: I = I - 1
167: END DO
168: *
169: * Compute (U' \ B) -> B [ U' \ (D \ (U \P' * B) ) ]
170: *
171: CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N)
172: *
173: * P * B [ P * (U' \ (D \ (U \P' * B) )) ]
174: *
175: K=1
176: DO WHILE ( K .LE. N )
177: IF( IPIV( K ).GT.0 ) THEN
178: * 1 x 1 diagonal block
179: * Interchange rows K and IPIV(K).
180: KP = IPIV( K )
181: IF( KP.NE.K )
182: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
183: K=K+1
184: ELSE
185: * 2 x 2 diagonal block
186: * Interchange rows K-1 and -IPIV(K).
187: KP = -IPIV( K )
188: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
189: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
190: K=K+2
191: ENDIF
192: END DO
193: *
194: ELSE
195: *
196: * Solve A*X = B, where A = L*D*L'.
197: *
198: * P' * B
199: K=1
200: DO WHILE ( K .LE. N )
201: IF( IPIV( K ).GT.0 ) THEN
202: * 1 x 1 diagonal block
203: * Interchange rows K and IPIV(K).
204: KP = IPIV( K )
205: IF( KP.NE.K )
206: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
207: K=K+1
208: ELSE
209: * 2 x 2 diagonal block
210: * Interchange rows K and -IPIV(K+1).
211: KP = -IPIV( K+1 )
212: IF( KP.EQ.-IPIV( K ) )
213: $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
214: K=K+2
215: ENDIF
216: END DO
217: *
218: * Compute (L \P' * B) -> B [ (L \P' * B) ]
219: *
220: CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
221: *
222: * Compute D \ B -> B [ D \ (L \P' * B) ]
223: *
224: I=1
225: DO WHILE ( I .LE. N )
226: IF( IPIV(I) .GT. 0 ) THEN
227: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
228: ELSE
229: AKM1K = WORK(I)
230: AKM1 = A( I, I ) / AKM1K
231: AK = A( I+1, I+1 ) / AKM1K
232: DENOM = AKM1*AK - ONE
233: DO 25 J = 1, NRHS
234: BKM1 = B( I, J ) / AKM1K
235: BK = B( I+1, J ) / AKM1K
236: B( I, J ) = ( AK*BKM1-BK ) / DENOM
237: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
238: 25 CONTINUE
239: I = I + 1
240: ENDIF
241: I = I + 1
242: END DO
243: *
244: * Compute (L' \ B) -> B [ L' \ (D \ (L \P' * B) ) ]
245: *
246: CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
247: *
248: * P * B [ P * (L' \ (D \ (L \P' * B) )) ]
249: *
250: K=N
251: DO WHILE ( K .GE. 1 )
252: IF( IPIV( K ).GT.0 ) THEN
253: * 1 x 1 diagonal block
254: * Interchange rows K and IPIV(K).
255: KP = IPIV( K )
256: IF( KP.NE.K )
257: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
258: K=K-1
259: ELSE
260: * 2 x 2 diagonal block
261: * Interchange rows K-1 and -IPIV(K).
262: KP = -IPIV( K )
263: IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
264: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
265: K=K-2
266: ENDIF
267: END DO
268: *
269: END IF
270: *
271: * Revert A
272: *
273: CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
274: *
275: RETURN
276: *
277: * End of DSYTRS2
278: *
279: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrs2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack