OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
clavsy_rook.f
Go to the documentation of this file.
1*> \brief \b CLAVSY_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 CLAVSY_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*> CLAVSY_ROOK performs one of the matrix-vector operations
30*> x := A*x or x := A'*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 CSYTRF_ROOK.
33*>
34*> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
35*> If TRANS = 'T', 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*> = 'T': x := A'*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 CSYTRF_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[in] 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 CSYTRF_ROOK.
100*>
101*> If UPLO = 'U':
102*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
103*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
104*> (If IPIV( k ) = k, no interchange was done).
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*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
113*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
114*> (If IPIV( k ) = k, no interchange was done).
115*>
116*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
117*> columns k and -IPIV(k) were interchanged and rows and
118*> columns k+1 and -IPIV(k+1) were inerchaged,
119*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
120*> \endverbatim
121*>
122*> \param[in,out] B
123*> \verbatim
124*> B is COMPLEX array, dimension (LDB,NRHS)
125*> On entry, B contains NRHS vectors of length N.
126*> On exit, B is overwritten with the product A * B.
127*> \endverbatim
128*>
129*> \param[in] LDB
130*> \verbatim
131*> LDB is INTEGER
132*> The leading dimension of the array B. LDB >= max(1,N).
133*> \endverbatim
134*>
135*> \param[out] INFO
136*> \verbatim
137*> INFO is INTEGER
138*> = 0: successful exit
139*> < 0: if INFO = -k, the k-th argument had an illegal value
140*> \endverbatim
141*
142* Authors:
143* ========
144*
145*> \author Univ. of Tennessee
146*> \author Univ. of California Berkeley
147*> \author Univ. of Colorado Denver
148*> \author NAG Ltd.
149*
150*> \ingroup complex_lin
151*
152* =====================================================================
153 SUBROUTINE clavsy_rook( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
154 $ B, LDB, INFO )
155*
156* -- LAPACK test routine --
157* -- LAPACK is a software package provided by Univ. of Tennessee, --
158* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159*
160* .. Scalar Arguments ..
161 CHARACTER DIAG, TRANS, UPLO
162 INTEGER INFO, LDA, LDB, N, NRHS
163* ..
164* .. Array Arguments ..
165 INTEGER IPIV( * )
166 COMPLEX A( LDA, * ), B( LDB, * )
167* ..
168*
169* =====================================================================
170*
171* .. Parameters ..
172 COMPLEX CONE
173 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
174* ..
175* .. Local Scalars ..
176 LOGICAL NOUNIT
177 INTEGER J, K, KP
178 COMPLEX D11, D12, D21, D22, T1, T2
179* ..
180* .. External Functions ..
181 LOGICAL LSAME
182 EXTERNAL lsame
183* ..
184* .. External Subroutines ..
185 EXTERNAL cgemv, cgeru, cscal, cswap, xerbla
186* ..
187* .. Intrinsic Functions ..
188 INTRINSIC abs, max
189* ..
190* .. Executable Statements ..
191*
192* Test the input parameters.
193*
194 info = 0
195 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'l' ) ) THEN
196 INFO = -1
197.NOT. ELSE IF( LSAME( TRANS, 'n.AND..NOT.' ) LSAME( TRANS, 't' ) )
198 $ THEN
199 INFO = -2
200.NOT. ELSE IF( LSAME( DIAG, 'u.AND..NOT.' ) LSAME( DIAG, 'n' ) )
201 $ THEN
202 INFO = -3
203.LT. ELSE IF( N0 ) THEN
204 INFO = -4
205.LT. ELSE IF( LDAMAX( 1, N ) ) THEN
206 INFO = -6
207.LT. ELSE IF( LDBMAX( 1, N ) ) THEN
208 INFO = -9
209 END IF
210.NE. IF( INFO0 ) THEN
211 CALL XERBLA( 'clavsy_rook ', -INFO )
212 RETURN
213 END IF
214*
215* Quick return if possible.
216*
217.EQ. IF( N0 )
218 $ RETURN
219*
220 NOUNIT = LSAME( DIAG, 'n' )
221*------------------------------------------
222*
223* Compute B := A * B (No transpose)
224*
225*------------------------------------------
226 IF( LSAME( TRANS, 'n' ) ) THEN
227*
228* Compute B := U*B
229* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
230*
231 IF( LSAME( UPLO, 'u' ) ) THEN
232*
233* Loop forward applying the transformations.
234*
235 K = 1
236 10 CONTINUE
237.GT. IF( KN )
238 $ GO TO 30
239.GT. IF( IPIV( K )0 ) THEN
240*
241* 1 x 1 pivot block
242*
243* Multiply by the diagonal element if forming U * D.
244*
245 IF( NOUNIT )
246 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
247*
248* Multiply by P(K) * inv(U(K)) if K > 1.
249*
250.GT. IF( K1 ) THEN
251*
252* Apply the transformation.
253*
254 CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
255 $ LDB, B( 1, 1 ), LDB )
256*
257* Interchange if P(K) != I.
258*
259 KP = IPIV( K )
260.NE. IF( KPK )
261 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262 END IF
263 K = K + 1
264 ELSE
265*
266* 2 x 2 pivot block
267*
268* Multiply by the diagonal block if forming U * D.
269*
270 IF( NOUNIT ) THEN
271 D11 = A( K, K )
272 D22 = A( K+1, K+1 )
273 D12 = A( K, K+1 )
274 D21 = D12
275 DO 20 J = 1, NRHS
276 T1 = B( K, J )
277 T2 = B( K+1, J )
278 B( K, J ) = D11*T1 + D12*T2
279 B( K+1, J ) = D21*T1 + D22*T2
280 20 CONTINUE
281 END IF
282*
283* Multiply by P(K) * inv(U(K)) if K > 1.
284*
285.GT. IF( K1 ) THEN
286*
287* Apply the transformations.
288*
289 CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
290 $ LDB, B( 1, 1 ), LDB )
291 CALL CGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
292 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
293*
294* Interchange if a permutation was applied at the
295* K-th step of the factorization.
296*
297* Swap the first of pair with IMAXth
298*
299 KP = ABS( IPIV( K ) )
300.NE. IF( KPK )
301 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
302*
303* NOW swap the first of pair with Pth
304*
305 KP = ABS( IPIV( K+1 ) )
306.NE. IF( KPK+1 )
307 $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
308 $ LDB )
309 END IF
310 K = K + 2
311 END IF
312 GO TO 10
313 30 CONTINUE
314*
315* Compute B := L*B
316* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
317*
318 ELSE
319*
320* Loop backward applying the transformations to B.
321*
322 K = N
323 40 CONTINUE
324.LT. IF( K1 )
325 $ GO TO 60
326*
327* Test the pivot index. If greater than zero, a 1 x 1
328* pivot was used, otherwise a 2 x 2 pivot was used.
329*
330.GT. IF( IPIV( K )0 ) THEN
331*
332* 1 x 1 pivot block:
333*
334* Multiply by the diagonal element if forming L * D.
335*
336 IF( NOUNIT )
337 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
338*
339* Multiply by P(K) * inv(L(K)) if K < N.
340*
341.NE. IF( KN ) THEN
342 KP = IPIV( K )
343*
344* Apply the transformation.
345*
346 CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
347 $ B( K, 1 ), LDB, B( K+1, 1 ), LDB )
348*
349* Interchange if a permutation was applied at the
350* K-th step of the factorization.
351*
352.NE. IF( KPK )
353 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
354 END IF
355 K = K - 1
356*
357 ELSE
358*
359* 2 x 2 pivot block:
360*
361* Multiply by the diagonal block if forming L * D.
362*
363 IF( NOUNIT ) THEN
364 D11 = A( K-1, K-1 )
365 D22 = A( K, K )
366 D21 = A( K, K-1 )
367 D12 = D21
368 DO 50 J = 1, NRHS
369 T1 = B( K-1, J )
370 T2 = B( K, J )
371 B( K-1, J ) = D11*T1 + D12*T2
372 B( K, J ) = D21*T1 + D22*T2
373 50 CONTINUE
374 END IF
375*
376* Multiply by P(K) * inv(L(K)) if K < N.
377*
378.NE. IF( KN ) THEN
379*
380* Apply the transformation.
381*
382 CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
383 $ B( K, 1 ), LDB, B( K+1, 1 ), LDB )
384 CALL CGERU( N-K, NRHS, CONE, A( K+1, K-1 ), 1,
385 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
386*
387* Interchange if a permutation was applied at the
388* K-th step of the factorization.
389*
390* Swap the second of pair with IMAXth
391*
392 KP = ABS( IPIV( K ) )
393.NE. IF( KPK )
394 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
395*
396* NOW swap the first of pair with Pth
397*
398 KP = ABS( IPIV( K-1 ) )
399.NE. IF( KPK-1 )
400 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
401 $ LDB )
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' * B (transpose)
411*
412*----------------------------------------
413 ELSE IF( LSAME( TRANS, 't' ) ) THEN
414*
415* Form B := U'*B
416* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
417* and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
418*
419 IF( LSAME( UPLO, 'u' ) ) THEN
420*
421* Loop backward applying the transformations.
422*
423 K = N
424.LT. 70 IF( K1 )
425 $ GO TO 90
426*
427* 1 x 1 pivot block.
428*
429.GT. IF( IPIV( K )0 ) THEN
430.GT. IF( K1 ) THEN
431*
432* Interchange if P(K) != I.
433*
434 KP = IPIV( K )
435.NE. IF( KPK )
436 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
437*
438* Apply the transformation
439*
440 CALL CGEMV( 'transpose', K-1, NRHS, CONE, B, LDB,
441 $ A( 1, K ), 1, CONE, B( K, 1 ), LDB )
442 END IF
443 IF( NOUNIT )
444 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
445 K = K - 1
446*
447* 2 x 2 pivot block.
448*
449 ELSE
450.GT. IF( K2 ) THEN
451*
452* Swap the second of pair with Pth
453*
454 KP = ABS( IPIV( K ) )
455.NE. IF( KPK )
456 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
457*
458* Now swap the first of pair with IMAX(r)th
459*
460 KP = ABS( IPIV( K-1 ) )
461.NE. IF( KPK-1 )
462 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
463 $ LDB )
464*
465* Apply the transformations
466*
467 CALL CGEMV( 'transpose', K-2, NRHS, CONE, B, LDB,
468 $ A( 1, K ), 1, CONE, B( K, 1 ), LDB )
469 CALL CGEMV( 'transpose', K-2, NRHS, CONE, B, LDB,
470 $ A( 1, K-1 ), 1, CONE, B( K-1, 1 ), LDB )
471 END IF
472*
473* Multiply by the diagonal block if non-unit.
474*
475 IF( NOUNIT ) THEN
476 D11 = A( K-1, K-1 )
477 D22 = A( K, K )
478 D12 = A( K-1, K )
479 D21 = D12
480 DO 80 J = 1, NRHS
481 T1 = B( K-1, J )
482 T2 = B( K, J )
483 B( K-1, J ) = D11*T1 + D12*T2
484 B( K, J ) = D21*T1 + D22*T2
485 80 CONTINUE
486 END IF
487 K = K - 2
488 END IF
489 GO TO 70
490 90 CONTINUE
491*
492* Form B := L'*B
493* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
494* and L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
495*
496 ELSE
497*
498* Loop forward applying the L-transformations.
499*
500 K = 1
501 100 CONTINUE
502.GT. IF( KN )
503 $ GO TO 120
504*
505* 1 x 1 pivot block
506*
507.GT. IF( IPIV( K )0 ) THEN
508.LT. IF( KN ) THEN
509*
510* Interchange if P(K) != I.
511*
512 KP = IPIV( K )
513.NE. IF( KPK )
514 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
515*
516* Apply the transformation
517*
518 CALL CGEMV( 'transpose', n-k, nrhs, cone, b( k+1, 1 ),
519 $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
520 END IF
521 IF( nounit )
522 $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
523 k = k + 1
524*
525* 2 x 2 pivot block.
526*
527 ELSE
528 IF( k.LT.n-1 ) THEN
529*
530* Swap the first of pair with Pth
531*
532 kp = abs( ipiv( k ) )
533 IF( kp.NE.k )
534 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
535*
536* Now swap the second of pair with IMAX(r)th
537*
538 kp = abs( ipiv( k+1 ) )
539 IF( kp.NE.k+1 )
540 $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
541 $ ldb )
542*
543* Apply the transformation
544*
545 CALL cgemv( 'Transpose', n-k-1, nrhs, cone,
546 $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, cone,
547 $ b( k+1, 1 ), ldb )
548 CALL cgemv( 'Transpose', n-k-1, nrhs, cone,
549 $ b( k+2, 1 ), ldb, a( k+2, k ), 1, cone,
550 $ b( k, 1 ), ldb )
551 END IF
552*
553* Multiply by the diagonal block if non-unit.
554*
555 IF( nounit ) THEN
556 d11 = a( k, k )
557 d22 = a( k+1, k+1 )
558 d21 = a( k+1, k )
559 d12 = d21
560 DO 110 j = 1, nrhs
561 t1 = b( k, j )
562 t2 = b( k+1, j )
563 b( k, j ) = d11*t1 + d12*t2
564 b( k+1, j ) = d21*t1 + d22*t2
565 110 CONTINUE
566 END IF
567 k = k + 2
568 END IF
569 GO TO 100
570 120 CONTINUE
571 END IF
572 END IF
573 RETURN
574*
575* End of CLAVSY_ROOK
576*
577 END
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
cgeru
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
clavsy_rook
subroutine clavsy_rook(uplo, trans, diag, n, nrhs, a, lda, ipiv, b, ldb, info)
CLAVSY_ROOK
Definition clavsy_rook.f:155
max
#define max(a, b)
Definition macros.h:21