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