LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
127 *
128 *> \ingroup complex16_lin
129 *
130 * =====================================================================
131  SUBROUTINE zlavhp( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
132  $ info )
133 *
134 * -- LAPACK test routine (version 3.4.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 * November 2011
138 *
139 * .. Scalar Arguments ..
140  CHARACTER diag, trans, uplo
141  INTEGER info, ldb, n, nrhs
142 * ..
143 * .. Array Arguments ..
144  INTEGER ipiv( * )
145  COMPLEX*16 a( * ), b( ldb, * )
146 * ..
147 *
148 * =====================================================================
149 *
150 * .. Parameters ..
151  COMPLEX*16 one
152  parameter( one = ( 1.0d+0, 0.0d+0 ) )
153 * ..
154 * .. Local Scalars ..
155  LOGICAL nounit
156  INTEGER j, k, kc, kcnext, kp
157  COMPLEX*16 d11, d12, d21, d22, t1, t2
158 * ..
159 * .. External Functions ..
160  LOGICAL lsame
161  EXTERNAL lsame
162 * ..
163 * .. External Subroutines ..
164  EXTERNAL xerbla, zgemv, zgeru, zlacgv, zscal, zswap
165 * ..
166 * .. Intrinsic Functions ..
167  INTRINSIC abs, dconjg, 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( 'ZLAVHP ', -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 zscal( 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 zgeru( 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 zswap( 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 = dconjg( 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 zgeru( k-1, nrhs, one, a( kc ), 1, b( k, 1 ),
272  $ ldb, b( 1, 1 ), ldb )
273  CALL zgeru( 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 zswap( 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 zscal( 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 zgeru( 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 zswap( 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 = dconjg( 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 zgeru( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
360  $ ldb, b( k+1, 1 ), ldb )
361  CALL zgeru( 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 zswap( 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 continue
395  IF( k.LT.1 )
396  $ go to 90
397  kc = kc - k
398 *
399 * 1 x 1 pivot block.
400 *
401  IF( ipiv( k ).GT.0 ) THEN
402  IF( k.GT.1 ) THEN
403 *
404 * Interchange if P(K) != I.
405 *
406  kp = ipiv( k )
407  IF( kp.NE.k )
408  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
409 *
410 * Apply the transformation:
411 * y := y - B' * conjg(x)
412 * where x is a column of A and y is a row of B.
413 *
414  CALL zlacgv( nrhs, b( k, 1 ), ldb )
415  CALL zgemv( 'Conjugate', k-1, nrhs, one, b, ldb,
416  $ a( kc ), 1, one, b( k, 1 ), ldb )
417  CALL zlacgv( nrhs, b( k, 1 ), ldb )
418  END IF
419  IF( nounit )
420  $ CALL zscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
421  k = k - 1
422 *
423 * 2 x 2 pivot block.
424 *
425  ELSE
426  kcnext = kc - ( k-1 )
427  IF( k.GT.2 ) THEN
428 *
429 * Interchange if P(K) != I.
430 *
431  kp = abs( ipiv( k ) )
432  IF( kp.NE.k-1 )
433  $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
434  $ ldb )
435 *
436 * Apply the transformations.
437 *
438  CALL zlacgv( nrhs, b( k, 1 ), ldb )
439  CALL zgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
440  $ a( kc ), 1, one, b( k, 1 ), ldb )
441  CALL zlacgv( nrhs, b( k, 1 ), ldb )
442 *
443  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
444  CALL zgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
445  $ a( kcnext ), 1, one, b( k-1, 1 ), ldb )
446  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
447  END IF
448 *
449 * Multiply by the diagonal block if non-unit.
450 *
451  IF( nounit ) THEN
452  d11 = a( kc-1 )
453  d22 = a( kc+k-1 )
454  d12 = a( kc+k-2 )
455  d21 = dconjg( d12 )
456  DO 80 j = 1, nrhs
457  t1 = b( k-1, j )
458  t2 = b( k, j )
459  b( k-1, j ) = d11*t1 + d12*t2
460  b( k, j ) = d21*t1 + d22*t2
461  80 continue
462  END IF
463  kc = kcnext
464  k = k - 2
465  END IF
466  go to 70
467  90 continue
468 *
469 * Form B := L^H*B
470 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
471 * and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
472 *
473  ELSE
474 *
475 * Loop forward applying the L-transformations.
476 *
477  k = 1
478  kc = 1
479  100 continue
480  IF( k.GT.n )
481  $ go to 120
482 *
483 * 1 x 1 pivot block
484 *
485  IF( ipiv( k ).GT.0 ) THEN
486  IF( k.LT.n ) THEN
487 *
488 * Interchange if P(K) != I.
489 *
490  kp = ipiv( k )
491  IF( kp.NE.k )
492  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
493 *
494 * Apply the transformation
495 *
496  CALL zlacgv( nrhs, b( k, 1 ), ldb )
497  CALL zgemv( 'Conjugate', n-k, nrhs, one, b( k+1, 1 ),
498  $ ldb, a( kc+1 ), 1, one, b( k, 1 ), ldb )
499  CALL zlacgv( nrhs, b( k, 1 ), ldb )
500  END IF
501  IF( nounit )
502  $ CALL zscal( nrhs, a( kc ), b( k, 1 ), ldb )
503  kc = kc + n - k + 1
504  k = k + 1
505 *
506 * 2 x 2 pivot block.
507 *
508  ELSE
509  kcnext = kc + n - k + 1
510  IF( k.LT.n-1 ) THEN
511 *
512 * Interchange if P(K) != I.
513 *
514  kp = abs( ipiv( k ) )
515  IF( kp.NE.k+1 )
516  $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
517  $ ldb )
518 *
519 * Apply the transformation
520 *
521  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
522  CALL zgemv( 'Conjugate', n-k-1, nrhs, one,
523  $ b( k+2, 1 ), ldb, a( kcnext+1 ), 1, one,
524  $ b( k+1, 1 ), ldb )
525  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
526 *
527  CALL zlacgv( nrhs, b( k, 1 ), ldb )
528  CALL zgemv( 'Conjugate', n-k-1, nrhs, one,
529  $ b( k+2, 1 ), ldb, a( kc+2 ), 1, one,
530  $ b( k, 1 ), ldb )
531  CALL zlacgv( nrhs, b( k, 1 ), ldb )
532  END IF
533 *
534 * Multiply by the diagonal block if non-unit.
535 *
536  IF( nounit ) THEN
537  d11 = a( kc )
538  d22 = a( kcnext )
539  d21 = a( kc+1 )
540  d12 = dconjg( d21 )
541  DO 110 j = 1, nrhs
542  t1 = b( k, j )
543  t2 = b( k+1, j )
544  b( k, j ) = d11*t1 + d12*t2
545  b( k+1, j ) = d21*t1 + d22*t2
546  110 continue
547  END IF
548  kc = kcnext + ( n-k )
549  k = k + 2
550  END IF
551  go to 100
552  120 continue
553  END IF
554 *
555  END IF
556  return
557 *
558 * End of ZLAVHP
559 *
560  END