Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Browse files Browse the repository at this point in the history 
Refactoring for clarity. Replaced old Gaussian elimination code with newer one, and incorporated more error checking
  • Loading branch information
@daminton
daminton committed May 13, 2021
1 parent a84a01c commit f89f2e2
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 121 deletions.
9 changes: 5 additions & 4 deletions src/modules/module_interfaces.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
14 changes: 9 additions & 5 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
132 changes: 25 additions & 107 deletions src/util/util_bfgs.f90 → src/util/util_minimize_bfgs.f90
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
end function util_minimize_bfgs
15 changes: 10 additions & 5 deletions src/util/util_solve_linear_system.f90
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand All @@ -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))
Expand All @@ -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.
Expand Down

0 comments on commit f89f2e2

Please sign in to comment.