Skip to content

Commit 9b96f89

Browse files
awvwgkawvwgk
authored
Merge pull request #34 from awvwgk/example
Include additional examples from certik/minpack
2 parents 25a689c + 2c5bdc0 commit 9b96f89

File tree

7 files changed

+461
-0
lines changed

7 files changed

+461
-0
lines changed

‎.github/workflows/CI.yml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,9 @@ jobs:
5050
- name: Run test
5151
run: fpm test
5252

53+
- name: Run examples
54+
run: fpm run --example --all
55+
5356
Docs:
5457
runs-on: ${{ matrix.os }}
5558
strategy:

‎examples/example_hybrd.f90‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
!> The problem is to determine the values of x(1), x(2), ..., x(9)
2+
!> which solve the system of tridiagonal equations
3+
!> (3-2*x(1))*x(1) -2*x(2) = -1
4+
!> -x(i-1) + (3-2*x(i))*x(i) -2*x(i+1) = -1, i=2-8
5+
!> -x(8) + (3-2*x(9))*x(9) = -1
6+
program example_hybrd
7+
8+
use minpack_module, only: hybrd, enorm, dpmpar
9+
implicit none
10+
integer j, n, maxfev, ml, mu, mode, nprint, info, nfev, ldfjac, lr, nwrite
11+
double precision xtol, epsfcn, factor, fnorm
12+
double precision x(9), fvec(9), diag(9), fjac(9, 9), r(45), qtf(9), &
13+
wa1(9), wa2(9), wa3(9), wa4(9)
14+
15+
!> Logical output unit is assumed to be number 6.
16+
data nwrite/6/
17+
18+
n = 9
19+
20+
!> The following starting values provide a rough solution.
21+
do j = 1, 9
22+
x(j) = -1.0d0
23+
end do
24+
25+
ldfjac = 9
26+
lr = 45
27+
28+
!> Set xtol to the square root of the machine precision.
29+
!> unless high precision solutions are required,
30+
!> this is the recommended setting.
31+
xtol = dsqrt(dpmpar(1))
32+
33+
maxfev = 2000
34+
ml = 1
35+
mu = 1
36+
epsfcn = 0.0d0
37+
mode = 2
38+
do j = 1, 9
39+
diag(j) = 1.0d0
40+
end do
41+
factor = 1.0d2
42+
nprint = 0
43+
44+
call hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, &
45+
mode, factor, nprint, info, nfev, fjac, ldfjac, &
46+
r, lr, qtf, wa1, wa2, wa3, wa4)
47+
fnorm = enorm(n, fvec)
48+
write (nwrite, 1000) fnorm, nfev, info, (x(j), j=1, n)
49+
50+
1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
51+
5x, "NUMBER OF FUNCTION EVALUATIONS", i10// &
52+
5x, "EXIT PARAMETER", 16x, i10// &
53+
5x, "FINAL APPROXIMATE SOLUTION"//(5x, 3d15.7))
54+
55+
!> Results obtained with different compilers or machines
56+
!> may be slightly different.
57+
!>
58+
!>> FINAL L2 NORM OF THE RESIDUALS 0.1192636D-07
59+
!>>
60+
!>> NUMBER OF FUNCTION EVALUATIONS 14
61+
!>>
62+
!>> EXIT PARAMETER 1
63+
!>>
64+
!>> FINAL APPROXIMATE SOLUTION
65+
!>>
66+
!>> -0.5706545D+00 -0.6816283D+00 -0.7017325D+00
67+
!>> -0.7042129D+00 -0.7013690D+00 -0.6918656D+00
68+
!>> -0.6657920D+00 -0.5960342D+00 -0.4164121D+00
69+
70+
contains
71+
72+
!> Subroutine fcn for hybrd example.
73+
subroutine fcn(n, x, fvec, iflag)
74+
75+
implicit none
76+
integer, intent(in) :: n
77+
integer, intent(inout) :: iflag
78+
double precision, intent(in) :: x(n)
79+
double precision, intent(out) :: fvec(n)
80+
81+
integer k
82+
double precision one, temp, temp1, temp2, three, two, zero
83+
data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
84+
85+
if (iflag /= 0) go to 5
86+
87+
!! Insert print statements here when nprint is positive.
88+
89+
return
90+
5 continue
91+
do k = 1, n
92+
temp = (three - two*x(k))*x(k)
93+
temp1 = zero
94+
if (k /= 1) temp1 = x(k - 1)
95+
temp2 = zero
96+
if (k /= n) temp2 = x(k + 1)
97+
fvec(k) = temp - temp1 - two*temp2 + one
98+
end do
99+
return
100+
101+
end subroutine fcn
102+
103+
end program example_hybrd

