diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index cae1120a4..2878321d8 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -1568,8 +1568,9 @@ FUNCTION util_kahan_sum(xsum_current, xi, xerror) END FUNCTION END INTERFACE + interface - function util_solve_linear_system(A,b,n,lerr) result(x) + function util_solve_linear_system_d(A,b,n,lerr) result(x) use swiftest_globals implicit none real(DP), dimension(:,:), intent(in) :: A @@ -1577,18 +1578,50 @@ function util_solve_linear_system(A,b,n,lerr) result(x) integer(I4B), intent(in) :: n logical, intent(out) :: lerr real(DP), dimension(n) :: x - end function util_solve_linear_system + end function util_solve_linear_system_d - function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) + function util_solve_linear_system_q(A,b,n,lerr) result(x) + use swiftest_globals + implicit none + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + integer(I4B), intent(in) :: n + logical, intent(out) :: lerr + real(QP), dimension(n) :: x + end function util_solve_linear_system_q + + function solve_wbs(u, lerr) result(x) + use swiftest_globals + implicit none + real(QP), intent(in), dimension(:,:), allocatable :: u + logical, intent(out) :: lerr + real(QP), dimension(:), allocatable :: x + end function solve_wbs + + function ge_wpp(A, b) result(u) + use swiftest_globals + implicit none + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + real(QP), dimension(:,:), allocatable :: u + end function ge_wpp + end interface + + interface util_solve_linear_system + procedure util_solve_linear_system_d, util_solve_linear_system_q + end interface + + interface + function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d) use swiftest_globals use lambda_function implicit none integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0 - real(DP), intent(in) :: eps + real(DP), dimension(:), intent(in) :: x0_d + real(DP), intent(in) :: eps_d logical, intent(out) :: lerr - real(DP), dimension(:), allocatable :: x1 + real(DP), dimension(:), allocatable :: x1_d end function util_minimize_bfgs end interface diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 165995f54..247af61fb 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -543,7 +543,7 @@ subroutine set_fragment_radial_velocities(lmerge) ! Arguments logical, intent(out) :: lmerge ! Internals - real(DP), parameter :: TOL = epsilon(1._DP) + real(DP), parameter :: TOL = epsilon(1.0_DP) real(DP), dimension(:), allocatable :: vflat logical :: lerr integer(I4B) :: i diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index ddb715918..b8b0edc17 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -1,4 +1,4 @@ -function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) +function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d) !! author: David A. Minton !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. @@ -20,64 +20,66 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) ! Arguments integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0 - real(DP), intent(in) :: eps + real(DP), dimension(:), intent(in) :: x0_d + real(DP), intent(in) :: eps_d logical, intent(out) :: lerr ! Result - real(DP), dimension(:), allocatable :: x1 + real(DP), dimension(:), allocatable :: x1_d ! Internals - integer(I4B) :: i, j, k, l, conv, num, fnum + integer(I4B) :: i, j, k, l, conv, num integer(I4B), parameter :: MAXLOOP = 500 !! Maximum number of loops before method is determined to have failed - real(DP), parameter :: gradeps = 1e-6_DP !! Tolerance for gradient calculations - real(DP), dimension(N) :: S !! Direction vectors - real(DP), dimension(N) :: Snorm !! normalized direction - real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix - real(DP), dimension(N) :: grad1 !! gradient of f - real(DP), dimension(N) :: grad0 !! old value of gradient - real(DP) :: astar !! 1D minimized value - real(DP), dimension(N) :: y, P - real(DP), dimension(N,N) :: PP, PyH, HyP - real(DP) :: yHy, Py + real(QP), parameter :: gradeps = 1e-8_QP !! Tolerance for gradient calculations + real(QP), dimension(N) :: S !! Direction vectors + real(QP), dimension(N) :: Snorm !! normalized direction + real(QP), dimension(N,N) :: H !! Approximated inverse Hessian matrix + real(QP), dimension(N) :: grad1 !! gradient of f + real(QP), dimension(N) :: grad0 !! old value of gradient + real(QP) :: astar !! 1D minimized value + real(QP), dimension(N) :: y, P, x0, x1 + real(QP), dimension(N,N) :: PP, PyH, HyP + real(QP) :: yHy, Py, eps - fnum = 0 lerr = .false. + eps = real(eps_d, kind=QP) + x0 = real(x0_d, kind=QP) + allocate(x1_d, source=x0_d) + x1 = x0_d ! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent) - H(:,:) = reshape([((0._DP, i=1, j-1), 1._DP, (0._DP, i=j+1, N), j=1, N)], [N,N]) ! Get initial gradient and initialize arrays for updated values of gradient and x - fnum = fnum + gradf(f, N, x0(:), grad0, gradeps) - allocate(x1, source=x0) + H(:,:) = reshape([((0._QP, i=1, j-1), 1._QP, (0._QP, i=j+1, N), j=1, N)], [N,N]) + grad0 = gradf(f, N, x0(:), gradeps) grad1(:) = grad0(:) do i = 1, MAXLOOP !check for convergence conv = 0 - S(:) = 0._DP + S(:) = 0._QP do k = 1, N if (abs(grad1(k)) > eps) conv = conv + 1 S(k) = -sum(H(:,k) * grad1(:)) end do if (conv == 0) then write(*,*) "Converged on gradient after ",i," iterations" + x1_d = x1 return end if ! normalize gradient Snorm(:) = S(:) / norm2(S) - num = fnum + minimize1D(f, x1, Snorm, N, gradeps, astar) - if (num == fnum) then + astar = minimize1D(f, x1, Snorm, N, gradeps, lerr) + if (lerr) then write(*,*) "Exiting BFGS with error in minimize1D step" - lerr = .true. + x1_d = x1 return end if - fnum = num ! Get new x values P(:) = astar * Snorm(:) x1(:) = x1(:) + P(:) ! Calculate new gradient grad0(:) = grad1(:) - fnum = fnum + gradf(f, N, x1, grad1, gradeps) + grad1 = gradf(f, N, x1, gradeps) y(:) = grad1(:) - grad0(:) Py = sum(P(:) * y(:)) ! set up factors for H matrix update - yHy = 0._DP + yHy = 0._QP do k = 1, N do j = 1, N yHy = yHy + y(j) * H(j,k) * y(k) @@ -86,11 +88,12 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) ! prevent divide by zero (convergence) if (abs(Py) < tiny(Py)) then write(*,*) "Converged on tiny Py after ",i," iterations" + x1_d = x1 return end if ! set up update - PyH(:,:) = 0._DP - HyP(:,:) = 0._DP + PyH(:,:) = 0._QP + HyP(:,:) = 0._QP do k = 1, N do j = 1, N PP(j, k) = P(j) * P(k) @@ -101,18 +104,20 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) end do end do ! update H matrix - H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py - if (any(H(:,:) > sqrt(huge(1._DP)) / N)) then + H(:,:) = H(:,:) + ((1._QP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py + if (any(H(:,:) > sqrt(huge(1._QP)) / N)) then write(*,*) 'Did not converge after ',i,'iterations: H too big' + x1_d = x1 return end if end do write(*,*) "Did not converge! Ran out of loops" + x1_d = x1 return contains - function gradf(f, N, x1, grad, dx) result(fnum) + function gradf(f, N, x1, dx) result(grad) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! Purpose: Estimates the gradient of a function using a central difference !! approximation @@ -121,45 +126,44 @@ function gradf(f, N, x1, grad, dx) result(fnum) !! N : number of variables N !! x1 : x value array !! dx : step size to use when calculating derivatives - !! Output: - !! grad : N sized array containing estimated gradient of f at x1 !! Returns - !! Number of function calls performed + !! grad : N sized array containing estimated gradient of f at x1 !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x1 - real(DP), dimension(:), intent(out) :: grad - real(DP), intent(in) :: dx + real(QP), dimension(:), intent(in) :: x1 + real(QP), intent(in) :: dx ! Result - integer(I4B) :: fnum + real(QP), dimension(N) :: grad ! Internals - integer(I4B) :: i, j, k + integer(I4B) :: i, j real(DP), dimension(N) :: xp, xm - real(DP) :: fp, fm + real(DP) :: fp_d, fm_d, dx_d + real(QP) :: fp, fm - fnum = 0 + dx_d = dx do i = 1, N do j = 1, N if (j == i) then - xp(j) = x1(j) + dx - xm(j) = x1(j) - dx + xp(j) = x1(j) + dx_d + xm(j) = x1(j) - dx_d else xp(j) = x1(j) xm(j) = x1(j) end if end do - fp = f%eval(xp) - fm = f%eval(xm) + fp_d = f%eval(xp) + fm_d = f%eval(xm) + fp = real(fp_d, kind=QP) + fm = real(fm_d, kind=QP) grad(i) = (fp - fm) / (2 * dx) - fnum = fnum + 2 end do return end function gradf - function minimize1D(f, x0, S, N, eps, astar) result(fnum) + function minimize1D(f, x0, S, N, eps, lerr) result(astar) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This program find the minimum of a function of N variables in a single direction !! S using in sequence: @@ -173,61 +177,62 @@ function minimize1D(f, x0, S, N, eps, astar) result(fnum) !! N : Number of variables of function f !! eps : Accuracy of 1 - dimensional minimization at each step !! Output - !! astar : Final minimum along direction S + !! lerr : .true. if an error occurred. Otherwise returns .false. !! Returns - !! Number of function calls performed or - !! 0 = No miniumum found + !! astar : Final minimum along direction S !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0, S - real(DP), intent(in) :: eps - real(DP), intent(out) :: astar + real(QP), dimension(:), intent(in) :: x0, S + real(QP), intent(in) :: eps + logical, intent(out) :: lerr ! Result - integer(I4B) :: fnum + real(QP) :: astar ! Internals integer(I4B) :: num = 0 - real(DP), parameter :: step = 0.7_DP !! Bracketing method step size - real(DP), parameter :: gam = 1.2_DP !! Bracketing method expansion parameter - real(DP), parameter :: greduce = 0.2_DP !! Golden section method reduction factor - real(DP), parameter :: greduce2 = 0.0_DP ! Secondary golden section method reduction factor - real(DP) :: alo, ahi !! High and low values for 1 - D minimization routines - real(DP), parameter :: a0 = 0.0_DP !! Initial guess of alpha + real(QP), parameter :: step = 0.7_QP !! Bracketing method step size + real(QP), parameter :: gam = 1.2_QP !! Bracketing method expansion parameter + real(QP), parameter :: greduce = 0.2_QP !! Golden section method reduction factor + real(QP), parameter :: greduce2 = 0.1_QP ! Secondary golden section method reduction factor + real(QP) :: alo, ahi !! High and low values for 1 - D minimization routines + real(QP), parameter :: a0 = epsilon(1.0_QP) !! Initial guess of alpha - fnum = 0 alo = a0 - num = fnum + bracket(f, x0, S, N, alo, ahi, gam, step) - if (num == fnum) then + call bracket(f, x0, S, N, gam, step, alo, ahi, lerr) + if (lerr) then write(*,*) "Bracketing step failed!" - fnum = 0 return end if - fnum = num if (abs(alo - ahi) < eps) then astar = alo + lerr = .false. return end if - num = fnum + golden(f, x0, S, N, alo, ahi, greduce) - if (num == fnum) then + call golden(f, x0, S, N, greduce, alo, ahi, lerr) + if (lerr) then write(*,*) "Golden section step failed!" - fnum = 0 return end if - fnum = num if (abs(alo - ahi) < eps) then astar = alo + lerr = .false. + return + end if + call quadfit(f, x0, S, N, eps, alo, ahi, lerr) + if (lerr) then + write(*,*) "quadfit failed!" return end if - fnum = fnum + quadfit(f, x0, S, N, alo, ahi, eps) if (abs(alo - ahi) < eps) then astar = alo + lerr = .false. return end if ! Quadratic fit method won't converge, so finish off with another golden section - fnum = fnum + golden(f, x0, S, N, alo, ahi, greduce2) - astar = (alo + ahi) / 2.0_DP + call golden(f, x0, S, N, greduce2, alo, ahi, lerr) + if (.not. lerr) astar = (alo + ahi) / 2.0_QP return end function minimize1D @@ -236,56 +241,55 @@ function n2one(f, x0, S, N, a) result(fnew) ! Arguments integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0, S - real(DP), intent(in) :: a + real(QP), dimension(:), intent(in) :: x0, S + real(QP), intent(in) :: a ! Return - real(DP) :: fnew + real(QP) :: fnew ! Internals real(DP), dimension(N) :: xnew integer(I4B) :: i - xnew(:) = x0(:) + a * S(:) - fnew = f%eval(xnew(:)) + xnew(:) = real(x0(:) + a * S(:), kind=DP) + fnew = real(f%eval(xnew(:)), kind=QP) return end function n2one ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - function bracket(f, x0, S, N, lo, hi, gam, step) result(fnum) + subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function brackets the minimum. It recieves as input: + !! This subroutine brackets the minimum. It recieves as input: !! f%eval(x) : lambda function object containing the objective function as the eval metho !! x0 : Array of size N of initial x values !! S : Array of size N that determines the direction of minimization - !! lo : initial guess of lo bracket value !! gam : expansion parameter !! step : step size + !! lo : initial guess of lo bracket value !! The outputs include - !! lo : lo bracket + !! lo : lo bracket !! hi : hi bracket - !! Returns - !! Number of function calls performed or - !! 0 = No miniumum found + !! lerr : .true. if an error occurred. Otherwise returns .false. !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments - integer(I4B), intent(in) :: N + integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0, S - real(DP), intent(inout) :: lo, hi - real(DP), intent(in) :: gam, step - ! Result - integer(I4B) :: fnum + real(QP), dimension(:), intent(in) :: x0, S + real(QP), intent(in) :: gam, step + real(QP), intent(inout) :: lo + real(QP), intent(out) :: hi + logical, intent(out) :: lerr ! Internals - real(DP) :: a0, a1, a2, a3, da - real(DP) :: f0, f1, f2, fcon + real(QP) :: a0, a1, a2, a3, da + real(QP) :: f0, f1, f2, fcon integer(I4B) :: i integer(I4B), parameter :: MAXLOOP = 1000 ! maximum number of loops before method is determined to have failed - real(DP), parameter :: eps = tiny(lo) ! small number precision to test floating point equality - real(DP), parameter :: dela = 12.4423_DP ! arbitrary number to test if function is constant + real(QP), parameter :: eps = epsilon(lo) ! small number precision to test floating point equality + real(QP), parameter :: dela = 12.4423_QP ! arbitrary number to test if function is constant ! set up initial bracket points + lerr = .false. a0 = lo da = step a1 = a0 + da @@ -293,7 +297,6 @@ function bracket(f, x0, S, N, lo, hi, gam, step) result(fnum) f0 = n2one(f, x0, S, N, a0) f1 = n2one(f, x0, S, N, a1) f2 = n2one(f, x0, S, N, a2) - fnum = 3 ! set number off function calls ! loop over bracket method until either min is bracketed method fails do i = 1, MAXLOOP if ((f0 > f1) .and. (f2 > f1)) then ! minimum was found @@ -309,7 +312,6 @@ function bracket(f, x0, S, N, lo, hi, gam, step) result(fnum) f0 = f1 f1 = f2 f2 = n2one(f, x0, S, N, a2) - fnum = fnum + 1 else if ((f0 < f1) .and. (f1 < f2)) then ! function appears to increase da = da * gam a3 = a0 - da @@ -318,7 +320,6 @@ function bracket(f, x0, S, N, lo, hi, gam, step) result(fnum) a0 = a3 f2 = f1 f0 = n2one(f, x0, S, N, a0) - fnum = fnum + 1 else if ((f0 < f1) .and. (f1 > f2)) then ! more than one minimum present, so in this case we arbitrarily choose the RHS min da = da * gam a3 = a2 + da @@ -328,90 +329,82 @@ function bracket(f, x0, S, N, lo, hi, gam, step) result(fnum) f0 = f1 f1 = f2 f2 = n2one(f, x0, S, N, a2) - fnum = fnum + 1 else if ((f0 > f1) .and. (abs(f2 - f1) <= eps)) then ! RHS equal da = da * gam a3 = a2 + da a2 = a3 f2 = n2one(f, x0, S, N, a2) - fnum = fnum + 1 else if ((abs(f0 - f1) < eps) .and. (f2 > f1)) then ! LHS equal da = da * gam a3 = a0 - da a0 = a3 f0 = n2one(f, x0, S, N, a0) - fnum = fnum + 1 else ! all values equal stops if there is no minimum or takes RHS min if it exists ! test if function itself is constant - fcon = n2one(f, x0, S, N, a2 + dela) !multiply be an arbitrary number to see if constant - fnum = fnum + 1 + fcon = n2one(f, x0, S, N, a2 + dela) !add by an arbitrary number to see if constant if (abs(f2 - fcon) < eps) then - fnum = 0 + lerr = .true. return ! function is constant end if - a3 = a0 + 0.5_DP * (a1 - a0) + a3 = a0 + 0.5_QP * (a1 - a0) a0 = a1 a1 = a2 a2 = a3 f0 = f1 f1 = f2 - a3 = a0 + 0.5_DP * (a1 - a0) + a3 = a0 + 0.5_QP * (a1 - a0) a1 = a2 a2 = a3 f1 = f2 f2 = n2one(f, x0, S, N, a2) - fnum = fnum + 1 end if end do - fnum = 0 + lerr = .true. return ! no minimum found - end function bracket + end subroutine bracket ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - function golden(f, x0, S, N, lo, hi, eps) result(fnum) + subroutine golden(f, x0, S, N, eps, lo, hi, lerr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function uses the golden section method to reduce the starting interval lo, hi by some amount sigma. !! It recieves as input: !! f%eval(x) : lambda function object containing the objective function as the eval metho !! x0 : Array of size N of initial x values !! S : Array of size N that determines the direction of minimization - !! lo : initial guess of lo bracket value !! gam : expansion parameter !! eps : reduction interval in range (0 < sigma < 1) such that: !! hi(new) - lo(new) = eps * (hi(old) - lo(old)) + !! lo : initial guess of lo bracket value !! The outputs include !! lo : lo bracket !! hi : hi bracket - !! Returns - !! Number of function calls performed or - !! 0 = No miniumum found + !! lerr : .true. if an error occurred. Otherwise returns .false. !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments - integer(I4B), intent(in) :: N + integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0, S - real(DP), intent(inout) :: lo, hi - real(DP), intent(in) :: eps - ! Result - integer(I4B) :: fnum + real(QP), dimension(:), intent(in) :: x0, S + real(QP), intent(in) :: eps + real(QP), intent(inout) :: lo + real(QP), intent(out) :: hi + logical, intent(out) :: lerr ! Internals - real(DP), parameter :: tau = 0.5_DP * (sqrt(5.0_DP) - 1.0_DP) ! Golden section constant + real(QP), parameter :: tau = 0.5_QP * (sqrt(5.0_QP) - 1.0_QP) ! Golden section constant integer(I4B), parameter :: MAXLOOP = 40 ! maximum number of loops before method is determined to have failed (unlikely, but could occur if no minimum exists between lo and hi) - real(DP) :: i0 ! Initial interval value - real(DP) :: a1, a2 - real(DP) :: f1, f2 + real(QP) :: i0 ! Initial interval value + real(QP) :: a1, a2 + real(QP) :: f1, f2 integer(I4B) :: i, j i0 = hi - lo - fnum = 0 a1 = hi - tau * i0 a2 = lo + tau * i0 f1 = n2one(f, x0, S, N, a1) f2 = n2one(f, x0, S, N, a2) - fnum = fnum + 2 + lerr = .false. do i = 1, MAXLOOP if (abs((hi - lo) / i0) <= eps) return ! interval reduced to input amount if (f2 > f1) then @@ -420,24 +413,22 @@ function golden(f, x0, S, N, lo, hi, eps) result(fnum) f2 = f1 a1 = hi - tau * (hi - lo) f1 = n2one(f, x0, S, N, a1) - fnum = fnum + 1 else lo = a1 a1 = a2 f2 = f1 - a2 = hi - (1.0_DP - tau) * (hi - lo) + a2 = hi - (1.0_QP - tau) * (hi - lo) f2 = n2one(f, x0, S, N, a2) - fnum = fnum + 1 end if end do - fnum = 0 + lerr = .true. return ! search took too many iterations - no minimum found - end function golden + end subroutine golden ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) + subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function uses a quadratic polynomial fit to locate the minimum of a function !! to some accuracy eps. It recieves as input: @@ -454,35 +445,33 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments - integer(I4B), intent(in) :: N + integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x0, S - real(DP), intent(inout) :: lo, hi - real(DP), intent(in) :: eps - ! Result - integer(I4B) :: fnum + real(QP), dimension(:), intent(in) :: x0, S + real(QP), intent(in) :: eps + real(QP), intent(inout) :: lo + real(QP), intent(out) :: hi + logical, intent(out) :: lerr ! Internals integer(I4B), parameter :: MAXLOOP = 20 ! maximum number of loops before method is determined to have failed. - real(DP), parameter :: small = tiny(eps) ! small number to prevent divide by zero errors - real(DP) :: a1, a2, a3, astar ! three points for the polynomial fit and polynomial minimum - real(DP) :: f1, f2, f3, fstar ! three function values for the polynomial and polynomial minimum - real(DP), dimension(3) :: row_1, row_2, row_3, rhs, soln ! matrix for 3 equation solver (gaussian elimination) - real(DP), dimension(3,3) :: lhs - real(DP) :: d1, d2, d3, aold, denom + real(QP), parameter :: small = epsilon(eps) ! small number to prevent divide by zero errors + real(QP) :: a1, a2, a3, astar ! three points for the polynomial fit and polynomial minimum + real(QP) :: f1, f2, f3, fstar ! three function values for the polynomial and polynomial minimum + real(QP), dimension(3) :: row_1, row_2, row_3, rhs, soln ! matrix for 3 equation solver (gaussian elimination) + real(QP), dimension(3,3) :: lhs + real(QP) :: d1, d2, d3, aold, denom integer(I4B) :: i - logical :: lerr + lerr = .false. ! Get initial a1, a2, a3 values a1 = lo - a2 = lo + 0.5_DP * (hi - lo) + a2 = lo + 0.5_QP * (hi - lo) a3 = hi - fnum = 0 aold = a1 astar = a2 f1 = n2one(f, x0, S, N, a1) f2 = n2one(f, x0, S, N, a2) f3 = n2one(f, x0, S, N, a3) - fnum = fnum + 3 do i = 1, MAXLOOP ! check to see if convergence is reached and exit if (abs(astar) < small) then @@ -496,9 +485,9 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) return end if ! Set up system for gaussian elimination equation solver - row_1 = [1.0_DP, a1, a1**2] - row_2 = [1.0_DP, a2, a2**2] - row_3 = [1.0_DP, a3, a3**2] + row_1 = [1.0_QP, a1, a1**2] + row_2 = [1.0_QP, a2, a2**2] + row_3 = [1.0_QP, a3, a3**2] rhs = [f1, f2, f3] lhs(1, :) = row_1 lhs(2, :) = row_2 @@ -506,19 +495,12 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) ! Solve system of equations soln(:) = util_solve_linear_system(lhs, rhs, 3, lerr) if (lerr) then - !write(*,*) "Could not solve polynomial on loop ", i - !write(*,'("a1 = ",f9.6," f1 = ",f9.6)') a1, f1 - !write(*,'("a2 = ",f9.6," f2 = ",f9.6)') a2, f2 - !write(*,'("a3 = ",f9.6," f3 = ",f9.6)') a3, f3 - !write(*,'("aold = ",f7.4)') aold - fnum = 0 return end if aold = astar if (abs(soln(3)) < small) soln(3) = small astar = -soln(2) / (2 * soln(3)) fstar = n2one(f, x0, S, N, astar) - fnum = fnum + 1 ! keep the three closest a values to astar and discard the fourth d1 = abs(a1 - astar) d2 = abs(a2 - astar) @@ -545,7 +527,6 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) lo = a1 hi = a3 return - end function quadfit - + end subroutine quadfit end function util_minimize_bfgs \ No newline at end of file diff --git a/src/util/util_solve_linear_system.f90 b/src/util/util_solve_linear_system.f90 index 3933d1f68..dd9218ec8 100644 --- a/src/util/util_solve_linear_system.f90 +++ b/src/util/util_solve_linear_system.f90 @@ -1,4 +1,4 @@ -function util_solve_linear_system(A,b,n,lerr) result(x) +function util_solve_linear_system_d(A,b,n,lerr) result(x) !! Author: David A. Minton !! !! Solves the linear equation of the form A*x = b for x. @@ -7,81 +7,106 @@ function util_solve_linear_system(A,b,n,lerr) result(x) !! Uses Gaussian elimination, so will have issues if system is ill-conditioned. !! Uses quad precision intermidiate values, so works best on small arrays. use swiftest - use module_interfaces, EXCEPT_THIS_ONE => util_solve_linear_system + use module_interfaces, EXCEPT_THIS_ONE => util_solve_linear_system_d implicit none ! Arguments - integer(I4B), intent(in) :: n - real(DP), dimension(:,:), intent(in) :: A - real(DP), dimension(:), intent(in) :: b + integer(I4B), intent(in) :: n + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:), intent(in) :: b logical, intent(out) :: lerr ! Result - real(DP), dimension(n) :: x + real(DP), dimension(n) :: x ! Internals real(QP), dimension(:), allocatable :: qx qx = solve_wbs(ge_wpp(real(A, kind=QP), real(b, kind=QP)),lerr) where(abs(qx(:)) < tiny(1._DP)) qx(:) = 0._QP x = real(qx, kind=DP) + return +end function util_solve_linear_system_d - contains +function util_solve_linear_system_q(A,b,n,lerr) result(x) + !! Author: David A. Minton + !! + !! Solves the linear equation of the form A*x = b for x. + !! A is an (n,n) arrays + !! x and b are (n) arrays + !! Uses Gaussian elimination, so will have issues if system is ill-conditioned. + !! Uses quad precision intermidiate values, so works best on small arrays. + use swiftest + use module_interfaces, EXCEPT_THIS_ONE => util_solve_linear_system_q + implicit none + ! Arguments + integer(I4B), intent(in) :: n + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + ! Result + real(QP), dimension(n) :: x - function solve_wbs(u, lerr) result(x) ! solve with backward substitution - !! Based on code available on Rosetta Code: https://rosettacode.org/wiki/Gaussian_elimination#Fortran - implicit none - ! Arguments - real(QP), intent(in), dimension(:,:), allocatable :: u - logical, intent(out) :: lerr - ! Result - real(QP), dimension(:), allocatable :: x - ! Internals - integer(I4B) :: i,n - real(DP), parameter :: epsilon = 10 * tiny(1._DP) + x = solve_wbs(ge_wpp(A, b),lerr) + where(abs(x(:)) < tiny(1._QP)) x(:) = 0._QP - n = size(u, 1) - if (allocated(x) .and. (size(x) /= n)) deallocate(x) - if (.not.allocated(x)) allocate(x(n)) - lerr = any(abs(u(:,:)) < epsilon) - if (lerr) then - x(:) = 0._DP - return - end if + return +end function util_solve_linear_system_q + +function solve_wbs(u, lerr) result(x) ! solve with backward substitution + !! Based on code available on Rosetta Code: https://rosettacode.org/wiki/Gaussian_elimination#Fortran + use swiftest + implicit none + ! Arguments + real(QP), intent(in), dimension(:,:), allocatable :: u + logical, intent(out) :: lerr + ! Result + real(QP), dimension(:), allocatable :: x + ! Internals + integer(I4B) :: i,n + real(QP), parameter :: epsilon = tiny(1._QP) - do i = n,1,-1 - x(i) = (u(i, n + 1) - sum(u(i, i + 1:n) * x(i + 1:n))) / u(i, i) - end do - return - end function solve_wbs + n = size(u, 1) + if (allocated(x) .and. (size(x) /= n)) deallocate(x) + if (.not.allocated(x)) allocate(x(n)) + lerr = any(abs(u(:,:)) <= epsilon) + if (lerr) then + x(:) = 0._QP + return + end if - function ge_wpp(A, b) result(u) ! gaussian eliminate with partial pivoting - !! Solve Ax=b using Gaussian elimination then backwards substitution. - !! A being an n by n matrix. - !! x and b are n by 1 vectors. - !! Based on code available on Rosetta Code: https://rosettacode.org/wiki/Gaussian_elimination#Fortran - implicit none - ! Arguments - real(QP), dimension(:,:), intent(in) :: A - real(QP), dimension(:), intent(in) :: b - ! Result - real(QP), dimension(:,:), allocatable :: u - ! Internals - integer(I4B) :: i,j,n,p - real(QP) :: upi + do i = n, 1, -1 + x(i) = (u(i, n + 1) - sum(u(i, i + 1:n) * x(i + 1:n))) / u(i, i) + end do + return +end function solve_wbs - n = size(a, 1) - allocate(u(n, (n + 1))) - u = reshape([A, b], [n, n + 1]) - do j = 1, n - p = maxloc(abs(u(j:n, j)), 1) + j - 1 ! maxloc returns indices between (1, n - j + 1) - if (p /= j) u([p, j], j) = u([j, p], j) - u(j + 1:, j) = u(j + 1:, j) / u(j, j) - do i = j + 1, n + 1 - upi = u(p, i) - if (p /= j) u([p, j], i) = u([j, p], i) - u(j + 1:n, i) = u(j + 1:n, i) - upi * u(j + 1:n, j) - end do - end do - return - end function ge_wpp +function ge_wpp(A, b) result(u) ! gaussian eliminate with partial pivoting + !! Solve Ax=b using Gaussian elimination then backwards substitution. + !! A being an n by n matrix. + !! x and b are n by 1 vectors. + !! Based on code available on Rosetta Code: https://rosettacode.org/wiki/Gaussian_elimination#Fortran + use swiftest + implicit none + ! Arguments + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + ! Result + real(QP), dimension(:,:), allocatable :: u + ! Internals + integer(I4B) :: i,j,n,p + real(QP) :: upi - end function util_solve_linear_system \ No newline at end of file + n = size(a, 1) + allocate(u(n, (n + 1))) + u = reshape([A, b], [n, n + 1]) + do j = 1, n + p = maxloc(abs(u(j:n, j)), 1) + j - 1 ! maxloc returns indices between (1, n - j + 1) + if (p /= j) u([p, j], j) = u([j, p], j) + u(j + 1:, j) = u(j + 1:, j) / u(j, j) + do i = j + 1, n + 1 + upi = u(p, i) + if (p /= j) u([p, j], i) = u([j, p], i) + u(j + 1:n, i) = u(j + 1:n, i) - upi * u(j + 1:n, j) + end do + end do + return +end function ge_wpp