LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
clavhe.f
Go to the documentation of this file.
1 *> \brief \b CLAVHE
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 CLAVHE( 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 A( LDA, * ), B( LDB, * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> CLAVHE 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 CHETRF.
33 *> CHETRF 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', CLAVHE multiplies either by U or U * D
44 *> (or L or L * D).
45 *> If TRANS = 'C' or 'c', CLAVHE 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 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 CSYTRF or CHETRF.
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 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 complex_lin
136 *
137 * =====================================================================
138  SUBROUTINE clavhe( 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 a( lda, * ), b( ldb, * )
153 * ..
154 *
155 * =====================================================================
156 *
157 * .. Parameters ..
158  COMPLEX one
159  parameter( one = ( 1.0e+0, 0.0e+0 ) )
160 * ..
161 * .. Local Scalars ..
162  LOGICAL nounit
163  INTEGER j, k, kp
164  COMPLEX d11, d12, d21, d22, t1, t2
165 * ..
166 * .. External Functions ..
167  LOGICAL lsame
168  EXTERNAL lsame
169 * ..
170 * .. External Subroutines ..
171  EXTERNAL cgemv, cgeru, clacgv, cscal, cswap, xerbla
172 * ..
173 * .. Intrinsic Functions ..
174  INTRINSIC abs, conjg, 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( 'CLAVHE ', -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 cscal( 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 cgeru( 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 cswap( 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 = conjg( 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 cgeru( k-1, nrhs, one, a( 1, k ), 1, b( k, 1 ),
276  $ ldb, b( 1, 1 ), ldb )
277  CALL cgeru( 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 cswap( 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 cscal( 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 cgeru( n-k, nrhs, one, a( k+1, k ), 1,
323  $ b( k, 1 ), 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 cswap( 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 = conjg( 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 cgeru( n-k, nrhs, one, a( k+1, k ), 1,
359  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
360  CALL cgeru( 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 cswap( 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 IF( k.LT.1 )
392  $ go to 90
393 *
394 * 1 x 1 pivot block.
395 *
396  IF( ipiv( k ).GT.0 ) THEN
397  IF( k.GT.1 ) THEN
398 *
399 * Interchange if P(K) != I.
400 *
401  kp = ipiv( k )
402  IF( kp.NE.k )
403  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
404 *
405 * Apply the transformation
406 * y = y - B' conjg(x),
407 * where x is a column of A and y is a row of B.
408 *
409  CALL clacgv( nrhs, b( k, 1 ), ldb )
410  CALL cgemv( 'Conjugate', k-1, nrhs, one, b, ldb,
411  $ a( 1, k ), 1, one, b( k, 1 ), ldb )
412  CALL clacgv( nrhs, b( k, 1 ), ldb )
413  END IF
414  IF( nounit )
415  $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
416  k = k - 1
417 *
418 * 2 x 2 pivot block.
419 *
420  ELSE
421  IF( k.GT.2 ) THEN
422 *
423 * Interchange if P(K) != I.
424 *
425  kp = abs( ipiv( k ) )
426  IF( kp.NE.k-1 )
427  $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
428  $ ldb )
429 *
430 * Apply the transformations
431 * y = y - B' conjg(x),
432 * where x is a block column of A and y is a block
433 * row of B.
434 *
435  CALL clacgv( nrhs, b( k, 1 ), ldb )
436  CALL cgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
437  $ a( 1, k ), 1, one, b( k, 1 ), ldb )
438  CALL clacgv( nrhs, b( k, 1 ), ldb )
439 *
440  CALL clacgv( nrhs, b( k-1, 1 ), ldb )
441  CALL cgemv( 'Conjugate', k-2, nrhs, one, b, ldb,
442  $ a( 1, k-1 ), 1, one, b( k-1, 1 ), ldb )
443  CALL clacgv( 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( k-1, k-1 )
450  d22 = a( k, k )
451  d12 = a( k-1, k )
452  d21 = conjg( 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  k = k - 2
461  END IF
462  go to 70
463  90 continue
464 *
465 * Form B := L^H*B
466 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
467 * and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
468 *
469  ELSE
470 *
471 * Loop forward applying the L-transformations.
472 *
473  k = 1
474  100 continue
475  IF( k.GT.n )
476  $ go to 120
477 *
478 * 1 x 1 pivot block
479 *
480  IF( ipiv( k ).GT.0 ) THEN
481  IF( k.LT.n ) THEN
482 *
483 * Interchange if P(K) != I.
484 *
485  kp = ipiv( k )
486  IF( kp.NE.k )
487  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
488 *
489 * Apply the transformation
490 *
491  CALL clacgv( nrhs, b( k, 1 ), ldb )
492  CALL cgemv( 'Conjugate', n-k, nrhs, one, b( k+1, 1 ),
493  $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
494  CALL clacgv( nrhs, b( k, 1 ), ldb )
495  END IF
496  IF( nounit )
497  $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
498  k = k + 1
499 *
500 * 2 x 2 pivot block.
501 *
502  ELSE
503  IF( k.LT.n-1 ) THEN
504 *
505 * Interchange if P(K) != I.
506 *
507  kp = abs( ipiv( k ) )
508  IF( kp.NE.k+1 )
509  $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
510  $ ldb )
511 *
512 * Apply the transformation
513 *
514  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
515  CALL cgemv( 'Conjugate', n-k-1, nrhs, one,
516  $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, one,
517  $ b( k+1, 1 ), ldb )
518  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
519 *
520  CALL clacgv( nrhs, b( k, 1 ), ldb )
521  CALL cgemv( 'Conjugate', n-k-1, nrhs, one,
522  $ b( k+2, 1 ), ldb, a( k+2, k ), 1, one,
523  $ b( k, 1 ), ldb )
524  CALL clacgv( nrhs, b( k, 1 ), ldb )
525  END IF
526 *
527 * Multiply by the diagonal block if non-unit.
528 *
529  IF( nounit ) THEN
530  d11 = a( k, k )
531  d22 = a( k+1, k+1 )
532  d21 = a( k+1, k )
533  d12 = conjg( d21 )
534  DO 110 j = 1, nrhs
535  t1 = b( k, j )
536  t2 = b( k+1, j )
537  b( k, j ) = d11*t1 + d12*t2
538  b( k+1, j ) = d21*t1 + d22*t2
539  110 continue
540  END IF
541  k = k + 2
542  END IF
543  go to 100
544  120 continue
545  END IF
546 *
547  END IF
548  return
549 *
550 * End of CLAVHE
551 *
552  END