LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
zgeru
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
Definition zgeru.f:130
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:72
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
zlavhe_rook
subroutine zlavhe_rook(uplo, trans, diag, n, nrhs, a, lda, ipiv, b, ldb, info)
ZLAVHE_ROOK
Definition zlavhe_rook.f:153