gcc/testsuite/gfortran.dg/graphite/run-id-2.f90

   1   IMPLICIT NONE
   2   INTEGER, PARAMETER :: dp=KIND(0.0D0)
   3   REAL(KIND=dp)      :: res
   4
   5   res=exp_radius_very_extended(  0    ,      1   ,       0      ,    1, &
   6                                (/0.0D0,0.0D0,0.0D0/),&
   7                                (/1.0D0,0.0D0,0.0D0/),&
   8                                (/1.0D0,0.0D0,0.0D0/),&
   9                                  1.0D0,1.0D0,1.0D0,1.0D0)
  10   if (res.ne.1.0d0) call abort()
  11
  12 CONTAINS
  13
  14  FUNCTION exp_radius_very_extended(la_min,la_max,lb_min,lb_max,ra,rb,rp,&
  15                           zetp,eps,prefactor,cutoff) RESULT(radius)
  16
  17     INTEGER, INTENT(IN)                      :: la_min, la_max, lb_min, lb_max
  18     REAL(KIND=dp), INTENT(IN)                :: ra(3), rb(3), rp(3), zetp, &
  19                                                 eps, prefactor, cutoff
  20     REAL(KIND=dp)                            :: radius
  21
  22     INTEGER                                  :: i, ico, j, jco, la(3), lb(3), &
  23                                                 lxa, lxb, lya, lyb, lza, lzb
  24     REAL(KIND=dp)                            :: bini, binj, coef(0:20), &
  25                                                 epsin_local, polycoef(0:60), &
  26                                                 prefactor_local, rad_a, &
  27                                                 rad_b, s1, s2
  28
  29     epsin_local=1.0E-2_dp
  30
  31     prefactor_local=prefactor*MAX(1.0_dp,cutoff)
  32     rad_a=SQRT(SUM((ra-rp)**2))
  33     rad_b=SQRT(SUM((rb-rp)**2))
  34
  35     polycoef(0:la_max+lb_max)=0.0_dp
  36     DO lxa=0,la_max
  37     DO lxb=0,lb_max
  38        coef(0:la_max+lb_max)=0.0_dp
  39        bini=1.0_dp
  40        s1=1.0_dp
  41        DO i=0,lxa
  42           binj=1.0_dp
  43           s2=1.0_dp
  44           DO j=0,lxb
  45              coef(lxa+lxb-i-j)=coef(lxa+lxb-i-j) + bini*binj*s1*s2
  46              binj=(binj*(lxb-j))/(j+1)
  47              s2=s2*(rad_b)
  48           ENDDO
  49           bini=(bini*(lxa-i))/(i+1)
  50           s1=s1*(rad_a)
  51        ENDDO
  52        DO i=0,lxa+lxb
  53           polycoef(i)=MAX(polycoef(i),coef(i))
  54        ENDDO
  55     ENDDO
  56     ENDDO
  57
  58     polycoef(0:la_max+lb_max)=polycoef(0:la_max+lb_max)*prefactor_local
  59     radius=0.0_dp
  60     DO i=0,la_max+lb_max
  61           radius=MAX(radius,polycoef(i)**(i+1))
  62     ENDDO
  63
  64   END FUNCTION exp_radius_very_extended
  65
  66 END