1 SUBROUTINE ZHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, LDA, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( LDA, * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHETRS solves a system of linear equations A*X = B with a complex
21 * Hermitian matrix A using the factorization A = U*D*U**H or
22 * A = L*D*L**H computed by ZHETRF.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the details of the factorization are stored
29 * as an upper or lower triangular matrix.
30 * = 'U': Upper triangular, form is A = U*D*U**H;
31 * = 'L': Lower triangular, form is A = L*D*L**H.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * NRHS (input) INTEGER
37 * The number of right hand sides, i.e., the number of columns
38 * of the matrix B. NRHS >= 0.
39 *
40 * A (input) COMPLEX*16 array, dimension (LDA,N)
41 * The block diagonal matrix D and the multipliers used to
42 * obtain the factor U or L as computed by ZHETRF.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,N).
46 *
47 * IPIV (input) INTEGER array, dimension (N)
48 * Details of the interchanges and the block structure of D
49 * as determined by ZHETRF.
50 *
51 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
52 * On entry, the right hand side matrix B.
53 * On exit, the solution matrix X.
54 *
55 * LDB (input) INTEGER
56 * The leading dimension of the array B. LDB >= max(1,N).
57 *
58 * INFO (output) INTEGER
59 * = 0: successful exit
60 * < 0: if INFO = -i, the i-th argument had an illegal value
61 *
62 * =====================================================================
63 *
64 * .. Parameters ..
65 COMPLEX*16 ONE
66 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
67 * ..
68 * .. Local Scalars ..
69 LOGICAL UPPER
70 INTEGER J, K, KP
71 DOUBLE PRECISION S
72 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
73 * ..
74 * .. External Functions ..
75 LOGICAL LSAME
76 EXTERNAL LSAME
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
80 * ..
81 * .. Intrinsic Functions ..
82 INTRINSIC DBLE, DCONJG, MAX
83 * ..
84 * .. Executable Statements ..
85 *
86 INFO = 0
87 UPPER = LSAME( UPLO, 'U' )
88 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89 INFO = -1
90 ELSE IF( N.LT.0 ) THEN
91 INFO = -2
92 ELSE IF( NRHS.LT.0 ) THEN
93 INFO = -3
94 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
95 INFO = -5
96 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
97 INFO = -8
98 END IF
99 IF( INFO.NE.0 ) THEN
100 CALL XERBLA( 'ZHETRS', -INFO )
101 RETURN
102 END IF
103 *
104 * Quick return if possible
105 *
106 IF( N.EQ.0 .OR. NRHS.EQ.0 )
107 $ RETURN
108 *
109 IF( UPPER ) THEN
110 *
111 * Solve A*X = B, where A = U*D*U**H.
112 *
113 * First solve U*D*X = B, overwriting B with X.
114 *
115 * K is the main loop index, decreasing from N to 1 in steps of
116 * 1 or 2, depending on the size of the diagonal blocks.
117 *
118 K = N
119 10 CONTINUE
120 *
121 * If K < 1, exit from loop.
122 *
123 IF( K.LT.1 )
124 $ GO TO 30
125 *
126 IF( IPIV( K ).GT.0 ) THEN
127 *
128 * 1 x 1 diagonal block
129 *
130 * Interchange rows K and IPIV(K).
131 *
132 KP = IPIV( K )
133 IF( KP.NE.K )
134 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
135 *
136 * Multiply by inv(U(K)), where U(K) is the transformation
137 * stored in column K of A.
138 *
139 CALL ZGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
140 $ B( 1, 1 ), LDB )
141 *
142 * Multiply by the inverse of the diagonal block.
143 *
144 S = DBLE( ONE ) / DBLE( A( K, K ) )
145 CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
146 K = K - 1
147 ELSE
148 *
149 * 2 x 2 diagonal block
150 *
151 * Interchange rows K-1 and -IPIV(K).
152 *
153 KP = -IPIV( K )
154 IF( KP.NE.K-1 )
155 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
156 *
157 * Multiply by inv(U(K)), where U(K) is the transformation
158 * stored in columns K-1 and K of A.
159 *
160 CALL ZGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
161 $ B( 1, 1 ), LDB )
162 CALL ZGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
163 $ LDB, B( 1, 1 ), LDB )
164 *
165 * Multiply by the inverse of the diagonal block.
166 *
167 AKM1K = A( K-1, K )
168 AKM1 = A( K-1, K-1 ) / AKM1K
169 AK = A( K, K ) / DCONJG( AKM1K )
170 DENOM = AKM1*AK - ONE
171 DO 20 J = 1, NRHS
172 BKM1 = B( K-1, J ) / AKM1K
173 BK = B( K, J ) / DCONJG( AKM1K )
174 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
175 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
176 20 CONTINUE
177 K = K - 2
178 END IF
179 *
180 GO TO 10
181 30 CONTINUE
182 *
183 * Next solve U**H *X = B, overwriting B with X.
184 *
185 * K is the main loop index, increasing from 1 to N in steps of
186 * 1 or 2, depending on the size of the diagonal blocks.
187 *
188 K = 1
189 40 CONTINUE
190 *
191 * If K > N, exit from loop.
192 *
193 IF( K.GT.N )
194 $ GO TO 50
195 *
196 IF( IPIV( K ).GT.0 ) THEN
197 *
198 * 1 x 1 diagonal block
199 *
200 * Multiply by inv(U**H(K)), where U(K) is the transformation
201 * stored in column K of A.
202 *
203 IF( K.GT.1 ) THEN
204 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
205 CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
206 $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
207 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
208 END IF
209 *
210 * Interchange rows K and IPIV(K).
211 *
212 KP = IPIV( K )
213 IF( KP.NE.K )
214 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
215 K = K + 1
216 ELSE
217 *
218 * 2 x 2 diagonal block
219 *
220 * Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
221 * stored in columns K and K+1 of A.
222 *
223 IF( K.GT.1 ) THEN
224 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
225 CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
226 $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
227 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
228 *
229 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
230 CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
231 $ LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
232 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
233 END IF
234 *
235 * Interchange rows K and -IPIV(K).
236 *
237 KP = -IPIV( K )
238 IF( KP.NE.K )
239 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
240 K = K + 2
241 END IF
242 *
243 GO TO 40
244 50 CONTINUE
245 *
246 ELSE
247 *
248 * Solve A*X = B, where A = L*D*L**H.
249 *
250 * First solve L*D*X = B, overwriting B with X.
251 *
252 * K is the main loop index, increasing from 1 to N in steps of
253 * 1 or 2, depending on the size of the diagonal blocks.
254 *
255 K = 1
256 60 CONTINUE
257 *
258 * If K > N, exit from loop.
259 *
260 IF( K.GT.N )
261 $ GO TO 80
262 *
263 IF( IPIV( K ).GT.0 ) THEN
264 *
265 * 1 x 1 diagonal block
266 *
267 * Interchange rows K and IPIV(K).
268 *
269 KP = IPIV( K )
270 IF( KP.NE.K )
271 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
272 *
273 * Multiply by inv(L(K)), where L(K) is the transformation
274 * stored in column K of A.
275 *
276 IF( K.LT.N )
277 $ CALL ZGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
278 $ LDB, B( K+1, 1 ), LDB )
279 *
280 * Multiply by the inverse of the diagonal block.
281 *
282 S = DBLE( ONE ) / DBLE( A( K, K ) )
283 CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
284 K = K + 1
285 ELSE
286 *
287 * 2 x 2 diagonal block
288 *
289 * Interchange rows K+1 and -IPIV(K).
290 *
291 KP = -IPIV( K )
292 IF( KP.NE.K+1 )
293 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
294 *
295 * Multiply by inv(L(K)), where L(K) is the transformation
296 * stored in columns K and K+1 of A.
297 *
298 IF( K.LT.N-1 ) THEN
299 CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
300 $ LDB, B( K+2, 1 ), LDB )
301 CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
302 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
303 END IF
304 *
305 * Multiply by the inverse of the diagonal block.
306 *
307 AKM1K = A( K+1, K )
308 AKM1 = A( K, K ) / DCONJG( AKM1K )
309 AK = A( K+1, K+1 ) / AKM1K
310 DENOM = AKM1*AK - ONE
311 DO 70 J = 1, NRHS
312 BKM1 = B( K, J ) / DCONJG( AKM1K )
313 BK = B( K+1, J ) / AKM1K
314 B( K, J ) = ( AK*BKM1-BK ) / DENOM
315 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
316 70 CONTINUE
317 K = K + 2
318 END IF
319 *
320 GO TO 60
321 80 CONTINUE
322 *
323 * Next solve L**H *X = B, overwriting B with X.
324 *
325 * K is the main loop index, decreasing from N to 1 in steps of
326 * 1 or 2, depending on the size of the diagonal blocks.
327 *
328 K = N
329 90 CONTINUE
330 *
331 * If K < 1, exit from loop.
332 *
333 IF( K.LT.1 )
334 $ GO TO 100
335 *
336 IF( IPIV( K ).GT.0 ) THEN
337 *
338 * 1 x 1 diagonal block
339 *
340 * Multiply by inv(L**H(K)), where L(K) is the transformation
341 * stored in column K of A.
342 *
343 IF( K.LT.N ) THEN
344 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
345 CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
346 $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
347 $ B( K, 1 ), LDB )
348 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
349 END IF
350 *
351 * Interchange rows K and IPIV(K).
352 *
353 KP = IPIV( K )
354 IF( KP.NE.K )
355 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
356 K = K - 1
357 ELSE
358 *
359 * 2 x 2 diagonal block
360 *
361 * Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
362 * stored in columns K-1 and K of A.
363 *
364 IF( K.LT.N ) THEN
365 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
366 CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
367 $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
368 $ B( K, 1 ), LDB )
369 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
370 *
371 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
372 CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
373 $ B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
374 $ B( K-1, 1 ), LDB )
375 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
376 END IF
377 *
378 * Interchange rows K and -IPIV(K).
379 *
380 KP = -IPIV( K )
381 IF( KP.NE.K )
382 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
383 K = K - 2
384 END IF
385 *
386 GO TO 90
387 100 CONTINUE
388 END IF
389 *
390 RETURN
391 *
392 * End of ZHETRS
393 *
394 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, LDA, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( LDA, * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHETRS solves a system of linear equations A*X = B with a complex
21 * Hermitian matrix A using the factorization A = U*D*U**H or
22 * A = L*D*L**H computed by ZHETRF.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the details of the factorization are stored
29 * as an upper or lower triangular matrix.
30 * = 'U': Upper triangular, form is A = U*D*U**H;
31 * = 'L': Lower triangular, form is A = L*D*L**H.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * NRHS (input) INTEGER
37 * The number of right hand sides, i.e., the number of columns
38 * of the matrix B. NRHS >= 0.
39 *
40 * A (input) COMPLEX*16 array, dimension (LDA,N)
41 * The block diagonal matrix D and the multipliers used to
42 * obtain the factor U or L as computed by ZHETRF.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,N).
46 *
47 * IPIV (input) INTEGER array, dimension (N)
48 * Details of the interchanges and the block structure of D
49 * as determined by ZHETRF.
50 *
51 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
52 * On entry, the right hand side matrix B.
53 * On exit, the solution matrix X.
54 *
55 * LDB (input) INTEGER
56 * The leading dimension of the array B. LDB >= max(1,N).
57 *
58 * INFO (output) INTEGER
59 * = 0: successful exit
60 * < 0: if INFO = -i, the i-th argument had an illegal value
61 *
62 * =====================================================================
63 *
64 * .. Parameters ..
65 COMPLEX*16 ONE
66 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
67 * ..
68 * .. Local Scalars ..
69 LOGICAL UPPER
70 INTEGER J, K, KP
71 DOUBLE PRECISION S
72 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
73 * ..
74 * .. External Functions ..
75 LOGICAL LSAME
76 EXTERNAL LSAME
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
80 * ..
81 * .. Intrinsic Functions ..
82 INTRINSIC DBLE, DCONJG, MAX
83 * ..
84 * .. Executable Statements ..
85 *
86 INFO = 0
87 UPPER = LSAME( UPLO, 'U' )
88 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89 INFO = -1
90 ELSE IF( N.LT.0 ) THEN
91 INFO = -2
92 ELSE IF( NRHS.LT.0 ) THEN
93 INFO = -3
94 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
95 INFO = -5
96 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
97 INFO = -8
98 END IF
99 IF( INFO.NE.0 ) THEN
100 CALL XERBLA( 'ZHETRS', -INFO )
101 RETURN
102 END IF
103 *
104 * Quick return if possible
105 *
106 IF( N.EQ.0 .OR. NRHS.EQ.0 )
107 $ RETURN
108 *
109 IF( UPPER ) THEN
110 *
111 * Solve A*X = B, where A = U*D*U**H.
112 *
113 * First solve U*D*X = B, overwriting B with X.
114 *
115 * K is the main loop index, decreasing from N to 1 in steps of
116 * 1 or 2, depending on the size of the diagonal blocks.
117 *
118 K = N
119 10 CONTINUE
120 *
121 * If K < 1, exit from loop.
122 *
123 IF( K.LT.1 )
124 $ GO TO 30
125 *
126 IF( IPIV( K ).GT.0 ) THEN
127 *
128 * 1 x 1 diagonal block
129 *
130 * Interchange rows K and IPIV(K).
131 *
132 KP = IPIV( K )
133 IF( KP.NE.K )
134 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
135 *
136 * Multiply by inv(U(K)), where U(K) is the transformation
137 * stored in column K of A.
138 *
139 CALL ZGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
140 $ B( 1, 1 ), LDB )
141 *
142 * Multiply by the inverse of the diagonal block.
143 *
144 S = DBLE( ONE ) / DBLE( A( K, K ) )
145 CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
146 K = K - 1
147 ELSE
148 *
149 * 2 x 2 diagonal block
150 *
151 * Interchange rows K-1 and -IPIV(K).
152 *
153 KP = -IPIV( K )
154 IF( KP.NE.K-1 )
155 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
156 *
157 * Multiply by inv(U(K)), where U(K) is the transformation
158 * stored in columns K-1 and K of A.
159 *
160 CALL ZGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
161 $ B( 1, 1 ), LDB )
162 CALL ZGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
163 $ LDB, B( 1, 1 ), LDB )
164 *
165 * Multiply by the inverse of the diagonal block.
166 *
167 AKM1K = A( K-1, K )
168 AKM1 = A( K-1, K-1 ) / AKM1K
169 AK = A( K, K ) / DCONJG( AKM1K )
170 DENOM = AKM1*AK - ONE
171 DO 20 J = 1, NRHS
172 BKM1 = B( K-1, J ) / AKM1K
173 BK = B( K, J ) / DCONJG( AKM1K )
174 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
175 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
176 20 CONTINUE
177 K = K - 2
178 END IF
179 *
180 GO TO 10
181 30 CONTINUE
182 *
183 * Next solve U**H *X = B, overwriting B with X.
184 *
185 * K is the main loop index, increasing from 1 to N in steps of
186 * 1 or 2, depending on the size of the diagonal blocks.
187 *
188 K = 1
189 40 CONTINUE
190 *
191 * If K > N, exit from loop.
192 *
193 IF( K.GT.N )
194 $ GO TO 50
195 *
196 IF( IPIV( K ).GT.0 ) THEN
197 *
198 * 1 x 1 diagonal block
199 *
200 * Multiply by inv(U**H(K)), where U(K) is the transformation
201 * stored in column K of A.
202 *
203 IF( K.GT.1 ) THEN
204 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
205 CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
206 $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
207 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
208 END IF
209 *
210 * Interchange rows K and IPIV(K).
211 *
212 KP = IPIV( K )
213 IF( KP.NE.K )
214 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
215 K = K + 1
216 ELSE
217 *
218 * 2 x 2 diagonal block
219 *
220 * Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
221 * stored in columns K and K+1 of A.
222 *
223 IF( K.GT.1 ) THEN
224 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
225 CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
226 $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
227 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
228 *
229 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
230 CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
231 $ LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
232 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
233 END IF
234 *
235 * Interchange rows K and -IPIV(K).
236 *
237 KP = -IPIV( K )
238 IF( KP.NE.K )
239 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
240 K = K + 2
241 END IF
242 *
243 GO TO 40
244 50 CONTINUE
245 *
246 ELSE
247 *
248 * Solve A*X = B, where A = L*D*L**H.
249 *
250 * First solve L*D*X = B, overwriting B with X.
251 *
252 * K is the main loop index, increasing from 1 to N in steps of
253 * 1 or 2, depending on the size of the diagonal blocks.
254 *
255 K = 1
256 60 CONTINUE
257 *
258 * If K > N, exit from loop.
259 *
260 IF( K.GT.N )
261 $ GO TO 80
262 *
263 IF( IPIV( K ).GT.0 ) THEN
264 *
265 * 1 x 1 diagonal block
266 *
267 * Interchange rows K and IPIV(K).
268 *
269 KP = IPIV( K )
270 IF( KP.NE.K )
271 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
272 *
273 * Multiply by inv(L(K)), where L(K) is the transformation
274 * stored in column K of A.
275 *
276 IF( K.LT.N )
277 $ CALL ZGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
278 $ LDB, B( K+1, 1 ), LDB )
279 *
280 * Multiply by the inverse of the diagonal block.
281 *
282 S = DBLE( ONE ) / DBLE( A( K, K ) )
283 CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
284 K = K + 1
285 ELSE
286 *
287 * 2 x 2 diagonal block
288 *
289 * Interchange rows K+1 and -IPIV(K).
290 *
291 KP = -IPIV( K )
292 IF( KP.NE.K+1 )
293 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
294 *
295 * Multiply by inv(L(K)), where L(K) is the transformation
296 * stored in columns K and K+1 of A.
297 *
298 IF( K.LT.N-1 ) THEN
299 CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
300 $ LDB, B( K+2, 1 ), LDB )
301 CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
302 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
303 END IF
304 *
305 * Multiply by the inverse of the diagonal block.
306 *
307 AKM1K = A( K+1, K )
308 AKM1 = A( K, K ) / DCONJG( AKM1K )
309 AK = A( K+1, K+1 ) / AKM1K
310 DENOM = AKM1*AK - ONE
311 DO 70 J = 1, NRHS
312 BKM1 = B( K, J ) / DCONJG( AKM1K )
313 BK = B( K+1, J ) / AKM1K
314 B( K, J ) = ( AK*BKM1-BK ) / DENOM
315 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
316 70 CONTINUE
317 K = K + 2
318 END IF
319 *
320 GO TO 60
321 80 CONTINUE
322 *
323 * Next solve L**H *X = B, overwriting B with X.
324 *
325 * K is the main loop index, decreasing from N to 1 in steps of
326 * 1 or 2, depending on the size of the diagonal blocks.
327 *
328 K = N
329 90 CONTINUE
330 *
331 * If K < 1, exit from loop.
332 *
333 IF( K.LT.1 )
334 $ GO TO 100
335 *
336 IF( IPIV( K ).GT.0 ) THEN
337 *
338 * 1 x 1 diagonal block
339 *
340 * Multiply by inv(L**H(K)), where L(K) is the transformation
341 * stored in column K of A.
342 *
343 IF( K.LT.N ) THEN
344 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
345 CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
346 $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
347 $ B( K, 1 ), LDB )
348 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
349 END IF
350 *
351 * Interchange rows K and IPIV(K).
352 *
353 KP = IPIV( K )
354 IF( KP.NE.K )
355 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
356 K = K - 1
357 ELSE
358 *
359 * 2 x 2 diagonal block
360 *
361 * Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
362 * stored in columns K-1 and K of A.
363 *
364 IF( K.LT.N ) THEN
365 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
366 CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
367 $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
368 $ B( K, 1 ), LDB )
369 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
370 *
371 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
372 CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
373 $ B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
374 $ B( K-1, 1 ), LDB )
375 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
376 END IF
377 *
378 * Interchange rows K and -IPIV(K).
379 *
380 KP = -IPIV( K )
381 IF( KP.NE.K )
382 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
383 K = K - 2
384 END IF
385 *
386 GO TO 90
387 100 CONTINUE
388 END IF
389 *
390 RETURN
391 *
392 * End of ZHETRS
393 *
394 END