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