Improved floating point exception handling in both the BFGS and linea…

…r solver routines using the ieee_exceptions module.

src/modules/module_interfaces.f90

@@ -1590,11 +1590,10 @@ function util_solve_linear_system_q(A,b,n,lerr) result(x)
             real(QP), dimension(n)                :: x
          end function util_solve_linear_system_q
 
-         function solve_wbs(u, lerr) result(x) 
+         function solve_wbs(u) 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
 

src/util/util_minimize_bfgs.f90

@@ -13,6 +13,7 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
    !!   x1   :  Final minimum (all 0 if none found)
    !!   0 = No miniumum found
    !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  - 
+   use, intrinsic :: ieee_exceptions
    use swiftest
    use swiftest_globals
    use module_interfaces, EXCEPT_THIS_ONE => util_minimize_bfgs
@@ -38,6 +39,12 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
    real(QP), dimension(N) :: y, P, x0, x1
    real(QP), dimension(N,N) :: PP, PyH, HyP
    real(QP) :: yHy, Py, eps
+   type(ieee_status_type) :: original_fpe_status
+   logical, dimension(:), allocatable :: fpe_flag 
+
+   call ieee_get_status(original_fpe_status) ! Save the original floating point exception status
+   call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet
+   allocate(fpe_flag(size(ieee_usual)))
 
    lerr = .false.
    eps = real(eps_d, kind=QP)
@@ -59,16 +66,14 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
       end do
       if (conv == 0) then
          write(*,*) "Converged on gradient after ",i," iterations" 
-         x1_d = x1
-         return 
+         exit 
       end if
       ! normalize gradient 
       Snorm(:) = S(:) / norm2(S)
       astar = minimize1D(f, x1, Snorm, N, gradeps, lerr)
       if (lerr) then
          write(*,*) "Exiting BFGS with error in minimize1D step"
-         x1_d = x1
-         return 
+         exit
       end if
       ! Get new x values 
       P(:) = astar * Snorm(:) 
@@ -88,8 +93,7 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
       ! prevent divide by zero (convergence) 
       if (abs(Py) < tiny(Py)) then
          write(*,*) "Converged on tiny Py after ",i," iterations"
-         x1_d = x1
-         return
+         exit
       end if
       ! set up update 
       PyH(:,:) = 0._QP
@@ -107,12 +111,16 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
       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
+      if (i == MAXLOOP) write(*,*) "BFGS ran out of loops!"
    end do
-   write(*,*) "Did not converge! Ran out of loops"
    x1_d = x1
+   call ieee_get_flag(ieee_usual, fpe_flag)
+   lerr = lerr .or. any(fpe_flag)  
+   if (lerr) write(*,*) "BFGS did not converge!"
+   call ieee_set_status(original_fpe_status)
+
    return 
 
    contains
@@ -439,9 +447,7 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr)
          !! The outputs include
          !!   lo   :  final minimum location
          !!   hi   :  final minimum location
-         !! Returns
-         !!   Number of function calls performed or
-         !!   0 = No miniumum found
+         !! Notes: Uses the ieee_exceptions intrinsic module to allow for graceful failure due to floating point exceptions, which won't terminate the run.
          !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  -   
          implicit none
          ! Arguments
@@ -454,12 +460,11 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr)
          logical,                intent(out)   :: lerr
          ! Internals 
          integer(I4B), parameter :: MAXLOOP = 20 ! maximum number of loops before method is determined to have failed.   
-         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
+         real(QP) :: d1, d2, d3, aold, denom, errval
          integer(I4B) :: i
 
          lerr = .false.
@@ -474,15 +479,16 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr)
          f3 = n2one(f, x0, S, N, a3)
          do i = 1, MAXLOOP 
             ! check to see if convergence is reached and exit   
-            if (abs(astar) < small) then
-               denom = small 
-            else 
-               denom = astar
+            errval = abs((astar - aold) / astar)
+            call ieee_get_flag(ieee_usual, fpe_flag)
+            if (any(fpe_flag)) then
+               lerr = .true.
+               exit
             end if
