LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ clavhp()

subroutine clavhp ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
complex, dimension( * )  A,
integer, dimension( * )  IPIV,
complex, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

CLAVHP

Purpose:
    CLAVHP  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 CHPTRF.
    CHPTRF 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', CLAVHP multiplies either by U or U * D
    (or L or L * D).
    If TRANS = 'C' or 'c', CLAVHP 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 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 CSPTRF or CHPTRF.
           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 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.
Date
December 2016

Definition at line 133 of file clavhp.f.

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