LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clavhe_rook.f
Go to the documentation of this file.
1*> \brief \b CLAVHE_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 CLAVHE_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 A( LDA, * ), B( LDB, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> CLAVHE_ROOK 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 block U*D*U' or L*D*L' factorization computed by CHETRF_ROOK.
33*>
34*> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
35*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] UPLO
42*> \verbatim
43*> UPLO is CHARACTER*1
44*> Specifies whether the factor stored in A is upper or lower
45*> triangular.
46*> = 'U': Upper triangular
47*> = 'L': Lower triangular
48*> \endverbatim
49*>
50*> \param[in] TRANS
51*> \verbatim
52*> TRANS is CHARACTER*1
53*> Specifies the operation to be performed:
54*> = 'N': x := A*x
55*> = 'C': x := A^H*x
56*> \endverbatim
57*>
58*> \param[in] DIAG
59*> \verbatim
60*> DIAG is CHARACTER*1
61*> Specifies whether or not the diagonal blocks are unit
62*> matrices. If the diagonal blocks are assumed to be unit,
63*> then A = U or A = L, otherwise A = U*D or A = L*D.
64*> = 'U': Diagonal blocks are assumed to be unit matrices.
65*> = 'N': Diagonal blocks are assumed to be non-unit matrices.
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*> N is INTEGER
71*> The number of rows and columns of the matrix A. N >= 0.
72*> \endverbatim
73*>
74*> \param[in] NRHS
75*> \verbatim
76*> NRHS is INTEGER
77*> The number of right hand sides, i.e., the number of vectors
78*> x to be multiplied by A. NRHS >= 0.
79*> \endverbatim
80*>
81*> \param[in] A
82*> \verbatim
83*> A is COMPLEX array, dimension (LDA,N)
84*> The block diagonal matrix D and the multipliers used to
85*> obtain the factor U or L as computed by CHETRF_ROOK.
86*> Stored as a 2-D triangular matrix.
87*> \endverbatim
88*>
89*> \param[in] LDA
90*> \verbatim
91*> LDA is INTEGER
92*> The leading dimension of the array A. LDA >= max(1,N).
93*> \endverbatim
94*>
95*> \param[out] IPIV
96*> \verbatim
97*> IPIV is INTEGER array, dimension (N)
98*> Details of the interchanges and the block structure of D,
99*> as determined by CHETRF_ROOK.
100*> If UPLO = 'U':
101*> Only the last KB elements of IPIV are set.
102*>
103*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
104*> interchanged and D(k,k) is a 1-by-1 diagonal block.
105*>
106*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
107*> columns k and -IPIV(k) were interchanged and rows and
108*> columns k-1 and -IPIV(k-1) were inerchaged,
109*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
110*>
111*> If UPLO = 'L':
112*> Only the first KB elements of IPIV are set.
113*>
114*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
115*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
116*>
117*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
118*> columns k and -IPIV(k) were interchanged and rows and
119*> columns k+1 and -IPIV(k+1) were inerchaged,
120*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
121*> \endverbatim
122*>
123*> \param[in,out] B
124*> \verbatim
125*> B is COMPLEX array, dimension (LDB,NRHS)
126*> On entry, B contains NRHS vectors of length N.
127*> On exit, B is overwritten with the product A * B.
128*> \endverbatim
129*>
130*> \param[in] LDB
131*> \verbatim
132*> LDB is INTEGER
133*> The leading dimension of the array B. LDB >= max(1,N).
134*> \endverbatim
135*>
136*> \param[out] INFO
137*> \verbatim
138*> INFO is INTEGER
139*> = 0: successful exit
140*> < 0: if INFO = -k, the k-th argument had an illegal value
141*> \endverbatim
142*
143* Authors:
144* ========
145*
146*> \author Univ. of Tennessee
147*> \author Univ. of California Berkeley
148*> \author Univ. of Colorado Denver
149*> \author NAG Ltd.
150*
151*> \ingroup complex_lin
152*
153* =====================================================================
154 SUBROUTINE clavhe_rook( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
155 $ B, LDB, INFO )
156*
157* -- LAPACK test routine --
158* -- LAPACK is a software package provided by Univ. of Tennessee, --
159* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
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 A( LDA, * ), B( LDB, * )
168* ..
169*
170* =====================================================================
171*
172* .. Parameters ..
173 COMPLEX CONE
174 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
175* ..
176* .. Local Scalars ..
177 LOGICAL NOUNIT
178 INTEGER J, K, KP
179 COMPLEX D11, D12, D21, D22, T1, T2
180* ..
181* .. External Functions ..
182 LOGICAL LSAME
183 EXTERNAL lsame
184* ..
185* .. External Subroutines ..
186 EXTERNAL cgemv, cgeru, clacgv, cscal, cswap, xerbla
187* ..
188* .. Intrinsic Functions ..
189 INTRINSIC abs, conjg, 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( 'CLAVHE_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 cscal( 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 cgeru( 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 cswap( 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 = conjg( 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 cgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
291 $ ldb, b( 1, 1 ), ldb )
292 CALL cgeru( 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 cswap( 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 cswap( 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 cscal( 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 cgeru( 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 cswap( 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 = conjg( 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 cgeru( n-k, nrhs, cone, a( k+1, k ), 1,
384 $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
385 CALL cgeru( 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 cswap( 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 cswap( 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
440*
441* Apply the transformation
442* y = y - B' conjg(x),
443* where x is a column of A and y is a row of B.
444*
445 CALL clacgv( nrhs, b( k, 1 ), ldb )
446 CALL cgemv( 'Conjugate', k-1, nrhs, cone, b, ldb,
447 $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
448 CALL clacgv( nrhs, b( k, 1 ), ldb )
449 END IF
450 IF( nounit )
451 $ CALL cscal( 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 cswap( 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 cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
470 $ ldb )
471*
472* Apply the transformations
473* y = y - B' conjg(x),
474* where x is a block column of A and y is a block
475* row of B.
476*
477 CALL clacgv( nrhs, b( k, 1 ), ldb )
478 CALL cgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
479 $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
480 CALL clacgv( nrhs, b( k, 1 ), ldb )
481*
482 CALL clacgv( nrhs, b( k-1, 1 ), ldb )
483 CALL cgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
484 $ a( 1, k-1 ), 1, cone, b( k-1, 1 ), ldb )
485 CALL clacgv( 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 = conjg( 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
530*
531* Apply the transformation
532*
533 CALL clacgv( nrhs, b( k, 1 ), ldb )
534 CALL cgemv( 'Conjugate', n-k, nrhs, cone, b( k+1, 1 ),
535 $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
536 CALL clacgv( nrhs, b( k, 1 ), ldb )
537 END IF
538 IF( nounit )
539 $ CALL cscal( 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 cswap( 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 cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
558 $ ldb )
559*
560* Apply the transformation
561*
562 CALL clacgv( nrhs, b( k+1, 1 ), ldb )
563 CALL cgemv( '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 clacgv( nrhs, b( k+1, 1 ), ldb )
567*
568 CALL clacgv( nrhs, b( k, 1 ), ldb )
569 CALL cgemv( '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 clacgv( 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 = conjg( 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 CLAVHE_ROOK
599*
600 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
clavhe_rook
subroutine clavhe_rook(uplo, trans, diag, n, nrhs, a, lda, ipiv, b, ldb, info)
CLAVHE_ROOK
Definition clavhe_rook.f:156
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
cgeru
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:72
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81