From e8a8d776eccb8ee07ea565e2c2997c43ee6a65fc Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 14 May 2021 17:36:00 -0400 Subject: [PATCH] Improved array handling and tolerance --- src/symba/symba_frag_pos.f90 | 2 +- src/util/util_minimize_bfgs.f90 | 5 +++-- src/util/util_solve_linear_system.f90 | 6 +++--- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 39d9157ec..9fdc8014c 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -347,7 +347,7 @@ subroutine symba_frag_pos_kinetic_energy(m_frag, ke_target, L_target, x_frag, v_ ! Internals real(DP) :: mtot !! Total mass of fragments integer(I4B) :: i, nfrag, neval - real(DP), parameter :: TOL = 1e-9_DP + real(DP), parameter :: TOL = 1e-12_DP real(DP), dimension(:), allocatable :: vflat logical :: lerr diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 8971cf6e3..a5425dd44 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -28,6 +28,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) ! Internals 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), parameter :: gradeps = 1e-4_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 @@ -43,7 +44,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) ! 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, eps) + fnum = fnum + gradf(f, N, x0(:), grad0, gradeps) allocate(x1, source=x0) grad1(:) = grad0(:) do i = 1, MAXLOOP @@ -69,7 +70,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) x1(:) = x1(:) + P(:) ! Calculate new gradient grad0(:) = grad1(:) - fnum = fnum + gradf(f, N, x1, grad1, eps) + fnum = fnum + gradf(f, N, x1, grad1, gradeps) y(:) = grad1(:) - grad0(:) Py = sum(P(:) * y(:)) ! set up factors for H matrix update diff --git a/src/util/util_solve_linear_system.f90 b/src/util/util_solve_linear_system.f90 index b52f08796..3933d1f68 100644 --- a/src/util/util_solve_linear_system.f90 +++ b/src/util/util_solve_linear_system.f90 @@ -38,9 +38,9 @@ function solve_wbs(u, lerr) result(x) ! solve with backward substitution integer(I4B) :: i,n real(DP), parameter :: epsilon = 10 * tiny(1._DP) - if (allocated(x)) deallocate(x) - allocate(x(n)) - + 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