4172
 4173
 4174
 4175
 4176
 4177
 4178
 4179      CHARACTER*1        TRANS
 4180      INTEGER            IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
 4181     $                   JY, M, N
 4182      DOUBLE PRECISION   ERR
 4183      COMPLEX*16         ALPHA, BETA
 4184
 4185
 4186      INTEGER            DESCA( * ), DESCX( * ), DESCY( * )
 4187      DOUBLE PRECISION   G( * )
 4188      COMPLEX*16         A( * ), PY( * ), X( * ), Y( * )
 4189
 4190
 4191
 4192
 4193
 4194
 4195
 4196
 4197
 4198
 4199
 4200
 4201
 4202
 4203
 4204
 4205
 4206
 4207
 4208
 4209
 4210
 4211
 4212
 4213
 4214
 4215
 4216
 4217
 4218
 4219
 4220
 4221
 4222
 4223
 4224
 4225
 4226
 4227
 4228
 4229
 4230
 4231
 4232
 4233
 4234
 4235
 4236
 4237
 4238
 4239
 4240
 4241
 4242
 4243
 4244
 4245
 4246
 4247
 4248
 4249
 4250
 4251
 4252
 4253
 4254
 4255
 4256
 4257
 4258
 4259
 4260
 4261
 4262
 4263
 4264
 4265
 4266
 4267
 4268
 4269
 4270
 4271
 4272
 4273
 4274
 4275
 4276
 4277
 4278
 4279
 4280
 4281
 4282
 4283
 4284
 4285
 4286
 4287
 4288
 4289
 4290
 4291
 4292
 4293
 4294
 4295
 4296
 4297
 4298
 4299
 4300
 4301
 4302
 4303
 4304
 4305
 4306
 4307
 4308
 4309
 4310
 4311
 4312
 4313
 4314
 4315
 4316
 4317
 4318
 4319
 4320
 4321
 4322
 4323
 4324
 4325
 4326
 4327
 4328
 4329
 4330
 4331
 4332
 4333
 4334
 4335
 4336
 4337
 4338
 4339
 4340
 4341
 4342
 4343
 4344
 4345
 4346
 4347
 4348
 4349
 4350
 4351
 4352
 4353
 4354
 4355
 4356
 4357
 4358
 4359
 4360
 4361
 4362
 4363
 4364
 4365
 4366
 4367
 4368
 4369      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
 4370     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
 4371     $                   RSRC_
 4372      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
 4373     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
 4374     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
 4375     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
 4376      DOUBLE PRECISION   RZERO, RONE
 4377      parameter( rzero = 0.0d+0, rone = 1.0d+0 )
 4378      COMPLEX*16         ZERO, ONE
 4379      parameter( zero = ( 0.0d+0, 0.0d+0 ),
 4380     $                   one = ( 1.0d+0, 0.0d+0 ) )
 4381
 4382
 4383      LOGICAL            COLREP, CTRAN, ROWREP, TRAN
 4384      INTEGER            I, IB, ICURCOL, ICURROW, IIY, IN, IOFFA, IOFFX,
 4385     $                   IOFFY, IYCOL, IYROW, J, JB, JJY, JN, KK, LDA,
 4386     $                   LDPY, LDX, LDY, ML, MYCOL, MYROW, NL, NPCOL,
 4387     $                   NPROW
 4388      DOUBLE PRECISION   EPS, ERRI, GTMP
 4389      COMPLEX*16         C, TBETA, YTMP
 4390
 4391
 4392      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, 
pb_infog2l 
 4393
 4394
 4395      LOGICAL            LSAME
 4396      DOUBLE PRECISION   PDLAMCH
 4398
 4399
 4400      INTRINSIC          abs, dble, dconjg, dimag, 
max, 
min, mod, sqrt
 
 4401
 4402
 4403      DOUBLE PRECISION   ABS1
 4404      abs1( c ) = abs( dble( c ) ) + abs( dimag( c ) )
 4405
 4406
 4407
 4408      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
 4409
 4411
 4412      IF( m.EQ.0 .OR. n.EQ.0 ) THEN
 4413         tbeta = one
 4414      ELSE
 4415         tbeta = beta
 4416      END IF
 4417
 4418      tran = 
lsame( trans, 
'T' )
 
 4419      ctran = 
lsame( trans, 
'C' )
 
 4420      IF( tran.OR.ctran ) THEN
 4421         ml = n
 4422         nl = m
 4423      ELSE
 4424         ml = m
 4425         nl = n
 4426      END IF
 4427
 4428      lda = 
max( 1, desca( m_ ) )
 
 4429      ldx = 
max( 1, descx( m_ ) )
 
 4430      ldy = 
