Skip to content

Commit 5937f90

Browse files
hsnyderhsnyder
committed
Some changes in accorance with preliminary review comments
1 parent 56822b5 commit 5937f90

File tree

3 files changed

+89
-83
lines changed

3 files changed

+89
-83
lines changed

‎src/stdlib_functions.f90

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -13,21 +13,21 @@ module stdlib_functions
1313
!! version: experimental
1414
!!
1515
!! Legendre polynomial
16-
pure elemental module function legendre_fp64(N,x) result(L)
17-
integer, intent(in) :: N
16+
pure elemental module function legendre_fp64(n,x) result(leg)
17+
integer, intent(in) :: n
1818
real(dp), intent(in) :: x
19-
real(dp) :: L
19+
real(dp) :: leg
2020
end function
2121
end interface
2222

2323
interface dlegendre
2424
!! version: experimental
2525
!!
2626
!! First derivative Legendre polynomial
27-
pure elemental module function dlegendre_fp64(N,x) result(DL)
28-
integer, intent(in) :: N
27+
pure elemental module function dlegendre_fp64(n,x) result(dleg)
28+
integer, intent(in) :: n
2929
real(dp), intent(in) :: x
30-
real(dp) :: DL
30+
real(dp) :: dleg
3131
end function
3232
end interface
3333

‎src/stdlib_functions_legendre.f90

Lines changed: 36 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -4,65 +4,65 @@
44

55
contains
66

7-
! derivatives of Legendre polynomials
8-
! unspecified behaviour if N is negative
9-
pure elemental module function dlegendre_fp64(N,x) result(DL)
10-
integer, intent(in) :: N
7+
! derivatives of legegendre polynomials
8+
! unspecified behaviour if n is negative
9+
pure elemental module function dlegendre_fp64(n,x) result(dleg)
10+
integer, intent(in) :: n
1111
real(dp), intent(in) :: x
12-
real(dp) :: DL
12+
real(dp) :: dleg
1313

14-
select case(N)
14+
select case(n)
1515
case(0)
16-
DL = 0
16+
dleg = 0
1717
case(1)
18-
DL = 1
18+
dleg = 1
1919
case default
2020
block
21-
real(dp) :: L_down1, L_down2, L
22-
real(dp) :: DL_down1, DL_down2
21+
real(dp) :: leg_down1, leg_down2, leg
22+
real(dp) :: dleg_down1, dleg_down2
2323
integer :: i
2424

25-
L_down1 = x
26-
DL_down1 = 1
25+
leg_down1 = x
26+
dleg_down1 = 1
2727

28-
L_down2 = 1
29-
DL_down2 = 0
28+
leg_down2 = 1
29+
dleg_down2 = 0
3030

31-
do i = 2, N
32-
L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
33-
DL = DL_down2 + (2*i-1)*L_down1
34-
L_down2 = L_down1
35-
L_down1 = L
36-
DL_down2 = DL_down1
37-
DL_down1 = DL
31+
do i = 2, n
32+
leg = (2*i-1)*x*leg_down1/i - (i-1)*leg_down2/i
33+
dleg = dleg_down2 + (2*i-1)*leg_down1
34+
leg_down2 = leg_down1
35+
leg_down1 = leg
36+
dleg_down2 = dleg_down1
37+
dleg_down1 = dleg
3838
end do
3939
end block
4040
end select
4141
end function
4242

43-
! Legendre polynomials
44-
! unspecified behaviour if N is negative
45-
pure elemental module function legendre_fp64(N,x) result(L)
46-
integer, intent(in) :: N
43+
! legegendre polynomials
44+
! unspecified behaviour if n is negative
45+
pure elemental module function legendre_fp64(n,x) result(leg)
46+
integer, intent(in) :: n
4747
real(dp), intent(in) :: x
48-
real(dp) :: L
49-
select case(N)
48+
real(dp) :: leg
49+
select case(n)
5050
case(0)
51-
L = 1
51+
leg = 1
5252
case(1)
53-
L = x
53+
leg = x
5454
case default
5555
block
56-
real(dp) :: L_down1, L_down2
56+
real(dp) :: leg_down1, leg_down2
5757
integer :: i
5858

59-
L_down1 = x
60-
L_down2 = 1
59+
leg_down1 = x
60+
leg_down2 = 1
6161

62-
do i = 2, N
63-
L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
64-
L_down2 = L_down1
65-
L_down1 = L
62+
do i = 2, n
63+
leg = (2*i-1)*x*leg_down1/i - (i-1)*leg_down2/i
64+
leg_down2 = leg_down1
65+
leg_down1 = leg
6666
end do
6767
end block
6868
end select

‎src/stdlib_quadrature_gauss.f90

Lines changed: 47 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
use stdlib_functions, only: legendre, dlegendre
44
implicit none
55

6-
real(dp), parameter :: PI = acos(-1._dp)
6+
real(dp), parameter :: pi = acos(-1._dp)
77
real(dp), parameter :: tolerance = 4._dp * epsilon(1._dp)
88
integer, parameter :: newton_iters = 100
99

@@ -13,49 +13,52 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
1313
real(dp), intent(out) :: x(:), w(:)
1414
real(dp), intent(in), optional :: interval(2)
1515

16-
associate (N => size(x)-1 )
17-
select case (N)
16+
associate (n => size(x)-1 )
17+
select case (n)
1818
case (0)
19-
x = 0._dp
20-
w = 2._dp
19+
x = 0
20+
w = 2
2121
case (1)
22-
x = [-sqrt(1._dp/3._dp), sqrt(1._dp/3._dp)]
23-
w = [1._dp, 1._dp]
22+
x(1) = -sqrt(1._dp/3._dp)
23+
x(2) = -x(1)
24+
w = 1
2425
case default
2526
block
2627
integer :: i,j
2728
real(dp) :: leg, dleg, delta
2829

