|
1 | 1 | module hompack_core |
2 | 2 | !! Core routines used by two or more of the `hompack` solvers. |
3 | 3 |
|
| 4 | + use iso_c_binding, only: c_ptr, c_null_ptr |
4 | 5 | use hompack_kinds, only: dp |
5 | 6 | implicit none |
| 7 | + private |
6 | 8 |
|
7 | | - type root_state |
| 9 | + public :: hompack_status, hompack_f_callbacks |
| 10 | + public :: root_state, root, qofs |
| 11 | + |
| 12 | + ! Problem "enum" |
| 13 | + integer, parameter, public :: problem_fix_point = 1 |
| 14 | + integer, parameter, public :: problem_zero_find = 2 |
| 15 | + integer, parameter, public :: problem_curve_track = 3 |
| 16 | + |
| 17 | + ! Status "enum" |
| 18 | + integer, parameter, public :: hompack_success = 1 |
| 19 | + |
| 20 | + abstract interface |
| 21 | + subroutine f_t(x, v, data) |
| 22 | + !! User routine to evaluate the function \(f(x)\). |
| 23 | + import :: dp, c_ptr |
| 24 | + real(dp), intent(in) :: x(:) |
| 25 | + real(dp), intent(out) :: v(:) |
| 26 | + type(c_ptr), value :: data |
| 27 | + end subroutine f_t |
| 28 | + |
| 29 | + subroutine fjac_t(x, v, k, data) |
| 30 | + !! User routine to evaluate the Jacobian of \(f(x)\). |
| 31 | + import :: dp, c_ptr |
| 32 | + real(dp), intent(in) :: x(:) |
| 33 | + real(dp), intent(out) :: v(:) |
| 34 | + integer, intent(in) :: k |
| 35 | + type(c_ptr), value :: data |
| 36 | + end subroutine fjac_t |
| 37 | + |
| 38 | + subroutine rho_t(a, lambda, x, v, data) |
| 39 | + !! User routine to evaluate the homotopy function \(\rho(a, \lambda, x)\). |
| 40 | + import :: dp, c_ptr |
| 41 | + real(dp), intent(in) :: a(:) |
| 42 | + real(dp), intent(in) :: lambda |
| 43 | + real(dp), intent(in) :: x(:) |
| 44 | + real(dp), intent(out) :: v(:) |
| 45 | + type(c_ptr), value :: data |
| 46 | + end subroutine rho_t |
| 47 | + |
| 48 | + subroutine rhojac_t(a, lambda, x, v, k, data) |
| 49 | + !! User routine to evaluate the Jacobian of \(\rho(a, \lambda, x)\). |
| 50 | + import :: dp, c_ptr |
| 51 | + real(dp), intent(in) :: a(:) |
| 52 | + real(dp), intent(in) :: lambda |
| 53 | + real(dp), intent(in) :: x(:) |
| 54 | + real(dp), intent(out) :: v(:) |
| 55 | + integer, intent(in) :: k |
| 56 | + type(c_ptr), value :: data |
| 57 | + end subroutine rhojac_t |
| 58 | + end interface |
| 59 | + |
| 60 | + type :: hompack_f_callbacks |
| 61 | + !! Container to hold user routines for f-suffix (dense Jacobian) methods. |
| 62 | + procedure(f_t), nopass, pointer :: f => null() |
| 63 | + procedure(fjac_t), nopass, pointer :: fjac => null() |
| 64 | + procedure(rho_t), nopass, pointer :: rho => null() |
| 65 | + procedure(rhojac_t), nopass, pointer :: rhojac => null() |
| 66 | + type(c_ptr) :: data = c_null_ptr |
| 67 | + end type |
| 68 | + |
| 69 | + type :: hompack_s_callbacks |
| 70 | + !! Container to hold user routines for s-suffix (sparse Jacobian) methods. |
| 71 | + !! To be implemented in the future. |
| 72 | + procedure(f_t), nopass, pointer :: f => null() |
| 73 | + ! procedure(fjacs_t), nopass, pointer :: fjacs => null() |
| 74 | + procedure(rho_t), nopass, pointer :: rho => null() |
| 75 | + ! procedure(rhojacs_t), nopass, pointer :: rhojacs => null() |
| 76 | + type(c_ptr) :: data = c_null_ptr |
| 77 | + end type |
| 78 | + |
| 79 | + type :: hompack_status |
| 80 | + integer :: code = hompack_success |
| 81 | + character(len=256) :: message = "" |
| 82 | + contains |
| 83 | + procedure :: ok => hompack_status_ok |
| 84 | + end type |
| 85 | + |
| 86 | + type :: root_state |
8 | 87 | !! State variables for [[root]]. |
9 | 88 | real(dp) :: a |
10 | 89 | !! Previous iterate point used for secant step calculation. |
@@ -36,10 +115,20 @@ module hompack_core |
36 | 115 | !! Counter of sequential secant steps used to force bisection if convergence stalls. |
37 | 116 | integer :: fcount |
38 | 117 | !! Total number of function evaluations completed so far. |
39 | | - end type root_state |
| 118 | + end type |
40 | 119 |
|
41 | 120 | contains |
42 | 121 |
|
| 122 | + logical pure function hompack_status_ok(self) result(res) |
| 123 | + !! This function checks if the status code indicates success. |
| 124 | + |
| 125 | + class(hompack_status), intent(in) :: self |
| 126 | + !! Status object. |
| 127 | + |
| 128 | + res = (self%code == hompack_success) |
| 129 | + |
| 130 | + end function |
| 131 | + |
43 | 132 | subroutine root(t, ft, b, c, relerr, abserr, iflag, state) |
44 | 133 | !! This subroutine computes a root of the nonlinear equation `f(x) = 0` where `f(x)` |
45 | 134 | !! is a continuous real function of a single real variable `x`. The method used is a |
@@ -235,7 +324,7 @@ subroutine root(t, ft, b, c, relerr, abserr, iflag, state) |
235 | 324 | end subroutine root |
236 | 325 |
|
237 | 326 | elemental pure function qofs(f0, fp0, f1, fp1, dels, s) result(res) |
238 | | - !! Computes the Hermite cubic interpolant at a point `s`. |
| 327 | + !! This function computes the Hermite cubic interpolant at a point `s`. |
239 | 328 |
|
240 | 329 | real(dp), intent(in) :: f0 |
241 | 330 | !! Function value at the start of the interval. |
|
0 commit comments