5789
5790
5791
5792
5793
5794
5795
5796 CHARACTER*1 TRANS, UPLO
5797 INTEGER IA, IC, ICTXT, INFO, JA, JC, K, N
5798 REAL ERR
5799 COMPLEX ALPHA, BETA
5800
5801
5802 INTEGER DESCA( * ), DESCC( * )
5803 REAL G( * )
5804 COMPLEX A( * ), C( * ), CT( * ), PC( * )
5805
5806
5807
5808
5809
5810
5811
5812
5813
5814
5815
5816
5817
5818
5819
5820
5821
5822
5823
5824
5825
5826
5827
5828
5829
5830
5831
5832
5833
5834
5835
5836
5837
5838
5839
5840
5841
5842
5843
5844
5845
5846
5847
5848
5849
5850
5851
5852
5853
5854
5855
5856
5857
5858
5859
5860
5861
5862
5863
5864
5865
5866
5867
5868
5869
5870
5871
5872
5873
5874
5875
5876
5877
5878
5879
5880
5881
5882
5883
5884
5885
5886
5887
5888
5889
5890
5891
5892
5893
5894
5895
5896
5897
5898
5899
5900
5901
5902
5903
5904
5905
5906
5907
5908
5909
5910
5911
5912
5913
5914
5915
5916
5917
5918
5919
5920
5921
5922
5923
5924
5925
5926
5927
5928
5929
5930
5931
5932
5933
5934
5935
5936
5937
5938
5939
5940
5941
5942
5943
5944
5945
5946
5947
5948
5949
5950
5951
5952
5953
5954
5955
5956
5957
5958
5959
5960
5961
5962 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
5963 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
5964 $ RSRC_
5965 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
5966 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
5967 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
5968 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
5969 REAL RZERO, RONE
5970 parameter( rzero = 0.0e+0, rone = 1.0e+0 )
5971 COMPLEX ZERO
5972 parameter( zero = ( 0.0e+0, 0.0e+0 ) )
5973
5974
5975 LOGICAL COLREP, HTRAN, NOTRAN, ROWREP, TRAN, UPPER
5976 INTEGER I, IBB, IBEG, ICCOL, ICROW, ICURROW, IEND, IIC,
5977 $ IN, IOFFAK, IOFFAN, IOFFC, J, JJC, KK, LDA,
5978 $ LDC, LDPC, MYCOL, MYROW, NPCOL, NPROW
5979 REAL EPS, ERRI
5980 COMPLEX Z
5981
5982
5983 EXTERNAL blacs_gridinfo, igsum2d,
pb_infog2l, sgamx2d
5984
5985
5986 LOGICAL LSAME
5987 REAL PSLAMCH
5989
5990
5991 INTRINSIC abs, aimag, conjg,
max,
min, mod, real, sqrt
5992
5993
5994 REAL ABS1
5995 abs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
5996
5997
5998
5999 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
6000
6002
6003 upper =
lsame( uplo,
'U' )
6004 notran =
lsame( trans,
'N' )
6005 tran =
lsame( trans,
'T' )
6006 htran =
lsame( trans,
'H' )
6007
6008 lda =
max( 1, desca( m_ ) )
6009 ldc =
max( 1, descc( m_ ) )
6010
6011
6012
6013
6014
6015 DO 140 j = 1, n
6016
6017 IF( upper ) THEN
6018 ibeg = 1
6019 iend = j
6020 ELSE
6021 ibeg = j
6022 iend = n
6023 END IF
6024
6025 DO 10 i = 1, n
6026 ct( i ) = zero
6027 g( i ) = rzero
6028 10 CONTINUE
6029
6030 IF( notran ) THEN
6031 DO 30 kk = 1, k
6032 ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
6033 DO 20 i = ibeg, iend
6034 ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
6035 ct( i ) = ct( i ) + a( ioffak ) * a( ioffan )
6036 g( i ) = g( i ) + abs1( a( ioffak ) ) *
6037 $ abs1( a( ioffan ) )
6038 20 CONTINUE
6039 30 CONTINUE
6040 ELSE IF( tran ) THEN
6041 DO 50 kk = 1, k
6042 ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
6043 DO 40 i = ibeg, iend
6044 ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
6045 ct( i ) = ct( i ) + a( ioffak ) * a( ioffan )
6046 g( i ) = g( i ) + abs1( a( ioffak ) ) *
6047 $ abs1( a( ioffan ) )
6048 40 CONTINUE
6049 50 CONTINUE
6050 ELSE IF( htran ) THEN
6051 DO 70 kk = 1, k
6052 ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
6053 DO 60 i = ibeg, iend
6054 ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
6055 ct( i ) = ct( i ) + a( ioffan ) *
6056 $ conjg( a( ioffak ) )
6057 g( i ) = g( i ) + abs1( a( ioffak ) ) *
6058 $ abs1( a( ioffan ) )
6059 60 CONTINUE
6060 70 CONTINUE
6061 ELSE
6062 DO 90 kk = 1, k
6063 ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
6064 DO 80 i = ibeg, iend
6065 ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
6066 ct( i ) = ct( i ) + conjg( a( ioffan ) ) * a( ioffak )
6067 g( i ) = g( i ) + abs1( conjg( a( ioffan ) ) ) *
6068 $ abs1( a( ioffak ) )
6069 80 CONTINUE
6070 90 CONTINUE
6071 END IF
6072
6073 ioffc = ic + ibeg - 1 + ( jc + j - 2 ) * ldc
6074
6075 DO 100 i = ibeg, iend
6076 ct( i ) = alpha*ct( i ) + beta * c( ioffc )
6077 g( i ) = abs1( alpha )*g( i ) +
6078 $ abs1( beta )*abs1( c( ioffc ) )
6079 c( ioffc ) = ct( i )
6080 ioffc = ioffc + 1
6081 100 CONTINUE
6082
6083
6084
6085 err = rzero
6086 info = 0
6087 ldpc = descc( lld_ )
6088 ioffc = ic + ( jc + j - 2 ) * ldc
6089 CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
6090 $ iic, jjc, icrow, iccol )
6091 icurrow = icrow
6092 rowrep = ( icrow.EQ.-1 )
6093 colrep = ( iccol.EQ.-1 )
6094
6095 IF( mycol.EQ.iccol .OR. colrep ) THEN
6096
6097 ibb = descc( imb_ ) - ic + 1
6098 IF( ibb.LE.0 )
6099 $ ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
6101 in = ic + ibb - 1
6102
6103 DO 110 i = ic, in
6104
6105 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
6106 erri = abs( pc( iic+(jjc-1)*ldpc ) -
6107 $ c( ioffc ) ) / eps
6108 IF( g( i-ic+1 ).NE.rzero )
6109 $ erri = erri / g( i-ic+1 )
6110 err =
max( err, erri )
6111 IF( err*sqrt( eps ).GE.rone )
6112 $ info = 1
6113 iic = iic + 1
6114 END IF
6115
6116 ioffc = ioffc + 1
6117
6118 110 CONTINUE
6119
6120 icurrow = mod( icurrow+1, nprow )
6121
6122 DO 130 i = in+1, ic+n-1, descc( mb_ )
6123 ibb =
min( ic+n-i, descc( mb_ ) )
6124
6125 DO 120 kk = 0, ibb-1
6126
6127 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
6128 erri = abs( pc( iic+(jjc-1)*ldpc ) -
6129 $ c( ioffc ) )/eps
6130 IF( g( i+kk-ic+1 ).NE.rzero )
6131 $ erri = erri / g( i+kk-ic+1 )
6132 err =
max( err, erri )
6133 IF( err*sqrt( eps ).GE.rone )
6134 $ info = 1
6135 iic = iic + 1
6136 END IF
6137
6138 ioffc = ioffc + 1
6139
6140 120 CONTINUE
6141
6142 icurrow = mod( icurrow+1, nprow )
6143
6144 130 CONTINUE
6145
6146 END IF
6147
6148
6149
6150 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
6151 CALL sgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
6152 $ mycol )
6153 IF( info.NE.0 )
6154 $ GO TO 150
6155
6156 140 CONTINUE
6157
6158 150 CONTINUE
6159
6160 RETURN
6161
6162
6163
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
real function pslamch(ictxt, cmach)