LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlavhp()

subroutine zlavhp ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  nrhs,
complex*16, dimension( * )  a,
integer, dimension( * )  ipiv,
complex*16, dimension( ldb, * )  b,
integer  ldb,
integer  info 
)

ZLAVHP

Purpose:
    ZLAVHP  performs one of the matrix-vector operations
       x := A*x  or  x := A^H*x,
    where x is an N element vector and  A is one of the factors
    from the symmetric factorization computed by ZHPTRF.
    ZHPTRF produces a factorization of the form
         U * D * U^H     or     L * D * L^H,
    where U (or L) is a product of permutation and unit upper (lower)
    triangular matrices, U^H (or L^H) is the conjugate transpose of
    U (or L), and D is Hermitian and block diagonal with 1 x 1 and
    2 x 2 diagonal blocks.  The multipliers for the transformations
    and the upper or lower triangular parts of the diagonal blocks
    are stored columnwise in packed format in the linear array A.

    If TRANS = 'N' or 'n', ZLAVHP multiplies either by U or U * D
    (or L or L * D).
    If TRANS = 'C' or 'c', ZLAVHP multiplies either by U^H or D * U^H
    (or L^H or D * L^H ).
  UPLO   - CHARACTER*1
           On entry, UPLO specifies whether the triangular matrix
           stored in A is upper or lower triangular.
              UPLO = 'U' or 'u'   The matrix is upper triangular.
              UPLO = 'L' or 'l'   The matrix is lower triangular.
           Unchanged on exit.

  TRANS  - CHARACTER*1
           On entry, TRANS specifies the operation to be performed as
           follows:
              TRANS = 'N' or 'n'   x := A*x.
              TRANS = 'C' or 'c'   x := A^H*x.
           Unchanged on exit.

  DIAG   - CHARACTER*1
           On entry, DIAG specifies whether the diagonal blocks are
           assumed to be unit matrices, as follows:
              DIAG = 'U' or 'u'   Diagonal blocks are unit matrices.
              DIAG = 'N' or 'n'   Diagonal blocks are non-unit.
           Unchanged on exit.

  N      - INTEGER
           On entry, N specifies the order of the matrix A.
           N must be at least zero.
           Unchanged on exit.

  NRHS   - INTEGER
           On entry, NRHS specifies the number of right hand sides,
           i.e., the number of vectors x to be multiplied by A.
           NRHS must be at least zero.
           Unchanged on exit.

  A      - COMPLEX*16 array, dimension( N*(N+1)/2 )
           On entry, A contains a block diagonal matrix and the
           multipliers of the transformations used to obtain it,
           stored as a packed triangular matrix.
           Unchanged on exit.

  IPIV   - INTEGER array, dimension( N )
           On entry, IPIV contains the vector of pivot indices as
           determined by ZSPTRF or ZHPTRF.
           If IPIV( K ) = K, no interchange was done.
           If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
           changed with row IPIV( K ) and a 1 x 1 pivot block was used.
           If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
           If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
           with row | IPIV( K ) | and a 2 x 2 pivot block was used.

  B      - COMPLEX*16 array, dimension( LDB, NRHS )
           On entry, B contains NRHS vectors of length N.
           On exit, B is overwritten with the product A * B.

  LDB    - INTEGER
           On entry, LDB contains the leading dimension of B as
           declared in the calling program.  LDB must be at least
           max( 1, N ).
           Unchanged on exit.

  INFO   - INTEGER
           INFO is the error flag.
           On exit, a value of 0 indicates a successful exit.
           A negative value, say -K, indicates that the K-th argument
           has an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 129 of file zlavhp.f.

