LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ stfsm()

subroutine stfsm ( character transr,
character side,
character uplo,
character trans,
character diag,
integer m,
integer n,
real alpha,
real, dimension( 0: * ) a,
real, dimension( 0: ldb-1, 0: * ) b,
integer ldb )

STFSM solves a matrix equation (one operand is a triangular matrix in RFP format).

Download STFSM + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> Level 3 BLAS like routine for A in RFP Format.
!>
!> STFSM  solves the matrix equation
!>
!>    op( A )*X = alpha*B  or  X*op( A ) = alpha*B
!>
!> where alpha is a scalar, X and B are m by n matrices, A is a unit, or
!> non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
!>
!>    op( A ) = A   or   op( A ) = A**T.
!>
!> A is in Rectangular Full Packed (RFP) Format.
!>
!> The matrix X is overwritten on B.
!>
Parameters
[in]TRANSR
!>          TRANSR is CHARACTER*1
!>          = 'N':  The Normal Form of RFP A is stored;
!>          = 'T':  The Transpose Form of RFP A is stored.
!>
[in]SIDE
!>          SIDE is CHARACTER*1
!>           On entry, SIDE specifies whether op( A ) appears on the left
!>           or right of X as follows:
!>
!>              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
!>
!>              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
!>
!>           Unchanged on exit.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>           On entry, UPLO specifies whether the RFP matrix A came from
!>           an upper or lower triangular matrix as follows:
!>           UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
!>           UPLO = 'L' or 'l' RFP A came from a  lower triangular matrix
!>
!>           Unchanged on exit.
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>           On entry, TRANS  specifies the form of op( A ) to be used
!>           in the matrix multiplication as follows:
!>
!>              TRANS  = 'N' or 'n'   op( A ) = A.
!>
!>              TRANS  = 'T' or 't'   op( A ) = A'.
!>
!>           Unchanged on exit.
!>
[in]DIAG
!>          DIAG is CHARACTER*1
!>           On entry, DIAG specifies whether or not RFP A is unit
!>           triangular as follows:
!>
!>              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!>
!>              DIAG = 'N' or 'n'   A is not assumed to be unit
!>                                  triangular.
!>
!>           Unchanged on exit.
!>
[in]M
!>          M is INTEGER
!>           On entry, M specifies the number of rows of B. M must be at
!>           least zero.
!>           Unchanged on exit.
!>
[in]N
!>          N is INTEGER
!>           On entry, N specifies the number of columns of B.  N must be
!>           at least zero.
!>           Unchanged on exit.
!>
[in]ALPHA
!>          ALPHA is REAL
!>           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
!>           zero then  A is not referenced and  B need not be set before
!>           entry.
!>           Unchanged on exit.
!>
[in]A
!>          A is REAL array, dimension (NT)
!>           NT = N*(N+1)/2 if SIDE='R' and NT = M*(M+1)/2 otherwise.
!>           On entry, the matrix A in RFP Format.
!>           RFP Format is described by TRANSR, UPLO and N as follows:
!>           If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
!>           K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If
!>           TRANSR = 'T' then RFP is the transpose of RFP A as
!>           defined when TRANSR = 'N'. The contents of RFP A are defined
!>           by UPLO as follows: If UPLO = 'U' the RFP A contains the NT
!>           elements of upper packed A either in normal or
!>           transpose Format. If UPLO = 'L' the RFP A contains
!>           the NT elements of lower packed A either in normal or
!>           transpose Format. The LDA of RFP A is (N+1)/2 when
!>           TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
!>           even and is N when is odd.
!>           See the Note below for more details. Unchanged on exit.
!>
[in,out]B
!>          B is REAL array, dimension (LDB,N)
!>           Before entry,  the leading  m by n part of the array  B must
!>           contain  the  right-hand  side  matrix  B,  and  on exit  is
!>           overwritten by the solution matrix  X.
!>
[in]LDB
!>          LDB is INTEGER
!>           On entry, LDB specifies the first dimension of B as declared
!>           in  the  calling  (sub)  program.   LDB  must  be  at  least
!>           max( 1, m ).
!>           Unchanged on exit.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  We first consider Rectangular Full Packed (RFP) Format when N is
!>  even. We give an example where N = 6.
!>
!>      AP is Upper             AP is Lower
!>
!>   00 01 02 03 04 05       00
!>      11 12 13 14 15       10 11
!>         22 23 24 25       20 21 22
!>            33 34 35       30 31 32 33
!>               44 45       40 41 42 43 44
!>                  55       50 51 52 53 54 55
!>
!>
!>  Let TRANSR = 'N'. RFP holds AP as follows:
!>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
!>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
!>  the transpose of the first three columns of AP upper.
!>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
!>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
!>  the transpose of the last three columns of AP lower.
!>  This covers the case N even and TRANSR = 'N'.
!>
!>         RFP A                   RFP A
!>
!>        03 04 05                33 43 53
!>        13 14 15                00 44 54
!>        23 24 25                10 11 55
!>        33 34 35                20 21 22
!>        00 44 45                30 31 32
!>        01 11 55                40 41 42
!>        02 12 22                50 51 52
!>
!>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
!>  transpose of RFP A above. One therefore gets:
!>
!>
!>           RFP A                   RFP A
!>
!>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
!>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
!>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
!>
!>
!>  We then consider Rectangular Full Packed (RFP) Format when N is
!>  odd. We give an example where N = 5.
!>
!>     AP is Upper                 AP is Lower
!>
!>   00 01 02 03 04              00
!>      11 12 13 14              10 11
!>         22 23 24              20 21 22
!>            33 34              30 31 32 33
!>               44              40 41 42 43 44
!>
!>
!>  Let TRANSR = 'N'. RFP holds AP as follows:
!>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
!>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
!>  the transpose of the first two columns of AP upper.
!>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
!>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
!>  the transpose of the last two columns of AP lower.
!>  This covers the case N odd and TRANSR = 'N'.
!>
!>         RFP A                   RFP A
!>
!>        02 03 04                00 33 43
!>        12 13 14                10 11 44
!>        22 23 24                20 21 22
!>        00 33 34                30 31 32
!>        01 11 44                40 41 42
!>
!>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
!>  transpose of RFP A above. One therefore gets:
!>
!>           RFP A                   RFP A
!>
!>     02 12 22 00 01             00 10 20 30 40 50
!>     03 13 23 33 11             33 11 21 31 41 51
!>     04 14 24 34 44             43 44 22 32 42 52
!>

