From fd19724ca9d149e42e3cd5d8692dc0f1773d35b8 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 14 May 2021 09:38:32 -0400 Subject: [PATCH] Re-wrote lambda function class so that it can implicitly initialize using a structure constructor (I think this is super cool, even though the minimization method isn't working here yet) --- src/modules/lambda_function.f90 | 17 ++-- src/modules/module_interfaces.f90 | 13 +-- src/modules/module_symba.f90 | 80 ++++++++-------- src/symba/symba_frag_pos.f90 | 127 ++++++++++++------------- src/util/util_minimize_bfgs.f90 | 151 ++++++++++++++++-------------- 5 files changed, 196 insertions(+), 192 deletions(-) diff --git a/src/modules/lambda_function.f90 b/src/modules/lambda_function.f90 index 58cc79f60..bdd8da74e 100644 --- a/src/modules/lambda_function.f90 +++ b/src/modules/lambda_function.f90 @@ -4,33 +4,36 @@ module lambda_function implicit none type, public :: lambda_obj - private procedure(lambda0), pointer, nopass :: lambdaptr => null() + contains generic :: init => lambda_init_0 procedure :: eval => lambda_eval_0 - procedure :: lambda_init_0 + procedure, nopass :: lambda_init_0 final :: lambda_destroy end type + interface lambda_obj + module procedure lambda_init_0 + end interface abstract interface function lambda0(x) result(y) ! Template for a 0 argument function import DP real(DP), dimension(:), intent(in) :: x - real(DP) :: y + real(DP) :: y end function end interface contains - subroutine lambda_init_0(self, lambda) + type(lambda_obj) function lambda_init_0(lambda) implicit none ! Arguments - class(lambda_obj), intent(out) :: self procedure(lambda0) :: lambda - self%lambdaptr => lambda - end subroutine lambda_init_0 + lambda_init_0%lambdaptr => lambda + return + end function lambda_init_0 function lambda_eval_0(self, x) result(y) implicit none diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 29b4fb42e..43553ce56 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -1579,15 +1579,16 @@ function util_solve_linear_system(A,b,n,lerr) result(x) real(DP), dimension(n) :: x end function util_solve_linear_system - function util_minimize_bfgs(f, N, x1, eps) result(fnum) + function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) use swiftest_globals use lambda_function implicit none - integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(inout) :: x1 - real(DP), intent(in) :: eps - integer(I4B) :: fnum + integer(I4B), intent(in) :: N + class(lambda_obj), intent(in) :: f + real(DP), dimension(:), intent(in) :: x0 + real(DP), intent(in) :: eps + logical, intent(out) :: lerr + real(DP), dimension(:), allocatable :: x1 end function util_minimize_bfgs end interface diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90 index 536ef4ae1..30b29efc4 100644 --- a/src/modules/module_symba.f90 +++ b/src/modules/module_symba.f90 @@ -103,76 +103,76 @@ MODULE module_symba end type symba_merger - type, public, extends(lambda_obj) :: symba_vel_lambda_obj + type, public, extends(lambda_obj) :: symba_frag_lambda procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null() - real(DP), dimension(:), allocatable :: m_frag, v_r_mag - real(DP), dimension(:,:), allocatable :: v_r_unit - real(DP), dimension(NDIM) :: L_lin_tan - real(DP) :: T_rad + real(DP), dimension(:), allocatable :: m_frag + real(DP), dimension(:,:), allocatable :: x_frag + real(DP), dimension(:,:), allocatable :: v_frag + real(DP), dimension(NDIM) :: L_target + real(DP) :: ke_target contains generic :: init => ke_objective_func_init procedure :: eval => ke_objective_func_eval - procedure :: ke_objective_func_init + procedure, nopass :: ke_objective_func_init final :: ke_objective_func_destroy - end type symba_vel_lambda_obj + end type symba_frag_lambda + interface symba_frag_lambda + module procedure ke_objective_func_init + end interface abstract interface - function abstract_objective_func(v_r_mag_unknowns, v_r_mag_knowns, m_frag, v_r_unit, L_lin_tan, T_rad) result(fnorm) + function abstract_objective_func(vflat, v_frag, x_frag, m_frag, L_target, ke_target) result(fnorm) ! Template for the kinetic energy constraint function used for minimizing import DP - real(DP), dimension(:), intent(in) :: v_r_mag_unknowns !! Unknown radial velocity magnitudes - real(DP), dimension(:), intent(in) :: v_r_mag_knowns !! Known Radial velocity magnitude - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors - real(DP), dimension(:), intent(in) :: L_lin_tan !! Tangential component of linear momentum - real(DP), intent(in) :: T_rad !! Target radial kinetic energ - real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy + real(DP), dimension(:), intent(in) :: vflat !! Unrolled unknown velocity vectors + real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum + real(DP), intent(in) :: ke_target !! Target kinetic energ + real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy end function end interface contains - subroutine ke_objective_func_init(self, lambda, v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) + type(symba_frag_lambda) function ke_objective_func_init(lambda, v_frag, x_frag, m_frag, L_target, ke_target) implicit none ! Arguments - class(symba_vel_lambda_obj), intent(out) :: self procedure(abstract_objective_func) :: lambda !! The lambda function - real(DP), dimension(:), intent(in) :: v_r_mag !! Radial velocity magnitude - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors - real(DP), dimension(:), intent(in) :: L_lin_tan !! Tangential component of linear momentum - real(DP), intent(in) :: T_rad !! Target radial kinetic ener - - self%ke_objective_func_ptr => lambda - allocate(self%m_frag, source=m_frag) - allocate(self%v_r_mag, source=v_r_mag) - allocate(self%v_r_unit, source=v_r_unit) - self%L_lin_tan(:) = L_lin_tan(:) - self%T_rad = T_rad - end subroutine ke_objective_func_init + real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum + real(DP), intent(in) :: ke_target !! Target kinetic energ + ! Internals + associate(self => ke_objective_func_init) + self%ke_objective_func_ptr => lambda + allocate(self%m_frag, source=m_frag) + allocate(self%v_frag, source=v_frag) + allocate(self%x_frag, source=x_frag) + self%L_target(:) = L_target(:) + self%ke_target = ke_target + end associate + return + end function ke_objective_func_init subroutine ke_objective_func_destroy(self) implicit none - type(symba_vel_lambda_obj) :: self + type(symba_frag_lambda) :: self if (allocated(self%m_frag)) deallocate(self%m_frag) - if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) - if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) + if (allocated(self%v_frag)) deallocate(self%v_frag) + if (allocated(self%x_frag)) deallocate(self%x_frag) if (associated(self%ke_objective_func_ptr)) nullify(self%ke_objective_func_ptr) end subroutine ke_objective_func_destroy function ke_objective_func_eval(self, x) result(fnorm) implicit none ! Arguments - class(symba_vel_lambda_obj), intent(in) :: self - real(DP), dimension(:), intent(in) :: x + class(symba_frag_lambda), intent(in) :: self + real(DP), dimension(:), intent(in) :: x ! Result real(DP) :: fnorm - ! Internals - integer(I4B) :: nfrag, nunknown if (associated(self%ke_objective_func_ptr)) then - nfrag = size(self%v_r_mag) - nunknown = size(x) - fnorm = self%ke_objective_func_ptr(x, self%v_r_mag(nunknown + 1:nfrag), self%m_frag, self%v_r_unit, self%L_lin_tan, self%T_rad) + fnorm = self%ke_objective_func_ptr(x, self%v_frag, self%x_frag, self%m_frag, self%L_target, self%ke_target) else error stop "KE Objective function was not initialized." end if diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 6f44f21c4..39d9157ec 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -25,7 +25,7 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad real(DP), dimension(:,:), allocatable :: x_frag, v_frag ! Fragment positions and velocities in the collision center of mass frame real(DP), dimension(NDIM, 2) :: rot, L_orb integer(I4B) :: i, j, nfrag, fam_size, istart - real(DP), dimension(NDIM) :: xcom, vcom, Ltot_before, Ltot_after, L_residual, L_spin_frag + real(DP), dimension(NDIM) :: xcom, vcom, Ltot_before, Ltot_after, L_residual, L_spin_frag, L_target, L_frag real(DP) :: mtot, Lmag_before, Lmag_after real(DP) :: Etot_before, Etot_after, ke_before, pe_before real(DP) :: pe_after, ke_spin_before, ke_spin_after, ke_after, ke_family, ke_target, ke_frag @@ -115,7 +115,12 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad ! Set the "target" ke_after (the value of the orbital kinetic energy that the fragments ought to have) ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss - call symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_frag, v_frag, ke_target, lmerge) + L_target(:) = 0.0_DP + do i = 1, nfrag + call util_crossproduct(x_frag(:, i), v_frag(:, i), L_frag(:)) + L_target(:) = L_target(:) + L_frag(:) + end do + call symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_frag, lmerge) write(*, "(' ---------------------------------------------------------------------------')") write(*,fmtlabel) ' T_family |',ke_family / abs(Etot_before) @@ -320,7 +325,7 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_ return end subroutine symba_frag_pos_ang_mtm - subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_frag, v_frag, ke_target, lmerge) + subroutine symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_frag, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! @@ -332,96 +337,82 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr !! the target kinetic energy required to satisfy the constraints. implicit none ! Arguments - real(DP), dimension(:), intent(in) :: xcom, vcom !! Center of mass position and velocity in the system barycenter frame - real(DP), dimension(:,:), intent(in) :: L_orb, L_spin !! Pre-impact orbital and spin angular momentum - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(inout) :: x_frag, v_frag !! Fragment position and velocities in the center of mass frame - real(DP), intent(in) :: ke_target !! Target kinetic energy - logical, intent(out) :: lmerge + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:,:), intent(in) :: x_frag !! Fragment position vectors + real(DP), intent(in) :: ke_target !! Target kinetic energy + real(DP), dimension(:), intent(in) :: L_target !! Target kinetic energy + real(DP), dimension(:,:), intent(inout) :: v_frag !! Fragment position and velocities in the center of mass frame + logical, intent(out) :: lmerge ! Internals real(DP) :: mtot !! Total mass of fragments - real(DP) :: T_rad !! Sum of the radial kinetic energy of all fragments - real(DP), dimension(NDIM) :: L_lin_tan !! Sum of the tangential momentum vector of all fragments integer(I4B) :: i, nfrag, neval - real(DP), dimension(:,:), allocatable :: v_r_unit, v_t - real(DP), dimension(NDIM) :: v_t_unit, h_unit, L_orb_frag - real(DP), dimension(:), allocatable :: v_r_mag, v_r_mag_sol - integer(I4B), parameter :: MAXITER = 500 - real(DP), parameter :: TOL = 1e-12_DP - real(DP) :: vmult - type(symba_vel_lambda_obj) :: ke_objective_func + real(DP), parameter :: TOL = 1e-9_DP + real(DP), dimension(:), allocatable :: vflat + logical :: lerr + nfrag = size(m_frag) mtot = sum(m_frag) - allocate(v_r_unit, mold=v_frag) - allocate(v_t, mold=v_frag) - allocate(v_r_mag, mold=m_frag) - - ! Create the radial unit vectors pointing away from the collision center of mass, and subtract that off of the current - ! fragment velocities in order to create the tangential component - do i = 1, nfrag - v_r_unit(:, i) = x_frag(:,i) / norm2(x_frag(:, i)) - v_t(:,i) = v_frag(:, i) - dot_product(v_frag(:,i), v_r_unit(:, i)) * v_r_unit(:, i) - end do - - L_lin_tan = 0.0_DP - T_rad = ke_target - 0.5_DP * mtot * dot_product(vcom(:), vcom(:)) - do i = 1, nfrag - T_rad = T_rad - 0.5_DP * m_frag(i) * dot_product(v_t(:, i), v_t(:, i)) - L_lin_tan(:) = L_lin_tan(:) + m_frag(i) * v_t(:, i) - end do - if (T_rad > 0.0_DP) then - lmerge = .false. - - call random_number(v_r_mag(:)) - v_r_mag(:) = sqrt(2 * T_rad / nfrag / m_frag(:)) * (v_r_mag(:) + 0.1_DP) - - ! Initialize the lambda function with all the parameters that stay constant during the minimization - call ke_objective_func%init(symba_frag_pos_ke_objective_function, v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) - ! Minimize error using the BFGS optimizer - neval = util_minimize_bfgs(ke_objective_func, 4, v_r_mag(1:4), TOL) - - do i = 1, nfrag - v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i) - end do - else + ! Initialize the lambda function using a structure constructor that calls the init method + ! Minimize error using the BFGS optimizer + vflat = util_minimize_bfgs(symba_frag_lambda(lambda=symba_frag_pos_ke_objective_function, v_frag=v_frag, x_frag=x_frag, m_frag=m_frag, L_target=L_target, ke_target=ke_target), & + NDIM*nfrag, reshape(v_frag,[NDIM*nfrag]), TOL, lerr) + if (lerr) then ! No solution exists for this case, so we need to indicate that this should be a merge ! This may happen due to setting the tangential velocities too high when setting the angular momentum constraint lmerge = .true. - v_frag(:, :) = 0.0_DP + v_frag(:,:) = 0.0_DP + else + v_frag(:,:) = reshape(vflat,shape(v_frag)) end if return end subroutine symba_frag_pos_kinetic_energy - function symba_frag_pos_ke_objective_function(v_r_mag_unknowns, v_r_mag_knowns, m_frag, v_r_unit, L_lin_tan, T_rad) result(fnorm) + function symba_frag_pos_ke_objective_function(vflat, v_frag, x_frag, m_frag, L_target, ke_target) result(fnorm) ! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value implicit none ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag_unknowns !! Radial velocity magnitude - real(DP), dimension(:), intent(in) :: v_r_mag_knowns !! Radial velocity magnitude - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: L_lin_tan !! Tangential component of linear momentum - real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors - real(DP), intent(in) :: T_rad !! Target radial kinetic energy + real(DP), dimension(:), intent(in) :: vflat !! Unrolled unknown velocity vectors + real(DP), dimension(:,:), intent(in) :: v_frag, x_frag !! Velocity and position vectors + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum + real(DP), intent(in) :: ke_target !! Target kinetic energ ! Result real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy !! Minimizing this brings us closer to our objective ! Internals - real(DP), dimension(4) :: f_vec, f_vec_0 - - f_vec_0(1:3) = L_lin_tan(1:3) - f_vec_0(4) = -T_rad - f_vec(1) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] * v_r_unit(1,:)) - f_vec(2) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] * v_r_unit(2,:)) - f_vec(3) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] * v_r_unit(3,:)) - f_vec(4) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] **2) + real(DP), dimension(7) :: f_vec + integer(I4B) :: i, nfrag, nsol + real(DP), dimension(NDIM) :: L + real(DP), dimension(:,:), allocatable :: v_sol - f_vec(:) = f_vec(:) + f_vec_0(:) + nfrag = size(m_frag) + nsol = size(vflat) / NDIM + allocate(v_sol(NDIM,nsol)) + v_sol = reshape(vflat, shape(v_sol)) + + ! Linear momentum constraint + f_vec(1) = sum(m_frag(:) * [v_sol(1, 1:nsol), v_frag(1, nsol + 1:nfrag)]) + f_vec(2) = sum(m_frag(:) * [v_sol(2, 1:nsol), v_frag(2, nsol + 1:nfrag)]) + f_vec(3) = sum(m_frag(:) * [v_sol(3, 1:nsol), v_frag(3, nsol + 1:nfrag)]) + ! Angular momentum and kinetic energy constraints + f_vec(4:6) = -L_target(:) + f_vec(7) = -ke_target + do i = 1, nsol + call util_crossproduct(x_frag(:,i), v_sol(:, i), L(:)) + f_vec(4:6) = f_vec(4:6) + m_frag(i) * L(:) + f_vec(7) = f_vec(7) + 0.5_DP * m_frag(i) * dot_product(v_sol(:,i), v_sol(:, i)) + end do + do i = nsol + 1, nfrag + call util_crossproduct(x_frag(:,i), v_frag(:, i), L(:)) + f_vec(4:6) = f_vec(4:6) + m_frag(i) * L(:) + f_vec(7) = f_vec(7) + 0.5_DP * m_frag(i) * dot_product(v_frag(:,i), v_frag(:, i)) + end do - fnorm = norm2(f_vec(:) / f_vec_0(:)) + fnorm = norm2(f_vec(:)) return diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 2e1d01304..14e7b40f4 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -1,16 +1,16 @@ -function util_minimize_bfgs(f, N, x1, eps) result(fnum) +function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) !! author: David A. Minton !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. !! It recieves as input: !! f%eval(x) : lambda function object containing the objective function as the eval metho !! N : Number of variables of function f - !! x1 : Initial starting value of x + !! x0 : Initial starting value of x !! eps : Accuracy of 1 - dimensional minimization at each step !! The outputs include - !! x : Final minimum (all 0 if none found) + !! lerr : Returns .true. if it could not find the minimum !! Returns - !! Number of function calls performed or + !! x1 : Final minimum (all 0 if none found) !! 0 = No miniumum found !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - use swiftest @@ -18,79 +18,85 @@ function util_minimize_bfgs(f, N, x1, eps) result(fnum) use module_interfaces, EXCEPT_THIS_ONE => util_minimize_bfgs implicit none ! Arguments - integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(inout) :: x1 - real(DP), intent(in) :: eps + integer(I4B), intent(in) :: N + class(lambda_obj), intent(in) :: f + real(DP), dimension(:), intent(in) :: x0 + real(DP), intent(in) :: eps + logical, intent(out) :: lerr ! Result - integer(I4B) :: fnum + real(DP), dimension(:), allocatable :: x1 ! Internals - integer(I4B) :: i, j, k, l, conv, num + integer(I4B) :: i, j, k, l, conv, num, fnum integer(I4B), parameter :: MAXLOOP = 2000 !! Maximum number of loops before method is determined to have failed real(DP), dimension(N) :: S !! Direction vectors - real(DP), dimension(N) :: x0 real(DP), dimension(N) :: Snorm !! normalized direction real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix - real(DP), dimension(N) :: grad !! gradient of f + 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(DP) :: xmag fnum = 0 + lerr = .false. + ! 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]) - grad(:) = 0.0_DP + ! Get initial gradient and initialize arrays for updated values of gradient and x + fnum = fnum + gradf(f, N, x0(:), grad0, eps) + allocate(x1, source=x0) + grad1(:) = grad0(:) do i = 1, MAXLOOP - xmag = norm2(x1(:)) - grad0(:) = grad(:) - fnum = fnum + gradf(f, N, x1(1:N), grad, eps) - if (i > 1) then - ! set up factors for H matrix update - y(:) = grad(:) - grad0(:) - P(:) = x1(1:N) - x0(:) - Py = sum(P(:) * y(:)) - yHy = 0._DP - do k = 1, N - yHy = yHy + y(k) * sum(H(:,k) * y(:)) - end do - ! prevent divide by zero (convergence) - if (abs(Py) < tiny(Py)) return - ! set up update - PyH(:,:) = 0._DP - HyP(:,:) = 0._DP - do k = 1, N - do j = 1, N - PP(j, k) = P(j) * P(k) - PyH(j, k) = P(j) * sum(y(:) * H(:,k)) - HyP(j, k) = P(k) * sum(y(:) * H(j,:)) - end do - end do - ! update H matrix - H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:, :) - PyH(:, :) - HyP(:, :)) / Py - end if !check for convergence conv = 0 S(:) = 0._DP do k = 1, N - if (abs(grad(k)) > eps) conv = conv + 1 - S(k) = -sum(H(:,k) * grad(:)) + if (abs(grad1(k)) > eps) conv = conv + 1 + S(k) = -sum(H(:,k) * grad1(:)) end do if (conv == 0) return ! normalize gradient Snorm(:) = S(:) / norm2(S) - num = fnum + minimize1D(f, x1(1:N), Snorm, N, eps, astar) + num = fnum + minimize1D(f, x1, Snorm, N, eps, astar) if (num == fnum) then - write(*,*) "Exiting BFGS" - fnum = 0 + write(*,*) "Exiting BFGS with error in minimize1D step" + lerr = .true. return end if fnum = num ! Get new x values - x0(:) = x1(1:N) - x1(1:N) = x1(1:N) + astar * Snorm(:) + P(:) = astar * Snorm(:) + x1(:) = x1(:) + P(:) + ! Calculate new gradient + grad0(:) = grad1(:) + fnum = fnum + gradf(f, N, x1, grad1, eps) + y(:) = grad1(:) - grad0(:) + Py = sum(P(:) * y(:)) + ! set up factors for H matrix update + yHy = 0._DP + do k = 1, N + do j = 1, N + yHy = yHy + y(j) * H(j,k) * y(k) + end do + end do + ! prevent divide by zero (convergence) + if (abs(Py) < tiny(Py)) return + ! set up update + PyH(:,:) = 0._DP + HyP(:,:) = 0._DP + do k = 1, N + do j = 1, N + PP(j, k) = P(j) * P(k) + do l = 1, N + PyH(j, k) = PyH(j, k) + P(j) * y(l) * H(l,k) + HyP(j, k) = HyP(j, k) + P(k) * y(l) * H(j,l) + end do + end do + end do + ! update H matrix + H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py end do + return contains @@ -111,16 +117,17 @@ function gradf(f, N, x1, grad, dx) result(fnum) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments - integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f - real(DP), dimension(:), intent(in) :: x1 + 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(DP), intent(in) :: dx ! Result integer(I4B) :: fnum ! Internals integer(I4B) :: i, j, k real(DP), dimension(N) :: xp, xm + real(DP) :: fp, fm fnum = 0 do i = 1, N @@ -133,7 +140,9 @@ function gradf(f, N, x1, grad, dx) result(fnum) xm(j) = x1(j) end if end do - grad(i) = (f%eval(xp) - f%eval(xm)) / (2 * dx) + fp = f%eval(xp) + fm = f%eval(xm) + grad(i) = (fp - fm) / (2 * dx) fnum = fnum + 2 end do return @@ -160,10 +169,10 @@ function minimize1D(f, x0, S, N, eps, astar) result(fnum) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 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 + 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 ! Result integer(I4B) :: fnum @@ -215,9 +224,9 @@ function n2one(f, x0, S, N, a) result(fnew) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + class(lambda_obj), intent(in) :: f real(DP), dimension(:), intent(in) :: x0, S - real(DP), intent(in) :: a + real(DP), intent(in) :: a ! Return real(DP) :: fnew ! Internals @@ -251,7 +260,7 @@ function bracket(f, x0, S, N, lo, hi, gam, step) result(fnum) implicit none ! Arguments integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + 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 @@ -370,8 +379,8 @@ function golden(f, x0, S, N, lo, hi, eps) result(fnum) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments - integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + 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 @@ -419,7 +428,7 @@ end function golden function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function uses a quadratic polynomial fit to * locate the minimum of a function + !! This function uses a quadratic polynomial fit to locate the minimum of a function !! to some accuracy eps. It recieves as input: !! f%eval(x) : lambda function object containing the objective function as the eval metho !! lo : low bracket value @@ -434,8 +443,8 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments - integer(I4B), intent(in) :: N - class(lambda_obj), intent(in) :: f + 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 @@ -480,16 +489,16 @@ function quadfit(f, x0, S, N, lo, hi, eps) result(fnum) row_2 = [1.0_DP, a2, a2**2] row_3 = [1.0_DP, a3, a3**2] rhs = [f1, f2, f3] - lhs(1, :) = row_1 - lhs(2, :) = row_2 - lhs(3, :) = row_3 + lhs(:, 1) = row_1 + lhs(:, 2) = row_2 + lhs(:, 3) = row_3 ! 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.6f)') a1, f1 - write(*,'("a2 = ",f9.6," f2 = ",f9.6f)') a2, f2 - write(*,'("a3 = ",f9.6," f3 = ",f9.6f)') a3, f3 + 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