4570
4571
4572
4573
4574
4575
4576
4577 CHARACTER*1 UPLO
4578 INTEGER IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
4579 $ JY, M, N
4580 DOUBLE PRECISION ALPHA, ERR
4581
4582
4583 INTEGER DESCA( * ), DESCX( * ), DESCY( * )
4584 DOUBLE PRECISION A( * ), G( * ), PA( * ), X( * ), Y( * )
4585
4586
4587
4588
4589
4590
4591
4592
4593
4594
4595
4596
4597
4598
4599
4600
4601
4602
4603
4604
4605
4606
4607
4608
4609
4610
4611
4612
4613
4614
4615
4616
4617
4618
4619
4620
4621
4622
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 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
4760 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
4761 $ RSRC_
4762 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
4763 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
4764 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
4765 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
4766 DOUBLE PRECISION ZERO, ONE
4767 parameter( zero = 0.0d+0, one = 1.0d+0 )
4768
4769
4770 LOGICAL COLREP, LOWER, ROWREP, UPPER
4771 INTEGER I, IACOL, IAROW, IB, IBEG, ICURROW, IEND, IIA,
4772 $ IN, IOFFA, IOFFX, IOFFY, J, JJA, KK, LDA, LDPA,
4773 $ LDX, LDY, MYCOL, MYROW, NPCOL, NPROW
4774 DOUBLE PRECISION ATMP, EPS, ERRI, GTMP
4775
4776
4777 EXTERNAL blacs_gridinfo, dgamx2d, igsum2d,
pb_infog2l
4778
4779
4780 LOGICAL LSAME
4781 DOUBLE PRECISION PDLAMCH
4783
4784
4785 INTRINSIC abs,
max,
min, mod, sqrt
4786
4787
4788
4789 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
4790
4792
4793 upper =
lsame( uplo,
'U' )
4794 lower =
lsame( uplo,
'L' )
4795
4796 lda =
max( 1, desca( m_ ) )
4797 ldx =
max( 1, descx( m_ ) )
4798 ldy =
max( 1, descy( m_ ) )
4799
4800
4801
4802
4803
4804 DO 70 j = 1, n
4805
4806 ioffy = iy + ( jy - 1 ) * ldy + ( j - 1 ) * incy
4807
4808 IF( lower ) THEN
4809 ibeg = j
4810 iend = m
4811 DO 10 i = 1, j-1
4812 g( i ) = zero
4813 10 CONTINUE
4814 ELSE IF( upper ) THEN
4815 ibeg = 1
4816 iend = j
4817 DO 20 i = j+1, m
4818 g( i ) = zero
4819 20 CONTINUE
4820 ELSE
4821 ibeg = 1
4822 iend = m
4823 END IF
4824
4825 DO 30 i = ibeg, iend
4826
4827 ioffx = ix + ( jx - 1 ) * ldx + ( i - 1 ) * incx
4828 ioffa = ia + i - 1 + ( ja + j - 2 ) * lda
4829 atmp = x( ioffx ) * y( ioffy )
4830 gtmp = abs( x( ioffx ) * y( ioffy ) )
4831 g( i ) = abs( alpha ) * gtmp + abs( a( ioffa ) )
4832 a( ioffa ) = alpha * atmp + a( ioffa )
4833
4834 30 CONTINUE
4835
4836
4837
4838 info = 0
4839 err = zero
4840 ldpa = desca( lld_ )
4841 ioffa = ia + ( ja + j - 2 ) * lda
4842 CALL pb_infog2l( ia, ja+j-1, desca, nprow, npcol, myrow, mycol,
4843 $ iia, jja, iarow, iacol )
4844 rowrep = ( iarow.EQ.-1 )
4845 colrep = ( iacol.EQ.-1 )
4846
4847 IF( mycol.EQ.iacol .OR. colrep ) THEN
4848
4849 icurrow = iarow
4850 ib = desca( imb_ ) - ia + 1
4851 IF( ib.LE.0 )
4852 $ ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
4854 in = ia + ib - 1
4855
4856 DO 40 i = ia, in
4857
4858 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
4859 erri = abs( pa( iia+(jja-1)*ldpa ) - a( ioffa ) )/eps
4860 IF( g( i-ia+1 ).NE.zero )
4861 $ erri = erri / g( i-ia+1 )
4862 err =
max( err, erri )
4863 IF( err*sqrt( eps ).GE.one )
4864 $ info = 1
4865 iia = iia + 1
4866 END IF
4867
4868 ioffa = ioffa + 1
4869
4870 40 CONTINUE
4871
4872 icurrow = mod( icurrow+1, nprow )
4873
4874 DO 60 i = in+1, ia+m-1, desca( mb_ )
4875 ib =
min( ia+m-i, desca( mb_ ) )
4876
4877 DO 50 kk = 0, ib-1
4878
4879 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
4880 erri = abs( pa( iia+(jja-1)*ldpa )-a( ioffa ) )/eps
4881 IF( g( i+kk-ia+1 ).NE.zero )
4882 $ erri = erri / g( i+kk-ia+1 )
4883 err =
max( err, erri )
4884 IF( err*sqrt( eps ).GE.one )
4885 $ info = 1
4886 iia = iia + 1
4887 END IF
4888
4889 ioffa = ioffa + 1
4890
4891 50 CONTINUE
4892
4893 icurrow = mod( icurrow+1, nprow )
4894
4895 60 CONTINUE
4896
4897 END IF
4898
4899
4900
4901 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
4902 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
4903 $ mycol )
4904 IF( info.NE.0 )
4905 $ GO TO 80
4906
4907 70 CONTINUE
4908
4909 80 CONTINUE
4910
4911 RETURN
4912
4913
4914
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
double precision function pdlamch(ictxt, cmach)