4157
4158
4159
4160
4161
4162
4163
4164 CHARACTER*1 TRANS
4165 INTEGER IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
4166 $ JY, M, N
4167 DOUBLE PRECISION ALPHA, BETA, ERR
4168
4169
4170 INTEGER DESCA( * ), DESCX( * ), DESCY( * )
4171 DOUBLE PRECISION A( * ), G( * ), PY( * ), X( * ), Y( * )
4172
4173
4174
4175
4176
4177
4178
4179
4180
4181
4182
4183
4184
4185
4186
4187
4188
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 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
4351 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
4352 $ RSRC_
4353 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
4354 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
4355 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
4356 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
4357 DOUBLE PRECISION ZERO, ONE
4358 parameter( zero = 0.0d+0, one = 1.0d+0 )
4359
4360
4361 LOGICAL COLREP, ROWREP, TRAN
4362 INTEGER I, IB, ICURCOL, ICURROW, IIY, IN, IOFFA, IOFFX,
4363 $ IOFFY, IYCOL, IYROW, J, JB, JJY, JN, KK, LDA,
4364 $ LDPY, LDX, LDY, ML, MYCOL, MYROW, NL, NPCOL,
4365 $ NPROW
4366 DOUBLE PRECISION EPS, ERRI, GTMP, TBETA, YTMP
4367
4368
4369 EXTERNAL blacs_gridinfo, dgamx2d, igsum2d,
pb_infog2l
4370
4371
4372 LOGICAL LSAME
4373 DOUBLE PRECISION PDLAMCH
4375
4376
4377 INTRINSIC abs,
max,
min, mod, sqrt
4378
4379
4380
4381 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
4382
4384
4385 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
4386 tbeta = one
4387 ELSE
4388 tbeta = beta
4389 END IF
4390
4391 tran =
lsame( trans,
'T' ).OR.
lsame( trans,
'C' )
4392 IF( tran ) THEN
4393 ml = n
4394 nl = m
4395 ELSE
4396 ml = m
4397 nl = n
4398 END IF
4399
4400 lda =
max( 1, desca( m_ ) )
4401 ldx =
max( 1, descx( m_ ) )
4402 ldy =
max( 1, descy( m_ ) )
4403
4404
4405
4406
4407
4408 ioffy = iy + ( jy - 1 ) * ldy
4409 DO 30 i = 1, ml
4410 ytmp = zero
4411 gtmp = zero
4412 ioffx = ix + ( jx - 1 ) * ldx
4413 IF( tran )THEN
4414 ioffa = ia + ( ja + i - 2 ) * lda
4415 DO 10 j = 1, nl
4416 ytmp = ytmp + a( ioffa ) * x( ioffx )
4417 gtmp = gtmp + abs( a( ioffa ) * x( ioffx ) )
4418 ioffa = ioffa + 1
4419 ioffx = ioffx + incx
4420 10 CONTINUE
4421 ELSE
4422 ioffa = ia + i - 1 + ( ja - 1 ) * lda
4423 DO 20 j = 1, nl
4424 ytmp = ytmp + a( ioffa ) * x( ioffx )
4425 gtmp = gtmp + abs( a( ioffa ) * x( ioffx ) )
4426 ioffa = ioffa + lda
4427 ioffx = ioffx + incx
4428 20 CONTINUE
4429 END IF
4430 g( i ) = abs( alpha ) * gtmp + abs( tbeta * y( ioffy ) )
4431 y( ioffy ) = alpha * ytmp + tbeta * y( ioffy )
4432 ioffy = ioffy + incy
4433 30 CONTINUE
4434
4435
4436
4437 err = zero
4438 info = 0
4439 ldpy = descy( lld_ )
4440 ioffy = iy + ( jy - 1 ) * ldy
4441 CALL pb_infog2l( iy, jy, descy, nprow, npcol, myrow, mycol, iiy,
4442 $ jjy, iyrow, iycol )
4443 icurrow = iyrow
4444 icurcol = iycol
4445 rowrep = ( iyrow.EQ.-1 )
4446 colrep = ( iycol.EQ.-1 )
4447
4448 IF( incy.EQ.descy( m_ ) ) THEN
4449
4450
4451
4452 jb = descy( inb_ ) - jy + 1
4453 IF( jb.LE.0 )
4454 $ jb = ( ( -jb ) / descy( nb_ ) + 1 ) * descy( nb_ ) + jb
4456 jn = jy + jb - 1
4457
4458 DO 50 j = jy, jn
4459
4460 IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
4461 $ ( mycol.EQ.icurcol .OR. colrep ) ) THEN
4462 erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) ) / eps
4463 IF( g( j-jy+1 ).NE.zero )
4464 $ erri = erri / g( j-jy+1 )
4465 err =
max( err, erri )
4466 IF( err*sqrt( eps ).GE.one )
4467 $ info = 1
4468 jjy = jjy + 1
4469 END IF
4470
4471 ioffy = ioffy + incy
4472
4473 50 CONTINUE
4474
4475 icurcol = mod( icurcol+1, npcol )
4476
4477 DO 70 j = jn+1, jy+ml-1, descy( nb_ )
4478 jb =
min( jy+ml-j, descy( nb_ ) )
4479
4480 DO 60 kk = 0, jb-1
4481
4482 IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
4483 $ ( mycol.EQ.icurcol .OR. colrep ) ) THEN
4484 erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) )/eps
4485 IF( g( j+kk-jy+1 ).NE.zero )
4486 $ erri = erri / g( j+kk-jy+1 )
4487 err =
max( err, erri )
4488 IF( err*sqrt( eps ).GE.one )
4489 $ info = 1
4490 jjy = jjy + 1
4491 END IF
4492
4493 ioffy = ioffy + incy
4494
4495 60 CONTINUE
4496
4497 icurcol = mod( icurcol+1, npcol )
4498
4499 70 CONTINUE
4500
4501 ELSE
4502
4503
4504
4505 ib = descy( imb_ ) - iy + 1
4506 IF( ib.LE.0 )
4507 $ ib = ( ( -ib ) / descy( mb_ ) + 1 ) * descy( mb_ ) + ib
4509 in = iy + ib - 1
4510
4511 DO 80 i = iy, in
4512
4513 IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
4514 $ ( mycol.EQ.icurcol .OR. colrep ) ) THEN
4515 erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) ) / eps
4516 IF( g( i-iy+1 ).NE.zero )
4517 $ erri = erri / g( i-iy+1 )
4518 err =
max( err, erri )
4519 IF( err*sqrt( eps ).GE.one )
4520 $ info = 1
4521 iiy = iiy + 1
4522 END IF
4523
4524 ioffy = ioffy + incy
4525
4526 80 CONTINUE
4527
4528 icurrow = mod( icurrow+1, nprow )
4529
4530 DO 100 i = in+1, iy+ml-1, descy( mb_ )
4531 ib =
min( iy+ml-i, descy( mb_ ) )
4532
4533 DO 90 kk = 0, ib-1
4534
4535 IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
4536 $ ( mycol.EQ.icurcol .OR. colrep ) ) THEN
4537 erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) )/eps
4538 IF( g( i+kk-iy+1 ).NE.zero )
4539 $ erri = erri / g( i+kk-iy+1 )
4540 err =
max( err, erri )
4541 IF( err*sqrt( eps ).GE.one )
4542 $ info = 1
4543 iiy = iiy + 1
4544 END IF
4545
4546 ioffy = ioffy + incy
4547
4548 90 CONTINUE
4549
4550 icurrow = mod( icurrow+1, nprow )
4551
4552 100 CONTINUE
4553
4554 END IF
4555
4556
4557
4558 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
4559 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
4560 $ mycol )
4561
4562 RETURN
4563
4564
4565
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
double precision function pdlamch(ictxt, cmach)