1 SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
2 $ WORK, INFO )
3 *
4 * -- LAPACK PROTOTYPE routine (version 3.3.1) --
5 *
6 * -- Written by Julie Langou of the Univ. of TN --
7 * -- April 2011 --
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 complex
25 * Hermitian matrix A using the factorization A = U*D*U**H or
26 * A = L*D*L**H computed by ZHETRF 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**H.
122 *
123 * P**T * 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**T * B) -> B [ (U \P**T * B) ]
144 *
145 CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
146 *
147 * Compute D \ B -> B [ D \ (U \P**T * 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**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
173 *
174 CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,LDA,B,LDB)
175 *
176 * P * B [ P * (U**H \ (D \ (U \P**T * 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**H.
200 *
201 * P**T * 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**T * B) -> B [ (L \P**T * B) ]
222 *
223 CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
224 *
225 * Compute D \ B -> B [ D \ (L \P**T * 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**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
249 *
250 CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,LDA,B,LDB)
251 *
252 * P * B [ P * (L**H \ (D \ (L \P**T * 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
2 $ WORK, INFO )
3 *
4 * -- LAPACK PROTOTYPE routine (version 3.3.1) --
5 *
6 * -- Written by Julie Langou of the Univ. of TN --
7 * -- April 2011 --
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 complex
25 * Hermitian matrix A using the factorization A = U*D*U**H or
26 * A = L*D*L**H computed by ZHETRF 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**H.
122 *
123 * P**T * 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**T * B) -> B [ (U \P**T * B) ]
144 *
145 CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
146 *
147 * Compute D \ B -> B [ D \ (U \P**T * 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**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
173 *
174 CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,LDA,B,LDB)
175 *
176 * P * B [ P * (U**H \ (D \ (U \P**T * 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**H.
200 *
201 * P**T * 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**T * B) -> B [ (L \P**T * B) ]
222 *
223 CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
224 *
225 * Compute D \ B -> B [ D \ (L \P**T * 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**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
249 *
250 CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,LDA,B,LDB)
251 *
252 * P * B [ P * (L**H \ (D \ (L \P**T * 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