LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlavhe.f
Go to the documentation of this file.
1 *> \brief \b ZLAVHE
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 ZLAVHE( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
12 * LDB, INFO )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER DIAG, TRANS, UPLO
16 * INTEGER INFO, LDA, LDB, N, NRHS
17 * ..
18 * .. Array Arguments ..
19 * INTEGER IPIV( * )
20 * COMPLEX*16 A( LDA, * ), B( LDB, * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> ZLAVHE 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 ZHETRF.
33 *> ZHETRF 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 in the leading upper or lower triangle of the 2-D
41 *> array A.
42 *>
43 *> If TRANS = 'N' or 'n', ZLAVHE multiplies either by U or U * D
44 *> (or L or L * D).
45 *> If TRANS = 'C' or 'c', ZLAVHE multiplies either by U^H or D * U^H
46 *> (or L^H or D * L^H ).
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \verbatim
53 *> UPLO - CHARACTER*1
54 *> On entry, UPLO specifies whether the triangular matrix
55 *> stored in A is upper or lower triangular.
56 *> UPLO = 'U' or 'u' The matrix is upper triangular.
57 *> UPLO = 'L' or 'l' The matrix is lower triangular.
58 *> Unchanged on exit.
59 *>
60 *> TRANS - CHARACTER*1
61 *> On entry, TRANS specifies the operation to be performed as
62 *> follows:
63 *> TRANS = 'N' or 'n' x := A*x.
64 *> TRANS = 'C' or 'c' x := A^H*x.
65 *> Unchanged on exit.
66 *>
67 *> DIAG - CHARACTER*1
68 *> On entry, DIAG specifies whether the diagonal blocks are
69 *> assumed to be unit matrices:
70 *> DIAG = 'U' or 'u' Diagonal blocks are unit matrices.
71 *> DIAG = 'N' or 'n' Diagonal blocks are non-unit.
72 *> Unchanged on exit.
73 *>
74 *> N - INTEGER
75 *> On entry, N specifies the order of the matrix A.
76 *> N must be at least zero.
77 *> Unchanged on exit.
78 *>
79 *> NRHS - INTEGER
80 *> On entry, NRHS specifies the number of right hand sides,
81 *> i.e., the number of vectors x to be multiplied by A.
82 *> NRHS must be at least zero.
83 *> Unchanged on exit.
84 *>
85 *> A - COMPLEX*16 array, dimension( LDA, N )
86 *> On entry, A contains a block diagonal matrix and the
87 *> multipliers of the transformations used to obtain it,
88 *> stored as a 2-D triangular matrix.
89 *> Unchanged on exit.
90 *>
91 *> LDA - INTEGER
92 *> On entry, LDA specifies the first dimension of A as declared
93 *> in the calling ( sub ) program. LDA must be at least
94 *> max( 1, N ).
95 *> Unchanged on exit.
96 *>
97 *> IPIV - INTEGER array, dimension( N )
98 *> On entry, IPIV contains the vector of pivot indices as
99 *> determined by ZSYTRF or ZHETRF.
100 *> If IPIV( K ) = K, no interchange was done.
101 *> If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
102 *> changed with row IPIV( K ) and a 1 x 1 pivot block was used.
103 *> If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
104 *> with row | IPIV( K ) | and a 2 x 2 pivot block was used.
105 *> If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
106 *> with row | IPIV( K ) | and a 2 x 2 pivot block was used.
107 *>
108 *> B - COMPLEX*16 array, dimension( LDB, NRHS )
109 *> On entry, B contains NRHS vectors of length N.
110 *> On exit, B is overwritten with the product A * B.
111 *>
112 *> LDB - INTEGER
113 *> On entry, LDB contains the leading dimension of B as
114 *> declared in the calling program. LDB must be at least
115 *> max( 1, N ).
116 *> Unchanged on exit.
117 *>
118 *> INFO - INTEGER
119 *> INFO is the error flag.
120 *> On exit, a value of 0 indicates a successful exit.
121 *> A negative value, say -K, indicates that the K-th argument
122 *> has an illegal value.
123 *> \endverbatim
124 *
125 * Authors:
126 * ========
127 *
128 *> \author Univ. of Tennessee
129 *> \author Univ. of California Berkeley
130 *> \author Univ. of Colorado Denver
131 *> \author NAG Ltd.
132 *
133 *> \date November 2011
134 *
135 *> \ingroup complex16_lin
136 *
137 * =====================================================================
138  SUBROUTINE zlavhe( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
139  $ ldb, info )
140 *
141 * -- LAPACK test routine (version 3.4.0) --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 * November 2011
145 *
146 * .. Scalar Arguments ..
147  CHARACTER diag, trans, uplo
148  INTEGER info, lda, ldb, n, nrhs
149 * ..
150 * .. Array Arguments ..
151  INTEGER ipiv( * )
152  COMPLEX*16 a( lda, * ), b( ldb, * )
153 * ..
154 *
155 * =====================================================================
156 *
157 * .. Parameters ..
158  COMPLEX*16 one
159  parameter( one = ( 1.0d+0, 0.0d+0 ) )
160 * ..
161 * .. Local Scalars ..
162  LOGICAL nounit
163  INTEGER j, k, kp
164  COMPLEX*16 d11, d12, d21, d22, t1, t2
165 * ..
166 * .. External Functions ..
167  LOGICAL lsame
168  EXTERNAL lsame
169 * ..
170 * .. External Subroutines ..
171  EXTERNAL xerbla, zgemv, zgeru, zlacgv, zscal, zswap
172 * ..
173 * .. Intrinsic Functions ..
174  INTRINSIC abs, dconjg, max
175 * ..
176 * .. Executable Statements ..
177 *
178 * Test the input parameters.
179 *
180  info = 0
181  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
182  info = -1
183  ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'C' ) )
184  $ THEN
185  info = -2
186  ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
187  $ THEN
188  info = -3
189  ELSE IF( n.LT.0 ) THEN
190  info = -4
191  ELSE IF( lda.LT.max( 1, n ) ) THEN
192  info = -6
193  ELSE IF( ldb.LT.max( 1, n ) ) THEN
194  info = -9
195  END IF
196  IF( info.NE.0 ) THEN
197  CALL xerbla( 'ZLAVHE ', -info )
198  return
199  END IF
200 *
201 * Quick return if possible.
202 *
203  IF( n.EQ.0 )
204  $ return
205 *
206  nounit = lsame( diag, 'N' )
207 *------------------------------------------
208 *
209 * Compute B := A * B (No transpose)
210 *
211 *------------------------------------------
212  IF( lsame( trans, 'N' ) ) THEN
213 *
214 * Compute B := U*B
215 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
216 *
217  IF( lsame( uplo, 'U' ) ) THEN
218 *
219 * Loop forward applying the transformations.
220 *
221  k = 1
222  10 continue
223  IF( k.GT.n )
224  $ go to 30
225  IF( ipiv( k ).GT.0 ) THEN
226 *
227 * 1 x 1 pivot block
228 *
229 * Multiply by the diagonal element if forming U * D.
230 *
231  IF( nounit )
232  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
233 *
234 * Multiply by P(K) * inv(U(K)) if K > 1.
235 *
236  IF( k.GT.1 ) THEN
237 *
238 * Apply the transformation.
239 *
240  CALL zgeru( k-1, nrhs, one, a( 1, k ), 1, b( k, 1 ),
241  $ ldb, b( 1, 1 ), ldb )
242 *
243 * Interchange if P(K) != I.
244 *
245  kp = ipiv( k )
246  IF( kp.NE.k )
247  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
248  END IF
249  k = k + 1
250  ELSE
251 *
252 * 2 x 2 pivot block
253 *
254 * Multiply by the diagonal block if forming U * D.
255 *
256  IF( nounit ) THEN
257  d11 = a( k, k )
258  d22 = a( k+1, k+1 )
259  d12 = a( k, k+1 )
260  d21 = dconjg( d12 )
261  DO 20 j = 1, nrhs
262  t1 = b( k, j )
263  t2 = b( k+1, j )
264  b( k, j ) = d11*t1 + d12*t2
265  b( k+1, j ) = d21*t1 + d22*t2
266  20 continue
267  END IF
268 *
269 * Multiply by P(K) * inv(U(K)) if K > 1.
270 *
271  IF( k.GT.1 ) THEN
272 *
273 * Apply the transformations.
274 *
275  CALL zgeru( k-1, nrhs, one, a( 1, k ), 1, b( k, 1 ),
276  $ ldb, b( 1, 1 ), ldb )
277  CALL zgeru( k-1, nrhs, one, a( 1, k+1 ), 1,
278  $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
279 *
280 * Interchange if P(K) != I.
281 *
282  kp = abs( ipiv( k ) )
283  IF( kp.NE.k )
284  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
285  END IF
286  k = k + 2
287  END IF
288  go to 10
289  30 continue
290 *
291 * Compute B := L*B
292 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
293 *
294  ELSE
295 *
296 * Loop backward applying the transformations to B.
297 *
298  k = n
299  40 continue
300  IF( k.LT.1 )
301  $ go to 60
302 *
303 * Test the pivot index. If greater than zero, a 1 x 1
304 * pivot was used, otherwise a 2 x 2 pivot was used.
305 *
306  IF( ipiv( k ).GT.0 ) THEN
307 *
308 * 1 x 1 pivot block:
309 *
310 * Multiply by the diagonal element if forming L * D.
311 *
312  IF( nounit )
313  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
314 *
315 * Multiply by P(K) * inv(L(K)) if K < N.
316 *
317  IF( k.NE.n ) THEN
318  kp = ipiv( k )
319 *
320 * Apply the transformation.
321 *
322  CALL zgeru( n-k, nrhs, one, a( k+1, k ), 1, b( k, 1 ),
323  $ ldb, b( k+1, 1 ), ldb )
324 *
325 * Interchange if a permutation was applied at the
326 * K-th step of the factorization.
327 *
328  IF( kp.NE.k )
329  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
330  END IF
331  k = k - 1
332 *
333  ELSE
334 *
335 * 2 x 2 pivot block:
336 *
337 * Multiply by the diagonal block if forming L * D.
338 *
339  IF( nounit ) THEN
340  d11 = a( k-1, k-1 )
341  d22 = a( k, k )
342  d21 = a( k, k-1 )
343  d12 = dconjg( d21 )
344  DO 50 j = 1, nrhs
345  t1 = b( k-1, j )
346  t2 = b( k, j )
347  b( k-1, j ) = d11*t1 + d12*t2
348  b( k, j ) = d21*t1 + d22*t2
349  50 continue
350  END IF
351 *
352 * Multiply by P(K) * inv(L(K)) if K < N.
353 *
354  IF( k.NE.n ) THEN
355 *
356 * Apply the transformation.
357 *
358  CALL zgeru( n-k, nrhs, one, a( k+1, k ), 1, b( k, 1 ),
359  $ ldb, b( k+1, 1 ), ldb )
360  CALL zgeru( n-k, nrhs, one, a( k+1, k-1 ), 1,
361  $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
362 *
363 * Interchange if a permutation was applied at the
364 * K-th step of the factorization.
365 *
366  kp = abs( ipiv( k ) )
367  IF( kp.NE.k )
368  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
369  END IF
370  k = k - 2
371  END IF
372  go to 40
373  60 continue
374  END IF
375 *--------------------------------------------------
376 *
377 * Compute B := A^H * B (conjugate transpose)
378 *
379 *--------------------------------------------------
380  ELSE
381 *
382 * Form B := U^H*B
383 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
384 * and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
385 *
386  IF( lsame( uplo, 'U' ) ) THEN
387 *
388 * Loop backward applying the transformations.
389 *
390  k = n
391  70 continue
392  IF( k.LT.1 )
393  $ go to 90
394 *
395 * 1 x 1 pivot block.
396 *
397  IF( ipiv( k ).GT.0 ) THEN
398  IF( k.GT.1 ) THEN
399 *
400 * Interchange if P(K) != I.
401 *
402  kp = ipiv( k )
403  IF( kp.NE.k )
404  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
405 *
406 * Apply the transformation
407 * y = y - B' conjg(x),
408 * where x is a column of A and y is a row of B.
409 *
410  CALL zlacgv( nrhs, b( k, 1 ), ldb )
411  CALL zgemv( 'Conjugate', k-1, nrhs, one, b, ldb,
412  $ a( 1, k ), 1, one, b( k, 1 ), ldb )
413  CALL zlacgv( nrhs, b( k, 1 ), ldb )
414  END IF
415  IF( nounit )
416  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
417  k = k - 1
418 *
419 * 2 x 2 pivot block.
420 *
421  ELSE
422  IF( k.GT.2 ) THEN
423 *
424 * Interchange if P(K) != I.
425 *
426  kp = abs( ipiv( k ) )
427  IF( kp.NE.k-1 )
428  $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
429  $ ldb )
430 *
431 * Apply the transformations
432 * y = y - B' conjg(x),
433 * where x is a block column of A and y is a block
434 * row of B.
435 *
436  CALL zlacgv( nrhs, b( k, 1 ), ldb )
437  CALL zgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
438  $ a( 1, k ), 1, one, b( k, 1 ), ldb )
439  CALL zlacgv( nrhs, b( k, 1 ), ldb )
440 *
441  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
442  CALL zgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
443  $ a( 1, k-1 ), 1, one, b( k-1, 1 ), ldb )
444  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
445  END IF
446 *
447 * Multiply by the diagonal block if non-unit.
448 *
449  IF( nounit ) THEN
450  d11 = a( k-1, k-1 )
451  d22 = a( k, k )
452  d12 = a( k-1, k )
453  d21 = dconjg( d12 )
454  DO 80 j = 1, nrhs
455  t1 = b( k-1, j )
456  t2 = b( k, j )
457  b( k-1, j ) = d11*t1 + d12*t2
458  b( k, j ) = d21*t1 + d22*t2
459  80 continue
460  END IF
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^H(m))*P(m)* ... *inv(L^H(1))*P(1)
469 *
470  ELSE
471 *
472 * Loop forward applying the L-transformations.
473 *
474  k = 1
475  100 continue
476  IF( k.GT.n )
477  $ go to 120
478 *
479 * 1 x 1 pivot block
480 *
481  IF( ipiv( k ).GT.0 ) THEN
482  IF( k.LT.n ) THEN
483 *
484 * Interchange if P(K) != I.
485 *
486  kp = ipiv( k )
487  IF( kp.NE.k )
488  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
489 *
490 * Apply the transformation
491 *
492  CALL zlacgv( nrhs, b( k, 1 ), ldb )
493  CALL zgemv( 'Conjugate', n-k, nrhs, one, b( k+1, 1 ),
494  $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
495  CALL zlacgv( nrhs, b( k, 1 ), ldb )
496  END IF
497  IF( nounit )
498  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
499  k = k + 1
500 *
501 * 2 x 2 pivot block.
502 *
503  ELSE
504  IF( k.LT.n-1 ) THEN
505 *
506 * Interchange if P(K) != I.
507 *
508  kp = abs( ipiv( k ) )
509  IF( kp.NE.k+1 )
510  $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
511  $ ldb )
512 *
513 * Apply the transformation
514 *
515  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
516  CALL zgemv( 'Conjugate', n-k-1, nrhs, one,
517  $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, one,
518  $ b( k+1, 1 ), ldb )
519  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
520 *
521  CALL zlacgv( nrhs, b( k, 1 ), ldb )
522  CALL zgemv( 'Conjugate', n-k-1, nrhs, one,
523  $ b( k+2, 1 ), ldb, a( k+2, k ), 1, one,
524  $ b( k, 1 ), ldb )
525  CALL zlacgv( nrhs, b( k, 1 ), ldb )
526  END IF
527 *
528 * Multiply by the diagonal block if non-unit.
529 *
530  IF( nounit ) THEN
531  d11 = a( k, k )
532  d22 = a( k+1, k+1 )
533  d21 = a( k+1, k )
534  d12 = dconjg( d21 )
535  DO 110 j = 1, nrhs
536  t1 = b( k, j )
537  t2 = b( k+1, j )
538  b( k, j ) = d11*t1 + d12*t2
539  b( k+1, j ) = d21*t1 + d22*t2
540  110 continue
541  END IF
542  k = k + 2
543  END IF
544  go to 100
545  120 continue
546  END IF
547 *
548  END IF
549  return
550 *
551 * End of ZLAVHE
552 *
553  END