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