4919
4920
4921
4922
4923
4924
4925
4926 CHARACTER*1 UPLO
4927 INTEGER IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
4928 $ JY, M, N
4929 DOUBLE PRECISION ALPHA, ERR
4930
4931
4932 INTEGER DESCA( * ), DESCX( * ), DESCY( * )
4933 DOUBLE PRECISION A( * ), G( * ), PA( * ), X( * ), Y( * )
4934
4935
4936
4937
4938
4939
4940
4941
4942
4943
4944
4945
4946
4947
4948
4949
4950
4951
4952
4953
4954
4955
4956
4957
4958
4959
4960
4961
4962
4963
4964
4965
4966
4967
4968
4969
4970
4971
4972
4973
4974
4975
4976
4977
4978
4979
4980
4981
4982
4983
4984
4985
4986
4987
4988
4989
4990
4991
4992
4993
4994
4995
4996
4997
4998
4999
5000
5001
5002
5003
5004
5005
5006
5007
5008
5009
5010
5011
5012
5013
5014
5015
5016
5017
5018
5019
5020
5021
5022
5023
5024
5025
5026
5027
5028
5029
5030
5031
5032
5033
5034
5035
5036
5037
5038
5039
5040
5041
5042
5043
5044
5045
5046
5047
5048
5049
5050
5051
5052
5053
5054
5055
5056
5057
5058
5059
5060
5061
5062
5063
5064
5065
5066
5067
5068
5069
5070
5071
5072
5073
5074
5075
5076
5077
5078
5079
5080
5081
5082
5083
5084
5085
5086
5087
5088
5089
5090
5091
5092
5093
5094
5095
5096
5097
5098
5099
5100
5101
5102
5103
5104
5105
5106
5107
5108 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
5109 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
5110 $ RSRC_
5111 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
5112 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
5113 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
5114 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
5115 DOUBLE PRECISION ZERO, ONE
5116 parameter( zero = 0.0d+0, one = 1.0d+0 )
5117
5118
5119 LOGICAL COLREP, LOWER, ROWREP, UPPER
5120 INTEGER I, IACOL, IAROW, IB, IBEG, ICURROW, IEND, IIA,
5121 $ IN, IOFFA, IOFFXI, IOFFXJ, IOFFYI, IOFFYJ, J,
5122 $ JJA, KK, LDA, LDPA, LDX, LDY, MYCOL, MYROW,
5123 $ NPCOL, NPROW
5124 DOUBLE PRECISION EPS, ERRI, GTMP, ATMP
5125
5126
5127 EXTERNAL blacs_gridinfo, dgamx2d, igsum2d,
pb_infog2l
5128
5129
5130 LOGICAL LSAME
5131 DOUBLE PRECISION PDLAMCH
5133
5134
5135 INTRINSIC abs,
max,
min, mod, sqrt
5136
5137
5138
5139 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
5140
5142
5143 upper =
lsame( uplo,
'U' )
5144 lower =
lsame( uplo,
'L' )
5145
5146 lda =
max( 1, desca( m_ ) )
5147 ldx =
max( 1, descx( m_ ) )
5148 ldy =
max( 1, descy( m_ ) )
5149
5150
5151
5152
5153
5154 DO 70 j = 1, n
5155
5156 ioffxj = ix + ( jx - 1 ) * ldx + ( j - 1 ) * incx
5157 ioffyj = iy + ( jy - 1 ) * ldy + ( j - 1 ) * incy
5158
5159 IF( lower ) THEN
5160 ibeg = j
5161 iend = m
5162 DO 10 i = 1, j-1
5163 g( i ) = zero
5164 10 CONTINUE
5165 ELSE IF( upper ) THEN
5166 ibeg = 1
5167 iend = j
5168 DO 20 i = j+1, m
5169 g( i ) = zero
5170 20 CONTINUE
5171 ELSE
5172 ibeg = 1
5173 iend = m
5174 END IF
5175
5176 DO 30 i = ibeg, iend
5177 ioffa = ia + i - 1 + ( ja + j - 2 ) * lda
5178 ioffxi = ix + ( jx - 1 ) * ldx + ( i - 1 ) * incx
5179 ioffyi = iy + ( jy - 1 ) * ldy + ( i - 1 ) * incy
5180 atmp = x( ioffxi ) * y( ioffyj )
5181 atmp = atmp + y( ioffyi ) * x( ioffxj )
5182 gtmp = abs( x( ioffxi ) * y( ioffyj ) )
5183 gtmp = gtmp + abs( y( ioffyi ) * x( ioffxj ) )
5184 g( i ) = abs( alpha ) * gtmp + abs( a( ioffa ) )
5185 a( ioffa ) = alpha*atmp + a( ioffa )
5186
5187 30 CONTINUE
5188
5189
5190
5191 info = 0
5192 err = zero
5193 ldpa = desca( lld_ )
5194 ioffa = ia + ( ja + j - 2 ) * lda
5195 CALL pb_infog2l( ia, ja+j-1, desca, nprow, npcol, myrow, mycol,
5196 $ iia, jja, iarow, iacol )
5197 rowrep = ( iarow.EQ.-1 )
5198 colrep = ( iacol.EQ.-1 )
5199
5200 IF( mycol.EQ.iacol .OR. colrep ) THEN
5201
5202 icurrow = iarow
5203 ib = desca( imb_ ) - ia + 1
5204 IF( ib.LE.0 )
5205 $ ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
5207 in = ia + ib - 1
5208
5209 DO 40 i = ia, in
5210
5211 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5212 erri = abs( pa( iia+(jja-1)*ldpa ) - a( ioffa ) )/eps
5213 IF( g( i-ia+1 ).NE.zero )
5214 $ erri = erri / g( i-ia+1 )
5215 err =
max( err, erri )
5216 IF( err*sqrt( eps ).GE.one )
5217 $ info = 1
5218 iia = iia + 1
5219 END IF
5220
5221 ioffa = ioffa + 1
5222
5223 40 CONTINUE
5224
5225 icurrow = mod( icurrow+1, nprow )
5226
5227 DO 60 i = in+1, ia+m-1, desca( mb_ )
5228 ib =
min( ia+m-i, desca( mb_ ) )
5229
5230 DO 50 kk = 0, ib-1
5231
5232 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5233 erri = abs( pa( iia+(jja-1)*ldpa )-a( ioffa ) )/eps
5234 IF( g( i+kk-ia+1 ).NE.zero )
5235 $ erri = erri / g( i+kk-ia+1 )
5236 err =
max( err, erri )
5237 IF( err*sqrt( eps ).GE.one )
5238 $ info = 1
5239 iia = iia + 1
5240 END IF
5241
5242 ioffa = ioffa + 1
5243
5244 50 CONTINUE
5245
5246 icurrow = mod( icurrow+1, nprow )
5247
5248 60 CONTINUE
5249
5250 END IF
5251
5252
5253
5254 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
5255 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
5256 $ mycol )
5257 IF( info.NE.0 )
5258 $ GO TO 80
5259
5260 70 CONTINUE
5261
5262 80 CONTINUE
5263
5264 RETURN
5265
5266
5267
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
double precision function pdlamch(ictxt, cmach)