LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zlavhe_rook.f
Go to the documentation of this file.
1 *> \brief \b ZLAVHE_ROOK
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_ROOK( 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 *> ZLAVHE_ROOK performs one of the matrix-vector operations
28 *> x := A*x or x := A^H*x,
29 *> where x is an N element vector and A is one of the factors
30 *> from the block U*D*U' or L*D*L' factorization computed by ZHETRF_ROOK.
31 *>
32 *> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
33 *> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
34 *
35 * Arguments:
36 * ==========
37 *
38 *> \param[in] UPLO
39 *> \verbatim
40 *> UPLO is CHARACTER*1
41 *> Specifies whether the factor stored in A is upper or lower
42 *> triangular.
43 *> = 'U': Upper triangular
44 *> = 'L': Lower triangular
45 *> \endverbatim
46 *>
47 *> \param[in] TRANS
48 *> \verbatim
49 *> TRANS is CHARACTER*1
50 *> Specifies the operation to be performed:
51 *> = 'N': x := A*x
52 *> = 'C': x := A^H*x
53 *> \endverbatim
54 *>
55 *> \param[in] DIAG
56 *> \verbatim
57 *> DIAG is CHARACTER*1
58 *> Specifies whether or not the diagonal blocks are unit
59 *> matrices. If the diagonal blocks are assumed to be unit,
60 *> then A = U or A = L, otherwise A = U*D or A = L*D.
61 *> = 'U': Diagonal blocks are assumed to be unit matrices.
62 *> = 'N': Diagonal blocks are assumed to be non-unit matrices.
63 *> \endverbatim
64 *>
65 *> \param[in] N
66 *> \verbatim
67 *> N is INTEGER
68 *> The number of rows and columns of the matrix A. N >= 0.
69 *> \endverbatim
70 *>
71 *> \param[in] NRHS
72 *> \verbatim
73 *> NRHS is INTEGER
74 *> The number of right hand sides, i.e., the number of vectors
75 *> x to be multiplied by A. NRHS >= 0.
76 *> \endverbatim
77 *>
78 *> \param[in] A
79 *> \verbatim
80 *> A is COMPLEX*16 array, dimension (LDA,N)
81 *> The block diagonal matrix D and the multipliers used to
82 *> obtain the factor U or L as computed by ZHETRF_ROOK.
83 *> Stored as a 2-D triangular matrix.
84 *> \endverbatim
85 *>
86 *> \param[in] LDA
87 *> \verbatim
88 *> LDA is INTEGER
89 *> The leading dimension of the array A. LDA >= max(1,N).
90 *> \endverbatim
91 *>
92 *> \param[out] IPIV
93 *> \verbatim
94 *> IPIV is INTEGER array, dimension (N)
95 *> Details of the interchanges and the block structure of D,
96 *> as determined by ZHETRF_ROOK.
97 *> If UPLO = 'U':
98 *> Only the last KB elements of IPIV are set.
99 *>
100 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
101 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
102 *>
103 *> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
104 *> columns k and -IPIV(k) were interchanged and rows and
105 *> columns k-1 and -IPIV(k-1) were inerchaged,
106 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
107 *>
108 *> If UPLO = 'L':
109 *> Only the first KB elements of IPIV are set.
110 *>
111 *> If IPIV(k) > 0, then rows and columns k and IPIV(k)
112 *> were interchanged and D(k,k) is a 1-by-1 diagonal block.
113 *>
114 *> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
115 *> columns k and -IPIV(k) were interchanged and rows and
116 *> columns k+1 and -IPIV(k+1) were inerchaged,
117 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
118 *> \endverbatim
119 *>
120 *> \param[in,out] B
121 *> \verbatim
122 *> B is COMPLEX array, dimension (LDB,NRHS)
123 *> On entry, B contains NRHS vectors of length N.
124 *> On exit, B is overwritten with the product A * B.
125 *> \endverbatim
126 *>
127 *> \param[in] LDB
128 *> \verbatim
129 *> LDB is INTEGER
130 *> The leading dimension of the array B. LDB >= max(1,N).
131 *> \endverbatim
132 *>
133 *> \param[out] INFO
134 *> \verbatim
135 *> INFO is INTEGER
136 *> = 0: successful exit
137 *> < 0: if INFO = -k, the k-th argument had an illegal value
138 *> \endverbatim
139 *
140 * Authors:
141 * ========
142 *
143 *> \author Univ. of Tennessee
144 *> \author Univ. of California Berkeley
145 *> \author Univ. of Colorado Denver
146 *> \author NAG Ltd.
147 *
148 *> \date November 2013
149 *
150 *> \ingroup complex16_lin
151 *
152 * =====================================================================
153  SUBROUTINE zlavhe_rook( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
154  $ b, ldb, info )
155 *
156 * -- LAPACK test routine (version 3.5.0) --
157 * -- LAPACK is a software package provided by Univ. of Tennessee, --
158 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159 * November 2013
160 *
161 * .. Scalar Arguments ..
162  CHARACTER DIAG, TRANS, UPLO
163  INTEGER INFO, LDA, LDB, N, NRHS
164 * ..
165 * .. Array Arguments ..
166  INTEGER IPIV( * )
167  COMPLEX*16 A( lda, * ), B( ldb, * )
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..
173  COMPLEX*16 CONE
174  parameter ( cone = ( 1.0d+0, 0.0d+0 ) )
175 * ..
176 * .. Local Scalars ..
177  LOGICAL NOUNIT
178  INTEGER J, K, KP
179  COMPLEX*16 D11, D12, D21, D22, T1, T2
180 * ..
181 * .. External Functions ..
182  LOGICAL LSAME
183  EXTERNAL lsame
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL zgemv, zgeru, zlacgv, zscal, zswap, xerbla
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC abs, dconjg, max
190 * ..
191 * .. Executable Statements ..
192 *
193 * Test the input parameters.
194 *
195  info = 0
196  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
197  info = -1
198  ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'C' ) )
199  $ THEN
200  info = -2
201  ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
202  $ THEN
203  info = -3
204  ELSE IF( n.LT.0 ) THEN
205  info = -4
206  ELSE IF( lda.LT.max( 1, n ) ) THEN
207  info = -6
208  ELSE IF( ldb.LT.max( 1, n ) ) THEN
209  info = -9
210  END IF
211  IF( info.NE.0 ) THEN
212  CALL xerbla( 'ZLAVHE_ROOK ', -info )
213  RETURN
214  END IF
215 *
216 * Quick return if possible.
217 *
218  IF( n.EQ.0 )
219  $ RETURN
220 *
221  nounit = lsame( diag, 'N' )
222 *------------------------------------------
223 *
224 * Compute B := A * B (No transpose)
225 *
226 *------------------------------------------
227  IF( lsame( trans, 'N' ) ) THEN
228 *
229 * Compute B := U*B
230 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
231 *
232  IF( lsame( uplo, 'U' ) ) THEN
233 *
234 * Loop forward applying the transformations.
235 *
236  k = 1
237  10 CONTINUE
238  IF( k.GT.n )
239  $ GO TO 30
240  IF( ipiv( k ).GT.0 ) THEN
241 *
242 * 1 x 1 pivot block
243 *
244 * Multiply by the diagonal element if forming U * D.
245 *
246  IF( nounit )
247  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
248 *
249 * Multiply by P(K) * inv(U(K)) if K > 1.
250 *
251  IF( k.GT.1 ) THEN
252 *
253 * Apply the transformation.
254 *
255  CALL zgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
256  $ ldb, b( 1, 1 ), ldb )
257 *
258 * Interchange if P(K) != I.
259 *
260  kp = ipiv( k )
261  IF( kp.NE.k )
262  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
263  END IF
264  k = k + 1
265  ELSE
266 *
267 * 2 x 2 pivot block
268 *
269 * Multiply by the diagonal block if forming U * D.
270 *
271  IF( nounit ) THEN
272  d11 = a( k, k )
273  d22 = a( k+1, k+1 )
274  d12 = a( k, k+1 )
275  d21 = dconjg( d12 )
276  DO 20 j = 1, nrhs
277  t1 = b( k, j )
278  t2 = b( k+1, j )
279  b( k, j ) = d11*t1 + d12*t2
280  b( k+1, j ) = d21*t1 + d22*t2
281  20 CONTINUE
282  END IF
283 *
284 * Multiply by P(K) * inv(U(K)) if K > 1.
285 *
286  IF( k.GT.1 ) THEN
287 *
288 * Apply the transformations.
289 *
290  CALL zgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
291  $ ldb, b( 1, 1 ), ldb )
292  CALL zgeru( k-1, nrhs, cone, a( 1, k+1 ), 1,
293  $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
294 *
295 * Interchange if a permutation was applied at the
296 * K-th step of the factorization.
297 *
298 * Swap the first of pair with IMAXth
299 *
300  kp = abs( ipiv( k ) )
301  IF( kp.NE.k )
302  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
303 *
304 * NOW swap the first of pair with Pth
305 *
306  kp = abs( ipiv( k+1 ) )
307  IF( kp.NE.k+1 )
308  $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
309  $ ldb )
310  END IF
311  k = k + 2
312  END IF
313  GO TO 10
314  30 CONTINUE
315 *
316 * Compute B := L*B
317 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
318 *
319  ELSE
320 *
321 * Loop backward applying the transformations to B.
322 *
323  k = n
324  40 CONTINUE
325  IF( k.LT.1 )
326  $ GO TO 60
327 *
328 * Test the pivot index. If greater than zero, a 1 x 1
329 * pivot was used, otherwise a 2 x 2 pivot was used.
330 *
331  IF( ipiv( k ).GT.0 ) THEN
332 *
333 * 1 x 1 pivot block:
334 *
335 * Multiply by the diagonal element if forming L * D.
336 *
337  IF( nounit )
338  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
339 *
340 * Multiply by P(K) * inv(L(K)) if K < N.
341 *
342  IF( k.NE.n ) THEN
343  kp = ipiv( k )
344 *
345 * Apply the transformation.
346 *
347  CALL zgeru( n-k, nrhs, cone, a( k+1, k ), 1,
348  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
349 *
350 * Interchange if a permutation was applied at the
351 * K-th step of the factorization.
352 *
353  IF( kp.NE.k )
354  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
355  END IF
356  k = k - 1
357 *
358  ELSE
359 *
360 * 2 x 2 pivot block:
361 *
362 * Multiply by the diagonal block if forming L * D.
363 *
364  IF( nounit ) THEN
365  d11 = a( k-1, k-1 )
366  d22 = a( k, k )
367  d21 = a( k, k-1 )
368  d12 = dconjg( d21 )
369  DO 50 j = 1, nrhs
370  t1 = b( k-1, j )
371  t2 = b( k, j )
372  b( k-1, j ) = d11*t1 + d12*t2
373  b( k, j ) = d21*t1 + d22*t2
374  50 CONTINUE
375  END IF
376 *
377 * Multiply by P(K) * inv(L(K)) if K < N.
378 *
379  IF( k.NE.n ) THEN
380 *
381 * Apply the transformation.
382 *
383  CALL zgeru( n-k, nrhs, cone, a( k+1, k ), 1,
384  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
385  CALL zgeru( n-k, nrhs, cone, a( k+1, k-1 ), 1,
386  $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
387 *
388 * Interchange if a permutation was applied at the
389 * K-th step of the factorization.
390 *
391 *
392 * Swap the second of pair with IMAXth
393 *
394  kp = abs( ipiv( k ) )
395  IF( kp.NE.k )
396  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
397 *
398 * NOW swap the first of pair with Pth
399 *
400  kp = abs( ipiv( k-1 ) )
401  IF( kp.NE.k-1 )
402  $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
403  $ ldb )
404 *
405  END IF
406  k = k - 2
407  END IF
408  GO TO 40
409  60 CONTINUE
410  END IF
411 *--------------------------------------------------
412 *
413 * Compute B := A^H * B (conjugate transpose)
414 *
415 *--------------------------------------------------
416  ELSE
417 *
418 * Form B := U^H*B
419 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
420 * and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
421 *
422  IF( lsame( uplo, 'U' ) ) THEN
423 *
424 * Loop backward applying the transformations.
425 *
426  k = n
427  70 IF( k.LT.1 )
428  $ GO TO 90
429 *
430 * 1 x 1 pivot block.
431 *
432  IF( ipiv( k ).GT.0 ) THEN
433  IF( k.GT.1 ) THEN
434 *
435 * Interchange if P(K) != I.
436 *
437  kp = ipiv( k )
438  IF( kp.NE.k )
439  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
440 *
441 * Apply the transformation
442 * y = y - B' DCONJG(x),
443 * where x is a column of A and y is a row of B.
444 *
445  CALL zlacgv( nrhs, b( k, 1 ), ldb )
446  CALL zgemv( 'Conjugate', k-1, nrhs, cone, b, ldb,
447  $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
448  CALL zlacgv( nrhs, b( k, 1 ), ldb )
449  END IF
450  IF( nounit )
451  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
452  k = k - 1
453 *
454 * 2 x 2 pivot block.
455 *
456  ELSE
457  IF( k.GT.2 ) THEN
458 *
459 * Swap the second of pair with Pth
460 *
461  kp = abs( ipiv( k ) )
462  IF( kp.NE.k )
463  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
464 *
465 * Now swap the first of pair with IMAX(r)th
466 *
467  kp = abs( ipiv( k-1 ) )
468  IF( kp.NE.k-1 )
469  $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
470  $ ldb )
471 *
472 * Apply the transformations
473 * y = y - B' DCONJG(x),
474 * where x is a block column of A and y is a block
475 * row of B.
476 *
477  CALL zlacgv( nrhs, b( k, 1 ), ldb )
478  CALL zgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
479  $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
480  CALL zlacgv( nrhs, b( k, 1 ), ldb )
481 *
482  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
483  CALL zgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
484  $ a( 1, k-1 ), 1, cone, b( k-1, 1 ), ldb )
485  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
486  END IF
487 *
488 * Multiply by the diagonal block if non-unit.
489 *
490  IF( nounit ) THEN
491  d11 = a( k-1, k-1 )
492  d22 = a( k, k )
493  d12 = a( k-1, k )
494  d21 = dconjg( d12 )
495  DO 80 j = 1, nrhs
496  t1 = b( k-1, j )
497  t2 = b( k, j )
498  b( k-1, j ) = d11*t1 + d12*t2
499  b( k, j ) = d21*t1 + d22*t2
500  80 CONTINUE
501  END IF
502  k = k - 2
503  END IF
504  GO TO 70
505  90 CONTINUE
506 *
507 * Form B := L^H*B
508 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
509 * and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
510 *
511  ELSE
512 *
513 * Loop forward applying the L-transformations.
514 *
515  k = 1
516  100 CONTINUE
517  IF( k.GT.n )
518  $ GO TO 120
519 *
520 * 1 x 1 pivot block
521 *
522  IF( ipiv( k ).GT.0 ) THEN
523  IF( k.LT.n ) THEN
524 *
525 * Interchange if P(K) != I.
526 *
527  kp = ipiv( k )
528  IF( kp.NE.k )
529  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
530 *
531 * Apply the transformation
532 *
533  CALL zlacgv( nrhs, b( k, 1 ), ldb )
534  CALL zgemv( 'Conjugate', n-k, nrhs, cone, b( k+1, 1 ),
535  $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
536  CALL zlacgv( nrhs, b( k, 1 ), ldb )
537  END IF
538  IF( nounit )
539  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
540  k = k + 1
541 *
542 * 2 x 2 pivot block.
543 *
544  ELSE
545  IF( k.LT.n-1 ) THEN
546 *
547 * Swap the first of pair with Pth
548 *
549  kp = abs( ipiv( k ) )
550  IF( kp.NE.k )
551  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
552 *
553 * Now swap the second of pair with IMAX(r)th
554 *
555  kp = abs( ipiv( k+1 ) )
556  IF( kp.NE.k+1 )
557  $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
558  $ ldb )
559 *
560 * Apply the transformation
561 *
562  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
563  CALL zgemv( 'Conjugate', n-k-1, nrhs, cone,
564  $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, cone,
565  $ b( k+1, 1 ), ldb )
566  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
567 *
568  CALL zlacgv( nrhs, b( k, 1 ), ldb )
569  CALL zgemv( 'Conjugate', n-k-1, nrhs, cone,
570  $ b( k+2, 1 ), ldb, a( k+2, k ), 1, cone,
571  $ b( k, 1 ), ldb )
572  CALL zlacgv( nrhs, b( k, 1 ), ldb )
573  END IF
574 *
575 * Multiply by the diagonal block if non-unit.
576 *
577  IF( nounit ) THEN
578  d11 = a( k, k )
579  d22 = a( k+1, k+1 )
580  d21 = a( k+1, k )
581  d12 = dconjg( d21 )
582  DO 110 j = 1, nrhs
583  t1 = b( k, j )
584  t2 = b( k+1, j )
585  b( k, j ) = d11*t1 + d12*t2
586  b( k+1, j ) = d21*t1 + d22*t2
587  110 CONTINUE
588  END IF
589  k = k + 2
590  END IF
591  GO TO 100
592  120 CONTINUE
593  END IF
594 *
595  END IF
596  RETURN
597 *
598 * End of ZLAVHE_ROOK
599 *
600  END
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zlavhe_rook
subroutine zlavhe_rook(UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZLAVHE_ROOK
Definition: zlavhe_rook.f:155
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zgeru
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
Definition: zgeru.f:132
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
zlacgv
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76