5336
5337
5338
5339
5340
5341
5342
5343 CHARACTER*1 TRANSA, TRANSB
5344 INTEGER IA, IB, IC, ICTXT, INFO, JA, JB, JC, K, M, N
5345 REAL ERR
5346 COMPLEX ALPHA, BETA
5347
5348
5349 INTEGER DESCA( * ), DESCB( * ), DESCC( * )
5350 REAL G( * )
5351 COMPLEX A( * ), B( * ), C( * ), CT( * ), PC( * )
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
5463
5464
5465
5466
5467
5468
5469
5470
5471
5472
5473
5474
5475
5476
5477
5478
5479
5480
5481
5482
5483
5484
5485
5486
5487
5488
5489
5490
5491
5492
5493
5494
5495
5496
5497
5498
5499
5500
5501
5502
5503
5504
5505
5506
5507
5508
5509
5510
5511
5512
5513
5514
5515
5516
5517
5518
5519
5520
5521
5522
5523
5524
5525
5526
5527 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
5528 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
5529 $ RSRC_
5530 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
5531 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
5532 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
5533 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
5534 REAL RZERO, RONE
5535 parameter( rzero = 0.0e+0, rone = 1.0e+0 )
5536 COMPLEX ZERO
5537 parameter( zero = ( 0.0e+0, 0.0e+0 ) )
5538
5539
5540 LOGICAL COLREP, CTRANA, CTRANB, ROWREP, TRANA, TRANB
5541 INTEGER I, IBB, ICCOL, ICROW, ICURROW, IIC, IN, IOFFA,
5542 $ IOFFB, IOFFC, J, JJC, KK, LDA, LDB, LDC, LDPC,
5543 $ MYCOL, MYROW, NPCOL, NPROW
5544 REAL EPS, ERRI
5545 COMPLEX Z
5546
5547
5548 EXTERNAL blacs_gridinfo, igsum2d,
pb_infog2l, sgamx2d
5549
5550
5551 LOGICAL LSAME
5552 REAL PSLAMCH
5554
5555
5556 INTRINSIC abs, aimag, conjg,
max,
min, mod, real, sqrt
5557
5558
5559 REAL ABS1
5560 abs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
5561
5562
5563
5564 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
5565
5567
5568 trana =
lsame( transa,
'T' ).OR.
lsame( transa,
'C' )
5569 tranb =
lsame( transb,
'T' ).OR.
lsame( transb,
'C' )
5570 ctrana =
lsame( transa,
'C' )
5571 ctranb =
lsame( transb,
'C' )
5572
5573 lda =
max( 1, desca( m_ ) )
5574 ldb =
max( 1, descb( m_ ) )
5575 ldc =
max( 1, descc( m_ ) )
5576
5577
5578
5579
5580
5581 DO 240 j = 1, n
5582
5583 ioffc = ic + ( jc + j - 2 ) * ldc
5584 DO 10 i = 1, m
5585 ct( i ) = zero
5586 g( i ) = rzero
5587 10 CONTINUE
5588
5589 IF( .NOT.trana .AND. .NOT.tranb ) THEN
5590 DO 30 kk = 1, k
5591 ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
5592 DO 20 i = 1, m
5593 ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
5594 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5595 g( i ) = g( i ) + abs( a( ioffa ) ) *
5596 $ abs( b( ioffb ) )
5597 20 CONTINUE
5598 30 CONTINUE
5599 ELSE IF( trana .AND. .NOT.tranb ) THEN
5600 IF( ctrana ) THEN
5601 DO 50 kk = 1, k
5602 ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
5603 DO 40 i = 1, m
5604 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5605 ct( i ) = ct( i ) + conjg( a( ioffa ) ) *
5606 $ b( ioffb )
5607 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5608 $ abs1( b( ioffb ) )
5609 40 CONTINUE
5610 50 CONTINUE
5611 ELSE
5612 DO 70 kk = 1, k
5613 ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
5614 DO 60 i = 1, m
5615 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5616 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5617 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5618 $ abs1( b( ioffb ) )
5619 60 CONTINUE
5620 70 CONTINUE
5621 END IF
5622 ELSE IF( .NOT.trana .AND. tranb ) THEN
5623 IF( ctranb ) THEN
5624 DO 90 kk = 1, k
5625 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5626 DO 80 i = 1, m
5627 ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
5628 ct( i ) = ct( i ) + a( ioffa ) *
5629 $ conjg( b( ioffb ) )
5630 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5631 $ abs1( b( ioffb ) )
5632 80 CONTINUE
5633 90 CONTINUE
5634 ELSE
5635 DO 110 kk = 1, k
5636 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5637 DO 100 i = 1, m
5638 ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
5639 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5640 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5641 $ abs1( b( ioffb ) )
5642 100 CONTINUE
5643 110 CONTINUE
5644 END IF
5645 ELSE IF( trana .AND. tranb ) THEN
5646 IF( ctrana ) THEN
5647 IF( ctranb ) THEN
5648 DO 130 kk = 1, k
5649 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5650 DO 120 i = 1, m
5651 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5652 ct( i ) = ct( i ) + conjg( a( ioffa ) ) *
5653 $ conjg( b( ioffb ) )
5654 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5655 $ abs1( b( ioffb ) )
5656 120 CONTINUE
5657 130 CONTINUE
5658 ELSE
5659 DO 150 kk = 1, k
5660 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5661 DO 140 i = 1, m
5662 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5663 ct( i ) = ct( i ) + conjg( a( ioffa ) ) *
5664 $ b( ioffb )
5665 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5666 $ abs1( b( ioffb ) )
5667 140 CONTINUE
5668 150 CONTINUE
5669 END IF
5670 ELSE
5671 IF( ctranb ) THEN
5672 DO 170 kk = 1, k
5673 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5674 DO 160 i = 1, m
5675 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5676 ct( i ) = ct( i ) + a( ioffa ) *
5677 $ conjg( b( ioffb ) )
5678 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5679 $ abs1( b( ioffb ) )
5680 160 CONTINUE
5681 170 CONTINUE
5682 ELSE
5683 DO 190 kk = 1, k
5684 ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
5685 DO 180 i = 1, m
5686 ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
5687 ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
5688 g( i ) = g( i ) + abs1( a( ioffa ) ) *
5689 $ abs1( b( ioffb ) )
5690 180 CONTINUE
5691 190 CONTINUE
5692 END IF
5693 END IF
5694 END IF
5695
5696 DO 200 i = 1, m
5697 ct( i ) = alpha*ct( i ) + beta * c( ioffc )
5698 g( i ) = abs1( alpha )*g( i ) +
5699 $ abs1( beta )*abs1( c( ioffc ) )
5700 c( ioffc ) = ct( i )
5701 ioffc = ioffc + 1
5702 200 CONTINUE
5703
5704
5705
5706 err = rzero
5707 info = 0
5708 ldpc = descc( lld_ )
5709 ioffc = ic + ( jc + j - 2 ) * ldc
5710 CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
5711 $ iic, jjc, icrow, iccol )
5712 icurrow = icrow
5713 rowrep = ( icrow.EQ.-1 )
5714 colrep = ( iccol.EQ.-1 )
5715
5716 IF( mycol.EQ.iccol .OR. colrep ) THEN
5717
5718 ibb = descc( imb_ ) - ic + 1
5719 IF( ibb.LE.0 )
5720 $ ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
5722 in = ic + ibb - 1
5723
5724 DO 210 i = ic, in
5725
5726 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5727 erri = abs( pc( iic+(jjc-1)*ldpc ) -
5728 $ c( ioffc ) ) / eps
5729 IF( g( i-ic+1 ).NE.rzero )
5730 $ erri = erri / g( i-ic+1 )
5731 err =
max( err, erri )
5732 IF( err*sqrt( eps ).GE.rone )
5733 $ info = 1
5734 iic = iic + 1
5735 END IF
5736
5737 ioffc = ioffc + 1
5738
5739 210 CONTINUE
5740
5741 icurrow = mod( icurrow+1, nprow )
5742
5743 DO 230 i = in+1, ic+m-1, descc( mb_ )
5744 ibb =
min( ic+m-i, descc( mb_ ) )
5745
5746 DO 220 kk = 0, ibb-1
5747
5748 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
5749 erri = abs( pc( iic+(jjc-1)*ldpc ) -
5750 $ c( ioffc ) )/eps
5751 IF( g( i+kk-ic+1 ).NE.rzero )
5752 $ erri = erri / g( i+kk-ic+1 )
5753 err =
max( err, erri )
5754 IF( err*sqrt( eps ).GE.rone )
5755 $ info = 1
5756 iic = iic + 1
5757 END IF
5758
5759 ioffc = ioffc + 1
5760
5761 220 CONTINUE
5762
5763 icurrow = mod( icurrow+1, nprow )
5764
5765 230 CONTINUE
5766
5767 END IF
5768
5769
5770
5771 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
5772 CALL sgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
5773 $ mycol )
5774 IF( info.NE.0 )
5775 $ GO TO 250
5776
5777 240 CONTINUE
5778
5779 250 CONTINUE
5780
5781 RETURN
5782
5783
5784
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
real function pslamch(ictxt, cmach)