Skip to content

Commit e45e3e5

Browse files
scivisionscivision
committed
more-examples
1 parent e2cf068 commit e45e3e5

File tree

4 files changed

+197
-7
lines changed

4 files changed

+197
-7
lines changed

‎.github/workflows/cmake.yml‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,8 @@ jobs:
5151
FC: gfortran-${{ matrix.fc }}
5252

5353
steps:
54-
- uses: actions/checkout@v4
54+
- &checkout
55+
uses: actions/checkout@v4
5556

5657
- run: cmake --workflow debug
5758
- run: cmake --workflow release
@@ -71,12 +72,13 @@ jobs:
7172
FC: flang-new
7273

7374
steps:
74-
- uses: actions/checkout@v4
75-
7675
- name: install Flang
7776
run: sudo apt install --no-install-recommends flang
7877

79-
- run: cmake --workflow default
78+
- *checkout
79+
80+
- &default
81+
run: cmake --workflow default
8082

8183

8284
windows:
@@ -94,7 +96,6 @@ jobs:
9496
- name: Put MSYS2_MinGW64 on PATH
9597
run: echo "${{ steps.msys2.outputs.msys2-location }}/ucrt64/bin" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append
9698

97-
- uses: actions/checkout@v4
99+
- *checkout
98100

99-
- run: cmake --workflow debug
100-
- run: cmake --workflow release
101+
- *default

‎test/standard/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,8 +21,12 @@ endforeach()
2121

2222
if(f18abstract)
2323
add_executable(abstract_interface abstract.f90)
24+
add_executable(abstract_procedure abstract_procedure.f90)
2425
endif()
2526

27+
add_executable(lubounds bounds.f90)
28+
add_test(NAME LUbounds COMMAND lubounds)
29+
2630
add_executable(shortcircuit ${PROJECT_SOURCE_DIR}/app/standard/short_circuit.f90)
2731
# shaky result, don't CI test
2832

