4606
 4607
 4608
 4609
 4610
 4611
 4612
 4613      CHARACTER*1        TRANS, UPLO
 4614      INTEGER            IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
 4615     $                   JY, M, N
 4616      REAL               ERR
 4617      COMPLEX            ALPHA
 4618
 4619
 4620      INTEGER            DESCA( * ), DESCX( * ), DESCY( * )
 4621      REAL               G( * )
 4622      COMPLEX            A( * ), PA( * ), X( * ), Y( * )
 4623
 4624
 4625
 4626
 4627
 4628
 4629
 4630
 4631
 4632
 4633
 4634
 4635
 4636
 4637
 4638
 4639
 4640
 4641
 4642
 4643
 4644
 4645
 4646
 4647
 4648
 4649
 4650
 4651
 4652
 4653
 4654
 4655
 4656
 4657
 4658
 4659
 4660
 4661
 4662
 4663
 4664
 4665
 4666
 4667
 4668
 4669
 4670
 4671
 4672
 4673
 4674
 4675
 4676
 4677
 4678
 4679
 4680
 4681
 4682
 4683
 4684
 4685
 4686
 4687
 4688
 4689
 4690
 4691
 4692
 4693
 4694
 4695
 4696
 4697
 4698
 4699
 4700
 4701
 4702
 4703
 4704
 4705
 4706
 4707
 4708
 4709
 4710
 4711
 4712
 4713
 4714
 4715
 4716
 4717
 4718
 4719
 4720
 4721
 4722
 4723
 4724
 4725
 4726
 4727
 4728
 4729
 4730
 4731
 4732
 4733
 4734
 4735
 4736
 4737
 4738
 4739
 4740
 4741
 4742
 4743
 4744
 4745
 4746
 4747
 4748
 4749
 4750
 4751
 4752
 4753
 4754
 4755
 4756
 4757
 4758
 4759
 4760
 4761
 4762
 4763
 4764
 4765
 4766
 4767
 4768
 4769
 4770
 4771
 4772
 4773
 4774
 4775
 4776
 4777
 4778
 4779
 4780
 4781
 4782
 4783
 4784
 4785
 4786
 4787
 4788
 4789
 4790
 4791
 4792
 4793
 4794
 4795
 4796
 4797
 4798
 4799
 4800
 4801
 4802
 4803
 4804
 4805      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
 4806     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
 4807     $                   RSRC_
 4808      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
 4809     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
 4810     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
 4811     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
 4812      REAL               ZERO, ONE
 4813      parameter( zero = 0.0e+0, one = 1.0e+0 )
 4814
 4815
 4816      LOGICAL            COLREP, CTRAN, LOWER, ROWREP, UPPER
 4817      INTEGER            I, IACOL, IAROW, IB, IBEG, ICURROW, IEND, IIA,
 4818     $                   IN, IOFFA, IOFFX, IOFFY, J, JJA, KK, LDA, LDPA,
 4819     $                   LDX, LDY, MYCOL, MYROW, NPCOL, NPROW
 4820      REAL               EPS, ERRI, GTMP
 4821      COMPLEX            ATMP, C
 4822
 4823
 4824      EXTERNAL           blacs_gridinfo, igsum2d, 
pb_infog2l, sgamx2d
 
 4825
 4826
 4827      LOGICAL            LSAME
 4828      REAL               PSLAMCH
 4830
 4831
 4832      INTRINSIC          abs, aimag, conjg, 
max, 
min, mod, real, sqrt
 
 4833
 4834
 4835      REAL               ABS1
 4836      abs1( c ) = abs( real( c ) ) + abs( aimag( c ) )
 4837
 4838
 4839
 4840      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
 4841
 4843
 4844      ctran = 
lsame( trans, 
'C' )
 
 4845      upper = 
lsame( uplo, 
'U' )
 
 4846      lower = 
lsame( uplo, 
'L' )
 
 4847
 4848      lda = 
max( 1, desca( m_ ) )
 
 4849      ldx = 
max( 1, descx( m_ ) )
 
 4850      ldy = 
