5649
5650
5651
5652
5653
5654
5655
5656 CHARACTER*1 TRANS, UPLO
5657 INTEGER IA, IC, ICTXT, INFO, JA, JC, K, N
5658 REAL ALPHA, BETA, ERR
5659
5660
5661 INTEGER DESCA( * ), DESCC( * )
5662 REAL A( * ), C( * ), CT( * ), G( * ), PC( * )
5663
5664
5665
5666
5667
5668
5669
5670
5671
5672
5673
5674
5675
5676
5677
5678
5679
5680
5681
5682
5683
5684
5685
5686
5687
5688
5689
5690
5691
5692
5693
5694
5695
5696
5697
5698
5699
5700
5701
5702
5703
5704
5705
5706
5707
5708
5709
5710
5711
5712
5713
5714
5715
5716
5717
5718
5719
5720
5721
5722
5723
5724
5725
5726
5727
5728
5729
5730
5731
5732
5733
5734
5735
5736
5737
5738
5739
5740
5741
5742
5743
5744
5745
5746
5747
5748
5749
5750
5751
5752
5753
5754
5755
5756
5757
5758
5759
5760
5761
5762
5763
5764
5765
5766
5767
5768
5769
5770
5771
5772
5773
5774
5775
5776
5777
5778
5779
5780
5781
5782
5783
5784
5785
5786
5787
5788
5789
5790
5791
5792
5793
5794
5795
5796
5797
5798
5799
5800
5801
5802
5803
5804
5805
5806
5807
5808
5809
5810
5811
5812
5813
5814
5815
5816
5817
5818
5819
5820 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
5821 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
5822 $ RSRC_
5823 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
5824 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
5825 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
5826 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
5827 REAL ZERO, ONE
5828 parameter( zero = 0.0e+0, one = 1.0e+0 )
5829
5830
5831 LOGICAL COLREP, NOTRAN, ROWREP, TRAN, UPPER
5832 INTEGER I, IBB, IBEG, ICCOL, ICROW, ICURROW, IEND, IIC,
5833 $ IN, IOFFAK, IOFFAN, IOFFC, J, JJC, KK, LDA,
5834 $ LDC, LDPC, MYCOL, MYROW, NPCOL, NPROW
5835 REAL EPS, ERRI
5836
5837
5838 EXTERNAL blacs_gridinfo, igsum2d,
pb_infog2l, sgamx2d
5839
5840
5841 LOGICAL LSAME
5842 REAL PSLAMCH
5844
5845
5846 INTRINSIC abs,
max,
min, mod, sqrt
5847
5848
5849
5850 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
5851
5853
5854 upper =
lsame( uplo,
'U' )
5855 notran =
lsame( trans,
'N' )
5856 tran =
lsame( trans,
'T' )
5857
5858 lda =
max( 1, desca( m_ ) )
5859 ldc =
max( 1, descc( m_ ) )
5860
5861
5862
5863
5864
5865 DO 140 j = 1, n
5866
5867 IF( upper ) THEN
5868 ibeg = 1
5869 iend = j
5870 ELSE
5871 ibeg = j
5872 iend = n
5873 END IF
5874
5875 DO 10 i = 1, n
5876 ct( i ) = zero
5877 g( i ) = zero
5878 10 CONTINUE
5879
5880 IF( notran ) THEN
5881 DO 30 kk = 1, k
5882 ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
5883 DO 20 i = ibeg, iend
5884 ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
5885 ct( i ) = ct( i ) + a( ioffak ) * a( ioffan )
5886 g( i ) = g( i ) + abs( a( ioffak ) ) *
5887 $ abs( a( ioffan ) )
5888 20 CONTINUE
5889 30 CONTINUE
5890 ELSE IF( tran ) THEN
5891 DO 50 kk = 1, k
5892 ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
5893 DO 40 i = ibeg, iend
5894 ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
5895 ct( i ) = ct( i ) + a( ioffak ) * a( ioffan )
5896 g( i ) = g( i ) + abs( a( ioffak ) ) *
5897 $ abs( a( ioffan ) )
5898 40 CONTINUE
5899 50 CONTINUE
5900 END IF
5901
5902 ioffc = ic + ibeg - 1 + ( jc + j - 2 ) * ldc
5903
5904 DO 100 i = ibeg, iend
5905 ct( i ) = alpha*ct( i ) + beta * c( ioffc )
5906 g( i ) = abs( alpha )*g( i ) + abs( beta )*abs( c( ioffc ) )
5907 c( ioffc ) = ct( i )
5908 ioffc = ioffc + 1
5909 100 CONTINUE
5910
5911
5912
5913 err = zero
5914 info = 0
5915 ldpc = descc( lld_ )
5916 ioffc = ic + ( jc + j - 2 ) * ldc
5917 CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
5918 $ iic, jjc, icrow, iccol )
5919 icurrow = icrow
5920 rowrep = ( icrow.EQ.-1 )
5921 colrep = ( iccol.EQ.-1 )
5922
5923 IF( mycol.EQ.iccol .OR. colrep ) THEN
5924
5925 ibb = descc( imb_ ) - ic + 1
5926 IF( ibb.LE.0 )
5927 $ ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
5929 in = ic + ibb - 1
5930
5931 DO 110 i = ic, in
5932
5933 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5934 erri = abs( pc( iic+(jjc-1)*ldpc ) -
5935 $ c( ioffc ) ) / eps
5936 IF( g( i-ic+1 ).NE.zero )
5937 $ erri = erri / g( i-ic+1 )
5938 err =
max( err, erri )
5939 IF( err*sqrt( eps ).GE.one )
5940 $ info = 1
5941 iic = iic + 1
5942 END IF
5943
5944 ioffc = ioffc + 1
5945
5946 110 CONTINUE
5947
5948 icurrow = mod( icurrow+1, nprow )
5949
5950 DO 130 i = in+1, ic+n-1, descc( mb_ )
5951 ibb =
min( ic+n-i, descc( mb_ ) )
5952
5953 DO 120 kk = 0, ibb-1
5954
5955 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5956 erri = abs( pc( iic+(jjc-1)*ldpc ) -
5957 $ c( ioffc ) )/eps
5958 IF( g( i+kk-ic+1 ).NE.zero )
5959 $ erri = erri / g( i+kk-ic+1 )
5960 err =
max( err, erri )
5961 IF( err*sqrt( eps ).GE.one )
5962 $ info = 1
5963 iic = iic + 1
5964 END IF
5965
5966 ioffc = ioffc + 1
5967
5968 120 CONTINUE
5969
5970 icurrow = mod( icurrow+1, nprow )
5971
5972 130 CONTINUE
5973
5974 END IF
5975
5976
5977
5978 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
5979 CALL sgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
5980 $ mycol )
5981 IF( info.NE.0 )
5982 $ GO TO 150
5983
5984 140 CONTINUE
5985
5986 150 CONTINUE
5987
5988 RETURN
5989
5990
5991
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
real function pslamch(ictxt, cmach)