RECURSIVE SUBROUTINE powell(nn, bl, bu, yy, val, n_iuser, iuser, n_ruser, ruser, istat) !************************************************************* ! ! Minimisation routine in multiple dimensions without derivatives ! ! Based on Powell's method as recommended by Numerical Recipes, 1987, p. 299. ! ! INPUTS: ! nn = number of dimensions to optimise over ! bl = lower limit of optimisation space ! bu = upper limit of optimisation space ! n_iuser = dimension of vector of integer parameter inputs fed into routine ! iuser = vector of integer parameter inputs fed into routine ! n_ruser = dimension of vector of real parameter inputs fed into routine ! ruser = vector of real parameter inputs fed into routine ! OUTPUTS: ! yy = optimised co-ordinate ! val = optimised value of objective function ! istat = output indicator code (set to 0 here) ! ! soubroutine objfn(nn, BX0, FX0, nstate, n_iuser, iuser, n_ruser, ruser, inform) returns ! the value of the objective function at ordinate BX0, as variable FX0 ! ! subroutine brent optimises over a single dimension without reference to derivatives ! ! subroutine boundary_pnts(nn, AX, CX, xit, BX, AAX, CCX, bound_chk) evaluates the boundary ! points of the linear projection through BX on the search space defined by AX and CX ! INPUTS: ! nn = number of dimensions of optimisation space ! AX = lower co-ordinate of search space ! CX = upper co-ordinate of search space ! xit = gradient of linear segment ! BX = point on linear segment ! OUTPUTS: ! AAX = lower intersection of linear segment with search domain ! CCX = upper intersection of linear segment with search domain ! bound_chk = indicator of corner co-ordinate ! ! code developed for the SIDD microsimulation model ! ! Author: J. van de Ven ! Date: 02 Aug 2016 ! !************************************************************* USE global_ParamStore USE global_CommonAll IMPLICIT none integer(4), intent(in) :: nn, n_iuser, n_ruser integer(4) :: istat integer(4) :: iuser(n_iuser) real(8), intent(in) :: bl(nn), bu(nn) real(8) :: val real(8) :: yy(nn), ruser(n_ruser) ! local variables integer(4), parameter :: itmax=200, nn_max = 3 real(8), parameter :: ftol = (epsilon(1.0_8))**(0.5), eps = 1.0E-10, eps_v = epsilon(1.0_8) integer(4) :: nstate, inform integer(4) :: ii, jj, iter, ibig, finish, bound_chk real(8) :: xit(nn), AX(nn), BX(nn), CX(nn), AAX(nn), CCX(nn), FX, FE real(8) :: xi(nn,nn), pt(nn), ptt(nn) real(8) :: fp, fret, fptt, del, T_val, adjust_factor !************************************************************* ! begin code !************************************************************* ii = 0 AX = bl BX = yy CX = bu nstate = 1 call objfn(nn, BX, FX, nstate, n_iuser, iuser, n_ruser, ruser, inform) nstate = 0 fret = FX ! initialise xi to the unit matrix xi = 0.0 xit = 0.0 do ii = 1, nn xi(ii,ii) = 1.0 end do finish = 0 iter = 0 do while ( finish.eq.0 ) iter = iter+1 fp = fret pt = BX ibig = 0 del = 0.0 do ii = 1, nn do jj = 1, nn xit(jj) = xi(jj,ii) end do fptt = fret call boundary_pnts(nn, AX, CX, xit, BX, AAX, CCX, bound_chk) if ( bound_chk.eq.0 ) then call brent(nn, AAX, CCX, BX, FX, nstate, n_iuser, iuser, n_ruser, ruser, istat) else ! corner solution call objfn(nn, BX, FX, nstate, n_iuser, iuser, n_ruser, ruser, inform) end if fret = FX if ( (fptt-fret).gt.del ) then ! record dimension of largest decrease del = fptt-fret ibig = ii end if end do if ( 2.0*(fp-fret).le.ftol*(abs(fret)+abs(fp)+eps) ) then finish = 1 else adjust_factor = 1.0_8 do jj = 1, nn ptt(jj) = BX(jj) + BX(jj) - pt(jj) ! extrapolated point along direction of average shift xit(jj) = BX(jj) - pt(jj) ! average direction of shift if ( ptt(jj).lt.AX(jj) ) then adjust_factor = min( adjust_factor, (AX(jj)-BX(jj))/(BX(jj)-pt(jj)) ) elseif( ptt(jj).gt.CX(jj) ) then adjust_factor = min( adjust_factor, (CX(jj)-BX(jj))/(BX(jj)-pt(jj)) ) end if end do if ( adjust_factor.lt.0.0 ) then !$OMP CRITICAL(error_1) print *, 'problem in calculating adjustment factor' pause 'press enter to exit' stop continue !$OMP END CRITICAL(error_1) elseif ( adjust_factor.ne.1.0 ) then do jj = 1, nn ptt(jj) = BX(jj) + adjust_factor * (BX(jj) - pt(jj)) end do end if call objfn(nn, ptt, FE, nstate, n_iuser, iuser, n_ruser, ruser, inform) fptt = FE if ( fptt+eps_v.lt.fp ) then T_val = 2.0 * (fp-2.0*fret+fptt) * (fp-fret-del)**2.0 - del*(fp-fptt)**2 if ( T_val.lt.0.0 ) then call boundary_pnts(nn, AX, CX, xit, BX, AAX, CCX, bound_chk) if ( bound_chk.eq.0 ) then call brent(nn, AAX, CCX, BX, FX, nstate, n_iuser, iuser, n_ruser, ruser, istat) else ! corner solution call objfn(nn, BX, FX, nstate, n_iuser, iuser, n_ruser, ruser, inform) end if fret = FX do jj = 1, nn xi(jj,ibig) = xi(jj,nn) xi(jj,nn) = xit(jj) end do end if end if end if if ( (iter.eq.itmax).and.(finish.eq.0) ) then print *, 'Failure to converge in surface solution routine without gradients' finish = 1 end if end do yy = BX val = FX END SUBROUTINE powell