5272
5273
5274
5275
5276
5277
5278
5279 CHARACTER*1 TRANSA, TRANSB
5280 INTEGER IA, IB, IC, ICTXT, INFO, JA, JB, JC, K, M, N
5281 DOUBLE PRECISION ALPHA, BETA, ERR
5282
5283
5284 INTEGER DESCA( * ), DESCB( * ), DESCC( * )
5285 DOUBLE PRECISION A( * ), B( * ), C( * ), CT( * ), G( * ),
5286 $ PC( * )
5287
5288
5289
5290
5291
5292
5293
5294
5295
5296
5297
5298
5299
5300
5301
5302
5303
5304
5305
5306
5307
5308
5309
5310
5311
5312
5313
5314
5315
5316
5317
5318
5319
5320
5321
5322
5323
5324
5325
5326
5327
5328
5329
5330
5331
5332
5333
5334
5335
5336
5337
5338
5339
5340
5341
5342
5343
5344
5345
5346
5347
5348
5349
5350
5351
5352
5353
5354
5355
5356
5357
5358
5359
5360
5361
5362
5363
5364
5365
5366
5367
5368
5369
5370
5371
5372
5373
5374
5375
5376
5377
5378
5379
5380
5381
5382
5383
5384
5385
5386
5387
5388
5389
5390
5391
5392
5393
5394
5395
5396
5397
5398
5399
5400
5401
5402
5403
5404
5405
5406
5407
5408
5409
5410
5411
5412
5413
5414
5415
5416
5417
5418
5419
5420
5421
5422
5423
5424
5425
5426
5427
5428
5429
5430
5431
5432
5433
5434
5435
5436
5437
5438
5439
5440
5441
5442
5443
5444
5445
5446
5447
5448
5449
5450
5451
5452
5453
5454
5455
5456
5457
5458
5459
5460
5461
5462 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
5463 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
5464 $ RSRC_
5465 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
5466 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
5467 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
5468 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
5469 DOUBLE PRECISION ZERO, ONE
5470 parameter( zero = 0.0d+0, one = 1.0d+0 )
5471
5472
5473 LOGICAL COLREP, ROWREP, TRANA, TRANB
5474 INTEGER I, IBB, ICCOL, ICROW, ICURROW, IIC, IN, IOFFA,
5475 $ IOFFB, IOFFC, J, JJC, KK, LDA, LDB, LDC, LDPC,
5476 $ MYCOL, MYROW, NPCOL, NPROW
5477 DOUBLE PRECISION EPS, ERRI
5478
5479
5480 EXTERNAL blacs_gridinfo, dgamx2d, igsum2d,
pb_infog2l
5481
5482
5483 LOGICAL LSAME
5484 DOUBLE PRECISION PDLAMCH
5486
5487
5488 INTRINSIC abs,
max,
min, mod, sqrt
5489
5490
5491
5492 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
5493
5495
5496 trana =
lsame( transa,
'T' ).OR.
lsame( transa,
'C' )
5497 tranb =
lsame( transb,
'T' ).OR.
lsame( transb,
'C' )
5498
5499 lda =
max( 1, desca( m_ ) )
5500 ldb =
max( 1, descb( m_ ) )
5501 ldc =
max( 1, descc( m_ ) )
5502
5503
5504
5505
5506
5507 DO 240 j = 1, n
5508
5509 ioffc = ic + ( jc + j - 2 ) * ldc
5510 DO 10 i = 1, m
5511 ct( i ) = zero
5512 g( i ) = zero
5513 10 CONTINUE
5514
5515 IF( .NOT.trana .AND. .NOT.tranb ) THEN
5516 DO 30 kk = 1, k
5517 ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
5518 DO 20 i = 1, m
5519 ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
5520 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5521 g( i ) = g( i ) + abs( a( ioffa ) ) *
5522 $ abs( b( ioffb ) )
5523 20 CONTINUE
5524 30 CONTINUE
5525 ELSE IF( trana .AND. .NOT.tranb ) THEN
5526 DO 50 kk = 1, k
5527 ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
5528 DO 40 i = 1, m
5529 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5530 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5531 g( i ) = g( i ) + abs( a( ioffa ) ) *
5532 $ abs( b( ioffb ) )
5533 40 CONTINUE
5534 50 CONTINUE
5535 ELSE IF( .NOT.trana .AND. tranb ) THEN
5536 DO 70 kk = 1, k
5537 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5538 DO 60 i = 1, m
5539 ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
5540 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5541 g( i ) = g( i ) + abs( a( ioffa ) ) *
5542 $ abs( b( ioffb ) )
5543 60 CONTINUE
5544 70 CONTINUE
5545 ELSE IF( trana .AND. tranb ) THEN
5546 DO 90 kk = 1, k
5547 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5548 DO 80 i = 1, m
5549 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5550 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5551 g( i ) = g( i ) + abs( a( ioffa ) ) *
5552 $ abs( b( ioffb ) )
5553 80 CONTINUE
5554 90 CONTINUE
5555 END IF
5556
5557 DO 200 i = 1, m
5558 ct( i ) = alpha*ct( i ) + beta * c( ioffc )
5559 g( i ) = abs( alpha )*g( i ) + abs( beta )*abs( c( ioffc ) )
5560 c( ioffc ) = ct( i )
5561 ioffc = ioffc + 1
5562 200 CONTINUE
5563
5564
5565
5566 err = zero
5567 info = 0
5568 ldpc = descc( lld_ )
5569 ioffc = ic + ( jc + j - 2 ) * ldc
5570 CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
5571 $ iic, jjc, icrow, iccol )
5572 icurrow = icrow
5573 rowrep = ( icrow.EQ.-1 )
5574 colrep = ( iccol.EQ.-1 )
5575
5576 IF( mycol.EQ.iccol .OR. colrep ) THEN
5577
5578 ibb = descc( imb_ ) - ic + 1
5579 IF( ibb.LE.0 )
5580 $ ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
5582 in = ic + ibb - 1
5583
5584 DO 210 i = ic, in
5585
5586 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5587 erri = abs( pc( iic+(jjc-1)*ldpc ) -
5588 $ c( ioffc ) ) / eps
5589 IF( g( i-ic+1 ).NE.zero )
5590 $ erri = erri / g( i-ic+1 )
5591 err =
max( err, erri )
5592 IF( err*sqrt( eps ).GE.one )
5593 $ info = 1
5594 iic = iic + 1
5595 END IF
5596
5597 ioffc = ioffc + 1
5598
5599 210 CONTINUE
5600
5601 icurrow = mod( icurrow+1, nprow )
5602
5603 DO 230 i = in+1, ic+m-1, descc( mb_ )
5604 ibb =
min( ic+m-i, descc( mb_ ) )
5605
5606 DO 220 kk = 0, ibb-1
5607
5608 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5609 erri = abs( pc( iic+(jjc-1)*ldpc ) -
5610 $ c( ioffc ) )/eps
5611 IF( g( i+kk-ic+1 ).NE.zero )
5612 $ erri = erri / g( i+kk-ic+1 )
5613 err =
max( err, erri )
5614 IF( err*sqrt( eps ).GE.one )
5615 $ info = 1
5616 iic = iic + 1
5617 END IF
5618
5619 ioffc = ioffc + 1
5620
5621 220 CONTINUE
5622
5623 icurrow = mod( icurrow+1, nprow )
5624
5625 230 CONTINUE
5626
5627 END IF
5628
5629
5630
5631 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
5632 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
5633 $ mycol )
5634 IF( info.NE.0 )
5635 $ GO TO 250
5636
5637 240 CONTINUE
5638
5639 250 CONTINUE
5640
5641 RETURN
5642
5643
5644
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
double precision function pdlamch(ictxt, cmach)