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)