29-
do i = 0, int(floor((N+1)/2._dp)-1)
30-
x(i+1) = -cos((2*i+1)/(2._dp*N+2._dp) * PI)
30+
do i = 0, (n+1)/2 - 1
31+
x(i+1) = -cos((2*i+1)/(2._dp*n+2._dp) * pi)
3132
do j = 0, newton_iters-1
32-
leg = legendre(N+1,x(i+1))
33-
dleg = dlegendre(N+1,x(i+1))
33+
leg = legendre(n+1,x(i+1))
34+
dleg = dlegendre(n+1,x(i+1))
3435
delta = -leg/dleg
3536
x(i+1) = x(i+1) + delta
3637
if ( abs(delta) <= tolerance * abs(x(i+1)) ) exit
3738
end do
38-
x(N-i+1) = -x(i+1)
39+
x(n-i+1) = -x(i+1)
3940

40-
dleg = dlegendre(N+1,x(i+1))
41+
dleg = dlegendre(n+1,x(i+1))
4142
w(i+1) = 2._dp/((1-x(i+1)**2)*dleg**2)
42-
w(N-i+1) = w(i+1)
43+
w(n-i+1) = w(i+1)
4344
end do
4445

45-
if (mod(N,2) == 0) then
46-
x(N/2+1) = 0.0
46+
if (mod(n,2) == 0) then
47+
x(n/2+1) = 0
4748

48-
dleg = dlegendre(N+1, 0.0_dp)
49-
w(N/2+1) = 2._dp/(dleg**2)
49+
dleg = dlegendre(n+1, 0.0_dp)
50+
w(n/2+1) = 2._dp/(dleg**2)
5051
end if
5152
end block
5253
end select
5354
end associate
5455

5556
if (present(interval)) then
5657
associate ( a => interval(1) , b => interval(2) )
57-
x = 0.5*(b-a)*x+0.5*(b+a)
58-
w = 0.5*(b-a)*w
58+
x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
59+
x(1) = interval(1)
60+
x(size(x)) = interval(2)
61+
w = 0.5_dp*(b-a)*w
5962
end associate
6063
end if
6164
end subroutine
@@ -64,51 +67,54 @@ pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
6467
real(dp), intent(out) :: x(:), w(:)
6568
real(dp), intent(in), optional :: interval(2)
6669

67-
associate (N => size(x)-1)
68-
select case (N)
70+
associate (n => size(x)-1)
71+
select case (n)
6972
case (1)
70-
x = [-1._dp, 1._dp]
71-
w = [ 1._dp, 1._dp]
73+
x(1) = -1
74+
x(2) = 1
75+
w = 1
7276
case default
7377
block
7478
integer :: i,j
7579
real(dp) :: leg, dleg, delta
7680

7781
x(1) = -1._dp
78-
x(N+1) = 1._dp
79-
w(1) = 2._dp/(N*(N+1._dp))
80-
w(N+1) = 2._dp/(N*(N+1._dp))
82+
x(n+1) = 1._dp
83+
w(1) = 2._dp/(n*(n+1._dp))
84+
w(n+1) = 2._dp/(n*(n+1._dp))
8185

82-
do i = 1, int(floor((N+1)/2._dp)-1)
83-
x(i+1) = -cos( (i+0.25_dp)*PI/N - 3/(8*N*PI*(i+0.25_dp)))
86+
do i = 1, (n+1)/2 - 1
87+
x(i+1) = -cos( (i+0.25_dp)*pi/n - 3/(8*n*pi*(i+0.25_dp)))
8488
do j = 0, newton_iters-1
85-
leg = legendre(N+1,x(i+1)) - legendre(N-1,x(i+1))
86-
dleg = dlegendre(N+1,x(i+1)) - dlegendre(N-1,x(i+1))
89+
leg = legendre(n+1,x(i+1)) - legendre(n-1,x(i+1))
90+
dleg = dlegendre(n+1,x(i+1)) - dlegendre(n-1,x(i+1))
8791
delta = -leg/dleg
8892
x(i+1) = x(i+1) + delta
8993
if ( abs(delta) <= tolerance * abs(x(i+1)) ) exit
9094
end do
91-
x(N-i+1) = -x(i+1)
95+
x(n-i+1) = -x(i+1)
9296

93-
leg = legendre(N, x(i+1))
94-
w(i+1) = 2._dp/(N*(N+1._dp)*leg**2)
95-
w(N-i+1) = w(i+1)
97+
leg = legendre(n, x(i+1))
98+
w(i+1) = 2._dp/(n*(n+1._dp)*leg**2)
99+
w(n-i+1) = w(i+1)
96100
end do
97101

98-
if (mod(N,2) == 0) then
99-
x(N/2+1) = 0.0
102+
if (mod(n,2) == 0) then
103+
x(n/2+1) = 0
100104

101-
leg = legendre(N, 0.0_dp)
102-
w(N/2+1) = 2._dp/(N*(N+1._dp)*leg**2)
105+
leg = legendre(n, 0.0_dp)
106+
w(n/2+1) = 2._dp/(n*(n+1._dp)*leg**2)
103107
end if
104108
end block
105109
end select
106110
end associate
107111

108112
if (present(interval)) then
109113
associate ( a => interval(1) , b => interval(2) )
110-
x = 0.5*(b-a)*x+0.5*(b+a)
111-
w = 0.5*(b-a)*w
114+
x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
115+
x(1) = interval(1)
116+
x(size(x)) = interval(2)
117+
w = 0.5_dp*(b-a)*w
112118
end associate
113119
end if
114120
end subroutine

0 commit comments

Comments
 (0)