max( 1, descy( m_ ) )
 
 4431
 4432
 4433
 4434
 4435
 4436      ioffy = iy + ( jy - 1 ) * ldy
 4437      DO 40 i = 1, ml
 4438         ytmp = zero
 4439         gtmp = rzero
 4440         ioffx = ix + ( jx - 1 ) * ldx
 4441         IF( tran )THEN
 4442            ioffa = ia + ( ja + i - 2 ) * lda
 4443            DO 10 j = 1, nl
 4444               ytmp = ytmp + a( ioffa ) * x( ioffx )
 4445               gtmp = gtmp + abs1( a( ioffa ) ) * abs1( x( ioffx ) )
 4446               ioffa = ioffa + 1
 4447               ioffx = ioffx + incx
 4448   10       CONTINUE
 4449         ELSE IF( ctran )THEN
 4450            ioffa = ia + ( ja + i - 2 ) * lda
 4451            DO 20 j = 1, nl
 4452               ytmp = ytmp + dconjg( a( ioffa ) ) * x( ioffx )
 4453               gtmp = gtmp + abs1( a( ioffa ) ) * abs1( x( ioffx ) )
 4454               ioffa = ioffa + 1
 4455               ioffx = ioffx + incx
 4456   20       CONTINUE
 4457         ELSE
 4458            ioffa = ia + i - 1 + ( ja - 1 ) * lda
 4459            DO 30 j = 1, nl
 4460               ytmp = ytmp + a( ioffa ) * x( ioffx )
 4461               gtmp = gtmp + abs1( a( ioffa ) ) * abs1( x( ioffx ) )
 4462               ioffa = ioffa + lda
 4463               ioffx = ioffx + incx
 4464   30       CONTINUE
 4465         END IF
 4466         g( i ) = abs1( alpha )*gtmp + abs1( tbeta )*abs1( y( ioffy ) )
 4467         y( ioffy ) = alpha * ytmp + tbeta * y( ioffy )
 4468         ioffy = ioffy + incy
 4469   40 CONTINUE
 4470
 4471
 4472
 4473      err  = rzero
 4474      info = 0
 4475      ldpy = descy( lld_ )
 4476      ioffy = iy + ( jy - 1 ) * ldy
 4477      CALL pb_infog2l( iy, jy, descy, nprow, npcol, myrow, mycol, iiy,
 
 4478     $                 jjy, iyrow, iycol )
 4479      icurrow = iyrow
 4480      icurcol = iycol
 4481      rowrep  = ( iyrow.EQ.-1 )
 4482      colrep  = ( iycol.EQ.-1 )
 4483
 4484      IF( incy.EQ.descy( m_ ) ) THEN
 4485
 4486
 4487
 4488         jb = descy( inb_ ) - jy + 1
 4489         IF( jb.LE.0 )
 4490     $      jb = ( ( -jb ) / descy( nb_ ) + 1 ) * descy( nb_ ) + jb
 4492         jn = jy + jb - 1
 4493
 4494         DO 50 j = jy, jn
 4495
 4496            IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
 4497     $          ( mycol.EQ.icurcol .OR. colrep ) ) THEN
 4498               erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) ) / eps
 4499               IF( g( j-jy+1 ).NE.rzero )
 4500     $            erri = erri / g( j-jy+1 )
 4501               err = 
max( err, erri )
 
 4502               IF( err*sqrt( eps ).GE.rone )
 4503     $            info = 1
 4504               jjy = jjy + 1
 4505            END IF
 4506
 4507            ioffy = ioffy + incy
 4508
 4509   50    CONTINUE
 4510
 4511         icurcol = mod( icurcol+1, npcol )
 4512
 4513         DO 70 j = jn+1, jy+ml-1, descy( nb_ )
 4514            jb = 
min( jy+ml-j, descy( nb_ ) )
 
 4515
 4516            DO 60 kk = 0, jb-1
 4517
 4518               IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
 4519     $             ( mycol.EQ.icurcol .OR. colrep ) ) THEN
 4520                  erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) )/eps
 4521                  IF( g( j+kk-jy+1 ).NE.rzero )
 4522     $               erri = erri / g( j+kk-jy+1 )
 4523                  err = 
max( err, erri )
 
 4524                  IF( err*sqrt( eps ).GE.rone )
 4525     $               info = 1
 4526                  jjy = jjy + 1
 4527               END IF
 4528
 4529               ioffy = ioffy + incy
 4530
 4531   60       CONTINUE
 4532
 4533            icurcol = mod( icurcol+1, npcol )
 4534
 4535   70    CONTINUE
 4536
 4537      ELSE
 4538
 4539
 4540
 4541         ib = descy( imb_ ) - iy + 1
 4542         IF( ib.LE.0 )
 4543     $      ib = ( ( -ib ) / descy( mb_ ) + 1 ) * descy( mb_ ) + ib
 4545         in = iy + ib - 1
 4546
 4547         DO 80 i = iy, in
 4548
 4549            IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
 4550     $          ( mycol.EQ.icurcol .OR. colrep ) ) THEN
 4551               erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) ) / eps
 4552               IF( g( i-iy+1 ).NE.rzero )
 4553     $            erri = erri / g( i-iy+1 )
 4554               err = 
max( err, erri )
 
 4555               IF( err*sqrt( eps ).GE.rone )
 4556     $            info = 1
 4557               iiy = iiy + 1
 4558            END IF
 4559
 4560            ioffy = ioffy + incy
 4561
 4562   80    CONTINUE
 4563
 4564         icurrow = mod( icurrow+1, nprow )
 4565
 4566         DO 100 i = in+1, iy+ml-1, descy( mb_ )
 4567            ib = 
min( iy+ml-i, descy( mb_ ) )
 
 4568
 4569            DO 90 kk = 0, ib-1
 4570
 4571               IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
 4572     $             ( mycol.EQ.icurcol .OR. colrep ) ) THEN
 4573                  erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) )/eps
 4574                  IF( g( i+kk-iy+1 ).NE.rzero )
 4575     $               erri = erri / g( i+kk-iy+1 )
 4576                  err = 
max( err, erri )
 
 4577                  IF( err*sqrt( eps ).GE.rone )
 4578     $               info = 1
 4579                  iiy = iiy + 1
 4580               END IF
 4581
 4582               ioffy = ioffy + incy
 4583
 4584   90       CONTINUE
 4585
 4586            icurrow = mod( icurrow+1, nprow )
 4587
 4588  100    CONTINUE
 4589
 4590      END IF
 4591
 4592
 4593
 4594      CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
 4595      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
 4596     $              mycol )
 4597
 4598      RETURN
 4599
 4600
 4601
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
 
double precision function pdlamch(ictxt, cmach)