From f4f60b00c56c64215097ac235f087e30d5d81a15 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Sat, 26 Feb 2022 12:50:27 +0100 Subject: [PATCH 1/6] Replace enorm with intrinsic norm2 --- src/enorm.f90 | 71 ++++++++++++++++++++++ src/minpack.f90 | 158 +++++++++++++----------------------------------- 2 files changed, 114 insertions(+), 115 deletions(-) create mode 100644 src/enorm.f90 diff --git a/src/enorm.f90 b/src/enorm.f90 new file mode 100644 index 0000000..8bccd78 --- /dev/null +++ b/src/enorm.f90 @@ -0,0 +1,71 @@ +!***************************************************************************************** +!> +! 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) + + 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 + + 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 +!***************************************************************************************** \ No newline at end of file diff --git a/src/minpack.f90 b/src/minpack.f90 index 562d886..431a48f 100644 --- a/src/minpack.f90 +++ b/src/minpack.f90 @@ -279,7 +279,7 @@ subroutine dogleg(n, r, Lr, Diag, Qtb, Delta, x, Wa1, Wa2) Wa1(j) = zero Wa2(j) = Diag(j)*x(j) end do - qnorm = enorm(n, Wa2) + qnorm = norm2(Wa2) if (qnorm > Delta) then ! the gauss-newton direction is not acceptable. @@ -298,7 +298,7 @@ subroutine dogleg(n, r, Lr, Diag, Qtb, Delta, x, Wa1, Wa2) ! calculate the norm of the scaled gradient and test for ! the special case in which the scaled gradient is zero. - gnorm = enorm(n, Wa1) + gnorm = norm2(Wa1) sgnorm = zero alpha = Delta/qnorm if (gnorm /= zero) then @@ -318,7 +318,7 @@ subroutine dogleg(n, r, Lr, Diag, Qtb, Delta, x, Wa1, Wa2) end do Wa2(j) = sum end do - temp = enorm(n, Wa2) + temp = norm2(Wa2) sgnorm = (gnorm/temp)/temp ! test whether the scaled gradient direction is acceptable. @@ -330,7 +330,7 @@ subroutine dogleg(n, r, Lr, Diag, Qtb, Delta, x, Wa1, Wa2) ! finally, calculate the point along the dogleg ! at which the quadratic is minimized. - bnorm = enorm(n, Qtb) + bnorm = norm2(Qtb) temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/Delta) temp = temp - (Delta/qnorm)*(sgnorm/Delta)**2 + & sqrt((temp - (Delta/qnorm))**2 + & @@ -351,78 +351,6 @@ subroutine dogleg(n, r, Lr, Diag, Qtb, Delta, x, Wa1, Wa2) end subroutine dogleg !***************************************************************************************** -!***************************************************************************************** -!> -! 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) - - 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 - - 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 -!***************************************************************************************** - !***************************************************************************************** !> ! this subroutine computes a forward-difference approximation @@ -698,7 +626,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & call fcn(n, x, Fvec, iflag) Nfev = 1 if (iflag < 0) goto 300 - fnorm = enorm(n, Fvec) + fnorm = norm2(Fvec) ! determine the number of calls to fcn needed to compute ! the jacobian matrix. @@ -743,7 +671,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & do j = 1, n Wa3(j) = Diag(j)*x(j) end do - xnorm = enorm(n, Wa3) + xnorm = norm2(Wa3) delta = Factor*xnorm if (delta == zero) delta = Factor end if @@ -815,7 +743,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & Wa2(j) = x(j) + Wa1(j) Wa3(j) = Diag(j)*Wa1(j) end do - pnorm = enorm(n, Wa3) + pnorm = norm2(Wa3) ! on the first iteration, adjust the initial step bound. @@ -827,7 +755,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & call fcn(n, Wa2, Wa4, iflag) Nfev = Nfev + 1 if (iflag >= 0) then - fnorm1 = enorm(n, Wa4) + fnorm1 = norm2(Wa4) ! compute the scaled actual reduction. @@ -845,7 +773,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & end do Wa3(i) = Qtf(i) + sum end do - temp = enorm(n, Wa3) + temp = norm2(Wa3) prered = zero if (temp < fnorm) prered = one - (temp/fnorm)**2 @@ -877,7 +805,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & Wa2(j) = Diag(j)*x(j) Fvec(j) = Wa4(j) end do - xnorm = enorm(n, Wa2) + xnorm = norm2(Wa2) fnorm = fnorm1 iter = iter + 1 end if @@ -1140,7 +1068,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & call fcn(n, x, Fvec, Fjac, Ldfjac, iflag) Nfev = 1 if (iflag < 0) goto 300 - fnorm = enorm(n, Fvec) + fnorm = norm2(Fvec) ! initialize iteration counter and monitors. @@ -1182,7 +1110,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & do j = 1, n Wa3(j) = Diag(j)*x(j) end do - xnorm = enorm(n, Wa3) + xnorm = norm2(Wa3) delta = Factor*xnorm if (delta == zero) delta = Factor end if @@ -1255,7 +1183,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & Wa2(j) = x(j) + Wa1(j) Wa3(j) = Diag(j)*Wa1(j) end do - pnorm = enorm(n, Wa3) + pnorm = norm2(Wa3) ! on the first iteration, adjust the initial step bound. @@ -1267,7 +1195,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & call fcn(n, Wa2, Wa4, Fjac, Ldfjac, iflag) Nfev = Nfev + 1 if (iflag >= 0) then - fnorm1 = enorm(n, Wa4) + fnorm1 = norm2(Wa4) ! compute the scaled actual reduction. @@ -1285,7 +1213,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & end do Wa3(i) = Qtf(i) + sum end do - temp = enorm(n, Wa3) + temp = norm2(Wa3) prered = zero if (temp < fnorm) prered = one - (temp/fnorm)**2 @@ -1319,7 +1247,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & Wa2(j) = Diag(j)*x(j) Fvec(j) = Wa4(j) end do - xnorm = enorm(n, Wa2) + xnorm = norm2(Wa2) fnorm = fnorm1 iter = iter + 1 end if @@ -1611,7 +1539,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag) Nfev = 1 if (iflag >= 0) then - fnorm = enorm(m, Fvec) + fnorm = norm2(Fvec) ! initialize levenberg-marquardt parameter and iteration counter. @@ -1657,7 +1585,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & do j = 1, n Wa3(j) = Diag(j)*x(j) end do - xnorm = enorm(n, Wa3) + xnorm = norm2(Wa3) delta = Factor*xnorm if (delta == zero) delta = Factor end if @@ -1725,7 +1653,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & Wa2(j) = x(j) + Wa1(j) Wa3(j) = Diag(j)*Wa1(j) end do - pnorm = enorm(n, Wa3) + pnorm = norm2(Wa3) ! on the first iteration, adjust the initial step bound. @@ -1737,7 +1665,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & call fcn(m, n, Wa2, Wa4, Fjac, Ldfjac, iflag) Nfev = Nfev + 1 if (iflag >= 0) then - fnorm1 = enorm(m, Wa4) + fnorm1 = norm2(Wa4) ! compute the scaled actual reduction. @@ -1755,7 +1683,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & Wa3(i) = Wa3(i) + Fjac(i, j)*temp end do end do - temp1 = enorm(n, Wa3)/fnorm + temp1 = norm2(Wa3)/fnorm temp2 = (sqrt(par)*pnorm)/fnorm prered = temp1**2 + temp2**2/p5 dirder = -(temp1**2 + temp2**2) @@ -1790,7 +1718,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & do i = 1, m Fvec(i) = Wa4(i) end do - xnorm = enorm(n, Wa2) + xnorm = norm2(Wa2) fnorm = fnorm1 iter = iter + 1 end if @@ -2080,7 +2008,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & call fcn(m, n, x, Fvec, iflag) Nfev = 1 if (iflag >= 0) then - fnorm = enorm(m, Fvec) + fnorm = norm2(Fvec) ! initialize levenberg-marquardt parameter and iteration counter. @@ -2126,7 +2054,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & do j = 1, n Wa3(j) = Diag(j)*x(j) end do - xnorm = enorm(n, Wa3) + xnorm = norm2(Wa3) delta = Factor*xnorm if (delta == zero) delta = Factor end if @@ -2195,7 +2123,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & Wa2(j) = x(j) + Wa1(j) Wa3(j) = Diag(j)*Wa1(j) end do - pnorm = enorm(n, Wa3) + pnorm = norm2(Wa3) ! on the first iteration, adjust the initial step bound. @@ -2207,7 +2135,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & call fcn(m, n, Wa2, Wa4, iflag) Nfev = Nfev + 1 if (iflag >= 0) then - fnorm1 = enorm(m, Wa4) + fnorm1 = norm2(Wa4) ! compute the scaled actual reduction. @@ -2225,7 +2153,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & Wa3(i) = Wa3(i) + Fjac(i, j)*temp end do end do - temp1 = enorm(n, Wa3)/fnorm + temp1 = norm2(Wa3)/fnorm temp2 = (sqrt(par)*pnorm)/fnorm prered = temp1**2 + temp2**2/p5 dirder = -(temp1**2 + temp2**2) @@ -2263,7 +2191,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & do i = 1, m Fvec(i) = Wa4(i) end do - xnorm = enorm(n, Wa2) + xnorm = norm2(Wa2) fnorm = fnorm1 iter = iter + 1 end if @@ -2500,7 +2428,7 @@ subroutine lmpar(n, r, Ldr, Ipvt, Diag, Qtb, Delta, Par, x, Sdiag, Wa1, Wa2) do j = 1, n Wa2(j) = Diag(j)*x(j) end do - dxnorm = enorm(n, Wa2) + dxnorm = norm2(Wa2) fp = dxnorm - Delta if (fp <= p1*Delta) then ! termination. @@ -2527,7 +2455,7 @@ subroutine lmpar(n, r, Ldr, Ipvt, Diag, Qtb, Delta, Par, x, Sdiag, Wa1, Wa2) end if Wa1(j) = (Wa1(j) - sum)/r(j, j) end do - temp = enorm(n, Wa1) + temp = norm2(Wa1) parl = ((fp/Delta)/temp)/temp end if @@ -2541,7 +2469,7 @@ subroutine lmpar(n, r, Ldr, Ipvt, Diag, Qtb, Delta, Par, x, Sdiag, Wa1, Wa2) l = Ipvt(j) Wa1(j) = sum/Diag(l) end do - gnorm = enorm(n, Wa1) + gnorm = norm2(Wa1) paru = gnorm/Delta if (paru == zero) paru = dwarf/min(Delta, p1) @@ -2567,7 +2495,7 @@ subroutine lmpar(n, r, Ldr, Ipvt, Diag, Qtb, Delta, Par, x, Sdiag, Wa1, Wa2) do j = 1, n Wa2(j) = Diag(j)*x(j) end do - dxnorm = enorm(n, Wa2) + dxnorm = norm2(Wa2) temp = fp fp = dxnorm - Delta @@ -2596,7 +2524,7 @@ subroutine lmpar(n, r, Ldr, Ipvt, Diag, Qtb, Delta, Par, x, Sdiag, Wa1, Wa2) end do end if end do - temp = enorm(n, Wa1) + temp = norm2(Wa1) parc = ((fp/Delta)/temp)/temp ! depending on the sign of the function, update parl or paru. @@ -2765,7 +2693,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & call fcn(m, n, x, Fvec, Wa3, iflag) Nfev = 1 if (iflag < 0) goto 200 - fnorm = enorm(m, Fvec) + fnorm = norm2(Fvec) ! initialize levenberg-marquardt parameter and iteration counter. @@ -2810,7 +2738,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & do j = 1, n if (Fjac(j, j) == zero) sing = .true. Ipvt(j) = j - Wa2(j) = enorm(j, Fjac(1, j)) + Wa2(j) = norm2(Fjac(1:j, j)) end do if (sing) then call qrfac(n, n, Fjac, Ldfjac, .true., Ipvt, n, Wa1, Wa2, Wa3) @@ -2846,7 +2774,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & do j = 1, n Wa3(j) = Diag(j)*x(j) end do - xnorm = enorm(n, Wa3) + xnorm = norm2(Wa3) delta = Factor*xnorm if (delta == zero) delta = Factor end if @@ -2893,7 +2821,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & Wa2(j) = x(j) + Wa1(j) Wa3(j) = Diag(j)*Wa1(j) end do - pnorm = enorm(n, Wa3) + pnorm = norm2(Wa3) ! on the first iteration, adjust the initial step bound. @@ -2905,7 +2833,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & call fcn(m, n, Wa2, Wa4, Wa3, iflag) Nfev = Nfev + 1 if (iflag >= 0) then - fnorm1 = enorm(m, Wa4) + fnorm1 = norm2(Wa4) ! compute the scaled actual reduction. @@ -2923,7 +2851,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & Wa3(i) = Wa3(i) + Fjac(i, j)*temp end do end do - temp1 = enorm(n, Wa3)/fnorm + temp1 = norm2(Wa3)/fnorm temp2 = (sqrt(par)*pnorm)/fnorm prered = temp1**2 + temp2**2/p5 dirder = -(temp1**2 + temp2**2) @@ -2960,7 +2888,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & do i = 1, m Fvec(i) = Wa4(i) end do - xnorm = enorm(n, Wa2) + xnorm = norm2(Wa2) fnorm = fnorm1 iter = iter + 1 end if @@ -3232,7 +3160,7 @@ subroutine qrfac(m, n, a, Lda, Pivot, Ipvt, Lipvt, Rdiag, Acnorm, Wa) ! compute the initial column norms and initialize several arrays. do j = 1, n - Acnorm(j) = enorm(m, a(1, j)) + Acnorm(j) = norm2(a(1:m, j)) Rdiag(j) = Acnorm(j) Wa(j) = Rdiag(j) if (Pivot) Ipvt(j) = j @@ -3267,7 +3195,7 @@ subroutine qrfac(m, n, a, Lda, Pivot, Ipvt, Lipvt, Rdiag, Acnorm, Wa) ! compute the householder transformation to reduce the ! j-th column of a to a multiple of the j-th unit vector. - ajnorm = enorm(m - j + 1, a(j, j)) + ajnorm = norm2(a(j:m, j)) if (ajnorm /= zero) then if (a(j, j) < zero) ajnorm = -ajnorm do i = j, m @@ -3293,7 +3221,7 @@ subroutine qrfac(m, n, a, Lda, Pivot, Ipvt, Lipvt, Rdiag, Acnorm, Wa) temp = a(j, k)/Rdiag(k) Rdiag(k) = Rdiag(k)*sqrt(max(zero, one - temp**2)) if (p05*(Rdiag(k)/Wa(k))**2 <= epsmch) then - Rdiag(k) = enorm(m - j, a(jp1, k)) + Rdiag(k) = norm2(a(jp1:m, k)) Wa(k) = Rdiag(k) end if end if From 0e48dd5d22d41670c395e1873de78f4c694943ab Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Sat, 26 Feb 2022 12:57:38 +0100 Subject: [PATCH 2/6] Import precision in enorm --- src/enorm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/enorm.f90 b/src/enorm.f90 index 8bccd78..cc0c8b5 100644 --- a/src/enorm.f90 +++ b/src/enorm.f90 @@ -16,7 +16,7 @@ ! given here are suitable for every known computer. pure real(wp) function enorm(n, x) - + use iso_fortran_env, only: wp => real64 implicit none integer, intent(in) :: n !! a positive integer input variable. From 01bad721fa6c2e7beb8d3488b7069ccd660d5cb8 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Sat, 26 Feb 2022 13:02:27 +0100 Subject: [PATCH 3/6] Add missing parameters in enorm --- src/enorm.f90 | 90 ++++++++++++++++++++++++++------------------------- 1 file changed, 46 insertions(+), 44 deletions(-) diff --git a/src/enorm.f90 b/src/enorm.f90 index cc0c8b5..f17f02e 100644 --- a/src/enorm.f90 +++ b/src/enorm.f90 @@ -15,57 +15,59 @@ ! underflow and rgiant**2 not overflow. the constants ! given here are suitable for every known computer. - pure real(wp) function enorm(n, x) - use iso_fortran_env, only: wp => real64 - implicit none +pure real(wp) function enorm(n, x) + 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, 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 + 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 :: 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 + 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 - s1 = one + s1*(x1max/xabs)**2 - x1max = xabs + s3 = one + s3*(x3max/xabs)**2 + x3max = 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) + ! sum for large components. + elseif (xabs <= x1max) then + s1 = s1 + (xabs/x1max)**2 else - if (s2 >= x3max) enorm = sqrt(s2*(one + (x3max/s2)*(x3max*s3))) - if (s2 < x3max) enorm = sqrt(x3max*((s2/x3max) + (x3max*s3))) + 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 +end function enorm !***************************************************************************************** \ No newline at end of file From b4b72aebaadd59697ebf6860f8d1bd5efd7916a7 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Sat, 26 Feb 2022 13:06:55 +0100 Subject: [PATCH 4/6] Replace enorm in examples --- examples/example_hybrd.f90 | 4 ++-- examples/example_hybrd1.f90 | 4 ++-- examples/example_lmder1.f90 | 4 ++-- examples/example_lmdif1.f90 | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/example_hybrd.f90 b/examples/example_hybrd.f90 index b6bdd81..4fa461a 100644 --- a/examples/example_hybrd.f90 +++ b/examples/example_hybrd.f90 @@ -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 @@ -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// & diff --git a/examples/example_hybrd1.f90 b/examples/example_hybrd1.f90 index 21ae354..7d93ce5 100644 --- a/examples/example_hybrd1.f90 +++ b/examples/example_hybrd1.f90 @@ -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 @@ -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// & diff --git a/examples/example_lmder1.f90 b/examples/example_lmder1.f90 index 8b235bb..6635cb8 100644 --- a/examples/example_lmder1.f90 +++ b/examples/example_lmder1.f90 @@ -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 @@ -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' // & diff --git a/examples/example_lmdif1.f90 b/examples/example_lmdif1.f90 index 8b508b0..0bd1d5a 100644 --- a/examples/example_lmdif1.f90 +++ b/examples/example_lmdif1.f90 @@ -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 @@ -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' // & From 008400c1bd56bd7b43b58c50b46f2223496d46c8 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Sat, 26 Feb 2022 13:13:12 +0100 Subject: [PATCH 5/6] Replace enorm with norm2 in tests --- test/test_hybrd.f90 | 8 ++++---- test/test_hybrj.f90 | 8 ++++---- test/test_lmder.f90 | 4 ++-- test/test_lmdif.f90 | 8 ++++---- test/test_lmstr.f90 | 8 ++++---- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/test/test_hybrd.f90 b/test/test_hybrd.f90 index 5e017d4..9841845 100644 --- a/test/test_hybrd.f90 +++ b/test/test_hybrd.f90 @@ -19,9 +19,9 @@ program test ! ! USER-SUPPLIED ...... FCN ! -! MINPACK-SUPPLIED ... DPMPAR,ENORM,HYBRD1,INITPT,VECFCN +! MINPACK-SUPPLIED ... DPMPAR,HYBRD1,INITPT,VECFCN ! -! FORTRAN-SUPPLIED ... sqrt +! FORTRAN-SUPPLIED ... sqrt,norm2 ! ! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. ! BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE @@ -63,12 +63,12 @@ program test ic = ic + 1 call initpt(n,x,NPRob,factor) call vecfcn(n,x,fvec,NPRob) - fnorm1 = enorm(n,fvec) + fnorm1 = norm2(fvec(1:n)) write (nwrite,99005) NPRob , n 99005 format (////5x,' PROBLEM',i5,5x,' DIMENSION',i5,5x//) NFEv = 0 call hybrd1(fcn,n,x,fvec,tol,info,wa,lwa) - fnorm2 = enorm(n,fvec) + fnorm2 = norm2(fvec(1:n)) np(ic) = NPRob na(ic) = n nf(ic) = NFEv diff --git a/test/test_hybrj.f90 b/test/test_hybrj.f90 index a1e3213..a868b27 100644 --- a/test/test_hybrj.f90 +++ b/test/test_hybrj.f90 @@ -19,9 +19,9 @@ program test ! ! USER-SUPPLIED ...... FCN ! -! MINPACK-SUPPLIED ... DPMPAR,ENORM,HYBRJ1,INITPT,VECFCN +! MINPACK-SUPPLIED ... DPMPAR,HYBRJ1,INITPT,VECFCN ! -! FORTRAN-SUPPLIED ... sqrt +! FORTRAN-SUPPLIED ... sqrt,norm2 ! ! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. ! BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE @@ -66,13 +66,13 @@ program test ic = ic + 1 call initpt(n,x,NPRob,factor) call vecfcn(n,x,fvec,NPRob) - fnorm1 = enorm(n,fvec) + fnorm1 = norm2(fvec(1:n)) write (nwrite,99005) NPRob , n 99005 format (////5x,' PROBLEM',i5,5x,' DIMENSION',i5,5x//) NFEv = 0 NJEv = 0 call hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa) - fnorm2 = enorm(n,fvec) + fnorm2 = norm2(fvec(1:n)) np(ic) = NPRob na(ic) = n nf(ic) = NFEv diff --git a/test/test_lmder.f90 b/test/test_lmder.f90 index 8a11340..714df9b 100644 --- a/test/test_lmder.f90 +++ b/test/test_lmder.f90 @@ -55,14 +55,14 @@ program test ic = ic + 1 call initpt(n, x, NPRob, factor) call ssqfcn(m, n, x, fvec, NPRob) - fnorm1 = enorm(m, fvec) + fnorm1 = norm2(fvec(1:m)) write (nwrite, 99005) NPRob, n, m 99005 format(////5x, ' PROBLEM', i5, 5x, ' DIMENSIONS', 2i5, 5x//) NFEv = 0 NJEv = 0 call lmder1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info, iwa, wa, lwa) call ssqfcn(m, n, x, fvec, NPRob) - fnorm2 = enorm(m, fvec) + fnorm2 = norm2(fvec(1:m)) np(ic) = NPRob na(ic) = n ma(ic) = m diff --git a/test/test_lmdif.f90 b/test/test_lmdif.f90 index 486ad4c..9e95a70 100644 --- a/test/test_lmdif.f90 +++ b/test/test_lmdif.f90 @@ -19,9 +19,9 @@ program test ! ! USER-SUPPLIED ...... FCN ! -! MINPACK-SUPPLIED ... DPMPAR,ENORM,INITPT,LMDIF1,SSQFCN +! MINPACK-SUPPLIED ... DPMPAR,INITPT,LMDIF1,SSQFCN ! -! FORTRAN-SUPPLIED ... sqrt +! FORTRAN-SUPPLIED ... sqrt,norm2 ! ! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. ! BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE @@ -66,14 +66,14 @@ program test ic = ic + 1 call initpt(n,x,NPRob,factor) call ssqfcn(m,n,x,fvec,NPRob) - fnorm1 = enorm(m,fvec) + fnorm1 = norm2(fvec(1:m)) write (nwrite,99005) NPRob , n , m 99005 format (////5x,' PROBLEM',i5,5x,' DIMENSIONS',2i5,5x//) NFEv = 0 NJEv = 0 call lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa) call ssqfcn(m,n,x,fvec,NPRob) - fnorm2 = enorm(m,fvec) + fnorm2 = norm2(fvec(1:m)) np(ic) = NPRob na(ic) = n ma(ic) = m diff --git a/test/test_lmstr.f90 b/test/test_lmstr.f90 index 85bbefc..4c0aaa7 100644 --- a/test/test_lmstr.f90 +++ b/test/test_lmstr.f90 @@ -17,9 +17,9 @@ program test ! ! USER-SUPPLIED ...... FCN ! -! MINPACK-SUPPLIED ... DPMPAR,ENORM,INITPT,LMSTR1,SSQFCN +! MINPACK-SUPPLIED ... DPMPAR,INITPT,LMSTR1,SSQFCN ! -! FORTRAN-SUPPLIED ... sqrt +! FORTRAN-SUPPLIED ... sqrt,norm2 ! ! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. ! BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE @@ -72,14 +72,14 @@ program test ic = ic + 1 call initpt(n,x,NPRob,factor) call ssqfcn(m,n,x,fvec,NPRob) - fnorm1 = enorm(m,fvec) + fnorm1 = norm2(fvec(1:m)) write (nwrite,99005) NPRob , n , m 99005 format (////5x,' PROBLEM',i5,5x,' DIMENSIONS',2i5,5x//) NFEv = 0 NJEv = 0 call lmstr1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,iwa,wa,lwa) call ssqfcn(m,n,x,fvec,NPRob) - fnorm2 = enorm(m,fvec) + fnorm2 = norm2(fvec(1:m)) np(ic) = NPRob na(ic) = n ma(ic) = m From 9abf8be8c0497b52455aedd2beacd7137fa86557 Mon Sep 17 00:00:00 2001 From: Ivan Pribec Date: Sat, 26 Feb 2022 13:50:14 +0100 Subject: [PATCH 6/6] Remove enorm function --- src/enorm.f90 | 73 --------------------------------------------------- 1 file changed, 73 deletions(-) delete mode 100644 src/enorm.f90 diff --git a/src/enorm.f90 b/src/enorm.f90 deleted file mode 100644 index f17f02e..0000000 --- a/src/enorm.f90 +++ /dev/null @@ -1,73 +0,0 @@ -!***************************************************************************************** -!> -! 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) - 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 -!***************************************************************************************** \ No newline at end of file