Definition at line 274 of file stfsm.f.

277*
278* -- LAPACK computational routine --
279* -- LAPACK is a software package provided by Univ. of Tennessee, --
280* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
281*
282* .. Scalar Arguments ..
283 CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
284 INTEGER LDB, M, N
285 REAL ALPHA
286* ..
287* .. Array Arguments ..
288 REAL A( 0: * ), B( 0: LDB-1, 0: * )
289* ..
290*
291* =====================================================================
292*
293* ..
294* .. Parameters ..
295 REAL ONE, ZERO
296 parameter( one = 1.0e+0, zero = 0.0e+0 )
297* ..
298* .. Local Scalars ..
299 LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
300 $ NOTRANS
301 INTEGER M1, M2, N1, N2, K, INFO, I, J
302* ..
303* .. External Functions ..
304 LOGICAL LSAME
305 EXTERNAL lsame
306* ..
307* .. External Subroutines ..
308 EXTERNAL sgemm, strsm, xerbla
309* ..
310* .. Intrinsic Functions ..
311 INTRINSIC max, mod
312* ..
313* .. Executable Statements ..
314*
315* Test the input parameters.
316*
317 info = 0
318 normaltransr = lsame( transr, 'N' )
319 lside = lsame( side, 'L' )
320 lower = lsame( uplo, 'L' )
321 notrans = lsame( trans, 'N' )
322 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
323 info = -1
324 ELSE IF( .NOT.lside .AND. .NOT.lsame( side, 'R' ) ) THEN
325 info = -2
326 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
327 info = -3
328 ELSE IF( .NOT.notrans .AND. .NOT.lsame( trans, 'T' ) ) THEN
329 info = -4
330 ELSE IF( .NOT.lsame( diag, 'N' ) .AND.
331 $ .NOT.lsame( diag, 'U' ) )
332 $ THEN
333 info = -5
334 ELSE IF( m.LT.0 ) THEN
335 info = -6
336 ELSE IF( n.LT.0 ) THEN
337 info = -7
338 ELSE IF( ldb.LT.max( 1, m ) ) THEN
339 info = -11
340 END IF
341 IF( info.NE.0 ) THEN
342 CALL xerbla( 'STFSM ', -info )
343 RETURN
344 END IF
345*
346* Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
347*
348 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )
349 $ RETURN
350*
351* Quick return when ALPHA.EQ.(0D+0)
352*
353 IF( alpha.EQ.zero ) THEN
354 DO 20 j = 0, n - 1
355 DO 10 i = 0, m - 1
356 b( i, j ) = zero
357 10 CONTINUE
358 20 CONTINUE
359 RETURN
360 END IF
361*
362 IF( lside ) THEN
363*
364* SIDE = 'L'
365*
366* A is M-by-M.
367* If M is odd, set NISODD = .TRUE., and M1 and M2.
368* If M is even, NISODD = .FALSE., and M.
369*
370 IF( mod( m, 2 ).EQ.0 ) THEN
371 misodd = .false.
372 k = m / 2
373 ELSE
374 misodd = .true.
375 IF( lower ) THEN
376 m2 = m / 2
377 m1 = m - m2
378 ELSE
379 m1 = m / 2
380 m2 = m - m1
381 END IF
382 END IF
383*
384 IF( misodd ) THEN
385*
386* SIDE = 'L' and N is odd
387*
388 IF( normaltransr ) THEN
389*
390* SIDE = 'L', N is odd, and TRANSR = 'N'
391*
392 IF( lower ) THEN
393*
394* SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
395*
396 IF( notrans ) THEN
397*
398* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
399* TRANS = 'N'
400*
401 IF( m.EQ.1 ) THEN
402 CALL strsm( 'L', 'L', 'N', diag, m1, n,
403 $ alpha,
404 $ a, m, b, ldb )
405 ELSE
406 CALL strsm( 'L', 'L', 'N', diag, m1, n,
407 $ alpha,
408 $ a( 0 ), m, b, ldb )
409 CALL sgemm( 'N', 'N', m2, n, m1, -one,
410 $ a( m1 ),
411 $ m, b, ldb, alpha, b( m1, 0 ), ldb )
412 CALL strsm( 'L', 'U', 'T', diag, m2, n, one,
413 $ a( m ), m, b( m1, 0 ), ldb )
414 END IF
415*
416 ELSE
417*
418* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
419* TRANS = 'T'
420*
421 IF( m.EQ.1 ) THEN
422 CALL strsm( 'L', 'L', 'T', diag, m1, n,
423 $ alpha,
424 $ a( 0 ), m, b, ldb )
425 ELSE
426 CALL strsm( 'L', 'U', 'N', diag, m2, n,
427 $ alpha,
428 $ a( m ), m, b( m1, 0 ), ldb )
429 CALL sgemm( 'T', 'N', m1, n, m2, -one,
430 $ a( m1 ),
431 $ m, b( m1, 0 ), ldb, alpha, b, ldb )
432 CALL strsm( 'L', 'L', 'T', diag, m1, n, one,
433 $ a( 0 ), m, b, ldb )
434 END IF
435*
436 END IF
437*
438 ELSE
439*
440* SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
441*
442 IF( .NOT.notrans ) THEN
443*
444* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
445* TRANS = 'N'
446*
447 CALL strsm( 'L', 'L', 'N', diag, m1, n, alpha,
448 $ a( m2 ), m, b, ldb )
449 CALL sgemm( 'T', 'N', m2, n, m1, -one, a( 0 ),
450 $ m,
451 $ b, ldb, alpha, b( m1, 0 ), ldb )
452 CALL strsm( 'L', 'U', 'T', diag, m2, n, one,
453 $ a( m1 ), m, b( m1, 0 ), ldb )
454*
455 ELSE
456*
457* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
458* TRANS = 'T'
459*
460 CALL strsm( 'L', 'U', 'N', diag, m2, n, alpha,
461 $ a( m1 ), m, b( m1, 0 ), ldb )
462 CALL sgemm( 'N', 'N', m1, n, m2, -one, a( 0 ),
463 $ m,
464 $ b( m1, 0 ), ldb, alpha, b, ldb )
465 CALL strsm( 'L', 'L', 'T', diag, m1, n, one,
466 $ a( m2 ), m, b, ldb )
467*
468 END IF
469*
470 END IF
471*
472 ELSE
473*
474* SIDE = 'L', N is odd, and TRANSR = 'T'
475*
476 IF( lower ) THEN
477*
478* SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
479*
480 IF( notrans ) THEN
481*
482* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
483* TRANS = 'N'
484*
485 IF( m.EQ.1 ) THEN
486 CALL strsm( 'L', 'U', 'T', diag, m1, n,
487 $ alpha,
488 $ a( 0 ), m1, b, ldb )
489 ELSE
490 CALL strsm( 'L', 'U', 'T', diag, m1, n,
491 $ alpha,
492 $ a( 0 ), m1, b, ldb )
493 CALL sgemm( 'T', 'N', m2, n, m1, -one,
494 $ a( m1*m1 ), m1, b, ldb, alpha,
495 $ b( m1, 0 ), ldb )
496 CALL strsm( 'L', 'L', 'N', diag, m2, n, one,
497 $ a( 1 ), m1, b( m1, 0 ), ldb )
498 END IF
499*
500 ELSE
501*
502* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
503* TRANS = 'T'
504*
505 IF( m.EQ.1 ) THEN
506 CALL strsm( 'L', 'U', 'N', diag, m1, n,
507 $ alpha,
508 $ a( 0 ), m1, b, ldb )
509 ELSE
510 CALL strsm( 'L', 'L', 'T', diag, m2, n,
511 $ alpha,
512 $ a( 1 ), m1, b( m1, 0 ), ldb )
513 CALL sgemm( 'N', 'N', m1, n, m2, -one,
514 $ a( m1*m1 ), m1, b( m1, 0 ), ldb,
515 $ alpha, b, ldb )
516 CALL strsm( 'L', 'U', 'N', diag, m1, n, one,
517 $ a( 0 ), m1, b, ldb )
518 END IF
519*
520 END IF
521*
522 ELSE
523*
524* SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
525*
526 IF( .NOT.notrans ) THEN
527*
528* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
529* TRANS = 'N'
530*
531 CALL strsm( 'L', 'U', 'T', diag, m1, n, alpha,
532 $ a( m2*m2 ), m2, b, ldb )
533 CALL sgemm( 'N', 'N', m2, n, m1, -one, a( 0 ),
534 $ m2,
535 $ b, ldb, alpha, b( m1, 0 ), ldb )
536 CALL strsm( 'L', 'L', 'N', diag, m2, n, one,
537 $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
538*
539 ELSE
540*
541* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
542* TRANS = 'T'
543*
544 CALL strsm( 'L', 'L', 'T', diag, m2, n, alpha,
545 $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
546 CALL sgemm( 'T', 'N', m1, n, m2, -one, a( 0 ),
547 $ m2,
548 $ b( m1, 0 ), ldb, alpha, b, ldb )
549 CALL strsm( 'L', 'U', 'N', diag, m1, n, one,
550 $ a( m2*m2 ), m2, b, ldb )
551*
552 END IF
553*
554 END IF
555*
556 END IF
557*
558 ELSE
559*
560* SIDE = 'L' and N is even
561*
562 IF( normaltransr ) THEN
563*
564* SIDE = 'L', N is even, and TRANSR = 'N'
565*
566 IF( lower ) THEN
567*
568* SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
569*
570 IF( notrans ) THEN
571*
572* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
573* and TRANS = 'N'
574*
575 CALL strsm( 'L', 'L', 'N', diag, k, n, alpha,
576 $ a( 1 ), m+1, b, ldb )
577 CALL sgemm( 'N', 'N', k, n, k, -one, a( k+1 ),
578 $ m+1, b, ldb, alpha, b( k, 0 ), ldb )
579 CALL strsm( 'L', 'U', 'T', diag, k, n, one,
580 $ a( 0 ), m+1, b( k, 0 ), ldb )
581*
582 ELSE
583*
584* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
585* and TRANS = 'T'
586*
587 CALL strsm( 'L', 'U', 'N', diag, k, n, alpha,
588 $ a( 0 ), m+1, b( k, 0 ), ldb )
589 CALL sgemm( 'T', 'N', k, n, k, -one, a( k+1 ),
590 $ m+1, b( k, 0 ), ldb, alpha, b, ldb )
591 CALL strsm( 'L', 'L', 'T', diag, k, n, one,
592 $ a( 1 ), m+1, b, ldb )
593*
594 END IF
595*
596 ELSE
597*
598* SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
599*
600 IF( .NOT.notrans ) THEN
601*
602* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
603* and TRANS = 'N'
604*
605 CALL strsm( 'L', 'L', 'N', diag, k, n, alpha,
606 $ a( k+1 ), m+1, b, ldb )
607 CALL sgemm( 'T', 'N', k, n, k, -one, a( 0 ),
608 $ m+1,
609 $ b, ldb, alpha, b( k, 0 ), ldb )
610 CALL strsm( 'L', 'U', 'T', diag, k, n, one,
611 $ a( k ), m+1, b( k, 0 ), ldb )
612*
613 ELSE
614*
615* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
616* and TRANS = 'T'
617 CALL strsm( 'L', 'U', 'N', diag, k, n, alpha,
618 $ a( k ), m+1, b( k, 0 ), ldb )
619 CALL sgemm( 'N', 'N', k, n, k, -one, a( 0 ),
620 $ m+1,
621 $ b( k, 0 ), ldb, alpha, b, ldb )
622 CALL strsm( 'L', 'L', 'T', diag, k, n, one,
623 $ a( k+1 ), m+1, b, ldb )
624*
625 END IF
626*
627 END IF
628*
629 ELSE
630*
631* SIDE = 'L', N is even, and TRANSR = 'T'
632*
633 IF( lower ) THEN
634*
635* SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
636*
637 IF( notrans ) THEN
638*
639* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
640* and TRANS = 'N'
641*
642 CALL strsm( 'L', 'U', 'T', diag, k, n, alpha,
643 $ a( k ), k, b, ldb )
644 CALL sgemm( 'T', 'N', k, n, k, -one,
645 $ a( k*( k+1 ) ), k, b, ldb, alpha,
646 $ b( k, 0 ), ldb )
647 CALL strsm( 'L', 'L', 'N', diag, k, n, one,
648 $ a( 0 ), k, b( k, 0 ), ldb )
649*
650 ELSE
651*
652* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
653* and TRANS = 'T'
654*
655 CALL strsm( 'L', 'L', 'T', diag, k, n, alpha,
656 $ a( 0 ), k, b( k, 0 ), ldb )
657 CALL sgemm( 'N', 'N', k, n, k, -one,
658 $ a( k*( k+1 ) ), k, b( k, 0 ), ldb,
659 $ alpha, b, ldb )
660 CALL strsm( 'L', 'U', 'N', diag, k, n, one,
661 $ a( k ), k, b, ldb )
662*
663 END IF
664*
665 ELSE
666*
667* SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
668*
669 IF( .NOT.notrans ) THEN
670*
671* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
672* and TRANS = 'N'
673*
674 CALL strsm( 'L', 'U', 'T', diag, k, n, alpha,
675 $ a( k*( k+1 ) ), k, b, ldb )
676 CALL sgemm( 'N', 'N', k, n, k, -one, a( 0 ), k,
677 $ b,
678 $ ldb, alpha, b( k, 0 ), ldb )
679 CALL strsm( 'L', 'L', 'N', diag, k, n, one,
680 $ a( k*k ), k, b( k, 0 ), ldb )
681*
682 ELSE
683*
684* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
685* and TRANS = 'T'
686*
687 CALL strsm( 'L', 'L', 'T', diag, k, n, alpha,
688 $ a( k*k ), k, b( k, 0 ), ldb )
689 CALL sgemm( 'T', 'N', k, n, k, -one, a( 0 ), k,
690 $ b( k, 0 ), ldb, alpha, b, ldb )
691 CALL strsm( 'L', 'U', 'N', diag, k, n, one,
692 $ a( k*( k+1 ) ), k, b, ldb )
693*
694 END IF
695*
696 END IF
697*
698 END IF
699*
700 END IF
701*
702 ELSE
703*
704* SIDE = 'R'
705*
706* A is N-by-N.
707* If N is odd, set NISODD = .TRUE., and N1 and N2.
708* If N is even, NISODD = .FALSE., and K.
709*
710 IF( mod( n, 2 ).EQ.0 ) THEN
711 nisodd = .false.
712 k = n / 2
713 ELSE
714 nisodd = .true.
715 IF( lower ) THEN
716 n2 = n / 2
717 n1 = n - n2
718 ELSE
719 n1 = n / 2
720 n2 = n - n1
721 END IF
722 END IF
723*
724 IF( nisodd ) THEN
725*
726* SIDE = 'R' and N is odd
727*
728 IF( normaltransr ) THEN
729*
730* SIDE = 'R', N is odd, and TRANSR = 'N'
731*
732 IF( lower ) THEN
733*
734* SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
735*
736 IF( notrans ) THEN
737*
738* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
739* TRANS = 'N'
740*
741 CALL strsm( 'R', 'U', 'T', diag, m, n2, alpha,
742 $ a( n ), n, b( 0, n1 ), ldb )
743 CALL sgemm( 'N', 'N', m, n1, n2, -one, b( 0,
744 $ n1 ),
745 $ ldb, a( n1 ), n, alpha, b( 0, 0 ),
746 $ ldb )
747 CALL strsm( 'R', 'L', 'N', diag, m, n1, one,
748 $ a( 0 ), n, b( 0, 0 ), ldb )
749*
750 ELSE
751*
752* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
753* TRANS = 'T'
754*
755 CALL strsm( 'R', 'L', 'T', diag, m, n1, alpha,
756 $ a( 0 ), n, b( 0, 0 ), ldb )
757 CALL sgemm( 'N', 'T', m, n2, n1, -one, b( 0,
758 $ 0 ),
759 $ ldb, a( n1 ), n, alpha, b( 0, n1 ),
760 $ ldb )
761 CALL strsm( 'R', 'U', 'N', diag, m, n2, one,
762 $ a( n ), n, b( 0, n1 ), ldb )
763*
764 END IF
765*
766 ELSE
767*
768* SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
769*
770 IF( notrans ) THEN
771*
772* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
773* TRANS = 'N'
774*
775 CALL strsm( 'R', 'L', 'T', diag, m, n1, alpha,
776 $ a( n2 ), n, b( 0, 0 ), ldb )
777 CALL sgemm( 'N', 'N', m, n2, n1, -one, b( 0,
778 $ 0 ),
779 $ ldb, a( 0 ), n, alpha, b( 0, n1 ),
780 $ ldb )
781 CALL strsm( 'R', 'U', 'N', diag, m, n2, one,
782 $ a( n1 ), n, b( 0, n1 ), ldb )
783*
784 ELSE
785*
786* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
787* TRANS = 'T'
788*
789 CALL strsm( 'R', 'U', 'T', diag, m, n2, alpha,
790 $ a( n1 ), n, b( 0, n1 ), ldb )
791 CALL sgemm( 'N', 'T', m, n1, n2, -one, b( 0,
792 $ n1 ),
793 $ ldb, a( 0 ), n, alpha, b( 0, 0 ), ldb )
794 CALL strsm( 'R', 'L', 'N', diag, m, n1, one,
795 $ a( n2 ), n, b( 0, 0 ), ldb )
796*
797 END IF
798*
799 END IF
800*
801 ELSE
802*
803* SIDE = 'R', N is odd, and TRANSR = 'T'
804*
805 IF( lower ) THEN
806*
807* SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
808*
809 IF( notrans ) THEN
810*
811* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
812* TRANS = 'N'
813*
814 CALL strsm( 'R', 'L', 'N', diag, m, n2, alpha,
815 $ a( 1 ), n1, b( 0, n1 ), ldb )
816 CALL sgemm( 'N', 'T', m, n1, n2, -one, b( 0,
817 $ n1 ),
818 $ ldb, a( n1*n1 ), n1, alpha, b( 0, 0 ),
819 $ ldb )
820 CALL strsm( 'R', 'U', 'T', diag, m, n1, one,
821 $ a( 0 ), n1, b( 0, 0 ), ldb )
822*
823 ELSE
824*
825* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
826* TRANS = 'T'
827*
828 CALL strsm( 'R', 'U', 'N', diag, m, n1, alpha,
829 $ a( 0 ), n1, b( 0, 0 ), ldb )
830 CALL sgemm( 'N', 'N', m, n2, n1, -one, b( 0,
831 $ 0 ),
832 $ ldb, a( n1*n1 ), n1, alpha, b( 0, n1 ),
833 $ ldb )
834 CALL strsm( 'R', 'L', 'T', diag, m, n2, one,
835 $ a( 1 ), n1, b( 0, n1 ), ldb )
836*
837 END IF
838*
839 ELSE
840*
841* SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
842*
843 IF( notrans ) THEN
844*
845* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
846* TRANS = 'N'
847*
848 CALL strsm( 'R', 'U', 'N', diag, m, n1, alpha,
849 $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
850 CALL sgemm( 'N', 'T', m, n2, n1, -one, b( 0,
851 $ 0 ),
852 $ ldb, a( 0 ), n2, alpha, b( 0, n1 ),
853 $ ldb )
854 CALL strsm( 'R', 'L', 'T', diag, m, n2, one,
855 $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
856*
857 ELSE
858*
859* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
860* TRANS = 'T'
861*
862 CALL strsm( 'R', 'L', 'N', diag, m, n2, alpha,
863 $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
864 CALL sgemm( 'N', 'N', m, n1, n2, -one, b( 0,
865 $ n1 ),
866 $ ldb, a( 0 ), n2, alpha, b( 0, 0 ),
867 $ ldb )
868 CALL strsm( 'R', 'U', 'T', diag, m, n1, one,
869 $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
870*
871 END IF
872*
873 END IF
874*
875 END IF
876*
877 ELSE
878*
879* SIDE = 'R' and N is even
880*
881 IF( normaltransr ) THEN
882*
883* SIDE = 'R', N is even, and TRANSR = 'N'
884*
885 IF( lower ) THEN
886*
887* SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
888*
889 IF( notrans ) THEN
890*
891* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
892* and TRANS = 'N'
893*
894 CALL strsm( 'R', 'U', 'T', diag, m, k, alpha,
895 $ a( 0 ), n+1, b( 0, k ), ldb )
896 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
897 $ ldb, a( k+1 ), n+1, alpha, b( 0, 0 ),
898 $ ldb )
899 CALL strsm( 'R', 'L', 'N', diag, m, k, one,
900 $ a( 1 ), n+1, b( 0, 0 ), ldb )
901*
902 ELSE
903*
904* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
905* and TRANS = 'T'
906*
907 CALL strsm( 'R', 'L', 'T', diag, m, k, alpha,
908 $ a( 1 ), n+1, b( 0, 0 ), ldb )
909 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
910 $ ldb, a( k+1 ), n+1, alpha, b( 0, k ),
911 $ ldb )
912 CALL strsm( 'R', 'U', 'N', diag, m, k, one,
913 $ a( 0 ), n+1, b( 0, k ), ldb )
914*
915 END IF
916*
917 ELSE
918*
919* SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
920*
921 IF( notrans ) THEN
922*
923* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
924* and TRANS = 'N'
925*
926 CALL strsm( 'R', 'L', 'T', diag, m, k, alpha,
927 $ a( k+1 ), n+1, b( 0, 0 ), ldb )
928 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
929 $ ldb, a( 0 ), n+1, alpha, b( 0, k ),
930 $ ldb )
931 CALL strsm( 'R', 'U', 'N', diag, m, k, one,
932 $ a( k ), n+1, b( 0, k ), ldb )
933*
934 ELSE
935*
936* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
937* and TRANS = 'T'
938*
939 CALL strsm( 'R', 'U', 'T', diag, m, k, alpha,
940 $ a( k ), n+1, b( 0, k ), ldb )
941 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
942 $ ldb, a( 0 ), n+1, alpha, b( 0, 0 ),
943 $ ldb )
944 CALL strsm( 'R', 'L', 'N', diag, m, k, one,
945 $ a( k+1 ), n+1, b( 0, 0 ), ldb )
946*
947 END IF
948*
949 END IF
950*
951 ELSE
952*
953* SIDE = 'R', N is even, and TRANSR = 'T'
954*
955 IF( lower ) THEN
956*
957* SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
958*
959 IF( notrans ) THEN
960*
961* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
962* and TRANS = 'N'
963*
964 CALL strsm( 'R', 'L', 'N', diag, m, k, alpha,
965 $ a( 0 ), k, b( 0, k ), ldb )
966 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
967 $ ldb, a( ( k+1 )*k ), k, alpha,
968 $ b( 0, 0 ), ldb )
969 CALL strsm( 'R', 'U', 'T', diag, m, k, one,
970 $ a( k ), k, b( 0, 0 ), ldb )
971*
972 ELSE
973*
974* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
975* and TRANS = 'T'
976*
977 CALL strsm( 'R', 'U', 'N', diag, m, k, alpha,
978 $ a( k ), k, b( 0, 0 ), ldb )
979 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
980 $ ldb, a( ( k+1 )*k ), k, alpha,
981 $ b( 0, k ), ldb )
982 CALL strsm( 'R', 'L', 'T', diag, m, k, one,
983 $ a( 0 ), k, b( 0, k ), ldb )
984*
985 END IF
986*
987 ELSE
988*
989* SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
990*
991 IF( notrans ) THEN
992*
993* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
994* and TRANS = 'N'
995*
996 CALL strsm( 'R', 'U', 'N', diag, m, k, alpha,
997 $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
998 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
999 $ ldb, a( 0 ), k, alpha, b( 0, k ), ldb )
1000 CALL strsm( 'R', 'L', 'T', diag, m, k, one,
1001 $ a( k*k ), k, b( 0, k ), ldb )
1002*
1003 ELSE
1004*
1005* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
1006* and TRANS = 'T'
1007*
1008 CALL strsm( 'R', 'L', 'N', diag, m, k, alpha,
1009 $ a( k*k ), k, b( 0, k ), ldb )
1010 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
1011 $ ldb, a( 0 ), k, alpha, b( 0, 0 ), ldb )
1012 CALL strsm( 'R', 'U', 'T', diag, m, k, one,
1013 $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
1014*
1015 END IF
1016*
1017 END IF
1018*
1019 END IF
1020*
1021 END IF
1022 END IF
1023*
1024 RETURN
1025*
1026* End of STFSM
1027*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
Here is the call graph for this function:
Here is the caller graph for this function: