LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zlavhp.f
Go to the documentation of this file.
1 *> \brief \b ZLAVHP
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZLAVHP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
12 * INFO )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER DIAG, TRANS, UPLO
16 * INTEGER INFO, LDB, N, NRHS
17 * ..
18 * .. Array Arguments ..
19 * INTEGER IPIV( * )
20 * COMPLEX*16 A( * ), B( LDB, * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> ZLAVHP performs one of the matrix-vector operations
30 *> x := A*x or x := A^H*x,
31 *> where x is an N element vector and A is one of the factors
32 *> from the symmetric factorization computed by ZHPTRF.
33 *> ZHPTRF produces a factorization of the form
34 *> U * D * U^H or L * D * L^H,
35 *> where U (or L) is a product of permutation and unit upper (lower)
36 *> triangular matrices, U^H (or L^H) is the conjugate transpose of
37 *> U (or L), and D is Hermitian and block diagonal with 1 x 1 and
38 *> 2 x 2 diagonal blocks. The multipliers for the transformations
39 *> and the upper or lower triangular parts of the diagonal blocks
40 *> are stored columnwise in packed format in the linear array A.
41 *>
42 *> If TRANS = 'N' or 'n', ZLAVHP multiplies either by U or U * D
43 *> (or L or L * D).
44 *> If TRANS = 'C' or 'c', ZLAVHP multiplies either by U^H or D * U^H
45 *> (or L^H or D * L^H ).
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \verbatim
52 *> UPLO - CHARACTER*1
53 *> On entry, UPLO specifies whether the triangular matrix
54 *> stored in A is upper or lower triangular.
55 *> UPLO = 'U' or 'u' The matrix is upper triangular.
56 *> UPLO = 'L' or 'l' The matrix is lower triangular.
57 *> Unchanged on exit.
58 *>
59 *> TRANS - CHARACTER*1
60 *> On entry, TRANS specifies the operation to be performed as
61 *> follows:
62 *> TRANS = 'N' or 'n' x := A*x.
63 *> TRANS = 'C' or 'c' x := A^H*x.
64 *> Unchanged on exit.
65 *>
66 *> DIAG - CHARACTER*1
67 *> On entry, DIAG specifies whether the diagonal blocks are
68 *> assumed to be unit matrices, as follows:
69 *> DIAG = 'U' or 'u' Diagonal blocks are unit matrices.
70 *> DIAG = 'N' or 'n' Diagonal blocks are non-unit.
71 *> Unchanged on exit.
72 *>
73 *> N - INTEGER
74 *> On entry, N specifies the order of the matrix A.
75 *> N must be at least zero.
76 *> Unchanged on exit.
77 *>
78 *> NRHS - INTEGER
79 *> On entry, NRHS specifies the number of right hand sides,
80 *> i.e., the number of vectors x to be multiplied by A.
81 *> NRHS must be at least zero.
82 *> Unchanged on exit.
83 *>
84 *> A - COMPLEX*16 array, dimension( N*(N+1)/2 )
85 *> On entry, A contains a block diagonal matrix and the
86 *> multipliers of the transformations used to obtain it,
87 *> stored as a packed triangular matrix.
88 *> Unchanged on exit.
89 *>
90 *> IPIV - INTEGER array, dimension( N )
91 *> On entry, IPIV contains the vector of pivot indices as
92 *> determined by ZSPTRF or ZHPTRF.
93 *> If IPIV( K ) = K, no interchange was done.
94 *> If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
95 *> changed with row IPIV( K ) and a 1 x 1 pivot block was used.
96 *> If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
97 *> with row | IPIV( K ) | and a 2 x 2 pivot block was used.
98 *> If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
99 *> with row | IPIV( K ) | and a 2 x 2 pivot block was used.
100 *>
101 *> B - COMPLEX*16 array, dimension( LDB, NRHS )
102 *> On entry, B contains NRHS vectors of length N.
103 *> On exit, B is overwritten with the product A * B.
104 *>
105 *> LDB - INTEGER
106 *> On entry, LDB contains the leading dimension of B as
107 *> declared in the calling program. LDB must be at least
108 *> max( 1, N ).
109 *> Unchanged on exit.
110 *>
111 *> INFO - INTEGER
112 *> INFO is the error flag.
113 *> On exit, a value of 0 indicates a successful exit.
114 *> A negative value, say -K, indicates that the K-th argument
115 *> has an illegal value.
116 *> \endverbatim
117 *
118 * Authors:
119 * ========
120 *
121 *> \author Univ. of Tennessee
122 *> \author Univ. of California Berkeley
123 *> \author Univ. of Colorado Denver
124 *> \author NAG Ltd.
125 *
126 *> \ingroup complex16_lin
127 *
128 * =====================================================================
129  SUBROUTINE zlavhp( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
130  $ INFO )
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 *
557  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
Definition: zgeru.f:130
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zlavhp(UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB, INFO)
ZLAVHP
Definition: zlavhp.f:131
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74