131*
132* -- LAPACK test routine --
133* -- LAPACK is a software package provided by Univ. of Tennessee, --
134* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135*
136* .. Scalar Arguments ..
137 CHARACTER DIAG, TRANS, UPLO
138 INTEGER INFO, LDB, N, NRHS
139* ..
140* .. Array Arguments ..
141 INTEGER IPIV( * )
142 COMPLEX*16 A( * ), B( LDB, * )
143* ..
144*
145* =====================================================================
146*
147* .. Parameters ..
148 COMPLEX*16 ONE
149 parameter( one = ( 1.0d+0, 0.0d+0 ) )
150* ..
151* .. Local Scalars ..
152 LOGICAL NOUNIT
153 INTEGER J, K, KC, KCNEXT, KP
154 COMPLEX*16 D11, D12, D21, D22, T1, T2
155* ..
156* .. External Functions ..
157 LOGICAL LSAME
158 EXTERNAL lsame
159* ..
160* .. External Subroutines ..
161 EXTERNAL xerbla, zgemv, zgeru, zlacgv, zscal, zswap
162* ..
163* .. Intrinsic Functions ..
164 INTRINSIC abs, dconjg, max
165* ..
166* .. Executable Statements ..
167*
168* Test the input parameters.
169*
170 info = 0
171 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
172 info = -1
173 ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'C' ) )
174 $ THEN
175 info = -2
176 ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
177 $ THEN
178 info = -3
179 ELSE IF( n.LT.0 ) THEN
180 info = -4
181 ELSE IF( ldb.LT.max( 1, n ) ) THEN
182 info = -8
183 END IF
184 IF( info.NE.0 ) THEN
185 CALL xerbla( 'ZLAVHP ', -info )
186 RETURN
187 END IF
188*
189* Quick return if possible.
190*
191 IF( n.EQ.0 )
192 $ RETURN
193*
194 nounit = lsame( diag, 'N' )
195*------------------------------------------
196*
197* Compute B := A * B (No transpose)
198*
199*------------------------------------------
200 IF( lsame( trans, 'N' ) ) THEN
201*
202* Compute B := U*B
203* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
204*
205 IF( lsame( uplo, 'U' ) ) THEN
206*
207* Loop forward applying the transformations.
208*
209 k = 1
210 kc = 1
211 10 CONTINUE
212 IF( k.GT.n )
213 $ GO TO 30
214*
215* 1 x 1 pivot block
216*
217 IF( ipiv( k ).GT.0 ) THEN
218*
219* Multiply by the diagonal element if forming U * D.
220*
221 IF( nounit )
222 $ CALL zscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
223*
224* Multiply by P(K) * inv(U(K)) if K > 1.
225*
226 IF( k.GT.1 ) THEN
227*
228* Apply the transformation.
229*
230 CALL zgeru( k-1, nrhs, one, a( kc ), 1, b( k, 1 ),
231 $ ldb, b( 1, 1 ), ldb )
232*
233* Interchange if P(K) != I.
234*
235 kp = ipiv( k )
236 IF( kp.NE.k )
237 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
238 END IF
239 kc = kc + k
240 k = k + 1
241 ELSE
242*
243* 2 x 2 pivot block
244*
245 kcnext = kc + k
246*
247* Multiply by the diagonal block if forming U * D.
248*
249 IF( nounit ) THEN
250 d11 = a( kcnext-1 )
251 d22 = a( kcnext+k )
252 d12 = a( kcnext+k-1 )
253 d21 = dconjg( d12 )
254 DO 20 j = 1, nrhs
255 t1 = b( k, j )
256 t2 = b( k+1, j )
257 b( k, j ) = d11*t1 + d12*t2
258 b( k+1, j ) = d21*t1 + d22*t2
259 20 CONTINUE
260 END IF
261*
262* Multiply by P(K) * inv(U(K)) if K > 1.
263*
264 IF( k.GT.1 ) THEN
265*
266* Apply the transformations.
267*
268 CALL zgeru( k-1, nrhs, one, a( kc ), 1, b( k, 1 ),
269 $ ldb, b( 1, 1 ), ldb )
270 CALL zgeru( k-1, nrhs, one, a( kcnext ), 1,
271 $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
272*
273* Interchange if P(K) != I.
274*
275 kp = abs( ipiv( k ) )
276 IF( kp.NE.k )
277 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
278 END IF
279 kc = kcnext + k + 1
280 k = k + 2
281 END IF
282 GO TO 10
283 30 CONTINUE
284*
285* Compute B := L*B
286* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
287*
288 ELSE
289*
290* Loop backward applying the transformations to B.
291*
292 k = n
293 kc = n*( n+1 ) / 2 + 1
294 40 CONTINUE
295 IF( k.LT.1 )
296 $ GO TO 60
297 kc = kc - ( n-k+1 )
298*
299* Test the pivot index. If greater than zero, a 1 x 1
300* pivot was used, otherwise a 2 x 2 pivot was used.
301*
302 IF( ipiv( k ).GT.0 ) THEN
303*
304* 1 x 1 pivot block:
305*
306* Multiply by the diagonal element if forming L * D.
307*
308 IF( nounit )
309 $ CALL zscal( nrhs, a( kc ), b( k, 1 ), ldb )
310*
311* Multiply by P(K) * inv(L(K)) if K < N.
312*
313 IF( k.NE.n ) THEN
314 kp = ipiv( k )
315*
316* Apply the transformation.
317*
318 CALL zgeru( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
319 $ 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 IF( kp.NE.k )
325 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
326 END IF
327 k = k - 1
328*
329 ELSE
330*
331* 2 x 2 pivot block:
332*
333 kcnext = kc - ( n-k+2 )
334*
335* Multiply by the diagonal block if forming L * D.
336*
337 IF( nounit ) THEN
338 d11 = a( kcnext )
339 d22 = a( kc )
340 d21 = a( kcnext+1 )
341 d12 = dconjg( d21 )
342 DO 50 j = 1, nrhs
343 t1 = b( k-1, j )
344 t2 = b( k, j )
345 b( k-1, j ) = d11*t1 + d12*t2
346 b( k, j ) = d21*t1 + d22*t2
347 50 CONTINUE
348 END IF
349*
350* Multiply by P(K) * inv(L(K)) if K < N.
351*
352 IF( k.NE.n ) THEN
353*
354* Apply the transformation.
355*
356 CALL zgeru( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
357 $ ldb, b( k+1, 1 ), ldb )
358 CALL zgeru( n-k, nrhs, one, a( kcnext+2 ), 1,
359 $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
360*
361* Interchange if a permutation was applied at the
362* K-th step of the factorization.
363*
364 kp = abs( ipiv( k ) )
365 IF( kp.NE.k )
366 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
367 END IF
368 kc = kcnext
369 k = k - 2
370 END IF
371 GO TO 40
372 60 CONTINUE
373 END IF
374*-------------------------------------------------
375*
376* Compute B := A^H * B (conjugate transpose)
377*
378*-------------------------------------------------
379 ELSE
380*
381* Form B := U^H*B
382* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
383* and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
384*
385 IF( lsame( uplo, 'U' ) ) THEN
386*
387* Loop backward applying the transformations.
388*
389 k = n
390 kc = n*( n+1 ) / 2 + 1
391 70 CONTINUE
392 IF( k.LT.1 )
393 $ GO TO 90
394 kc = kc - k
395*
396* 1 x 1 pivot block.
397*
398 IF( ipiv( k ).GT.0 ) THEN
399 IF( k.GT.1 ) THEN
400*
401* Interchange if P(K) != I.
402*
403 kp = ipiv( k )
404 IF( kp.NE.k )
405 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
406*
407* Apply the transformation:
408* y := y - B' * conjg(x)
409* where x is a column of A and y is a row of B.
410*
411 CALL zlacgv( nrhs, b( k, 1 ), ldb )
412 CALL zgemv( 'Conjugate', k-1, nrhs, one, b, ldb,
413 $ a( kc ), 1, one, b( k, 1 ), ldb )
414 CALL zlacgv( nrhs, b( k, 1 ), ldb )
415 END IF
416 IF( nounit )
417 $ CALL zscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
418 k = k - 1
419*
420* 2 x 2 pivot block.
421*
422 ELSE
423 kcnext = kc - ( k-1 )
424 IF( k.GT.2 ) THEN
425*
426* Interchange if P(K) != I.
427*
428 kp = abs( ipiv( k ) )
429 IF( kp.NE.k-1 )
430 $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
431 $ ldb )
432*
433* Apply the transformations.
434*
435 CALL zlacgv( nrhs, b( k, 1 ), ldb )
436 CALL zgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
437 $ a( kc ), 1, one, b( k, 1 ), ldb )
438 CALL zlacgv( nrhs, b( k, 1 ), ldb )
439*
440 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
441 CALL zgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
442 $ a( kcnext ), 1, one, b( k-1, 1 ), ldb )
443 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
444 END IF
445*
446* Multiply by the diagonal block if non-unit.
447*
448 IF( nounit ) THEN
449 d11 = a( kc-1 )
450 d22 = a( kc+k-1 )
451 d12 = a( kc+k-2 )
452 d21 = dconjg( d12 )
453 DO 80 j = 1, nrhs
454 t1 = b( k-1, j )
455 t2 = b( k, j )
456 b( k-1, j ) = d11*t1 + d12*t2
457 b( k, j ) = d21*t1 + d22*t2
458 80 CONTINUE
459 END IF
460 kc = kcnext
461 k = k - 2
462 END IF
463 GO TO 70
464 90 CONTINUE
465*
466* Form B := L^H*B
467* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
468* and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
469*
470 ELSE
471*
472* Loop forward applying the L-transformations.
473*
474 k = 1
475 kc = 1
476 100 CONTINUE
477 IF( k.GT.n )
478 $ GO TO 120
479*
480* 1 x 1 pivot block
481*
482 IF( ipiv( k ).GT.0 ) THEN
483 IF( k.LT.n ) THEN
484*
485* Interchange if P(K) != I.
486*
487 kp = ipiv( k )
488 IF( kp.NE.k )
489 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
490*
491* Apply the transformation
492*
493 CALL zlacgv( nrhs, b( k, 1 ), ldb )
494 CALL zgemv( 'Conjugate', n-k, nrhs, one, b( k+1, 1 ),
495 $ ldb, a( kc+1 ), 1, one, b( k, 1 ), ldb )
496 CALL zlacgv( nrhs, b( k, 1 ), ldb )
497 END IF
498 IF( nounit )
499 $ CALL zscal( nrhs, a( kc ), b( k, 1 ), ldb )
500 kc = kc + n - k + 1
501 k = k + 1
502*
503* 2 x 2 pivot block.
504*
505 ELSE
506 kcnext = kc + n - k + 1
507 IF( k.LT.n-1 ) THEN
508*
509* Interchange if P(K) != I.
510*
511 kp = abs( ipiv( k ) )
512 IF( kp.NE.k+1 )
513 $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
514 $ ldb )
515*
516* Apply the transformation
517*
518 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
519 CALL zgemv( 'Conjugate', n-k-1, nrhs, one,
520 $ b( k+2, 1 ), ldb, a( kcnext+1 ), 1, one,
521 $ b( k+1, 1 ), ldb )
522 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
523*
524 CALL zlacgv( nrhs, b( k, 1 ), ldb )
525 CALL zgemv( 'Conjugate', n-k-1, nrhs, one,
526 $ b( k+2, 1 ), ldb, a( kc+2 ), 1, one,
527 $ b( k, 1 ), ldb )
528 CALL zlacgv( nrhs, b( k, 1 ), ldb )
529 END IF
530*
531* Multiply by the diagonal block if non-unit.
532*
533 IF( nounit ) THEN
534 d11 = a( kc )
535 d22 = a( kcnext )
536 d21 = a( kc+1 )
537 d12 = dconjg( d21 )
538 DO 110 j = 1, nrhs
539 t1 = b( k, j )
540 t2 = b( k+1, j )
541 b( k, j ) = d11*t1 + d12*t2
542 b( k+1, j ) = d21*t1 + d22*t2
543 110 CONTINUE
544 END IF
545 kc = kcnext + ( n-k )
546 k = k + 2
547 END IF
548 GO TO 100
549 120 CONTINUE
550 END IF
551*
552 END IF
553 RETURN
554*
555* End of ZLAVHP
556*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
Definition zgeru.f:130
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:74
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: