From f89f2e273526373bde3aa3b039d5517e35e0a9cc Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 13 May 2021 14:55:01 -0400 Subject: [PATCH] Refactoring for clarity. Replaced old Gaussian elimination code with newer one, and incorporated more error checking --- src/modules/module_interfaces.f90 | 9 +- src/symba/symba_frag_pos.f90 | 14 +- .../{util_bfgs.f90 => util_minimize_bfgs.f90} | 132 ++++-------------- src/util/util_solve_linear_system.f90 | 15 +- 4 files changed, 49 insertions(+), 121 deletions(-) rename src/util/{util_bfgs.f90 => util_minimize_bfgs.f90} (82%) diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 018508aad..41e2ed20d 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -1569,16 +1569,17 @@ FUNCTION util_kahan_sum(xsum_current, xi, xerror) END INTERFACE interface - function util_solve_linear_system(n,A,b) result(x) + function util_solve_linear_system(A,b,n,lerr) result(x) use swiftest_globals implicit none - integer(I4B), intent(in) :: n real(DP), dimension(:,:), intent(in) :: A real(DP), dimension(:), intent(in) :: b + integer(I4B), intent(in) :: n + logical, intent(out) :: lerr real(DP), dimension(n) :: x end function util_solve_linear_system - function util_bfgs(f, N, x1, eps) result(fnum) + function util_minimize_bfgs(f, N, x1, eps) result(fnum) use swiftest_globals implicit none integer(I4B), intent(in) :: N @@ -1592,7 +1593,7 @@ end function f real(DP), dimension(:), intent(inout) :: x1 real(DP), intent(in) :: eps integer(I4B) :: fnum - end function util_bfgs + end function util_minimize_bfgs end interface INTERFACE diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 3300b4fff..b22a112c9 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -79,7 +79,8 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad nfrag = size(m_frag) ! Initialize positions and velocities of fragments that conserve angular momentum - call symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag) + call symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) + if (lmerge) return ! Energy calculation requires the fragments to be in the system barcyentric frame, so we need to temporarily shift them do i = 1, nfrag @@ -182,7 +183,7 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad contains - subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag) + subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the @@ -197,6 +198,7 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Fragment masses and radii real(DP), dimension(:,:), intent(in) :: Ip_frag !! Fragment prinicpal moments of inertia real(DP), dimension(:,:), intent(out) :: x_frag, v_frag, rot_frag !! Fragment position, velocities, and spin states + logical, intent(out) :: lmerge ! Internals real(DP) :: mtot, theta, v_frag_norm, r_frag_norm, v_col_norm, r_col_norm real(DP), dimension(NDIM) :: Ltot, delta_r, delta_v @@ -236,12 +238,12 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L call symba_frag_pos_com_adjust(xcom, m_frag, x_frag) v_frag(:,:) = 0._DP - call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag) + call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) return end subroutine symba_frag_pos_initialize_fragments - subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag) + subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag, lmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum @@ -253,6 +255,7 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_ real(DP), dimension(:), intent(in) :: m_frag, rad_frag !! Fragment masses and radii real(DP), dimension(:,:), intent(in) :: Ip_frag, x_frag !! Fragment prinicpal moments of inertia and position vectors real(DP), dimension(:,:), intent(out) :: v_frag, rot_frag !! Fragment velocities and spin states + logical, intent(out) :: lmerge ! Internals integer(I4B) :: i real(DP) :: L_orb_mag @@ -308,7 +311,8 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_ end do b(1:3) = -L_lin_others(:) b(4:6) = L_orb_old(:) - L_orb_others(:) - v_t_mag(1:6) = util_solve_linear_system(6, A, b) + v_t_mag(1:6) = util_solve_linear_system(A, b, 6, lmerge) + if (lmerge) return do i = 1, 6 v_frag(:, i) = v_t_mag(i) * v_t_unit(:, i) end do diff --git a/src/util/util_bfgs.f90 b/src/util/util_minimize_bfgs.f90 similarity index 82% rename from src/util/util_bfgs.f90 rename to src/util/util_minimize_bfgs.f90 index 85275cd32..069fb5b35 100644 --- a/src/util/util_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -1,4 +1,4 @@ -function util_bfgs(f, N, x1, eps) result(fnum) +function util_minimize_bfgs(f, N, x1, eps) result(fnum) !! author: David A. Minton !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. @@ -15,7 +15,7 @@ function util_bfgs(f, N, x1, eps) result(fnum) !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - use swiftest use swiftest_globals - use module_interfaces, EXCEPT_THIS_ONE => util_bfgs + use module_interfaces, EXCEPT_THIS_ONE => util_minimize_bfgs implicit none ! Arguments integer(I4B), intent(in) :: N @@ -447,20 +447,20 @@ 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 - !! to some accuracy eps. It recieves as input: - !! f(x) : function of one real(DP) :: variable as input - !! lo : low bracket value - !! hi : high bracket value - !! eps : desired accuracy of final minimum location - !! The outputs include - !! lo : final minimum location - !! hi : final minimum location - !! Returns - !! Number of function calls performed or - !! 0 = No miniumum found - !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! This function uses a quadratic polynomial fit to * locate the minimum of a function + !! to some accuracy eps. It recieves as input: + !! f(x) : function of one real(DP) :: variable as input + !! lo : low bracket value + !! hi : high bracket value + !! eps : desired accuracy of final minimum location + !! The outputs include + !! lo : final minimum location + !! hi : final minimum location + !! Returns + !! Number of function calls performed or + !! 0 = No miniumum found + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! Arguments integer(I4B), intent(in) :: N @@ -485,6 +485,7 @@ end function f real(DP), dimension(3,3) :: lhs real(DP) :: d1, d2, d3, aold, denom integer(I4B) :: i + logical :: lerr ! Get initial a1, a2, a3 values a1 = lo @@ -518,12 +519,13 @@ end function f lhs(2, :) = row_2 lhs(3, :) = row_3 ! Solve system of equations - if (gauss(lhs, rhs, soln, 3) /= 0) then - !write(*,*) "Could not solve polynomial on loop %d", i) - !write(*,*) "a1 = %9.6f f1 = %9.6f", a1, f1) - !write(*,*) "a2 = %9.6f f2 = %9.6f", a2, f2) - !write(*,*) "a3 = %9.6f f3 = %9.6f", a3, f3) - !write(*,*) "aold = %7.4f", aold) + 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(*,'("aold = ",f7.4)') aold fnum = 0 return end if @@ -560,89 +562,5 @@ end function f return end function quadfit - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - function gauss(a1, b, soln, N) result(errcode) - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! Purpose: Subroutine to solve a set of n linear equations in n unknowns using - !! Gaussian elimination and the maximum pivot techniques. - !! - !! From: Chapman, Stephen J., _Fortran 90 / 95 for Scientists and Engineers_, Firs - !! Edition, WCB McGraw - Hill, 1998, pgs. 433 - 435 - !! - !! Converted to C by David Minton - !! - !! Input - !! a : NxN input array - !! b : RHS vector - !! n : number of equations - !! Output - !! soln : solution vector - !! Returns - !! 0 : No error - !! 1 : Singular matrix - !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none - ! Arguments - integer(I4B), intent(in) :: N - interface - pure function f(x) ! Objective function template - import DP - real(DP), dimension(:), intent(in) :: x - real(DP) :: f - end function f - end interface - real(DP), dimension(:,:), intent(inout) :: a1 - real(DP), dimension(:), intent(in) :: b - real(DP), dimension(:), intent(out) :: soln - ! Result - integer(I4B) :: errcode - ! Internals - real(DP) :: factor, temp - real(DP), parameter :: epsilon = tiny(temp) ! A "small" number for comparison when determining singular matrix - integer(I4B) :: irow, jrow, ipeak, kcol - - soln(:) = b(:) - do irow = 1, N - ! Find peak pivot for column irow in rows irow to n - ipeak = irow - do jrow = irow + 1, N - if (abs(a1(jrow, irow)) > abs(a1(ipeak, irow))) ipeak = jrow - end do - - ! Check for singular equations - if (abs(a1(ipeak, irow)) < epsilon) then - errcode = 1 - return - end if - - ! Otherwise, if ipeak ! = irow, swap rows irow & ipeak - if (ipeak /= irow) then - do kcol = 1, N - temp = a1(ipeak, kcol) - a1(ipeak, kcol) = a1(irow, kcol) - a1(irow, kcol) = temp - end do - temp = soln(ipeak) - soln(ipeak) = soln(irow) - soln(irow) = temp - end if - - do jrow = 1, N - if (jrow /= irow) then - factor = -(a1(jrow,irow)) / (a1(irow, irow)) - a1(jrow, :) = a1(irow, :) * factor + a1(jrow, :) - soln(jrow) = soln(irow) * factor + soln(jrow) - end if - end do - end do - do irow = 1, N - soln(irow) = soln(irow) / (a1(irow, irow)) - end do - errcode = 0 - return - - end function gauss -end function util_bfgs \ No newline at end of file +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 cd42a5578..1f12c0711 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(n,A,b) result(x) +function util_solve_linear_system(A,b,n,lerr) result(x) !! Author: David A. Minton !! !! Solves the linear equation of the form A*x = b for x. @@ -13,27 +13,32 @@ function util_solve_linear_system(n,A,b) result(x) 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 ! Internals real(QP), dimension(:), allocatable :: qx - qx = solve_wbs(ge_wpp(real(A, kind=QP), real(b, kind=QP))) + qx = solve_wbs(ge_wpp(real(A, kind=QP), real(b, kind=QP)),lerr) x = real(qx, kind=DP) return - contains - function solve_wbs(u) result(x) ! solve with backward substitution + 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) + + lerr = any(abs(u(:,:)) < epsilon) + if (lerr) return n = size(u,1) allocate(x(n)) @@ -43,7 +48,7 @@ function solve_wbs(u) result(x) ! solve with backward substitution return end function solve_wbs - function ge_wpp(A, b) result(u) ! gaussian eliminate with partial pivoting + 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.