-            if (abs((astar - aold) / denom) < eps) then
+            if (errval < eps) then
                lo = astar
                hi = astar
-               return 
+               exit
             end if
             ! Set up system for gaussian elimination equation solver   
             row_1 = [1.0_QP, a1, a1**2]
@@ -494,12 +500,14 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr)
             lhs(3, :) = row_3
             ! Solve system of equations   
             soln(:) = util_solve_linear_system(lhs, rhs, 3, lerr)
-            if (lerr) then
-               return 
-            end if
+            if (lerr) exit
             aold = astar
-            if (abs(soln(3)) < small) soln(3) = small
             astar =  -soln(2) / (2 * soln(3))
+            call ieee_get_flag(ieee_usual, fpe_flag)
+            if (any(fpe_flag)) then
+               lerr = .true.
+               exit
+            end if
             fstar = n2one(f, x0, S, N, astar)
             ! keep the three closest a values to astar and discard the fourth  
             d1 = abs(a1 - astar)
@@ -524,6 +532,7 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr)
                end if
             end if
          end do
+         if (lerr) return
          lo = a1
          hi = a3
          return 

src/util/util_solve_linear_system.f90

@@ -6,6 +6,7 @@ function util_solve_linear_system_d(A,b,n,lerr) result(x)
    !!   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, intrinsic :: ieee_exceptions
    use swiftest
    use module_interfaces, EXCEPT_THIS_ONE => util_solve_linear_system_d
    implicit none
@@ -18,10 +19,23 @@ function util_solve_linear_system_d(A,b,n,lerr) result(x)
    real(DP), dimension(n)                :: x
    ! Internals
    real(QP), dimension(:), allocatable :: qx
+   type(ieee_status_type) :: original_fpe_status
+   logical, dimension(:), allocatable :: fpe_flag 
 
-   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)
+   call ieee_get_status(original_fpe_status) ! Save the original floating point exception status
+   call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet
+   allocate(fpe_flag(size(ieee_usual)))
+
+   qx = solve_wbs(ge_wpp(real(A, kind=QP), real(b, kind=QP)))
+
+   call ieee_get_flag(ieee_usual, fpe_flag)
+   lerr = any(fpe_flag) 
+   if (lerr) then
+      x = 0.0_DP
+   else
+      x = real(qx, kind=DP)
+   end if
+   call ieee_set_status(original_fpe_status)
 
    return
 end function util_solve_linear_system_d
@@ -34,6 +48,7 @@ function util_solve_linear_system_q(A,b,n,lerr) result(x)
    !!   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, intrinsic :: ieee_exceptions
    use swiftest
    use module_interfaces, EXCEPT_THIS_ONE => util_solve_linear_system_q
    implicit none
@@ -44,35 +59,37 @@ function util_solve_linear_system_q(A,b,n,lerr) result(x)
    logical,                  intent(out) :: lerr
    ! Result
    real(QP), dimension(n)  :: x
+   type(ieee_status_type) :: original_fpe_status
+   logical, dimension(:), allocatable :: fpe_flag 
+
+   call ieee_get_status(original_fpe_status) ! Save the original floating point exception status
+   call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet
+   allocate(fpe_flag(size(ieee_usual)))
 
-   x = solve_wbs(ge_wpp(A, b),lerr)
-   where(abs(x(:)) < tiny(1._QP)) x(:) = 0._QP
+   x = solve_wbs(ge_wpp(A, b))
+
+   call ieee_get_flag(ieee_usual, fpe_flag)
+   lerr = any(fpe_flag) 
+   if (lerr) x = 0.0_DP
+   call ieee_set_status(original_fpe_status) 
 
    return
 end function util_solve_linear_system_q
 
-function solve_wbs(u, lerr) result(x) ! solve with backward substitution
+function solve_wbs(u) 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) 
 
    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
-
    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