4975
 4976
 4977
 4978
 4979
 4980
 4981
 4982      CHARACTER*1        UPLO
 4983      INTEGER            IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
 4984     $                   JY, M, N
 4985      DOUBLE PRECISION   ERR
 4986      COMPLEX*16         ALPHA
 4987
 4988
 4989      INTEGER            DESCA( * ), DESCX( * ), DESCY( * )
 4990      DOUBLE PRECISION   G( * )
 4991      COMPLEX*16         A( * ), PA( * ), X( * ), Y( * )
 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
 5109
 5110
 5111
 5112
 5113
 5114
 5115
 5116
 5117
 5118
 5119
 5120
 5121
 5122
 5123
 5124
 5125
 5126
 5127
 5128
 5129
 5130
 5131
 5132
 5133
 5134
 5135
 5136
 5137
 5138
 5139
 5140
 5141
 5142
 5143
 5144
 5145
 5146
 5147
 5148
 5149
 5150
 5151
 5152
 5153
 5154
 5155
 5156
 5157
 5158
 5159
 5160
 5161
 5162
 5163
 5164
 5165
 5166      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
 5167     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
 5168     $                   RSRC_
 5169      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
 5170     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
 5171     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
 5172     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
 5173      DOUBLE PRECISION   ZERO, ONE
 5174      parameter( zero = 0.0d+0, one = 1.0d+0 )
 5175
 5176
 5177      LOGICAL            COLREP, LOWER, ROWREP, UPPER
 5178      INTEGER            I, IACOL, IAROW, IB, IBEG, ICURROW, IEND, IIA,
 5179     $                   IN, IOFFA, IOFFXI, IOFFXJ, IOFFYI, IOFFYJ, J,
 5180     $                   JJA, KK, LDA, LDPA, LDX, LDY, MYCOL, MYROW,
 5181     $                   NPCOL, NPROW
 5182      DOUBLE PRECISION   EPS, ERRI, GTMP
 5183      COMPLEX*16         C, ATMP
 5184
 5185
 5186      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, 
pb_infog2l 
 5187
 5188
 5189      LOGICAL            LSAME
 5190      DOUBLE PRECISION   PDLAMCH
 5192
 5193
 5194      INTRINSIC          abs, dble, dconjg, dimag, 
max, 
min, mod, sqrt
 
 5195
 5196
 5197      DOUBLE PRECISION   ABS1
 5198      abs1( c ) = abs( dble( c ) ) + abs( dimag( c ) )
 5199
 5200
 5201
 5202      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
 5203
 5205
 5206      upper = 
lsame( uplo, 
'U' )
 
 5207      lower = 
lsame( uplo, 
'L' )
 
 5208
 5209      lda = 
max( 1, desca( m_ ) )
 
 5210      ldx = 
max( 1, descx( m_ ) )
 
 5211      ldy = 
max( 1, descy( m_ ) )
 
 5212
 5213
 5214
 5215
 5216
 5217      DO 70 j = 1, n
 5218
 5219         ioffxj = ix + ( jx - 1 ) * ldx + ( j - 1 ) * incx
 5220         ioffyj = iy + ( jy - 1 ) * ldy + ( j - 1 ) * incy
 5221
 5222         IF( lower ) THEN
 5223            ibeg = j
 5224            iend = m
 5225            DO 10 i = 1, j-1
 5226               g( i ) = zero
 5227   10       CONTINUE
 5228         ELSE IF( upper ) THEN
 5229            ibeg = 1
 5230            iend = j
 5231            DO 20 i = j+1, m
 5232               g( i ) = zero
 5233   20       CONTINUE
 5234         ELSE
 5235            ibeg = 1
 5236            iend = m
 5237         END IF
 5238
 5239         DO 30 i = ibeg, iend
 5240            ioffa = ia + i - 1 + ( ja + j - 2 ) * lda
 5241            ioffxi = ix + ( jx - 1 ) * ldx + ( i - 1 ) * incx
 5242            ioffyi = iy + ( jy - 1 ) * ldy + ( i - 1 ) * incy
 5243            atmp = alpha * x( ioffxi ) * dconjg( y( ioffyj ) )
 5244            atmp = atmp + y( ioffyi ) * dconjg( alpha * x( ioffxj ) )
 5245            gtmp = abs1( alpha * x( ioffxi ) ) * abs1( y( ioffyj ) )
 5246            gtmp = gtmp + abs1( y( ioffyi ) ) *
 5247     $                    abs1( dconjg( alpha * x( ioffxj ) ) )
 5248            g( i ) = gtmp + abs1( a( ioffa ) )
 5249            a( ioffa ) = a( ioffa ) + atmp
 5250
 5251   30    CONTINUE
 5252
 5253
 5254
 5255         info = 0
 5256         err  = zero
 5257         ldpa = desca( lld_ )
 5258         ioffa = ia + ( ja + j - 2 ) * lda
 5259         CALL pb_infog2l( ia, ja+j-1, desca, nprow, npcol, myrow, mycol,
 
 5260     $                    iia, jja, iarow, iacol )
 5261         rowrep = ( iarow.EQ.-1 )
 5262         colrep = ( iacol.EQ.-1 )
 5263
 5264         IF( mycol.EQ.iacol .OR. colrep ) THEN
 5265
 5266            icurrow = iarow
 5267            ib = desca( imb_ ) - ia + 1
 5268            IF( ib.LE.0 )
 5269     $         ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
 5271            in = ia + ib - 1
 5272
 5273            DO 40 i = ia, in
 5274
 5275               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
 5276                  erri = abs( pa( iia+(jja-1)*ldpa ) - a( ioffa ) )/eps
 5277                  IF( g( i-ia+1 ).NE.zero )
 5278     $               erri = erri / g( i-ia+1 )
 5279                  err = 
max( err, erri )
 
 5280                  IF( err*sqrt( eps ).GE.one )
 5281     $               info = 1
 5282                  iia = iia + 1
 5283               END IF
 5284
 5285               ioffa = ioffa + 1
 5286
 5287   40       CONTINUE
 5288
 5289            icurrow = mod( icurrow+1, nprow )
 5290
 5291            DO 60 i = in+1, ia+m-1, desca( mb_ )
 5292               ib = 
min( ia+m-i, desca( mb_ ) )
 
 5293
 5294               DO 50 kk = 0, ib-1
 5295
 5296                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
 5297                     erri = abs( pa( iia+(jja-1)*ldpa )-a( ioffa ) )/eps
 5298                     IF( g( i+kk-ia+1 ).NE.zero )
 5299     $                  erri = erri / g( i+kk-ia+1 )
 5300                     err = 
max( err, erri )
 
 5301                     IF( err*sqrt( eps ).GE.one )
 5302     $                  info = 1
 5303                     iia = iia + 1
 5304                  END IF
 5305
 5306                  ioffa = ioffa + 1
 5307
 5308   50          CONTINUE
 5309
 5310               icurrow = mod( icurrow+1, nprow )
 5311
 5312   60       CONTINUE
 5313
 5314         END IF
 5315
 5316
 5317
 5318         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
 5319         CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
 5320     $                 mycol )
 5321         IF( info.NE.0 )
 5322     $      GO TO 80
 5323
 5324   70 CONTINUE
 5325
 5326   80 CONTINUE
 5327
 5328      RETURN
 5329
 5330
 5331
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
 
double precision function pdlamch(ictxt, cmach)