1 SUBROUTINE CLAVHP( 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 A( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CLAVHP 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 CHPTRF.
24 * CHPTRF 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', CLAVHP multiplies either by U or U * D
34 * (or L or L * D).
35 * If TRANS = 'C' or 'c', CLAVHP 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 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 CSPTRF or CHPTRF.
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 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 ONE
110 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
111 * ..
112 * .. Local Scalars ..
113 LOGICAL NOUNIT
114 INTEGER J, K, KC, KCNEXT, KP
115 COMPLEX D11, D12, D21, D22, T1, T2
116 * ..
117 * .. External Functions ..
118 LOGICAL LSAME
119 EXTERNAL LSAME
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL CGEMV, CGERU, CLACGV, CSCAL, CSWAP, XERBLA
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, CONJG, 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( 'CLAVHP ', -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 CSCAL( 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 CGERU( 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 CSWAP( 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 = CONJG( 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 CGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
230 $ LDB, B( 1, 1 ), LDB )
231 CALL CGERU( 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 CSWAP( 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 CSCAL( 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 CGERU( 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 CSWAP( 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 = CONJG( 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 CGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
318 $ LDB, B( K+1, 1 ), LDB )
319 CALL CGERU( 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 CSWAP( 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 IF( K.LT.1 )
353 $ GO TO 90
354 KC = KC - K
355 *
356 * 1 x 1 pivot block.
357 *
358 IF( IPIV( K ).GT.0 ) THEN
359 IF( K.GT.1 ) THEN
360 *
361 * Interchange if P(K) != I.
362 *
363 KP = IPIV( K )
364 IF( KP.NE.K )
365 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
366 *
367 * Apply the transformation:
368 * y := y - B' * conjg(x)
369 * where x is a column of A and y is a row of B.
370 *
371 CALL CLACGV( NRHS, B( K, 1 ), LDB )
372 CALL CGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB,
373 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
374 CALL CLACGV( NRHS, B( K, 1 ), LDB )
375 END IF
376 IF( NOUNIT )
377 $ CALL CSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
378 K = K - 1
379 *
380 * 2 x 2 pivot block.
381 *
382 ELSE
383 KCNEXT = KC - ( K-1 )
384 IF( K.GT.2 ) THEN
385 *
386 * Interchange if P(K) != I.
387 *
388 KP = ABS( IPIV( K ) )
389 IF( KP.NE.K-1 )
390 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
391 $ LDB )
392 *
393 * Apply the transformations.
394 *
395 CALL CLACGV( NRHS, B( K, 1 ), LDB )
396 CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
397 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
398 CALL CLACGV( NRHS, B( K, 1 ), LDB )
399 *
400 CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
401 CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
402 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
403 CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
404 END IF
405 *
406 * Multiply by the diagonal block if non-unit.
407 *
408 IF( NOUNIT ) THEN
409 D11 = A( KC-1 )
410 D22 = A( KC+K-1 )
411 D12 = A( KC+K-2 )
412 D21 = CONJG( D12 )
413 DO 80 J = 1, NRHS
414 T1 = B( K-1, J )
415 T2 = B( K, J )
416 B( K-1, J ) = D11*T1 + D12*T2
417 B( K, J ) = D21*T1 + D22*T2
418 80 CONTINUE
419 END IF
420 KC = KCNEXT
421 K = K - 2
422 END IF
423 GO TO 70
424 90 CONTINUE
425 *
426 * Form B := L^H*B
427 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
428 * and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
429 *
430 ELSE
431 *
432 * Loop forward applying the L-transformations.
433 *
434 K = 1
435 KC = 1
436 100 CONTINUE
437 IF( K.GT.N )
438 $ GO TO 120
439 *
440 * 1 x 1 pivot block
441 *
442 IF( IPIV( K ).GT.0 ) THEN
443 IF( K.LT.N ) THEN
444 *
445 * Interchange if P(K) != I.
446 *
447 KP = IPIV( K )
448 IF( KP.NE.K )
449 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
450 *
451 * Apply the transformation
452 *
453 CALL CLACGV( NRHS, B( K, 1 ), LDB )
454 CALL CGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
455 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
456 CALL CLACGV( NRHS, B( K, 1 ), LDB )
457 END IF
458 IF( NOUNIT )
459 $ CALL CSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
460 KC = KC + N - K + 1
461 K = K + 1
462 *
463 * 2 x 2 pivot block.
464 *
465 ELSE
466 KCNEXT = KC + N - K + 1
467 IF( K.LT.N-1 ) THEN
468 *
469 * Interchange if P(K) != I.
470 *
471 KP = ABS( IPIV( K ) )
472 IF( KP.NE.K+1 )
473 $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
474 $ LDB )
475 *
476 * Apply the transformation
477 *
478 CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
479 CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE,
480 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
481 $ B( K+1, 1 ), LDB )
482 CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
483 *
484 CALL CLACGV( NRHS, B( K, 1 ), LDB )
485 CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE,
486 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
487 $ B( K, 1 ), LDB )
488 CALL CLACGV( NRHS, B( K, 1 ), LDB )
489 END IF
490 *
491 * Multiply by the diagonal block if non-unit.
492 *
493 IF( NOUNIT ) THEN
494 D11 = A( KC )
495 D22 = A( KCNEXT )
496 D21 = A( KC+1 )
497 D12 = CONJG( D21 )
498 DO 110 J = 1, NRHS
499 T1 = B( K, J )
500 T2 = B( K+1, J )
501 B( K, J ) = D11*T1 + D12*T2
502 B( K+1, J ) = D21*T1 + D22*T2
503 110 CONTINUE
504 END IF
505 KC = KCNEXT + ( N-K )
506 K = K + 2
507 END IF
508 GO TO 100
509 120 CONTINUE
510 END IF
511 *
512 END IF
513 RETURN
514 *
515 * End of CLAVHP
516 *
517 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 A( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CLAVHP 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 CHPTRF.
24 * CHPTRF 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', CLAVHP multiplies either by U or U * D
34 * (or L or L * D).
35 * If TRANS = 'C' or 'c', CLAVHP 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 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 CSPTRF or CHPTRF.
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 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 ONE
110 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
111 * ..
112 * .. Local Scalars ..
113 LOGICAL NOUNIT
114 INTEGER J, K, KC, KCNEXT, KP
115 COMPLEX D11, D12, D21, D22, T1, T2
116 * ..
117 * .. External Functions ..
118 LOGICAL LSAME
119 EXTERNAL LSAME
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL CGEMV, CGERU, CLACGV, CSCAL, CSWAP, XERBLA
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, CONJG, 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( 'CLAVHP ', -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 CSCAL( 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 CGERU( 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 CSWAP( 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 = CONJG( 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 CGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
230 $ LDB, B( 1, 1 ), LDB )
231 CALL CGERU( 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 CSWAP( 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 CSCAL( 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 CGERU( 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 CSWAP( 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 = CONJG( 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 CGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
318 $ LDB, B( K+1, 1 ), LDB )
319 CALL CGERU( 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 CSWAP( 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 IF( K.LT.1 )
353 $ GO TO 90
354 KC = KC - K
355 *
356 * 1 x 1 pivot block.
357 *
358 IF( IPIV( K ).GT.0 ) THEN
359 IF( K.GT.1 ) THEN
360 *
361 * Interchange if P(K) != I.
362 *
363 KP = IPIV( K )
364 IF( KP.NE.K )
365 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
366 *
367 * Apply the transformation:
368 * y := y - B' * conjg(x)
369 * where x is a column of A and y is a row of B.
370 *
371 CALL CLACGV( NRHS, B( K, 1 ), LDB )
372 CALL CGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB,
373 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
374 CALL CLACGV( NRHS, B( K, 1 ), LDB )
375 END IF
376 IF( NOUNIT )
377 $ CALL CSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
378 K = K - 1
379 *
380 * 2 x 2 pivot block.
381 *
382 ELSE
383 KCNEXT = KC - ( K-1 )
384 IF( K.GT.2 ) THEN
385 *
386 * Interchange if P(K) != I.
387 *
388 KP = ABS( IPIV( K ) )
389 IF( KP.NE.K-1 )
390 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
391 $ LDB )
392 *
393 * Apply the transformations.
394 *
395 CALL CLACGV( NRHS, B( K, 1 ), LDB )
396 CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
397 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
398 CALL CLACGV( NRHS, B( K, 1 ), LDB )
399 *
400 CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
401 CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
402 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
403 CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
404 END IF
405 *
406 * Multiply by the diagonal block if non-unit.
407 *
408 IF( NOUNIT ) THEN
409 D11 = A( KC-1 )
410 D22 = A( KC+K-1 )
411 D12 = A( KC+K-2 )
412 D21 = CONJG( D12 )
413 DO 80 J = 1, NRHS
414 T1 = B( K-1, J )
415 T2 = B( K, J )
416 B( K-1, J ) = D11*T1 + D12*T2
417 B( K, J ) = D21*T1 + D22*T2
418 80 CONTINUE
419 END IF
420 KC = KCNEXT
421 K = K - 2
422 END IF
423 GO TO 70
424 90 CONTINUE
425 *
426 * Form B := L^H*B
427 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
428 * and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
429 *
430 ELSE
431 *
432 * Loop forward applying the L-transformations.
433 *
434 K = 1
435 KC = 1
436 100 CONTINUE
437 IF( K.GT.N )
438 $ GO TO 120
439 *
440 * 1 x 1 pivot block
441 *
442 IF( IPIV( K ).GT.0 ) THEN
443 IF( K.LT.N ) THEN
444 *
445 * Interchange if P(K) != I.
446 *
447 KP = IPIV( K )
448 IF( KP.NE.K )
449 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
450 *
451 * Apply the transformation
452 *
453 CALL CLACGV( NRHS, B( K, 1 ), LDB )
454 CALL CGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
455 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
456 CALL CLACGV( NRHS, B( K, 1 ), LDB )
457 END IF
458 IF( NOUNIT )
459 $ CALL CSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
460 KC = KC + N - K + 1
461 K = K + 1
462 *
463 * 2 x 2 pivot block.
464 *
465 ELSE
466 KCNEXT = KC + N - K + 1
467 IF( K.LT.N-1 ) THEN
468 *
469 * Interchange if P(K) != I.
470 *
471 KP = ABS( IPIV( K ) )
472 IF( KP.NE.K+1 )
473 $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
474 $ LDB )
475 *
476 * Apply the transformation
477 *
478 CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
479 CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE,
480 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
481 $ B( K+1, 1 ), LDB )
482 CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
483 *
484 CALL CLACGV( NRHS, B( K, 1 ), LDB )
485 CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE,
486 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
487 $ B( K, 1 ), LDB )
488 CALL CLACGV( NRHS, B( K, 1 ), LDB )
489 END IF
490 *
491 * Multiply by the diagonal block if non-unit.
492 *
493 IF( NOUNIT ) THEN
494 D11 = A( KC )
495 D22 = A( KCNEXT )
496 D21 = A( KC+1 )
497 D12 = CONJG( D21 )
498 DO 110 J = 1, NRHS
499 T1 = B( K, J )
500 T2 = B( K+1, J )
501 B( K, J ) = D11*T1 + D12*T2
502 B( K+1, J ) = D21*T1 + D22*T2
503 110 CONTINUE
504 END IF
505 KC = KCNEXT + ( N-K )
506 K = K + 2
507 END IF
508 GO TO 100
509 120 CONTINUE
510 END IF
511 *
512 END IF
513 RETURN
514 *
515 * End of CLAVHP
516 *
517 END