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