-
Notifications
You must be signed in to change notification settings - Fork 20
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #63 from jacobwilliams/62-modernize-examples
Modernize the examples
- Loading branch information
Showing
5 changed files
with
243 additions
and
253 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,93 +1,95 @@ | ||
| module testmod_der1 | ||
| implicit none | ||
| private | ||
| public fcn, dp | ||
|
|
||
| integer, parameter :: dp=kind(0d0) | ||
|
|
||
| contains | ||
|
|
||
| subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag) | ||
| integer, intent(in) :: m, n, ldfjac | ||
| integer, intent(inout) :: iflag | ||
| real(dp), intent(in) :: x(n) | ||
| real(dp), intent(inout) :: fvec(m), fjac(ldfjac, n) | ||
|
|
||
| integer :: i | ||
| real(dp) :: tmp1, tmp2, tmp3, tmp4, y(15) | ||
| 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, & | ||
| 3.7D-1, 5.8D-1, 7.3D-1, 9.6D-1, 1.34D0, 2.1D0, 4.39D0] | ||
|
|
||
| if (iflag == 1) then | ||
| do i = 1, 15 | ||
| tmp1 = i | ||
| tmp2 = 16 - i | ||
| tmp3 = tmp1 | ||
| if (i > 8) tmp3 = tmp2 | ||
| fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3)) | ||
| end do | ||
| else | ||
| do i = 1, 15 | ||
| tmp1 = i | ||
| tmp2 = 16 - i | ||
| tmp3 = tmp1 | ||
| if (i > 8) tmp3 = tmp2 | ||
| tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2 | ||
| fjac(i,1) = -1.D0 | ||
| fjac(i,2) = tmp1*tmp2/tmp4 | ||
| fjac(i,3) = tmp1*tmp3/tmp4 | ||
| end do | ||
| end if | ||
| end subroutine | ||
|
|
||
| end module | ||
| program example_lmder1 | ||
|
|
||
| use minpack_module, only: wp, enorm, lmder1, chkder | ||
| use iso_fortran_env, only: nwrite => output_unit | ||
|
|
||
| program example_lmder1 | ||
| use minpack_module, only: enorm, lmder1, chkder | ||
| use testmod_der1, only: dp, fcn | ||
| implicit none | ||
|
|
||
| integer, parameter :: n = 3 | ||
| integer, parameter :: m = 15 | ||
| integer, parameter :: lwa = 5*n+m | ||
|
|
||
| integer :: info | ||
| real(dp) :: tol, x(3), fvec(15), fjac(size(fvec), size(x)) | ||
| integer :: ipvt(size(x)) | ||
| real(dp), allocatable :: wa(:) | ||
| real(wp) :: tol, x(n), fvec(m), fjac(m,n) | ||
| integer :: ipvt(n) | ||
| real(wp) :: wa(lwa) | ||
|
|
||
| ! The following starting values provide a rough fit. | ||
| x = [1._dp, 1._dp, 1._dp] | ||
| x = [1.0_wp, 1.0_wp, 1.0_wp] | ||
|
|
||
| call check_deriv() | ||
|
|
||
| ! Set tol to the square root of the machine precision. Unless high precision | ||
| ! solutions are required, this is the recommended setting. | ||
| tol = sqrt(epsilon(1._dp)) | ||
| tol = sqrt(epsilon(1._wp)) | ||
|
|
||
| allocate(wa(5*size(x) + size(fvec))) | ||
| call lmder1(fcn, size(fvec), size(x), x, fvec, fjac, size(fjac, 1), tol, & | ||
| info, ipvt, wa, size(wa)) | ||
| print 1000, enorm(size(fvec), fvec), info, x | ||
| 1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // & | ||
| 5x, 'EXIT PARAMETER', 16x, i10 // & | ||
| 5x, 'FINAL APPROXIMATE SOLUTION' // & | ||
| 5x, 3d15.7) | ||
| call lmder1(fcn, m, n, x, fvec, fjac, m, tol, info, ipvt, wa, lwa) | ||
|
|
||
| write(nwrite, '(5x,a,d15.7//,5x,a,16x,i10//,5x,a//(5x,3d15.7))') & | ||
| 'FINAL L2 NORM OF THE RESIDUALS', enorm(m, fvec), & | ||
| 'EXIT PARAMETER', info, & | ||
| 'FINAL APPROXIMATE SOLUTION', x | ||
|
|
||
| contains | ||
|
|
||
| subroutine check_deriv() | ||
| integer :: iflag | ||
| real(dp) :: xp(size(x)), fvecp(size(fvec)), err(size(fvec)) | ||
| call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, & | ||
| 1, err) | ||
| iflag = 1 | ||
| call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag) | ||
| iflag = 2 | ||
| call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag) | ||
| iflag = 1 | ||
| call fcn(size(fvec), size(x), xp, fvecp, fjac, size(fjac, 1), iflag) | ||
| call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, & | ||
| 2, err) | ||
| print *, "Derivatives check (1.0 is correct, 0.0 is incorrect):" | ||
| print *, err | ||
| end subroutine | ||
|
|
||
| integer :: iflag | ||
| real(wp) :: xp(n), fvecp(m), err(m) | ||
|
|
||
| call chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err) | ||
| iflag = 1 | ||
| call fcn(m, n, x, fvec, fjac, m, iflag) | ||
| iflag = 2 | ||
| call fcn(m, n, x, fvec, fjac, m, iflag) | ||
| iflag = 1 | ||
| call fcn(m, n, xp, fvecp, fjac, m, iflag) | ||
| call chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err) | ||
|
|
||
| write(nwrite, '(a)') 'Derivatives check (1.0 is correct, 0.0 is incorrect):' | ||
| write(nwrite,'(1p,(5x,3d15.7))') err | ||
| if (any(abs(err-1.0_wp)>epsilon(1.0_wp))) error stop 'Derivative check failed' | ||
|
|
||
| end subroutine check_deriv | ||
|
|
||
| subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag) | ||
|
|
||
| integer, intent(in) :: m | ||
| integer, intent(in) :: n | ||
| real(wp), intent(in) :: x(n) | ||
| real(wp), intent(inout) :: fvec(m) | ||
| real(wp), intent(inout) :: fjac(ldfjac, n) | ||
| integer, intent(in) :: ldfjac | ||
| integer, intent(inout) :: iflag | ||
|
|
||
| integer :: i | ||
| real(wp) :: tmp1, tmp2, tmp3, tmp4 | ||
|
|
||
| real(wp),parameter :: y(15) = [1.4e-1_wp, 1.8e-1_wp, 2.2e-1_wp, 2.5e-1_wp, 2.9e-1_wp, & | ||
| 3.2e-1_wp, 3.5e-1_wp, 3.9e-1_wp, 3.7e-1_wp, 5.8e-1_wp, & | ||
| 7.3e-1_wp, 9.6e-1_wp, 1.34e0_wp, 2.1e0_wp, 4.39e0_wp] | ||
|
|
||
| if (iflag == 1) then | ||
| do i = 1, 15 | ||
| tmp1 = i | ||
| tmp2 = 16 - i | ||
| tmp3 = tmp1 | ||
| if (i > 8) tmp3 = tmp2 | ||
| fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3)) | ||
| end do | ||
| else | ||
| do i = 1, 15 | ||
| tmp1 = i | ||
| tmp2 = 16 - i | ||
| tmp3 = tmp1 | ||
| if (i > 8) tmp3 = tmp2 | ||
| tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2 | ||
| fjac(i,1) = -1.0_wp | ||
| fjac(i,2) = tmp1*tmp2/tmp4 | ||
| fjac(i,3) = tmp1*tmp3/tmp4 | ||
| end do | ||
| end if | ||
|
|
||
| end subroutine fcn | ||
|
|
||
| end program |
Oops, something went wrong.