‎test/standard/abstract_procedure.f90‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
module newton
2+
3+
use, intrinsic :: iso_fortran_env, only: wp => real64, stderr => error_unit
4+
5+
implicit none (type, external)
6+
7+
private
8+
public :: wp, newtopts, newton_exact, objfun, objfun_deriv
9+
10+
!> derived type for containing the options for the newton method procedure, by default
11+
! these work okay with the dipole to spherical conversion problem but can be adjusted
12+
! by the user for other applications
13+
type :: newtopts
14+
real(wp) :: derivtol=1e-18
15+
integer :: maxit=100
16+
real(wp) :: tol=1e-9
17+
logical :: verbose=.false.
18+
end type newtopts
19+
20+
!> these interfaced define the types of functions needed to run Newton iterations
21+
! need to match custom functions defined in using modules; these are defined so
22+
! that there are multiple objective functions (conforming to these patterns) that
23+
! can be used with Newton's method. These are abstract because there are multiple
24+
! objective functions (in principle) that can be used with Newton's method
25+
abstract interface
26+
real(wp) function objfun(x,parms)
27+
import wp
28+
real(wp), intent(in) :: x
29+
real(wp), dimension(:), intent(in) :: parms
30+
end function
31+
32+
real(wp) function objfun_deriv(x,parms)
33+
import wp
34+
real(wp), intent(in) :: x
35+
real(wp), dimension(:), intent(in) :: parms
36+
end function
37+
end interface
38+
39+
contains
40+
41+
!> this implmements the exact Newton method for solving a nonlinear equation
42+
subroutine newton_exact(f,fprime,x0,parms,newtparms,root,it,converged)
43+
procedure(objfun), pointer, intent(in) :: f
44+
procedure(objfun_deriv), pointer, intent(in) :: fprime
45+
real(wp), intent(in) :: x0 ! starting point for newton iteration
46+
real(wp),dimension(:), intent(in) :: parms ! fixed parameters of the newton iteration, f,fprime must accommodate whatever size array is passed in
47+
type(newtopts), intent(in) :: newtparms ! options for the iteration that can be set by the user
48+
real(wp), intent(out) :: root
49+
integer, intent(out) :: it
50+
logical, intent(out) :: converged
51+
52+
real(wp) :: fval,derivative
53+
54+
! check starting point is not too close to inflection
55+
if (abs(fprime(x0,parms))<newtparms%derivtol) then
56+
write(stderr,*) 'Warning: starting near inflection point, please change initial guess!'
57+
it=0; converged=.false.; root=x0;
58+
return
59+
end if
60+
61+
! Newton iteration main loop
62+
it=1; root=x0; fval=f(root,parms); converged=.false.
63+
do while (.not. converged .and. it <= newtparms%maxit)
64+
derivative=fprime(root,parms)
65+
if (abs(derivative)<newtparms%derivtol) then
66+
write(stderr,*) 'Warning: Encountered inflection point during iteration: ',it
67+
return
68+
else
69+
root=root-fval/derivative
70+
fval=f(root,parms)
71+
if (newtparms%verbose) then
72+
print*, ' Iteration ',it,'; root ',root,' fval ',fval,' derivative ',derivative
73+
end if
74+
it=it+1
75+
converged=abs(fval)<newtparms%tol
76+
end if
77+
end do
78+
it=it-1
79+
80+
end subroutine newton_exact
81+
82+
end module newton
83+
84+
85+
program main
86+
87+
use newton, only: wp, newtopts, objfun, objfun_deriv, newton_exact
88+
89+
implicit none
90+
91+
real(wp) :: q,p,r
92+
93+
real(wp), parameter :: Re=6371e3_wp ! Earth radius in meters
94+
95+
type(newtopts) :: newtparms
96+
real(wp), dimension(2) :: parms
97+
real(wp) :: r0
98+
procedure(objfun), pointer :: f
99+
procedure(objfun_deriv), pointer :: fprime
100+
integer :: maxrestart, maxr, r0step
101+
integer :: it,ir0
102+
logical :: converged
103+
104+
! Set parameters of the restart and Newton iterations
105+
maxrestart=400
106+
maxr=100*Re
107+
r0step=0.25*Re
108+
newtparms%maxit=100
109+
newtparms%derivtol=1e-18
110+
newtparms%tol=1e-11
111+
newtparms%verbose=.false.
112+
f=>rpoly
113+
fprime=>rpoly_deriv
114+
parms=[q,p]
115+
116+
! Newton iterations with restarting (see parameters above for limits) until we get a satisfactory result
117+
r=0; converged=.false.; ir0=1;
118+
do while (.not. converged .and. ir0<maxrestart .and. (r<=0 .or. r>maxr))
119+
r0=(ir0-1)*(r0step) ! change starting point in increments of 0.25 Re until we get a "good" answer
120+
call newton_exact(f,fprime,r0,parms,newtparms,r,it,converged)
121+
ir0=ir0+1
122+
end do
123+
124+
print '(a,2f20.10)', ' r, Re = ', r, Re
125+
print '(a,i0)', ' iterations = ', it
126+
print '(a,l1)', ' converged = ', converged
127+
128+
contains
129+
!> objective function for newton iterations for solutions of roots for r
130+
function rpoly(x,parms) result(fval)
131+
real(wp), intent(in) :: x
132+
real(wp), dimension(:), intent(in) :: parms
133+
real(wp) :: fval
134+
135+
real(wp) :: q,p
136+
137+
q=parms(1); p=parms(2);
138+
fval=q**2*(x/Re)**4 + 1/p*(x/Re) - 1
139+
end function rpoly
140+
141+
142+
!> derivative objective function for newton iterations for roots of r
143+
function rpoly_deriv(x,parms) result(fval_deriv)
144+
real(wp), intent(in) :: x
145+
real(wp), dimension(:), intent(in) :: parms
146+
real(wp) :: fval_deriv
147+
148+
real(wp) :: q,p
149+
150+
q=parms(1); p=parms(2);
151+
fval_deriv=4/Re*q**2*(x/Re)**3 + 1/p/Re
152+
end function rpoly_deriv
153+
154+
155+
end program

‎test/standard/bounds.f90‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
program test_bounds
2+
use iso_fortran_env, only : error_unit
3+
implicit none
4+
5+
real, allocatable :: a(:)
6+
integer :: L1,L2, U1,U2
7+
allocate(a(1:2))
8+
9+
L1 = lbound(a,1)
10+
U1 = ubound(a,1)
11+
12+
call c_bounder(a)
13+
14+
L2 = lbound(a,1)
15+
U2 = ubound(a,1)
16+
17+
if (L1 /= L2 .or. U1 /= U2) then
18+
write(error_unit, '(a,2i2,a,2i2)') 'FAIL: bounds changed before/after lower:', L1,L2, " upper: ", U1,U2
19+
error stop
20+
endif
21+
22+
print '(a)', "bounds check OK"
23+
24+
contains
25+
26+
subroutine c_bounder(a) bind(c)
27+
real, intent(inout) :: a(:)
28+
end subroutine c_bounder
29+
30+
end program

0 commit comments

Comments
 (0)