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