LAPACK  3.10.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 *> \ingroup complex16_lin
149 *
150 * =====================================================================
151  SUBROUTINE zlavhe_rook( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
152  $ B, LDB, INFO )
153 *
154 * -- LAPACK test routine --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 *
158 * .. Scalar Arguments ..
159  CHARACTER DIAG, TRANS, UPLO
160  INTEGER INFO, LDA, LDB, N, NRHS
161 * ..
162 * .. Array Arguments ..
163  INTEGER IPIV( * )
164  COMPLEX*16 A( LDA, * ), B( LDB, * )
165 * ..
166 *
167 * =====================================================================
168 *
169 * .. Parameters ..
170  COMPLEX*16 CONE
171  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL NOUNIT
175  INTEGER J, K, KP
176  COMPLEX*16 D11, D12, D21, D22, T1, T2
177 * ..
178 * .. External Functions ..
179  LOGICAL LSAME
180  EXTERNAL lsame
181 * ..
182 * .. External Subroutines ..
183  EXTERNAL zgemv, zgeru, zlacgv, zscal, zswap, xerbla
184 * ..
185 * .. Intrinsic Functions ..
186  INTRINSIC abs, dconjg, max
187 * ..
188 * .. Executable Statements ..
189 *
190 * Test the input parameters.
191 *
192  info = 0
193  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
194  info = -1
195  ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'C' ) )
196  $ THEN
197  info = -2
198  ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
199  $ THEN
200  info = -3
201  ELSE IF( n.LT.0 ) THEN
202  info = -4
203  ELSE IF( lda.LT.max( 1, n ) ) THEN
204  info = -6
205  ELSE IF( ldb.LT.max( 1, n ) ) THEN
206  info = -9
207  END IF
208  IF( info.NE.0 ) THEN
209  CALL xerbla( 'ZLAVHE_ROOK ', -info )
210  RETURN
211  END IF
212 *
213 * Quick return if possible.
214 *
215  IF( n.EQ.0 )
216  $ RETURN
217 *
218  nounit = lsame( diag, 'N' )
219 *------------------------------------------
220 *
221 * Compute B := A * B (No transpose)
222 *
223 *------------------------------------------
224  IF( lsame( trans, 'N' ) ) THEN
225 *
226 * Compute B := U*B
227 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
228 *
229  IF( lsame( uplo, 'U' ) ) THEN
230 *
231 * Loop forward applying the transformations.
232 *
233  k = 1
234  10 CONTINUE
235  IF( k.GT.n )
236  $ GO TO 30
237  IF( ipiv( k ).GT.0 ) THEN
238 *
239 * 1 x 1 pivot block
240 *
241 * Multiply by the diagonal element if forming U * D.
242 *
243  IF( nounit )
244  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
245 *
246 * Multiply by P(K) * inv(U(K)) if K > 1.
247 *
248  IF( k.GT.1 ) THEN
249 *
250 * Apply the transformation.
251 *
252  CALL zgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
253  $ ldb, b( 1, 1 ), ldb )
254 *
255 * Interchange if P(K) != I.
256 *
257  kp = ipiv( k )
258  IF( kp.NE.k )
259  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
260  END IF
261  k = k + 1
262  ELSE
263 *
264 * 2 x 2 pivot block
265 *
266 * Multiply by the diagonal block if forming U * D.
267 *
268  IF( nounit ) THEN
269  d11 = a( k, k )
270  d22 = a( k+1, k+1 )
271  d12 = a( k, k+1 )
272  d21 = dconjg( d12 )
273  DO 20 j = 1, nrhs
274  t1 = b( k, j )
275  t2 = b( k+1, j )
276  b( k, j ) = d11*t1 + d12*t2
277  b( k+1, j ) = d21*t1 + d22*t2
278  20 CONTINUE
279  END IF
280 *
281 * Multiply by P(K) * inv(U(K)) if K > 1.
282 *
283  IF( k.GT.1 ) THEN
284 *
285 * Apply the transformations.
286 *
287  CALL zgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
288  $ ldb, b( 1, 1 ), ldb )
289  CALL zgeru( k-1, nrhs, cone, a( 1, k+1 ), 1,
290  $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
291 *
292 * Interchange if a permutation was applied at the
293 * K-th step of the factorization.
294 *
295 * Swap the first of pair with IMAXth
296 *
297  kp = abs( ipiv( k ) )
298  IF( kp.NE.k )
299  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
300 *
301 * NOW swap the first of pair with Pth
302 *
303  kp = abs( ipiv( k+1 ) )
304  IF( kp.NE.k+1 )
305  $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
306  $ ldb )
307  END IF
308  k = k + 2
309  END IF
310  GO TO 10
311  30 CONTINUE
312 *
313 * Compute B := L*B
314 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
315 *
316  ELSE
317 *
318 * Loop backward applying the transformations to B.
319 *
320  k = n
321  40 CONTINUE
322  IF( k.LT.1 )
323  $ GO TO 60
324 *
325 * Test the pivot index. If greater than zero, a 1 x 1
326 * pivot was used, otherwise a 2 x 2 pivot was used.
327 *
328  IF( ipiv( k ).GT.0 ) THEN
329 *
330 * 1 x 1 pivot block:
331 *
332 * Multiply by the diagonal element if forming L * D.
333 *
334  IF( nounit )
335  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
336 *
337 * Multiply by P(K) * inv(L(K)) if K < N.
338 *
339  IF( k.NE.n ) THEN
340  kp = ipiv( k )
341 *
342 * Apply the transformation.
343 *
344  CALL zgeru( n-k, nrhs, cone, a( k+1, k ), 1,
345  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
346 *
347 * Interchange if a permutation was applied at the
348 * K-th step of the factorization.
349 *
350  IF( kp.NE.k )
351  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
352  END IF
353  k = k - 1
354 *
355  ELSE
356 *
357 * 2 x 2 pivot block:
358 *
359 * Multiply by the diagonal block if forming L * D.
360 *
361  IF( nounit ) THEN
362  d11 = a( k-1, k-1 )
363  d22 = a( k, k )
364  d21 = a( k, k-1 )
365  d12 = dconjg( d21 )
366  DO 50 j = 1, nrhs
367  t1 = b( k-1, j )
368  t2 = b( k, j )
369  b( k-1, j ) = d11*t1 + d12*t2
370  b( k, j ) = d21*t1 + d22*t2
371  50 CONTINUE
372  END IF
373 *
374 * Multiply by P(K) * inv(L(K)) if K < N.
375 *
376  IF( k.NE.n ) THEN
377 *
378 * Apply the transformation.
379 *
380  CALL zgeru( n-k, nrhs, cone, a( k+1, k ), 1,
381  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
382  CALL zgeru( n-k, nrhs, cone, a( k+1, k-1 ), 1,
383  $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
384 *
385 * Interchange if a permutation was applied at the
386 * K-th step of the factorization.
387 *
388 *
389 * Swap the second of pair with IMAXth
390 *
391  kp = abs( ipiv( k ) )
392  IF( kp.NE.k )
393  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
394 *
395 * NOW swap the first of pair with Pth
396 *
397  kp = abs( ipiv( k-1 ) )
398  IF( kp.NE.k-1 )
399  $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
400  $ ldb )
401 *
402  END IF
403  k = k - 2
404  END IF
405  GO TO 40
406  60 CONTINUE
407  END IF
408 *--------------------------------------------------
409 *
410 * Compute B := A^H * B (conjugate transpose)
411 *
412 *--------------------------------------------------
413  ELSE
414 *
415 * Form B := U^H*B
416 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
417 * and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
418 *
419  IF( lsame( uplo, 'U' ) ) THEN
420 *
421 * Loop backward applying the transformations.
422 *
423  k = n
424  70 IF( k.LT.1 )
425  $ GO TO 90
426 *
427 * 1 x 1 pivot block.
428 *
429  IF( ipiv( k ).GT.0 ) THEN
430  IF( k.GT.1 ) THEN
431 *
432 * Interchange if P(K) != I.
433 *
434  kp = ipiv( k )
435  IF( kp.NE.k )
436  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
437 *
438 * Apply the transformation
439 * y = y - B' DCONJG(x),
440 * where x is a column of A and y is a row of B.
441 *
442  CALL zlacgv( nrhs, b( k, 1 ), ldb )
443  CALL zgemv( 'Conjugate', k-1, nrhs, cone, b, ldb,
444  $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
445  CALL zlacgv( nrhs, b( k, 1 ), ldb )
446  END IF
447  IF( nounit )
448  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
449  k = k - 1
450 *
451 * 2 x 2 pivot block.
452 *
453  ELSE
454  IF( k.GT.2 ) THEN
455 *
456 * Swap the second of pair with Pth
457 *
458  kp = abs( ipiv( k ) )
459  IF( kp.NE.k )
460  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
461 *
462 * Now swap the first of pair with IMAX(r)th
463 *
464  kp = abs( ipiv( k-1 ) )
465  IF( kp.NE.k-1 )
466  $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
467  $ ldb )
468 *
469 * Apply the transformations
470 * y = y - B' DCONJG(x),
471 * where x is a block column of A and y is a block
472 * row of B.
473 *
474  CALL zlacgv( nrhs, b( k, 1 ), ldb )
475  CALL zgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
476  $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
477  CALL zlacgv( nrhs, b( k, 1 ), ldb )
478 *
479  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
480  CALL zgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
481  $ a( 1, k-1 ), 1, cone, b( k-1, 1 ), ldb )
482  CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
483  END IF
484 *
485 * Multiply by the diagonal block if non-unit.
486 *
487  IF( nounit ) THEN
488  d11 = a( k-1, k-1 )
489  d22 = a( k, k )
490  d12 = a( k-1, k )
491  d21 = dconjg( d12 )
492  DO 80 j = 1, nrhs
493  t1 = b( k-1, j )
494  t2 = b( k, j )
495  b( k-1, j ) = d11*t1 + d12*t2
496  b( k, j ) = d21*t1 + d22*t2
497  80 CONTINUE
498  END IF
499  k = k - 2
500  END IF
501  GO TO 70
502  90 CONTINUE
503 *
504 * Form B := L^H*B
505 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
506 * and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
507 *
508  ELSE
509 *
510 * Loop forward applying the L-transformations.
511 *
512  k = 1
513  100 CONTINUE
514  IF( k.GT.n )
515  $ GO TO 120
516 *
517 * 1 x 1 pivot block
518 *
519  IF( ipiv( k ).GT.0 ) THEN
520  IF( k.LT.n ) THEN
521 *
522 * Interchange if P(K) != I.
523 *
524  kp = ipiv( k )
525  IF( kp.NE.k )
526  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
527 *
528 * Apply the transformation
529 *
530  CALL zlacgv( nrhs, b( k, 1 ), ldb )
531  CALL zgemv( 'Conjugate', n-k, nrhs, cone, b( k+1, 1 ),
532  $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
533  CALL zlacgv( nrhs, b( k, 1 ), ldb )
534  END IF
535  IF( nounit )
536  $ CALL zscal( nrhs, a( k, k ), b( k, 1 ), ldb )
537  k = k + 1
538 *
539 * 2 x 2 pivot block.
540 *
541  ELSE
542  IF( k.LT.n-1 ) THEN
543 *
544 * Swap the first of pair with Pth
545 *
546  kp = abs( ipiv( k ) )
547  IF( kp.NE.k )
548  $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
549 *
550 * Now swap the second of pair with IMAX(r)th
551 *
552  kp = abs( ipiv( k+1 ) )
553  IF( kp.NE.k+1 )
554  $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
555  $ ldb )
556 *
557 * Apply the transformation
558 *
559  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
560  CALL zgemv( 'Conjugate', n-k-1, nrhs, cone,
561  $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, cone,
562  $ b( k+1, 1 ), ldb )
563  CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
564 *
565  CALL zlacgv( nrhs, b( k, 1 ), ldb )
566  CALL zgemv( 'Conjugate', n-k-1, nrhs, cone,
567  $ b( k+2, 1 ), ldb, a( k+2, k ), 1, cone,
568  $ b( k, 1 ), ldb )
569  CALL zlacgv( nrhs, b( k, 1 ), ldb )
570  END IF
571 *
572 * Multiply by the diagonal block if non-unit.
573 *
574  IF( nounit ) THEN
575  d11 = a( k, k )
576  d22 = a( k+1, k+1 )
577  d21 = a( k+1, k )
578  d12 = dconjg( d21 )
579  DO 110 j = 1, nrhs
580  t1 = b( k, j )
581  t2 = b( k+1, j )
582  b( k, j ) = d11*t1 + d12*t2
583  b( k+1, j ) = d21*t1 + d22*t2
584  110 CONTINUE
585  END IF
586  k = k + 2
587  END IF
588  GO TO 100
589  120 CONTINUE
590  END IF
591 *
592  END IF
593  RETURN
594 *
595 * End of ZLAVHE_ROOK
596 *
597  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 zlavhe_rook(UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZLAVHE_ROOK
Definition: zlavhe_rook.f:153
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74