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