‎examples/example_hybrd1.f90‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
!> the problem is to determine the values of x(1), x(2), ..., x(9),
2+
!> which solve the system of tridiagonal equations.
3+
!>
4+
!> (3-2*x(1))*x(1) -2*x(2) = -1
5+
!> -x(i-1) + (3-2*x(i))*x(i) -2*x(i+1) = -1, i=2-8
6+
!> -x(8) + (3-2*x(9))*x(9) = -1
7+
program example_hybrd1
8+
9+
use minpack_module, only: hybrd1, dpmpar, enorm
10+
implicit none
11+
integer j, n, info, lwa, nwrite
12+
double precision tol, fnorm
13+
double precision x(9), fvec(9), wa(180)
14+
15+
!> Logical output unit is assumed to be number 6.
16+
data nwrite/6/
17+
18+
n = 9
19+
20+
!> The following starting values provide a rough solution.
21+
do j = 1, 9
22+
x(j) = -1.d0
23+
end do
24+
25+
lwa = 180
26+
27+
!> Set tol to the square root of the machine precision.
28+
!> unless high precision solutions are required,
29+
!> this is the recommended setting.
30+
tol = dsqrt(dpmpar(1))
31+
32+
call hybrd1(fcn, n, x, fvec, tol, info, wa, lwa)
33+
fnorm = enorm(n, fvec)
34+
write (nwrite, 1000) fnorm, info, (x(j), j=1, n)
35+
36+
1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
37+
5x, "EXIT PARAMETER", 16x, i10// &
38+
5x, "FINAL APPROXIMATE SOLUTION"// &
39+
(5x, 3d15.7))
40+
41+
!> Results obtained with different compilers or machines
42+
!> may be slightly different.
43+
!>
44+
!>> FINAL L2 NORM OF THE RESIDUALS 0.1192636D-07
45+
!>>
46+
!>> EXIT PARAMETER 1
47+
!>>
48+
!>> FINAL APPROXIMATE SOLUTION
49+
!>>
50+
!>> -0.5706545D+00 -0.6816283D+00 -0.7017325D+00
51+
!>> -0.7042129D+00 -0.7013690D+00 -0.6918656D+00
52+
!>> -0.6657920D+00 -0.5960342D+00 -0.4164121D+00
53+
54+
contains
55+
56+
!> Subroutine fcn for hybrd1 example.
57+
subroutine fcn(n, x, fvec, iflag)
58+
59+
implicit none
60+
integer, intent(in) :: n
61+
integer, intent(inout) :: iflag
62+
double precision, intent(in) :: x(n)
63+
double precision, intent(out) :: fvec(n)
64+
65+
integer k
66+
double precision one, temp, temp1, temp2, three, two, zero
67+
data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
68+
69+
do k = 1, n
70+
temp = (three - two*x(k))*x(k)
71+
temp1 = zero
72+
if (k /= 1) temp1 = x(k - 1)
73+
temp2 = zero
74+
if (k /= n) temp2 = x(k + 1)
75+
fvec(k) = temp - temp1 - two*temp2 + one
76+
end do
77+
78+
end subroutine fcn
79+
80+
end program example_hybrd1

