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 DOUBLE PRECISION ERR
4617 COMPLEX*16 ALPHA
4618
4619
4620 INTEGER DESCA( * ), DESCX( * ), DESCY( * )
4621 DOUBLE PRECISION G( * )
4622 COMPLEX*16 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 DOUBLE PRECISION ZERO, ONE
4813 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION EPS, ERRI, GTMP
4821 COMPLEX*16 ATMP, C
4822
4823
4824 EXTERNAL blacs_gridinfo, dgamx2d, igsum2d,
pb_infog2l
4825
4826
4827 LOGICAL LSAME
4828 DOUBLE PRECISION PDLAMCH
4830
4831
4832 INTRINSIC abs, dble, dconjg, dimag,
max,
min, mod, sqrt
4833
4834
4835 DOUBLE PRECISION ABS1
4836 abs1( c ) = abs( dble( c ) ) + abs( dimag( 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 ) * dconjg( 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 dgamx2d( 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)
double precision function pdlamch(ictxt, cmach)