max( 1, descy( m_ ) )
 
 4851
 4852
 4853
 4854
 4855
 4856      DO 70 j = 1, n
 4857
 4858         ioffy = iy + ( jy - 1 ) * ldy + ( j - 1 ) * incy
 4859
 4860         IF( lower ) THEN
 4861            ibeg = j
 4862            iend = m
 4863            DO 10 i = 1, j-1
 4864               g( i ) = zero
 4865   10       CONTINUE
 4866         ELSE IF( upper ) THEN
 4867            ibeg = 1
 4868            iend = j
 4869            DO 20 i = j+1, m
 4870               g( i ) = zero
 4871   20       CONTINUE
 4872         ELSE
 4873            ibeg = 1
 4874            iend = m
 4875         END IF
 4876
 4877         DO 30 i = ibeg, iend
 4878
 4879            ioffx = ix + ( jx - 1 ) * ldx + ( i - 1 ) * incx
 4880            ioffa = ia + i - 1 + ( ja + j - 2 ) * lda
 4881            IF( ctran ) THEN
 4882               atmp = x( ioffx ) * conjg( y( ioffy ) )
 4883            ELSE
 4884               atmp = x( ioffx ) * y( ioffy )
 4885            END IF
 4886            gtmp = abs1( x( ioffx ) ) * abs1( y( ioffy ) )
 4887            g( i ) = abs1( alpha ) * gtmp + abs1( a( ioffa ) )
 4888            a( ioffa ) = alpha * atmp + a( ioffa )
 4889
 4890   30    CONTINUE
 4891
 4892
 4893
 4894         info = 0
 4895         err  = zero
 4896         ldpa = desca( lld_ )
 4897         ioffa = ia + ( ja + j - 2 ) * lda
 4898         CALL pb_infog2l( ia, ja+j-1, desca, nprow, npcol, myrow, mycol,
 
 4899     $                    iia, jja, iarow, iacol )
 4900         rowrep = ( iarow.EQ.-1 )
 4901         colrep = ( iacol.EQ.-1 )
 4902
 4903         IF( mycol.EQ.iacol .OR. colrep ) THEN
 4904
 4905            icurrow = iarow
 4906            ib = desca( imb_ ) - ia + 1
 4907            IF( ib.LE.0 )
 4908     $         ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
 4910            in = ia + ib - 1
 4911
 4912            DO 40 i = ia, in
 4913
 4914               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
 4915                  erri = abs( pa( iia+(jja-1)*ldpa ) - a( ioffa ) )/eps
 4916                  IF( g( i-ia+1 ).NE.zero )
 4917     $               erri = erri / g( i-ia+1 )
 4918                  err = 
max( err, erri )
 
 4919                  IF( err*sqrt( eps ).GE.one )
 4920     $               info = 1
 4921                  iia = iia + 1
 4922               END IF
 4923
 4924               ioffa = ioffa + 1
 4925
 4926   40       CONTINUE
 4927
 4928            icurrow = mod( icurrow+1, nprow )
 4929
 4930            DO 60 i = in+1, ia+m-1, desca( mb_ )
 4931               ib = 
min( ia+m-i, desca( mb_ ) )
 
 4932
 4933               DO 50 kk = 0, ib-1
 4934
 4935                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
 4936                     erri = abs( pa( iia+(jja-1)*ldpa )-a( ioffa ) )/eps
 4937                     IF( g( i+kk-ia+1 ).NE.zero )
 4938     $                  erri = erri / g( i+kk-ia+1 )
 4939                     err = 
max( err, erri )
 
 4940                     IF( err*sqrt( eps ).GE.one )
 4941     $                  info = 1
 4942                     iia = iia + 1
 4943                  END IF
 4944
 4945                  ioffa = ioffa + 1
 4946
 4947   50          CONTINUE
 4948
 4949               icurrow = mod( icurrow+1, nprow )
 4950
 4951   60       CONTINUE
 4952
 4953         END IF
 4954
 4955
 4956
 4957         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
 4958         CALL sgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
 4959     $                 mycol )
 4960         IF( info.NE.0 )
 4961     $      GO TO 80
 4962
 4963   70 CONTINUE
 4964
 4965   80 CONTINUE
 4966
 4967      RETURN
 4968
 4969
 4970
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
 
real function pslamch(ictxt, cmach)