1 SUBROUTINE ZLAVHP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
2 $ INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, TRANS, UPLO
10 INTEGER INFO, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAVHP performs one of the matrix-vector operations
21 * x := A*x or x := A^H*x,
22 * where x is an N element vector and A is one of the factors
23 * from the symmetric factorization computed by ZHPTRF.
24 * ZHPTRF produces a factorization of the form
25 * U * D * U^H or L * D * L^H,
26 * where U (or L) is a product of permutation and unit upper (lower)
27 * triangular matrices, U^H (or L^H) is the conjugate transpose of
28 * U (or L), and D is Hermitian and block diagonal with 1 x 1 and
29 * 2 x 2 diagonal blocks. The multipliers for the transformations
30 * and the upper or lower triangular parts of the diagonal blocks
31 * are stored columnwise in packed format in the linear array A.
32 *
33 * If TRANS = 'N' or 'n', ZLAVHP multiplies either by U or U * D
34 * (or L or L * D).
35 * If TRANS = 'C' or 'c', ZLAVHP multiplies either by U^H or D * U^H
36 * (or L^H or D * L^H ).
37 *
38 * Arguments
39 * ==========
40 *
41 * UPLO - CHARACTER*1
42 * On entry, UPLO specifies whether the triangular matrix
43 * stored in A is upper or lower triangular.
44 * UPLO = 'U' or 'u' The matrix is upper triangular.
45 * UPLO = 'L' or 'l' The matrix is lower triangular.
46 * Unchanged on exit.
47 *
48 * TRANS - CHARACTER*1
49 * On entry, TRANS specifies the operation to be performed as
50 * follows:
51 * TRANS = 'N' or 'n' x := A*x.
52 * TRANS = 'C' or 'c' x := A^H*x.
53 * Unchanged on exit.
54 *
55 * DIAG - CHARACTER*1
56 * On entry, DIAG specifies whether the diagonal blocks are
57 * assumed to be unit matrices, as follows:
58 * DIAG = 'U' or 'u' Diagonal blocks are unit matrices.
59 * DIAG = 'N' or 'n' Diagonal blocks are non-unit.
60 * Unchanged on exit.
61 *
62 * N - INTEGER
63 * On entry, N specifies the order of the matrix A.
64 * N must be at least zero.
65 * Unchanged on exit.
66 *
67 * NRHS - INTEGER
68 * On entry, NRHS specifies the number of right hand sides,
69 * i.e., the number of vectors x to be multiplied by A.
70 * NRHS must be at least zero.
71 * Unchanged on exit.
72 *
73 * A - COMPLEX*16 array, dimension( N*(N+1)/2 )
74 * On entry, A contains a block diagonal matrix and the
75 * multipliers of the transformations used to obtain it,
76 * stored as a packed triangular matrix.
77 * Unchanged on exit.
78 *
79 * IPIV - INTEGER array, dimension( N )
80 * On entry, IPIV contains the vector of pivot indices as
81 * determined by ZSPTRF or ZHPTRF.
82 * If IPIV( K ) = K, no interchange was done.
83 * If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
84 * changed with row IPIV( K ) and a 1 x 1 pivot block was used.
85 * If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
86 * with row | IPIV( K ) | and a 2 x 2 pivot block was used.
87 * If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
88 * with row | IPIV( K ) | and a 2 x 2 pivot block was used.
89 *
90 * B - COMPLEX*16 array, dimension( LDB, NRHS )
91 * On entry, B contains NRHS vectors of length N.
92 * On exit, B is overwritten with the product A * B.
93 *
94 * LDB - INTEGER
95 * On entry, LDB contains the leading dimension of B as
96 * declared in the calling program. LDB must be at least
97 * max( 1, N ).
98 * Unchanged on exit.
99 *
100 * INFO - INTEGER
101 * INFO is the error flag.
102 * On exit, a value of 0 indicates a successful exit.
103 * A negative value, say -K, indicates that the K-th argument
104 * has an illegal value.
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109 COMPLEX*16 ONE
110 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
111 * ..
112 * .. Local Scalars ..
113 LOGICAL NOUNIT
114 INTEGER J, K, KC, KCNEXT, KP
115 COMPLEX*16 D11, D12, D21, D22, T1, T2
116 * ..
117 * .. External Functions ..
118 LOGICAL LSAME
119 EXTERNAL LSAME
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL XERBLA, ZGEMV, ZGERU, ZLACGV, ZSCAL, ZSWAP
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, DCONJG, MAX
126 * ..
127 * .. Executable Statements ..
128 *
129 * Test the input parameters.
130 *
131 INFO = 0
132 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
133 INFO = -1
134 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
135 $ THEN
136 INFO = -2
137 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
138 $ THEN
139 INFO = -3
140 ELSE IF( N.LT.0 ) THEN
141 INFO = -4
142 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
143 INFO = -8
144 END IF
145 IF( INFO.NE.0 ) THEN
146 CALL XERBLA( 'ZLAVHP ', -INFO )
147 RETURN
148 END IF
149 *
150 * Quick return if possible.
151 *
152 IF( N.EQ.0 )
153 $ RETURN
154 *
155 NOUNIT = LSAME( DIAG, 'N' )
156 *------------------------------------------
157 *
158 * Compute B := A * B (No transpose)
159 *
160 *------------------------------------------
161 IF( LSAME( TRANS, 'N' ) ) THEN
162 *
163 * Compute B := U*B
164 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
165 *
166 IF( LSAME( UPLO, 'U' ) ) THEN
167 *
168 * Loop forward applying the transformations.
169 *
170 K = 1
171 KC = 1
172 10 CONTINUE
173 IF( K.GT.N )
174 $ GO TO 30
175 *
176 * 1 x 1 pivot block
177 *
178 IF( IPIV( K ).GT.0 ) THEN
179 *
180 * Multiply by the diagonal element if forming U * D.
181 *
182 IF( NOUNIT )
183 $ CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
184 *
185 * Multiply by P(K) * inv(U(K)) if K > 1.
186 *
187 IF( K.GT.1 ) THEN
188 *
189 * Apply the transformation.
190 *
191 CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
192 $ LDB, B( 1, 1 ), LDB )
193 *
194 * Interchange if P(K) != I.
195 *
196 KP = IPIV( K )
197 IF( KP.NE.K )
198 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
199 END IF
200 KC = KC + K
201 K = K + 1
202 ELSE
203 *
204 * 2 x 2 pivot block
205 *
206 KCNEXT = KC + K
207 *
208 * Multiply by the diagonal block if forming U * D.
209 *
210 IF( NOUNIT ) THEN
211 D11 = A( KCNEXT-1 )
212 D22 = A( KCNEXT+K )
213 D12 = A( KCNEXT+K-1 )
214 D21 = DCONJG( D12 )
215 DO 20 J = 1, NRHS
216 T1 = B( K, J )
217 T2 = B( K+1, J )
218 B( K, J ) = D11*T1 + D12*T2
219 B( K+1, J ) = D21*T1 + D22*T2
220 20 CONTINUE
221 END IF
222 *
223 * Multiply by P(K) * inv(U(K)) if K > 1.
224 *
225 IF( K.GT.1 ) THEN
226 *
227 * Apply the transformations.
228 *
229 CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
230 $ LDB, B( 1, 1 ), LDB )
231 CALL ZGERU( K-1, NRHS, ONE, A( KCNEXT ), 1,
232 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
233 *
234 * Interchange if P(K) != I.
235 *
236 KP = ABS( IPIV( K ) )
237 IF( KP.NE.K )
238 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
239 END IF
240 KC = KCNEXT + K + 1
241 K = K + 2
242 END IF
243 GO TO 10
244 30 CONTINUE
245 *
246 * Compute B := L*B
247 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
248 *
249 ELSE
250 *
251 * Loop backward applying the transformations to B.
252 *
253 K = N
254 KC = N*( N+1 ) / 2 + 1
255 40 CONTINUE
256 IF( K.LT.1 )
257 $ GO TO 60
258 KC = KC - ( N-K+1 )
259 *
260 * Test the pivot index. If greater than zero, a 1 x 1
261 * pivot was used, otherwise a 2 x 2 pivot was used.
262 *
263 IF( IPIV( K ).GT.0 ) THEN
264 *
265 * 1 x 1 pivot block:
266 *
267 * Multiply by the diagonal element if forming L * D.
268 *
269 IF( NOUNIT )
270 $ CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
271 *
272 * Multiply by P(K) * inv(L(K)) if K < N.
273 *
274 IF( K.NE.N ) THEN
275 KP = IPIV( K )
276 *
277 * Apply the transformation.
278 *
279 CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
280 $ LDB, B( K+1, 1 ), LDB )
281 *
282 * Interchange if a permutation was applied at the
283 * K-th step of the factorization.
284 *
285 IF( KP.NE.K )
286 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
287 END IF
288 K = K - 1
289 *
290 ELSE
291 *
292 * 2 x 2 pivot block:
293 *
294 KCNEXT = KC - ( N-K+2 )
295 *
296 * Multiply by the diagonal block if forming L * D.
297 *
298 IF( NOUNIT ) THEN
299 D11 = A( KCNEXT )
300 D22 = A( KC )
301 D21 = A( KCNEXT+1 )
302 D12 = DCONJG( D21 )
303 DO 50 J = 1, NRHS
304 T1 = B( K-1, J )
305 T2 = B( K, J )
306 B( K-1, J ) = D11*T1 + D12*T2
307 B( K, J ) = D21*T1 + D22*T2
308 50 CONTINUE
309 END IF
310 *
311 * Multiply by P(K) * inv(L(K)) if K < N.
312 *
313 IF( K.NE.N ) THEN
314 *
315 * Apply the transformation.
316 *
317 CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
318 $ LDB, B( K+1, 1 ), LDB )
319 CALL ZGERU( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
320 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
321 *
322 * Interchange if a permutation was applied at the
323 * K-th step of the factorization.
324 *
325 KP = ABS( IPIV( K ) )
326 IF( KP.NE.K )
327 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
328 END IF
329 KC = KCNEXT
330 K = K - 2
331 END IF
332 GO TO 40
333 60 CONTINUE
334 END IF
335 *-------------------------------------------------
336 *
337 * Compute B := A^H * B (conjugate transpose)
338 *
339 *-------------------------------------------------
340 ELSE
341 *
342 * Form B := U^H*B
343 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
344 * and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
345 *
346 IF( LSAME( UPLO, 'U' ) ) THEN
347 *
348 * Loop backward applying the transformations.
349 *
350 K = N
351 KC = N*( N+1 ) / 2 + 1
352 70 CONTINUE
353 IF( K.LT.1 )
354 $ GO TO 90
355 KC = KC - K
356 *
357 * 1 x 1 pivot block.
358 *
359 IF( IPIV( K ).GT.0 ) THEN
360 IF( K.GT.1 ) THEN
361 *
362 * Interchange if P(K) != I.
363 *
364 KP = IPIV( K )
365 IF( KP.NE.K )
366 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
367 *
368 * Apply the transformation:
369 * y := y - B' * conjg(x)
370 * where x is a column of A and y is a row of B.
371 *
372 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
373 CALL ZGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB,
374 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
375 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
376 END IF
377 IF( NOUNIT )
378 $ CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
379 K = K - 1
380 *
381 * 2 x 2 pivot block.
382 *
383 ELSE
384 KCNEXT = KC - ( K-1 )
385 IF( K.GT.2 ) THEN
386 *
387 * Interchange if P(K) != I.
388 *
389 KP = ABS( IPIV( K ) )
390 IF( KP.NE.K-1 )
391 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
392 $ LDB )
393 *
394 * Apply the transformations.
395 *
396 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
397 CALL ZGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
398 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
399 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
400 *
401 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
402 CALL ZGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
403 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
404 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
405 END IF
406 *
407 * Multiply by the diagonal block if non-unit.
408 *
409 IF( NOUNIT ) THEN
410 D11 = A( KC-1 )
411 D22 = A( KC+K-1 )
412 D12 = A( KC+K-2 )
413 D21 = DCONJG( D12 )
414 DO 80 J = 1, NRHS
415 T1 = B( K-1, J )
416 T2 = B( K, J )
417 B( K-1, J ) = D11*T1 + D12*T2
418 B( K, J ) = D21*T1 + D22*T2
419 80 CONTINUE
420 END IF
421 KC = KCNEXT
422 K = K - 2
423 END IF
424 GO TO 70
425 90 CONTINUE
426 *
427 * Form B := L^H*B
428 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
429 * and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
430 *
431 ELSE
432 *
433 * Loop forward applying the L-transformations.
434 *
435 K = 1
436 KC = 1
437 100 CONTINUE
438 IF( K.GT.N )
439 $ GO TO 120
440 *
441 * 1 x 1 pivot block
442 *
443 IF( IPIV( K ).GT.0 ) THEN
444 IF( K.LT.N ) THEN
445 *
446 * Interchange if P(K) != I.
447 *
448 KP = IPIV( K )
449 IF( KP.NE.K )
450 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
451 *
452 * Apply the transformation
453 *
454 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
455 CALL ZGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
456 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
457 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
458 END IF
459 IF( NOUNIT )
460 $ CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
461 KC = KC + N - K + 1
462 K = K + 1
463 *
464 * 2 x 2 pivot block.
465 *
466 ELSE
467 KCNEXT = KC + N - K + 1
468 IF( K.LT.N-1 ) THEN
469 *
470 * Interchange if P(K) != I.
471 *
472 KP = ABS( IPIV( K ) )
473 IF( KP.NE.K+1 )
474 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
475 $ LDB )
476 *
477 * Apply the transformation
478 *
479 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
480 CALL ZGEMV( 'Conjugate', N-K-1, NRHS, ONE,
481 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
482 $ B( K+1, 1 ), LDB )
483 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
484 *
485 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
486 CALL ZGEMV( 'Conjugate', N-K-1, NRHS, ONE,
487 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
488 $ B( K, 1 ), LDB )
489 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
490 END IF
491 *
492 * Multiply by the diagonal block if non-unit.
493 *
494 IF( NOUNIT ) THEN
495 D11 = A( KC )
496 D22 = A( KCNEXT )
497 D21 = A( KC+1 )
498 D12 = DCONJG( D21 )
499 DO 110 J = 1, NRHS
500 T1 = B( K, J )
501 T2 = B( K+1, J )
502 B( K, J ) = D11*T1 + D12*T2
503 B( K+1, J ) = D21*T1 + D22*T2
504 110 CONTINUE
505 END IF
506 KC = KCNEXT + ( N-K )
507 K = K + 2
508 END IF
509 GO TO 100
510 120 CONTINUE
511 END IF
512 *
513 END IF
514 RETURN
515 *
516 * End of ZLAVHP
517 *
518 END
2 $ INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, TRANS, UPLO
10 INTEGER INFO, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAVHP performs one of the matrix-vector operations
21 * x := A*x or x := A^H*x,
22 * where x is an N element vector and A is one of the factors
23 * from the symmetric factorization computed by ZHPTRF.
24 * ZHPTRF produces a factorization of the form
25 * U * D * U^H or L * D * L^H,
26 * where U (or L) is a product of permutation and unit upper (lower)
27 * triangular matrices, U^H (or L^H) is the conjugate transpose of
28 * U (or L), and D is Hermitian and block diagonal with 1 x 1 and
29 * 2 x 2 diagonal blocks. The multipliers for the transformations
30 * and the upper or lower triangular parts of the diagonal blocks
31 * are stored columnwise in packed format in the linear array A.
32 *
33 * If TRANS = 'N' or 'n', ZLAVHP multiplies either by U or U * D
34 * (or L or L * D).
35 * If TRANS = 'C' or 'c', ZLAVHP multiplies either by U^H or D * U^H
36 * (or L^H or D * L^H ).
37 *
38 * Arguments
39 * ==========
40 *
41 * UPLO - CHARACTER*1
42 * On entry, UPLO specifies whether the triangular matrix
43 * stored in A is upper or lower triangular.
44 * UPLO = 'U' or 'u' The matrix is upper triangular.
45 * UPLO = 'L' or 'l' The matrix is lower triangular.
46 * Unchanged on exit.
47 *
48 * TRANS - CHARACTER*1
49 * On entry, TRANS specifies the operation to be performed as
50 * follows:
51 * TRANS = 'N' or 'n' x := A*x.
52 * TRANS = 'C' or 'c' x := A^H*x.
53 * Unchanged on exit.
54 *
55 * DIAG - CHARACTER*1
56 * On entry, DIAG specifies whether the diagonal blocks are
57 * assumed to be unit matrices, as follows:
58 * DIAG = 'U' or 'u' Diagonal blocks are unit matrices.
59 * DIAG = 'N' or 'n' Diagonal blocks are non-unit.
60 * Unchanged on exit.
61 *
62 * N - INTEGER
63 * On entry, N specifies the order of the matrix A.
64 * N must be at least zero.
65 * Unchanged on exit.
66 *
67 * NRHS - INTEGER
68 * On entry, NRHS specifies the number of right hand sides,
69 * i.e., the number of vectors x to be multiplied by A.
70 * NRHS must be at least zero.
71 * Unchanged on exit.
72 *
73 * A - COMPLEX*16 array, dimension( N*(N+1)/2 )
74 * On entry, A contains a block diagonal matrix and the
75 * multipliers of the transformations used to obtain it,
76 * stored as a packed triangular matrix.
77 * Unchanged on exit.
78 *
79 * IPIV - INTEGER array, dimension( N )
80 * On entry, IPIV contains the vector of pivot indices as
81 * determined by ZSPTRF or ZHPTRF.
82 * If IPIV( K ) = K, no interchange was done.
83 * If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
84 * changed with row IPIV( K ) and a 1 x 1 pivot block was used.
85 * If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
86 * with row | IPIV( K ) | and a 2 x 2 pivot block was used.
87 * If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
88 * with row | IPIV( K ) | and a 2 x 2 pivot block was used.
89 *
90 * B - COMPLEX*16 array, dimension( LDB, NRHS )
91 * On entry, B contains NRHS vectors of length N.
92 * On exit, B is overwritten with the product A * B.
93 *
94 * LDB - INTEGER
95 * On entry, LDB contains the leading dimension of B as
96 * declared in the calling program. LDB must be at least
97 * max( 1, N ).
98 * Unchanged on exit.
99 *
100 * INFO - INTEGER
101 * INFO is the error flag.
102 * On exit, a value of 0 indicates a successful exit.
103 * A negative value, say -K, indicates that the K-th argument
104 * has an illegal value.
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109 COMPLEX*16 ONE
110 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
111 * ..
112 * .. Local Scalars ..
113 LOGICAL NOUNIT
114 INTEGER J, K, KC, KCNEXT, KP
115 COMPLEX*16 D11, D12, D21, D22, T1, T2
116 * ..
117 * .. External Functions ..
118 LOGICAL LSAME
119 EXTERNAL LSAME
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL XERBLA, ZGEMV, ZGERU, ZLACGV, ZSCAL, ZSWAP
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, DCONJG, MAX
126 * ..
127 * .. Executable Statements ..
128 *
129 * Test the input parameters.
130 *
131 INFO = 0
132 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
133 INFO = -1
134 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
135 $ THEN
136 INFO = -2
137 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
138 $ THEN
139 INFO = -3
140 ELSE IF( N.LT.0 ) THEN
141 INFO = -4
142 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
143 INFO = -8
144 END IF
145 IF( INFO.NE.0 ) THEN
146 CALL XERBLA( 'ZLAVHP ', -INFO )
147 RETURN
148 END IF
149 *
150 * Quick return if possible.
151 *
152 IF( N.EQ.0 )
153 $ RETURN
154 *
155 NOUNIT = LSAME( DIAG, 'N' )
156 *------------------------------------------
157 *
158 * Compute B := A * B (No transpose)
159 *
160 *------------------------------------------
161 IF( LSAME( TRANS, 'N' ) ) THEN
162 *
163 * Compute B := U*B
164 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
165 *
166 IF( LSAME( UPLO, 'U' ) ) THEN
167 *
168 * Loop forward applying the transformations.
169 *
170 K = 1
171 KC = 1
172 10 CONTINUE
173 IF( K.GT.N )
174 $ GO TO 30
175 *
176 * 1 x 1 pivot block
177 *
178 IF( IPIV( K ).GT.0 ) THEN
179 *
180 * Multiply by the diagonal element if forming U * D.
181 *
182 IF( NOUNIT )
183 $ CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
184 *
185 * Multiply by P(K) * inv(U(K)) if K > 1.
186 *
187 IF( K.GT.1 ) THEN
188 *
189 * Apply the transformation.
190 *
191 CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
192 $ LDB, B( 1, 1 ), LDB )
193 *
194 * Interchange if P(K) != I.
195 *
196 KP = IPIV( K )
197 IF( KP.NE.K )
198 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
199 END IF
200 KC = KC + K
201 K = K + 1
202 ELSE
203 *
204 * 2 x 2 pivot block
205 *
206 KCNEXT = KC + K
207 *
208 * Multiply by the diagonal block if forming U * D.
209 *
210 IF( NOUNIT ) THEN
211 D11 = A( KCNEXT-1 )
212 D22 = A( KCNEXT+K )
213 D12 = A( KCNEXT+K-1 )
214 D21 = DCONJG( D12 )
215 DO 20 J = 1, NRHS
216 T1 = B( K, J )
217 T2 = B( K+1, J )
218 B( K, J ) = D11*T1 + D12*T2
219 B( K+1, J ) = D21*T1 + D22*T2
220 20 CONTINUE
221 END IF
222 *
223 * Multiply by P(K) * inv(U(K)) if K > 1.
224 *
225 IF( K.GT.1 ) THEN
226 *
227 * Apply the transformations.
228 *
229 CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
230 $ LDB, B( 1, 1 ), LDB )
231 CALL ZGERU( K-1, NRHS, ONE, A( KCNEXT ), 1,
232 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
233 *
234 * Interchange if P(K) != I.
235 *
236 KP = ABS( IPIV( K ) )
237 IF( KP.NE.K )
238 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
239 END IF
240 KC = KCNEXT + K + 1
241 K = K + 2
242 END IF
243 GO TO 10
244 30 CONTINUE
245 *
246 * Compute B := L*B
247 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
248 *
249 ELSE
250 *
251 * Loop backward applying the transformations to B.
252 *
253 K = N
254 KC = N*( N+1 ) / 2 + 1
255 40 CONTINUE
256 IF( K.LT.1 )
257 $ GO TO 60
258 KC = KC - ( N-K+1 )
259 *
260 * Test the pivot index. If greater than zero, a 1 x 1
261 * pivot was used, otherwise a 2 x 2 pivot was used.
262 *
263 IF( IPIV( K ).GT.0 ) THEN
264 *
265 * 1 x 1 pivot block:
266 *
267 * Multiply by the diagonal element if forming L * D.
268 *
269 IF( NOUNIT )
270 $ CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
271 *
272 * Multiply by P(K) * inv(L(K)) if K < N.
273 *
274 IF( K.NE.N ) THEN
275 KP = IPIV( K )
276 *
277 * Apply the transformation.
278 *
279 CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
280 $ LDB, B( K+1, 1 ), LDB )
281 *
282 * Interchange if a permutation was applied at the
283 * K-th step of the factorization.
284 *
285 IF( KP.NE.K )
286 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
287 END IF
288 K = K - 1
289 *
290 ELSE
291 *
292 * 2 x 2 pivot block:
293 *
294 KCNEXT = KC - ( N-K+2 )
295 *
296 * Multiply by the diagonal block if forming L * D.
297 *
298 IF( NOUNIT ) THEN
299 D11 = A( KCNEXT )
300 D22 = A( KC )
301 D21 = A( KCNEXT+1 )
302 D12 = DCONJG( D21 )
303 DO 50 J = 1, NRHS
304 T1 = B( K-1, J )
305 T2 = B( K, J )
306 B( K-1, J ) = D11*T1 + D12*T2
307 B( K, J ) = D21*T1 + D22*T2
308 50 CONTINUE
309 END IF
310 *
311 * Multiply by P(K) * inv(L(K)) if K < N.
312 *
313 IF( K.NE.N ) THEN
314 *
315 * Apply the transformation.
316 *
317 CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
318 $ LDB, B( K+1, 1 ), LDB )
319 CALL ZGERU( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
320 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
321 *
322 * Interchange if a permutation was applied at the
323 * K-th step of the factorization.
324 *
325 KP = ABS( IPIV( K ) )
326 IF( KP.NE.K )
327 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
328 END IF
329 KC = KCNEXT
330 K = K - 2
331 END IF
332 GO TO 40
333 60 CONTINUE
334 END IF
335 *-------------------------------------------------
336 *
337 * Compute B := A^H * B (conjugate transpose)
338 *
339 *-------------------------------------------------
340 ELSE
341 *
342 * Form B := U^H*B
343 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
344 * and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
345 *
346 IF( LSAME( UPLO, 'U' ) ) THEN
347 *
348 * Loop backward applying the transformations.
349 *
350 K = N
351 KC = N*( N+1 ) / 2 + 1
352 70 CONTINUE
353 IF( K.LT.1 )
354 $ GO TO 90
355 KC = KC - K
356 *
357 * 1 x 1 pivot block.
358 *
359 IF( IPIV( K ).GT.0 ) THEN
360 IF( K.GT.1 ) THEN
361 *
362 * Interchange if P(K) != I.
363 *
364 KP = IPIV( K )
365 IF( KP.NE.K )
366 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
367 *
368 * Apply the transformation:
369 * y := y - B' * conjg(x)
370 * where x is a column of A and y is a row of B.
371 *
372 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
373 CALL ZGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB,
374 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
375 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
376 END IF
377 IF( NOUNIT )
378 $ CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
379 K = K - 1
380 *
381 * 2 x 2 pivot block.
382 *
383 ELSE
384 KCNEXT = KC - ( K-1 )
385 IF( K.GT.2 ) THEN
386 *
387 * Interchange if P(K) != I.
388 *
389 KP = ABS( IPIV( K ) )
390 IF( KP.NE.K-1 )
391 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
392 $ LDB )
393 *
394 * Apply the transformations.
395 *
396 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
397 CALL ZGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
398 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
399 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
400 *
401 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
402 CALL ZGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
403 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
404 CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
405 END IF
406 *
407 * Multiply by the diagonal block if non-unit.
408 *
409 IF( NOUNIT ) THEN
410 D11 = A( KC-1 )
411 D22 = A( KC+K-1 )
412 D12 = A( KC+K-2 )
413 D21 = DCONJG( D12 )
414 DO 80 J = 1, NRHS
415 T1 = B( K-1, J )
416 T2 = B( K, J )
417 B( K-1, J ) = D11*T1 + D12*T2
418 B( K, J ) = D21*T1 + D22*T2
419 80 CONTINUE
420 END IF
421 KC = KCNEXT
422 K = K - 2
423 END IF
424 GO TO 70
425 90 CONTINUE
426 *
427 * Form B := L^H*B
428 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
429 * and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
430 *
431 ELSE
432 *
433 * Loop forward applying the L-transformations.
434 *
435 K = 1
436 KC = 1
437 100 CONTINUE
438 IF( K.GT.N )
439 $ GO TO 120
440 *
441 * 1 x 1 pivot block
442 *
443 IF( IPIV( K ).GT.0 ) THEN
444 IF( K.LT.N ) THEN
445 *
446 * Interchange if P(K) != I.
447 *
448 KP = IPIV( K )
449 IF( KP.NE.K )
450 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
451 *
452 * Apply the transformation
453 *
454 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
455 CALL ZGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
456 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
457 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
458 END IF
459 IF( NOUNIT )
460 $ CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
461 KC = KC + N - K + 1
462 K = K + 1
463 *
464 * 2 x 2 pivot block.
465 *
466 ELSE
467 KCNEXT = KC + N - K + 1
468 IF( K.LT.N-1 ) THEN
469 *
470 * Interchange if P(K) != I.
471 *
472 KP = ABS( IPIV( K ) )
473 IF( KP.NE.K+1 )
474 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
475 $ LDB )
476 *
477 * Apply the transformation
478 *
479 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
480 CALL ZGEMV( 'Conjugate', N-K-1, NRHS, ONE,
481 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
482 $ B( K+1, 1 ), LDB )
483 CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
484 *
485 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
486 CALL ZGEMV( 'Conjugate', N-K-1, NRHS, ONE,
487 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
488 $ B( K, 1 ), LDB )
489 CALL ZLACGV( NRHS, B( K, 1 ), LDB )
490 END IF
491 *
492 * Multiply by the diagonal block if non-unit.
493 *
494 IF( NOUNIT ) THEN
495 D11 = A( KC )
496 D22 = A( KCNEXT )
497 D21 = A( KC+1 )
498 D12 = DCONJG( D21 )
499 DO 110 J = 1, NRHS
500 T1 = B( K, J )
501 T2 = B( K+1, J )
502 B( K, J ) = D11*T1 + D12*T2
503 B( K+1, J ) = D21*T1 + D22*T2
504 110 CONTINUE
505 END IF
506 KC = KCNEXT + ( N-K )
507 K = K + 2
508 END IF
509 GO TO 100
510 120 CONTINUE
511 END IF
512 *
513 END IF
514 RETURN
515 *
516 * End of ZLAVHP
517 *
518 END