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