Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Replace enorm with intrinsic norm2 #45

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions examples/example_hybrd.f90
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
!> -x(8) + (3-2*x(9))*x(9) = -1
program example_hybrd

use minpack_module, only: hybrd, enorm, dpmpar
use minpack_module, only: hybrd, dpmpar
implicit none
integer j, n, maxfev, ml, mu, mode, nprint, info, nfev, ldfjac, lr, nwrite
double precision xtol, epsfcn, factor, fnorm
Expand Down Expand Up @@ -44,7 +44,7 @@ program example_hybrd
call hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, &
mode, factor, nprint, info, nfev, fjac, ldfjac, &
r, lr, qtf, wa1, wa2, wa3, wa4)
fnorm = enorm(n, fvec)
fnorm = norm2(fvec)
write (nwrite, 1000) fnorm, nfev, info, (x(j), j=1, n)

1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
Expand Down
4 changes: 2 additions & 2 deletions examples/example_hybrd1.f90
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
!> -x(8) + (3-2*x(9))*x(9) = -1
program example_hybrd1

use minpack_module, only: hybrd1, dpmpar, enorm
use minpack_module, only: hybrd1, dpmpar
implicit none
integer j, n, info, lwa, nwrite
double precision tol, fnorm
Expand All @@ -30,7 +30,7 @@ program example_hybrd1
tol = dsqrt(dpmpar(1))

call hybrd1(fcn, n, x, fvec, tol, info, wa, lwa)
fnorm = enorm(n, fvec)
fnorm = norm2(fvec)
write (nwrite, 1000) fnorm, info, (x(j), j=1, n)

1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
Expand Down
4 changes: 2 additions & 2 deletions examples/example_lmder1.f90
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag)


program example_lmder1
use minpack_module, only: enorm, lmder1, chkder
use minpack_module, only: lmder1, chkder
use testmod_der1, only: dp, fcn
implicit none

Expand All @@ -65,7 +65,7 @@ program example_lmder1
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
print 1000, norm2(fvec), info, x
1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
5x, 'EXIT PARAMETER', 16x, i10 // &
5x, 'FINAL APPROXIMATE SOLUTION' // &
Expand Down
4 changes: 2 additions & 2 deletions examples/example_lmdif1.f90
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ subroutine fcn(m, n, x, fvec, iflag)


program example_lmdif1
use minpack_module, only: enorm, lmdif1
use minpack_module, only: lmdif1
use testmod_dif1, only: dp, fcn
implicit none

Expand All @@ -53,7 +53,7 @@ program example_lmdif1
n = size(x)
allocate(wa(m*n + 5*n + m))
call lmdif1(fcn, size(fvec), size(x), x, fvec, tol, info, iwa, wa, size(wa))
print 1000, enorm(size(fvec), fvec), info, x
print 1000, norm2(fvec), info, x
1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
5x, 'EXIT PARAMETER', 16x, i10 // &
5x, 'FINAL APPROXIMATE SOLUTION' // &
Expand Down
73 changes: 73 additions & 0 deletions src/enorm.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
!*****************************************************************************************
!>
! given an n-vector x, this function calculates the
! euclidean norm of x.
!
! the euclidean norm is computed by accumulating the sum of
! squares in three different sums. the sums of squares for the
! small and large components are scaled so that no overflows
! occur. non-destructive underflows are permitted. underflows
! and overflows do not occur in the computation of the unscaled
! sum of squares for the intermediate components.
! the definitions of small, intermediate and large components
! depend on two constants, rdwarf and rgiant. the main
! restrictions on these constants are that rdwarf**2 not
! underflow and rgiant**2 not overflow. the constants
! given here are suitable for every known computer.

pure real(wp) function enorm(n, x)
awvwgk marked this conversation as resolved.
Show resolved Hide resolved
use, intrinsic :: iso_fortran_env, only: wp => real64
implicit none

integer, intent(in) :: n !! a positive integer input variable.
real(wp), intent(in) :: x(n) !! an input array of length n.

integer :: i
real(wp) :: agiant, s1, s2, s3, xabs, x1max, x3max

real(wp), parameter :: rdwarf = 3.834e-20_wp
real(wp), parameter :: rgiant = 1.304e19_wp
real(wp), parameter :: one = 1.0_wp
real(wp), parameter :: zero = 0.0_wp

s1 = zero
s2 = zero
s3 = zero
x1max = zero
x3max = zero
agiant = rgiant/real(n, wp)
do i = 1, n
xabs = abs(x(i))
if (xabs > rdwarf .and. xabs < agiant) then
! sum for intermediate components.
s2 = s2 + xabs**2
elseif (xabs <= rdwarf) then
! sum for small components.
if (xabs <= x3max) then
if (xabs /= zero) s3 = s3 + (xabs/x3max)**2
else
s3 = one + s3*(x3max/xabs)**2
x3max = xabs
end if
! sum for large components.
elseif (xabs <= x1max) then
s1 = s1 + (xabs/x1max)**2
else
s1 = one + s1*(x1max/xabs)**2
x1max = xabs
end if
end do

! calculation of norm.

if (s1 /= zero) then
enorm = x1max*sqrt(s1 + (s2/x1max)/x1max)
elseif (s2 == zero) then
enorm = x3max*sqrt(s3)
else
if (s2 >= x3max) enorm = sqrt(s2*(one + (x3max/s2)*(x3max*s3)))
if (s2 < x3max) enorm = sqrt(x3max*((s2/x3max) + (x3max*s3)))
end if

end function enorm
!*****************************************************************************************
Loading