‎examples/example_lmder1.f90‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
module testmod_der1
2+
implicit none
3+
private
4+
public fcn, dp
5+
6+
integer, parameter :: dp=kind(0d0)
7+
8+
contains
9+
10+
subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag)
11+
integer, intent(in) :: m, n, ldfjac
12+
integer, intent(inout) :: iflag
13+
real(dp), intent(in) :: x(n)
14+
real(dp), intent(inout) :: fvec(m), fjac(ldfjac, n)
15+
16+
integer :: i
17+
real(dp) :: tmp1, tmp2, tmp3, tmp4, y(15)
18+
y = [1.4D-1, 1.8D-1, 2.2D-1, 2.5D-1, 2.9D-1, 3.2D-1, 3.5D-1, 3.9D-1, &
19+
3.7D-1, 5.8D-1, 7.3D-1, 9.6D-1, 1.34D0, 2.1D0, 4.39D0]
20+
21+
if (iflag == 1) then
22+
do i = 1, 15
23+
tmp1 = i
24+
tmp2 = 16 - i
25+
tmp3 = tmp1
26+
if (i .gt. 8) tmp3 = tmp2
27+
fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
28+
end do
29+
else
30+
do i = 1, 15
31+
tmp1 = i
32+
tmp2 = 16 - i
33+
tmp3 = tmp1
34+
if (i .gt. 8) tmp3 = tmp2
35+
tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
36+
fjac(i,1) = -1.D0
37+
fjac(i,2) = tmp1*tmp2/tmp4
38+
fjac(i,3) = tmp1*tmp3/tmp4
39+
end do
40+
end if
41+
end subroutine
42+
43+
end module
44+
45+
46+
program example_lmder1
47+
use minpack_module, only: enorm, lmder1, chkder
48+
use testmod_der1, only: dp, fcn
49+
implicit none
50+
51+
integer :: info
52+
real(dp) :: tol, x(3), fvec(15), fjac(size(fvec), size(x))
53+
integer :: ipvt(size(x))
54+
real(dp), allocatable :: wa(:)
55+
56+
! The following starting values provide a rough fit.
57+
x = [1._dp, 1._dp, 1._dp]
58+
59+
call check_deriv()
60+
61+
! Set tol to the square root of the machine precision. Unless high precision
62+
! solutions are required, this is the recommended setting.
63+
tol = sqrt(epsilon(1._dp))
64+
65+
allocate(wa(5*size(x) + size(fvec)))
66+
call lmder1(fcn, size(fvec), size(x), x, fvec, fjac, size(fjac, 1), tol, &
67+
info, ipvt, wa, size(wa))
68+
print 1000, enorm(size(fvec), fvec), info, x
69+
1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
70+
5x, 'EXIT PARAMETER', 16x, i10 // &
71+
5x, 'FINAL APPROXIMATE SOLUTION' // &
72+
5x, 3d15.7)
73+
74+
contains
75+
76+
subroutine check_deriv()
77+
integer :: iflag
78+
real(dp) :: xp(size(x)), fvecp(size(fvec)), err(size(fvec))
79+
call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, &
80+
1, err)
81+
iflag = 1
82+
call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag)
83+
iflag = 2
84+
call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag)
85+
iflag = 1
86+
call fcn(size(fvec), size(x), xp, fvecp, fjac, size(fjac, 1), iflag)
87+
call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, &
88+
2, err)
89+
print *, "Derivatives check (1.0 is correct, 0.0 is incorrect):"
90+
print *, err
91+
end subroutine
92+
93+
end program

‎examples/example_lmdif1.f90‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
module testmod_dif1
2+
implicit none
3+
private
4+
public fcn, dp
5+
6+
integer, parameter :: dp=kind(0d0)
7+
8+
contains
9+
10+
subroutine fcn(m, n, x, fvec, iflag)
11+
integer, intent(in) :: m, n
12+
integer, intent(inout) :: iflag
13+
real(dp), intent(in) :: x(n)
14+
real(dp), intent(out) :: fvec(m)
15+
16+
integer :: i
17+
real(dp) :: tmp1, tmp2, tmp3, y(15)
18+
! Suppress compiler warning:
19+
y(1) = iflag
20+
y = [1.4D-1, 1.8D-1, 2.2D-1, 2.5D-1, 2.9D-1, 3.2D-1, 3.5D-1, 3.9D-1, &
21+
3.7D-1, 5.8D-1, 7.3D-1, 9.6D-1, 1.34D0, 2.1D0, 4.39D0]
22+
23+
do i = 1, 15
24+
tmp1 = i
25+
tmp2 = 16 - i
26+
tmp3 = tmp1
27+
if (i .gt. 8) tmp3 = tmp2
28+
fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
29+
end do
30+
end subroutine
31+
32+
end module
33+
34+
35+
program example_lmdif1
36+
use minpack_module, only: enorm, lmdif1
37+
use testmod_dif1, only: dp, fcn
38+
implicit none
39+
40+
integer :: info, m, n
41+
real(dp) :: tol, x(3), fvec(15)
42+
integer :: iwa(size(x))
43+
real(dp), allocatable :: wa(:)
44+
45+
! The following starting values provide a rough fit.
46+
x = [1._dp, 1._dp, 1._dp]
47+
48+
! Set tol to the square root of the machine precision. Unless high precision
49+
! solutions are required, this is the recommended setting.
50+
tol = sqrt(epsilon(1._dp))
51+
52+
m = size(fvec)
53+
n = size(x)
54+
allocate(wa(m*n + 5*n + m))
55+
call lmdif1(fcn, size(fvec), size(x), x, fvec, tol, info, iwa, wa, size(wa))
56+
print 1000, enorm(size(fvec), fvec), info, x
57+
1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
58+
5x, 'EXIT PARAMETER', 16x, i10 // &
59+
5x, 'FINAL APPROXIMATE SOLUTION' // &
60+
5x, 3d15.7)
61+
end program

0 commit comments

Comments
 (0)