|
| 1 | +!--------------------------------------------------------------------------------------------------! |
| 2 | +! The MEX gateway for COBYLA |
| 3 | +! This is a temporary MEX gateway to circumvent the MATLAB R2025a bug that it segfaults when |
| 4 | +! the Fortran MEX function contains an internal procedure that is passed as an actual argument. |
| 5 | +! This MEX uses a module variable FUN_PTR to store the function handle, which is essentially a |
| 6 | +! global variable and is not thread-safe or recursion-safe. |
| 7 | +! See MathWorks Technical Support Case 07931486 and |
| 8 | +! https://www.mathworks.com/matlabcentral/answers/2178414-bug-matlab-2025a-segfaults-on-ubuntu-when-handling-fortran-mex-files-with-internal-subroutines |
| 9 | +! https://stackoverflow.com/questions/79699706/matlab-2025a-vs-fortran-mex-files-with-internal-subroutines |
| 10 | +! https://fortran-lang.discourse.group/t/implementation-of-a-parametrized-objective-function-without-using-module-variables-or-internal-subroutines |
| 11 | +! https://stackoverflow.com/questions/79705107/fortran-implementating-a-parametrized-objective-function-without-using-module-v |
| 12 | +! |
| 13 | +! Authors: |
| 14 | +! Tom M. RAGONNEAU ([email protected]) |
| 15 | +! and Zaikun ZHANG ([email protected]) |
| 16 | +! Department of Applied Mathematics, |
| 17 | +! The Hong Kong Polytechnic University |
| 18 | +! |
| 19 | +! Dedicated to the late Professor M. J. D. Powell FRS (1936--2015) |
| 20 | +! |
| 21 | +! Started in July 2020 |
| 22 | +! |
| 23 | +! Last Modified: Sat 19 Jul 2025 11:40:45 PM PDT |
| 24 | +!--------------------------------------------------------------------------------------------------! |
| 25 | + |
| 26 | +#include "fintrf.h" |
| 27 | + |
| 28 | + |
| 29 | +module calcfc_mod |
| 30 | +implicit none |
| 31 | +private |
| 32 | +public :: funcon_ptr, calcfc |
| 33 | + |
| 34 | +mwPointer :: funcon_ptr ! Pointer to the objective function handle |
| 35 | + |
| 36 | +contains |
| 37 | + |
| 38 | +subroutine calcfc(x, f, nlconstr) |
| 39 | +use, non_intrinsic :: consts_mod, only : RP |
| 40 | +use, non_intrinsic :: cbfun_mod, only : evalcb |
| 41 | +implicit none |
| 42 | +real(RP), intent(in) :: x(:) |
| 43 | +real(RP), intent(out) :: f |
| 44 | +real(RP), intent(out) :: nlconstr(:) |
| 45 | +call evalcb(funcon_ptr, x, f, nlconstr) |
| 46 | +end subroutine calcfc |
| 47 | + |
| 48 | +end module calcfc_mod |
| 49 | + |
| 50 | + |
| 51 | +subroutine mexFunction(nargout, poutput, nargin, pinput) |
| 52 | +!--------------------------------------------------------------------------------------------------! |
| 53 | +! This is the entry point to the Fortran MEX function. If the compiled MEX file is named as |
| 54 | +! FUNCTION_NAME.mex*** (extension depends on the platform), then in MATLAB we can call: |
| 55 | +! [x, f, cstrv, nlconstr, info, nf, xhist, fhist, chist, nlchist] = ... |
| 56 | +! FUNCTION_NAME(funcon, x0, f0, nlconstr0, Aineq, bineq, Aeq, beq, lb, ub, rhobeg, rhoend, ... |
| 57 | +! eta1, eta2, gamma1, gamma2, ftarget, ctol, cweight, maxfun, iprint, maxhist, output_xhist, ... |
| 58 | +! output_nlchist, maxfilt) |
| 59 | +!--------------------------------------------------------------------------------------------------! |
| 60 | + |
| 61 | +! Generic modules |
| 62 | +use, non_intrinsic :: consts_mod, only : RP, IK |
| 63 | +use, non_intrinsic :: memory_mod, only : safealloc |
| 64 | + |
| 65 | +! Fortran MEX API modules |
| 66 | +use, non_intrinsic :: fmxapi_mod, only : fmxVerifyNArgin, fmxVerifyNArgout |
| 67 | +use, non_intrinsic :: fmxapi_mod, only : fmxVerifyClassShape |
| 68 | +use, non_intrinsic :: fmxapi_mod, only : fmxReadMPtr, fmxWriteMPtr |
| 69 | + |
| 70 | +! Solver-specific modules |
| 71 | +use, non_intrinsic :: calcfc_mod, only : funcon_ptr, calcfc |
| 72 | +use, non_intrinsic :: cobyla_mod, only : cobyla |
| 73 | + |
| 74 | +implicit none |
| 75 | + |
| 76 | +! mexFunction arguments nargout and nargin are of type INTEGER in MATLAB 2019a documents. |
| 77 | +integer, intent(in) :: nargout, nargin |
| 78 | +mwPointer, intent(in) :: pinput(nargin) |
| 79 | +mwPointer, intent(out) :: poutput(nargout) |
| 80 | + |
| 81 | +! Local variables |
| 82 | +integer(IK) :: info |
| 83 | +integer(IK) :: iprint |
| 84 | +integer(IK) :: m_nlcon |
| 85 | +integer(IK) :: maxfilt |
| 86 | +integer(IK) :: maxfun |
| 87 | +integer(IK) :: maxhist |
| 88 | +integer(IK) :: nf |
| 89 | +logical :: output_nlchist |
| 90 | +logical :: output_xhist |
| 91 | +real(RP) :: cstrv |
| 92 | +real(RP) :: ctol |
| 93 | +real(RP) :: cweight |
| 94 | +real(RP) :: eta1 |
| 95 | +real(RP) :: eta2 |
| 96 | +real(RP) :: f |
| 97 | +real(RP) :: f0 |
| 98 | +real(RP) :: ftarget |
| 99 | +real(RP) :: gamma1 |
| 100 | +real(RP) :: gamma2 |
| 101 | +real(RP) :: rhobeg |
| 102 | +real(RP) :: rhoend |
| 103 | +real(RP), allocatable :: Aeq(:, :) |
| 104 | +real(RP), allocatable :: Aineq(:, :) |
| 105 | +real(RP), allocatable :: beq(:) |
| 106 | +real(RP), allocatable :: bineq(:) |
| 107 | +real(RP), allocatable :: chist(:) |
| 108 | +real(RP), allocatable :: nlchist(:, :) |
| 109 | +real(RP), allocatable :: nlconstr(:) |
| 110 | +real(RP), allocatable :: nlconstr0(:) |
| 111 | +real(RP), allocatable :: fhist(:) |
| 112 | +real(RP), allocatable :: lb(:) |
| 113 | +real(RP), allocatable :: ub(:) |
| 114 | +real(RP), allocatable :: x(:) |
| 115 | +real(RP), allocatable :: xhist(:, :) |
| 116 | + |
| 117 | +! Validate the number of arguments |
| 118 | +call fmxVerifyNArgin(nargin, 25) |
| 119 | +call fmxVerifyNArgout(nargout, 10) |
| 120 | + |
| 121 | +! Verify that input 1 is a function handle; the other inputs will be verified when read. |
| 122 | +call fmxVerifyClassShape(pinput(1), 'function_handle', 'rank0') |
| 123 | + |
| 124 | +! Read the inputs |
| 125 | +funcon_ptr = pinput(1) ! FUNCON_PTR is a pointer to the function handle |
| 126 | +call fmxReadMPtr(pinput(2), x) |
| 127 | +call fmxReadMPtr(pinput(3), f0) |
| 128 | +call fmxReadMPtr(pinput(4), nlconstr0) |
| 129 | +call fmxReadMPtr(pinput(5), Aineq) |
| 130 | +call fmxReadMPtr(pinput(6), bineq) |
| 131 | +call fmxReadMPtr(pinput(7), Aeq) |
| 132 | +call fmxReadMPtr(pinput(8), beq) |
| 133 | +call fmxReadMPtr(pinput(9), lb) |
| 134 | +call fmxReadMPtr(pinput(10), ub) |
| 135 | +call fmxReadMPtr(pinput(11), rhobeg) |
| 136 | +call fmxReadMPtr(pinput(12), rhoend) |
| 137 | +call fmxReadMPtr(pinput(13), eta1) |
| 138 | +call fmxReadMPtr(pinput(14), eta2) |
| 139 | +call fmxReadMPtr(pinput(15), gamma1) |
| 140 | +call fmxReadMPtr(pinput(16), gamma2) |
| 141 | +call fmxReadMPtr(pinput(17), ftarget) |
| 142 | +call fmxReadMPtr(pinput(18), ctol) |
| 143 | +call fmxReadMPtr(pinput(19), cweight) |
| 144 | +call fmxReadMPtr(pinput(20), maxfun) |
| 145 | +call fmxReadMPtr(pinput(21), iprint) |
| 146 | +call fmxReadMPtr(pinput(22), maxhist) |
| 147 | +call fmxReadMPtr(pinput(23), output_xhist) |
| 148 | +call fmxReadMPtr(pinput(24), output_nlchist) |
| 149 | +call fmxReadMPtr(pinput(25), maxfilt) |
| 150 | + |
| 151 | +! Get the sizes |
| 152 | +m_nlcon = int(size(nlconstr0), kind(m_nlcon)) ! M_NLCON is a compulsory input of the Fortran code. |
| 153 | + |
| 154 | +! Allocate memory for nlconstr |
| 155 | +call safealloc(nlconstr, m_nlcon) |
| 156 | + |
| 157 | +! Call the Fortran code |
| 158 | +! There are different cases because XHIST/CONHIST may or may not be passed to the Fortran code. |
| 159 | +if (output_xhist .and. output_nlchist) then |
| 160 | + call cobyla(calcfc, m_nlcon, x, f, cstrv, nlconstr, Aineq, bineq, Aeq, beq, lb, ub, & |
| 161 | + & f0, nlconstr0, nf, rhobeg, rhoend, ftarget, ctol, cweight, maxfun, iprint, eta1, eta2, & |
| 162 | + & gamma1, gamma2, xhist = xhist, fhist = fhist, chist = chist, nlchist = nlchist, maxhist = maxhist, & |
| 163 | + & maxfilt = maxfilt, info = info) |
| 164 | +elseif (output_xhist) then |
| 165 | + call cobyla(calcfc, m_nlcon, x, f, cstrv, nlconstr, Aineq, bineq, Aeq, beq, lb, ub, & |
| 166 | + & f0, nlconstr0, nf, rhobeg, rhoend, ftarget, ctol, cweight, maxfun, iprint, eta1, eta2, & |
| 167 | + & gamma1, gamma2, xhist = xhist, fhist = fhist, chist = chist, maxhist = maxhist, maxfilt = maxfilt, & |
| 168 | + & info = info) |
| 169 | +elseif (output_nlchist) then |
| 170 | + call cobyla(calcfc, m_nlcon, x, f, cstrv, nlconstr, Aineq, bineq, Aeq, beq, lb, ub, & |
| 171 | + & f0, nlconstr0, nf, rhobeg, rhoend, ftarget, ctol, cweight, maxfun, iprint, eta1, eta2, & |
| 172 | + & gamma1, gamma2, fhist = fhist, chist = chist, nlchist = nlchist, maxhist = maxhist, & |
| 173 | + & maxfilt = maxfilt, info = info) |
| 174 | +else |
| 175 | + call cobyla(calcfc, m_nlcon, x, f, cstrv, nlconstr, Aineq, bineq, Aeq, beq, lb, ub, & |
| 176 | + & f0, nlconstr0, nf, rhobeg, rhoend, ftarget, ctol, cweight, maxfun, iprint, eta1, eta2, & |
| 177 | + & gamma1, gamma2, fhist = fhist, chist = chist, maxhist = maxhist, maxfilt = maxfilt, info = info) |
| 178 | +end if |
| 179 | + |
| 180 | +! After the Fortran code, XHIST or CONHIST may not be allocated, because it may not have been passed |
| 181 | +! to the Fortran code. We allocate it here. Otherwise, fmxWriteMPtr will fail. |
| 182 | +if (.not. allocated(xhist)) then |
| 183 | + call safealloc(xhist, int(size(x), IK), 0_IK) |
| 184 | +end if |
| 185 | +if (.not. allocated(nlchist)) then |
| 186 | + call safealloc(nlchist, int(size(nlconstr), IK), 0_IK) |
| 187 | +end if |
| 188 | + |
| 189 | +! Write the outputs |
| 190 | +call fmxWriteMPtr(x, poutput(1)) |
| 191 | +call fmxWriteMPtr(f, poutput(2)) |
| 192 | +call fmxWriteMPtr(cstrv, poutput(3)) |
| 193 | +call fmxWriteMPtr(nlconstr, poutput(4)) |
| 194 | +call fmxWriteMPtr(info, poutput(5)) |
| 195 | +call fmxWriteMPtr(nf, poutput(6)) |
| 196 | +call fmxWriteMPtr(xhist(:, 1:min(nf, int(size(xhist, 2), IK))), poutput(7)) |
| 197 | +call fmxWriteMPtr(fhist(1:min(nf, int(size(fhist), IK))), poutput(8), 'row') |
| 198 | +call fmxWriteMPtr(chist(1:min(nf, int(size(chist), IK))), poutput(9), 'row') |
| 199 | +call fmxWriteMPtr(nlchist(:, 1:min(nf, int(size(nlchist, 2), IK))), poutput(10)) |
| 200 | +! N.B.: |
| 201 | +! It can happen that 0 < SIZE(XHIST, 2) < MAXHIST or 0 < SIZE(FHIST) < MAXHIST due to the memory |
| 202 | +! limit in the Fortran code. Similar for CHIST and CONHIST. |
| 203 | + |
| 204 | +! Free memory. We prefer explicit deallocation to the automatic one. |
| 205 | +deallocate (x) ! Allocated by fmxReadMPtr. |
| 206 | +deallocate (nlconstr0) ! Allocated by fmxReadMPtr. |
| 207 | +deallocate (Aineq) ! Allocated by fmxReadMPtr. |
| 208 | +deallocate (bineq) ! Allocated by fmxReadMPtr. |
| 209 | +deallocate (Aeq) ! Allocated by fmxReadMPtr. |
| 210 | +deallocate (beq) ! Allocated by fmxReadMPtr. |
| 211 | +deallocate (lb) ! Allocated by fmxReadMPtr. |
| 212 | +deallocate (ub) ! Allocated by fmxReadMPtr. |
| 213 | +deallocate (nlconstr) ! Allocated manually |
| 214 | +deallocate (xhist) ! Allocated by the solver |
| 215 | +deallocate (fhist) ! Allocated by the solver |
| 216 | +deallocate (chist) ! Allocated by the solver |
| 217 | +deallocate (nlchist) ! Allocated by the solver |
| 218 | + |
| 219 | +end subroutine mexFunction |
0 commit comments