module odrpack_core !! Core mathematical routines, except drivers, and BLAS/LINPACK. use odrpack_kinds, only: wp implicit none abstract interface subroutine fcn_t(beta, xplusd, ifixb, ifixx, ideval, f, fjacb, fjacd, istop) !! User-supplied subroutine for evaluating the model. import :: wp implicit none real(wp), intent(in) :: beta(:) !! Current values of parameters. real(wp), intent(in) :: xplusd(:, :) !! Current value of explanatory variable, i.e., `x + delta`. Shape: `(n, m)`. integer, intent(in) :: ifixb(:) !! Indicators for "fixing" parameters (`beta`). integer, intent(in) :: ifixx(:, :) !! Indicators for "fixing" explanatory variable (`x`). Shape: `(ldifx, m)`. integer, intent(in) :: ideval !! Indicator for selecting computation to be performed. real(wp), intent(out) :: f(:, :) !! Predicted function values. Shape: `(n, q)`. real(wp), intent(out) :: fjacb(:, :, :) !! Jacobian with respect to `beta`. Shape: `(n, np, q)`. real(wp), intent(out) :: fjacd(:, :, :) !! Jacobian with respect to errors `delta`. Shape: `(n, m, q)`. integer, intent(out) :: istop !! Stopping condition, with meaning as follows. !! `0`: Current `beta` and `x + delta` were acceptable and values were computed !! successfully. !! `1`: Current `beta` and `x + delta` are not acceptable; 'odrpack' should select !! values closer to most recently used values if possible. !! `-1`: Current `beta` and `x + delta` are not acceptable; 'odrpack' should stop. end subroutine fcn_t end interface contains subroutine trust_region & (n, m, np, q, npp, & f, fjacb, fjacd, & wd, ldwd, ld2wd, ss, tt, ldtt, delta, & alpha2, tau, epsfcn, isodr, & tfjacb, omega, u, qraux, jpvt, & s, t, nlms, rcond, irank, & wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc) !! Compute Levenberg-Marquardt parameter and steps `s` and `t` using analog of the !! trust-region Levenberg-Marquardt algorithm. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: npp !! Number of function parameters being estimated. real(wp), intent(in) :: f(n, q) !! Weighted estimated values of `epsilon`. real(wp), intent(in) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. real(wp), intent(in) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. real(wp), intent(in) :: wd(ldwd, ld2wd, m) !! `delta` weights. integer, intent(in) :: ldwd !! Leading dimension of array `wd`. integer, intent(in) :: ld2wd !! Second dimension of array `wd`. real(wp), intent(in) :: ss(np) !! Scaling values used for the unfixed `beta`s. real(wp), intent(in) :: tt(ldtt, m) !! Scale used for the `delta`s. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. real(wp), intent(in) :: delta(n, m) !! Estimated errors in the explanatory variables. real(wp), intent(inout) :: alpha2 !! Current Levenberg-Marquardt parameter. real(wp), intent(inout) :: tau !! Trust region diameter. real(wp), intent(in) :: epsfcn !! Function's precision. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(out) :: tfjacb(n, q, np) !! Array `omega*fjacb`. real(wp), intent(out) :: omega(q, q) !! Array `(I-fjacd*inv(p)*trans(fjacd))**(-1/2)`. real(wp), intent(out) :: u(np) !! Approximate null vector for `tfjacb`. real(wp), intent(out) :: qraux(np) !! Array required to recover the orthogonal part of the QR decomposition. integer, intent(out) :: jpvt(np) !! Pivot vector. real(wp), intent(out) :: s(np) !! Step for `beta`. real(wp), intent(out) :: t(n, m) !! Step for `delta`. integer, intent(out) :: nlms !! Number of Levenberg-Marquardt steps taken. real(wp), intent(out) :: rcond !! Approximate reciprocal condition of `tfjacb`. integer, intent(out) :: irank !! Aank deficiency of the Jacobian wrt `beta`. real(wp), intent(out) :: wrk1(n, q, m) !! Work array of `(n, q, m)` elements. real(wp), intent(out) :: wrk2(n, q) !! Work array of `(n, q)` elements. real(wp), intent(out) :: wrk3(np) !! Work array of `(np)` elements. real(wp), intent(out) :: wrk4(m, m) !! Work array of `(m, m)` elements. real(wp), intent(out) :: wrk5(m) !! Work array of `(m)` elements. real(wp), intent(out) :: wrk(lwrk) !! Work array of `(lwrk)` elements, _equivalenced_ to `wrk1` and `wrk2`. integer, intent(in) :: lwrk !! Length of vector `wrk`. integer, intent(out) :: istopc !! Variable designating whether the computations were stopped due to some other !! numerical error detected within subroutine `dodstp`. ! Local scalars real(wp), parameter :: p001 = 0.001_wp, p1 = 0.1_wp real(wp) :: alpha1, alphan, bot, phi1, phi2, sa, top integer :: i, iwrk, j, k logical :: forvcv ! Variable Definitions (alphabetically) ! ALPHAN: The new Levenberg-Marquardt parameter. ! ALPHA1: The previous Levenberg-Marquardt parameter. ! BOT: The lower limit for setting ALPHA. ! FORVCV: The variable designating whether this subroutine was called to set up for the ! covariance matrix computations (TRUE) or not (FALSE). ! I: An indexing variable. ! IWRK: An indexing variable. ! J: An indexing variable. ! K: An indexing variable. ! PHI1: The previous difference between the norm of the scaled step and the trust ! region diameter. ! PHI2: The current difference between the norm of the scaled step and the trust region ! diameter. ! SA: The scalar PHI2*(ALPHA1-ALPHA2)/(PHI1-PHI2). ! TOP: The upper limit for setting ALPHA. forvcv = .false. istopc = 0 ! Compute full Gauss-Newton step (ALPHA=0) alpha1 = zero call lcstep(n, m, np, q, npp, & f, fjacb, fjacd, & wd, ldwd, ld2wd, ss, tt, ldtt, delta, & alpha1, epsfcn, isodr, & tfjacb, omega, u, qraux, jpvt, & s, t, phi1, irank, rcond, forvcv, & wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc) if (istopc /= 0) then return end if ! Initialize TAU if necessary if (tau < zero) then tau = abs(tau)*phi1 end if ! Check if full Gauss-Newton step is optimal if ((phi1 - tau) <= p1*tau) then nlms = 1 alpha2 = zero return end if ! Full Gauss-Newton step is outside trust region - find locally constrained optimal step phi1 = phi1 - tau ! Initialize upper and lower bounds for ALPHA bot = zero do k = 1, npp tfjacb(1:n, 1:q, k) = fjacb(1:n, k, 1:q) wrk(k) = sum(tfjacb(:, :, k)*f) end do call scale_vec(npp, 1, ss, npp, wrk, npp, wrk, npp) ! work is input (as t) and output (as sclt) if (isodr) then call scale_mat(n, m, wd, ldwd, ld2wd, delta, wrk(npp + 1:npp + 1 + n*m - 1)) iwrk = npp do j = 1, m do i = 1, n iwrk = iwrk + 1 wrk(iwrk) = wrk(iwrk) + dot_product(fjacd(i, j, :), f(i, :)) end do end do call scale_vec(n, m, tt, ldtt, wrk(npp + 1), n, wrk(npp + 1), n) top = norm2(wrk(1:npp + n*m))/tau else top = norm2(wrk(1:npp))/tau end if if (alpha2 > top .or. alpha2 == zero) then alpha2 = p001*top end if ! Main loop do i = 1, 10 ! Compute locally constrained steps S and T and PHI(ALPHA) for current value of ALPHA call lcstep(n, m, np, q, npp, & f, fjacb, fjacd, & wd, ldwd, ld2wd, ss, tt, ldtt, delta, & alpha2, epsfcn, isodr, & tfjacb, omega, u, qraux, jpvt, & s, t, phi2, irank, rcond, forvcv, & wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc) if (istopc /= 0) then return end if phi2 = phi2 - tau ! Check whether current step is optimal if (abs(phi2) <= p1*tau .or. (alpha2 == bot .and. phi2 < zero)) then nlms = i + 1 return end if ! Current step is not optimaL ! Update bounds for ALPHA and compute new ALPHA if (phi1 - phi2 == zero) then nlms = 12 return end if sa = phi2*(alpha1 - alpha2)/(phi1 - phi2) if (phi2 < zero) then top = min(top, alpha2) else bot = max(bot, alpha2) end if if (phi1*phi2 > zero) then bot = max(bot, alpha2 - sa) else top = min(top, alpha2 - sa) end if alphan = alpha2 - sa*(phi1 + tau)/tau if (alphan >= top .or. alphan <= bot) then alphan = max(p001*top, sqrt(top*bot)) end if ! Get ready for next iteration alpha1 = alpha2 alpha2 = alphan phi1 = phi2 end do ! Set NLMS to indicate an optimal step could not be found in 10 trys nlms = 12 end subroutine trust_region pure subroutine access_workspace & (n, m, np, q, ldwe, ld2we, & rwork, lrwork, iwork, liwork, & access, isodr, & jpvt, omega, u, qraux, sd, vcv, wrk1, wrk2, wrk3, wrk4, wrk5, wrk6, & nnzw, npp, & job, partol, sstol, maxit, taufac, eta, neta, & lunrpt, ipr1, ipr2, ipr2f, ipr3, & wss, rvar, idf, & tau, alpha, niter, nfev, njev, int2, olmavg, & rcond, irank, actrs, pnorm, prers, rnorms, istop) !! Access or store values in the work arrays. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: ldwe !! Leading dimension of array `we`. integer, intent(in) :: ld2we !! Second dimension of array `we`. real(wp), intent(inout) :: rwork(lrwork) !! Real work space. integer, intent(in) :: lrwork !! Length of vector `rwork`. integer, intent(inout) :: iwork(liwork) !! Integer work space. integer, intent(in) :: liwork !! Length of vector `iwork`. logical, intent(in) :: access !! Variable designating whether information is to be accessed from the work !! arrays (`.true.`) or stored in them (`.false.`). logical, intent(in) :: isodr !! Variable designating whether the solution is to be found by ODR (`.true.`) !! or by OLS (`.false.`). integer, intent(out) :: jpvt !! Pivot vector. integer, intent(out) :: omega !! Starting location in array `rwork` of array `omega(q**2)`. integer, intent(out) :: u !! Starting location in array `rwork` of array `u(np)`. integer, intent(out) :: qraux !! Starting location in array `rwork` of array `qraux(np)`. integer, intent(out) :: sd !! Starting location in array `rwork` of array `sd(np)`. integer, intent(out) :: vcv !! Starting location in array `rwork` of array `vcv(np**2)`. integer, intent(out) :: wrk1 !! Starting location in array `rwork` of array `wrk1(n, m, q)`. integer, intent(out) :: wrk2 !! Starting location in array `rwork` of array `wrk2(n, q)`. integer, intent(out) :: wrk3 !! Starting location in array `rwork` of array `wrk3(np)`. integer, intent(out) :: wrk4 !! Starting location in array `rwork` of array `wrk4(m, m)`. integer, intent(out) :: wrk5 !! Starting location in array `rwork` of array `wrk5(m)`. integer, intent(out) :: wrk6 !! Starting location in array `rwork` of array `wrk6(n, np, q)`. integer, intent(out) :: nnzw !! Number of nonzero weighted observations. integer, intent(out) :: npp !! Number of function parameters actually estimated. integer, intent(out) :: job !! Variable controlling problem initialization and computational method. real(wp), intent(inout) :: partol !! Parameter convergence stopping tolerance. real(wp), intent(inout) :: sstol !! Sum-of-squares convergence stopping tolerance. integer, intent(out) :: maxit !! Maximum number of iterations allowed. real(wp), intent(out) :: taufac !! Factor used to compute the initial trust region diameter. real(wp), intent(out) :: eta !! Relative noise in the function results. integer, intent(out) :: neta !! Number of accurate digits in the function results. integer, intent(out) :: lunrpt !! Logical unit number used for computation reports. integer, intent(out) :: ipr1 !! Value of the fourth digit (from the right) of `iprint`, which controls the !! initial summary report. integer, intent(out) :: ipr2 !! Value of the third digit (from the right) of `iprint`, which controls the !! iteration reports. integer, intent(out) :: ipr2f !! Value of the second digit (from the right) of `iprint`, which controls the !! frequency of the iteration reports. integer, intent(out) :: ipr3 !! Value of the first digit (from the right) of `iprint`, which controls the final !! summary report. real(wp), intent(inout) :: wss(3) !! Sum of the squares of the weighted `epsilons` and `deltas`, the sum of the squares !! of the weighted `deltas`, and the sum of the squares of the weighted `epsilons`. real(wp), intent(inout) :: rvar !! Residual variance, i.e. the standard deviation squared. integer, intent(inout) :: idf !! Degrees of freedom of the fit, equal to the number of observations with nonzero !! weighted derivatives minus the number of parameters being estimated. real(wp), intent(inout) :: tau !! Trust region diameter. real(wp), intent(inout) :: alpha !! Levenberg-Marquardt parameter. integer, intent(inout) :: niter !! Number of iterations taken. integer, intent(inout) :: nfev !! Number of function evaluations. integer, intent(inout) :: njev !! Number of Jacobian evaluations. integer, intent(inout) :: int2 !! Number of internal doubling steps. real(wp), intent(inout) :: olmavg !! Average number of Levenberg-Marquardt steps per iteration. real(wp), intent(inout) :: rcond !! Approximate reciprocal condition of `fjacb`. integer, intent(inout) :: irank !! Rank deficiency of the Jacobian wrt `beta`. real(wp), intent(inout) :: actrs !! Saved actual relative reduction in the sum-of-squares. real(wp), intent(inout) :: pnorm !! Norm of the scaled estimated parameters. real(wp), intent(inout) :: prers !! Saved predicted relative reduction in the sum-of-squares. real(wp), intent(inout) :: rnorms !! Norm of the saved weighted `epsilons` and `deltas`. integer, intent(inout) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. ! Local scalars integer :: actrsi, alphai, betaci, betani, betasi, beta0i, boundi, deltai, deltani, & deltasi, diffi, epsi, epsmaci, etai, fjacbi, fjacdi, fni, fsi, idfi, int2i, & iprinti, iprint, iranki, istopi, jobi, jpvti, ldtti, liwkmin, loweri, lunerri, & lunrpti, lrwkmin, maxiti, msgb, msgd, netai, nfevi, niteri, njevi, nnzwi, nppi, & nrowi, ntoli, olmavgi, omegai, partoli, pnormi, prersi, qrauxi, rcondi, rnormsi, & rvari, sdi, si, ssfi, ssi, sstoli, taufaci, taui, ti, tti, ui, upperi, vcvi, & we1i, wrk1i, wrk2i, wrk3i, wrk4i, wrk5i, wrk6i, wrk7i, wssi, wssdeli, wssepsi, & xplusdi ! Variable Definitions (alphabetically) ! ACTRSI: The location in array RWORK of variable ACTRS. ! ALPHAI: The location in array RWORK of variable ALPHA. ! BETACI: The starting location in array RWORK of array BETAC. ! BETANI: The starting location in array RWORK of array BETAN. ! BETASI: The starting location in array RWORK of array BETAS. ! BETA0I: The starting location in array RWORK of array BETA0. ! DELTAI: The starting location in array RWORK of array DELTA. ! DELTANI: The starting location in array RWORK of array DELTAN. ! DELTASI: The starting location in array RWORK of array DELTAS. ! DIFFI: The starting location in array RWORK of array DIFF. ! EPSI: The starting location in array RWORK of array EPS. ! EPSMACI: The location in array RWORK of variable EPSMAC. ! ETAI: The location in array RWORK of variable ETA. ! FJACBI: The starting location in array RWORK of array FJACB. ! FJACDI: The starting location in array RWORK of array FJACD. ! FNI: The starting location in array RWORK of array FN. ! FSI: The starting location in array RWORK of array FS. ! IDFI: The starting location in array IWORK of variable IDF. ! INT2I: The location in array IWORK of variable INT2. ! IPRINTI: The location in array IWORK of variable IPRINT. ! IPRINT: The print control variable. ! IRANKI: The location in array IWORK of variable IRANK. ! ISTOPI: The location in array IWORK of variable ISTOP. ! JOBI: The location in array IWORK of variable JOB. ! JPVTI: The starting location in array IWORK of variable JPVT. ! LDTTI: The starting location in array IWORK of variable LDTT. ! LUNERRI: The location in array IWORK of variable LUNERR. ! LUNRPTI: The location in array IWORK of variable LUNRPT. ! LRWKMIN: The minimum acceptable length of array RWORK. ! MAXITI: The location in array IWORK of variable MAXIT. ! MSGB: The starting location in array IWORK of array MSGB. ! MSGD: The starting location in array IWORK of array MSGD. ! NETAI: The location in array IWORK of variable NETA. ! NFEVI: The location in array IWORK of variable NFEV. ! NITERI: The location in array IWORK of variable NITER. ! NJEVI: The location in array IWORK of variable NJEV. ! NNZWI: The location in array IWORK of variable NNZW. ! NPPI: The location in array IWORK of variable NPP. ! NROWI: The location in array IWORK of variable NROW. ! NTOLI: The location in array IWORK of variable NTOL. ! OLMAVGI: The location in array RWORK of variable OLMAVG. ! OMEGAI: The starting location in array RWORK of array OMEGA. ! PARTOLI: The location in array work of variable PARTOL. ! PNORMI: The location in array RWORK of variable PNORM. ! PRERSI: The location in array RWORK of variable PRERS. ! QRAUXI: The starting location in array RWORK of array QRAUX. ! RCONDI: The location in array RWORK of variable RCOND. ! RNORMSI: The location in array RWORK of variable RNORMS. ! RVARI: The location in array RWORK of variable RVAR. ! SDI: The starting location in array RWORK of array SD. ! SSFI: The starting location in array RWORK of array SSF. ! SSI: The starting location in array RWORK of array SS. ! SSTOLI: The location in array RWORK of variable SSTOL. ! TAUFACI: The location in array RWORK of variable TAUFAC. ! TAUI: the location in array RWORK of variable TAU. ! TI: The starting location in array RWORK of array T. ! TTI: The starting location in array RWORK of array TT. ! UI: The starting location in array RWORK of array U. ! VCVI: The starting location in array RWORK of array VCV. ! WE1I: The starting location in array RWORK of array WE1. ! WRK1I: The starting location in array RWORK of array WRK1. ! WRK2I: The starting location in array RWORK of array WRK2. ! WRK3I: The starting location in array RWORK of array wrk3. ! WRK4I: The starting location in array RWORK of array wrk4. ! WRK5I: The starting location in array RWORK of array wrk5. ! WRK6I: The starting location in array RWORK of array wrk6. ! WRK7I: The starting location in array RWORK of array wrk7. ! WSSI: The starting location in array RWORK of variable WSS(1). ! WSSDELI: The starting location in array RWORK of variable WSS(2). ! WSSEPSI: The starting location in array RWORK of variable WSS(3). ! XPLUSDI: The starting location in array RWORK of array XPLUSD. ! Find starting locations within integer workspace call loc_iwork(m, q, np, & msgb, msgd, jpvti, istopi, & nnzwi, nppi, idfi, & jobi, iprinti, lunerri, lunrpti, & nrowi, ntoli, netai, & maxiti, niteri, nfevi, njevi, int2i, iranki, ldtti, & boundi, & liwkmin) ! Find starting locations within REAL work space call loc_rwork(n, m, q, np, ldwe, ld2we, isodr, & deltai, epsi, xplusdi, fni, sdi, vcvi, & rvari, wssi, wssdeli, wssepsi, rcondi, etai, & olmavgi, taui, alphai, actrsi, pnormi, rnormsi, prersi, & partoli, sstoli, taufaci, epsmaci, & beta0i, betaci, betasi, betani, si, ssi, ssfi, qrauxi, ui, & fsi, fjacbi, we1i, diffi, & deltasi, deltani, ti, tti, omegai, fjacdi, & wrk1i, wrk2i, wrk3i, wrk4i, wrk5i, wrk6i, wrk7i, & loweri, upperi, & lrwkmin) if (access) then ! Set starting locations for work vectors jpvt = jpvti omega = omegai qraux = qrauxi sd = sdi vcv = vcvi u = ui wrk1 = wrk1i wrk2 = wrk2i wrk3 = wrk3i wrk4 = wrk4i wrk5 = wrk5i wrk6 = wrk6i ! Access values from the work vectors actrs = rwork(actrsi) alpha = rwork(alphai) eta = rwork(etai) olmavg = rwork(olmavgi) partol = rwork(partoli) pnorm = rwork(pnormi) prers = rwork(prersi) rcond = rwork(rcondi) wss(1) = rwork(wssi) wss(2) = rwork(wssdeli) wss(3) = rwork(wssepsi) rvar = rwork(rvari) rnorms = rwork(rnormsi) sstol = rwork(sstoli) tau = rwork(taui) taufac = rwork(taufaci) neta = iwork(netai) irank = iwork(iranki) job = iwork(jobi) lunrpt = iwork(lunrpti) maxit = iwork(maxiti) nfev = iwork(nfevi) niter = iwork(niteri) njev = iwork(njevi) nnzw = iwork(nnzwi) npp = iwork(nppi) idf = iwork(idfi) int2 = iwork(int2i) ! Set up print control variables iprint = iwork(iprinti) ipr1 = mod(iprint, 10000)/1000 ipr2 = mod(iprint, 1000)/100 ipr2f = mod(iprint, 100)/10 ipr3 = mod(iprint, 10) else ! Store values into the work vectors rwork(actrsi) = actrs rwork(alphai) = alpha rwork(olmavgi) = olmavg rwork(partoli) = partol rwork(pnormi) = pnorm rwork(prersi) = prers rwork(rcondi) = rcond rwork(wssi) = wss(1) rwork(wssdeli) = wss(2) rwork(wssepsi) = wss(3) rwork(rvari) = rvar rwork(rnormsi) = rnorms rwork(sstoli) = sstol rwork(taui) = tau iwork(iranki) = irank iwork(istopi) = istop iwork(nfevi) = nfev iwork(niteri) = niter iwork(njevi) = njev iwork(idfi) = idf iwork(int2i) = int2 end if end subroutine access_workspace real(wp) pure function derstep(itype, k, betak, ssf, stpb, neta) result(res) !! Compute step size for center and forward difference calculations. use odrpack_kinds, only: zero, one integer, intent(in) :: itype !! Finite difference method being used. !! `0`: Forward finite differences. !! `1`: Central finite differences. integer, intent(in) :: k !! Index into `beta` where `betak` resides. real(wp), intent(in) :: betak !! `k`-th function parameter. real(wp), intent(in) :: ssf(k) !! Scale used for the `beta`s. real(wp), intent(in) :: stpb(k) !! Relative step used for computing finite difference derivatives with respect !! to `beta`. integer, intent(in) :: neta !! Number of good digits in the function results. ! Local scalars real(wp) :: typj ! Variable definitions (alphabetically) ! TYPJ: The typical size of the J-th unkonwn BETA. if (betak == zero) then if (ssf(1) < zero) then typj = 1/abs(ssf(1)) else typj = 1/ssf(k) end if else typj = abs(betak) end if res = sign(one, betak)*typj*hstep(itype, neta, 1, k, stpb, 1) end function derstep pure subroutine esubi(n, m, wd, ldwd, ld2wd, alpha, tt, ldtt, i, e) !! Compute `e = wd + alpha*tt**2`. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. real(wp), intent(in) :: wd(ldwd, ld2wd, m) !! Squared `delta` weights. integer, intent(in) :: ldwd !! Leading dimension of array `wd`. integer, intent(in) :: ld2wd !! Second dimension of array `wd`. real(wp), intent(in) :: alpha !! Levenberg-Marquardt parameter. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values used for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. integer, intent(in) :: i !! Indexing variable. real(wp), intent(out) :: e(m, m) !! Value of the array `e = wd + alpha*tt**2`. ! Local scalars integer :: j ! Variable Definitions (alphabetically) ! J: An indexing variable. ! N.B. the locations of WD and TT accessed depend on the value of the first element of ! each array and the leading dimensions of the multiply subscripted arrays. if (n == 0 .or. m == 0) return if (wd(1, 1, 1) >= zero) then if (ldwd >= n) then ! The elements of WD have been individually specified if (ld2wd == 1) then ! The arrays stored in WD are diagonal e = zero do j = 1, m e(j, j) = wd(i, 1, j) end do else ! The arrays stored in WD are full positive semidefinite matrices e = wd(i, :, :) end if if (tt(1, 1) > zero) then if (ldtt >= n) then do j = 1, m e(j, j) = e(j, j) + alpha*tt(i, j)**2 end do else do j = 1, m e(j, j) = e(j, j) + alpha*tt(1, j)**2 end do end if else do j = 1, m e(j, j) = e(j, j) + alpha*tt(1, 1)**2 end do end if else ! WD is an M by M matrix if (ld2wd == 1) then ! The array stored in WD is diagonal e = zero do j = 1, m e(j, j) = wd(1, 1, j) end do else ! The array stored in WD is a full positive semidefinite matrix e = wd(1, :, :) end if if (tt(1, 1) > zero) then if (ldtt >= n) then do j = 1, m e(j, j) = e(j, j) + alpha*tt(i, j)**2 end do else do j = 1, m e(j, j) = e(j, j) + alpha*tt(1, j)**2 end do end if else do j = 1, m e(j, j) = e(j, j) + alpha*tt(1, 1)**2 end do end if end if else ! WD is a diagonal matrix with elements ABS(WD(1,1,1)) e = zero if (tt(1, 1) > zero) then if (ldtt >= n) then do j = 1, m e(j, j) = abs(wd(1, 1, 1)) + alpha*tt(i, j)**2 end do else do j = 1, m e(j, j) = abs(wd(1, 1, 1)) + alpha*tt(1, j)**2 end do end if else do j = 1, m e(j, j) = abs(wd(1, 1, 1)) + alpha*tt(1, 1)**2 end do end if end if end subroutine esubi subroutine eta_fcn & (fcn, & n, m, np, q, & xplusd, beta, epsmac, nrow, & partmp, pv0, & ifixb, ifixx, ldifx, & istop, nfev, eta, neta, & fjacd, f, fjacb, wrk7, & info, & lower, upper) !! Compute noise and number of good digits in function results. ! Adapted from STARPAC subroutine ETAFUN. use odrpack_kinds, only: zero, one, two procedure(fcn_t) :: fcn !! User-supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(in) :: xplusd(n, m) !! Values of `x + delta`. real(wp), intent(in) :: beta(np) !! Function parameters. real(wp), intent(in) :: epsmac !! Value of machine precision. integer, intent(in) :: nrow !! Row number at which the derivative is to be checked. real(wp), intent(out) :: partmp(np) !! Model parameters. real(wp), intent(in) :: pv0(n, q) !! Original predicted values. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: eta !! Noise in the model results. integer, intent(out) :: neta !! Number of accurate digits in the model results. real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian wrt `delta`. real(wp), intent(out) :: f(n, q) !! Function values. real(wp), intent(out) :: fjacb(n, np, q) !! Jacobian wrt `beta`. real(wp), intent(out) :: wrk7(-2:2, q) !! Work array of `(5, q)` elements. integer, intent(out) :: info !! Variable indicating the status of the computation. real(wp), intent(in) :: lower(np) !! Lower bound of `beta`. real(wp), intent(in) :: upper(np) !! Upper bound of `beta`. ! Local scalars real(wp), parameter :: p1 = 0.1_wp, p2 = 0.2_wp, p5 = 0.5_wp real(wp) :: a, b, fac, shift, stp integer :: j, k, sbk ! Local arrays real(wp) :: parpts(-2:2, np) ! Variable Definitions (ALPHABETICALLY) ! A: Parameters of the local fit. ! B: Parameters of the local fit. ! FAC: A factor used in the computations. ! J: An index variable. ! K: An index variable. ! SHIFT: When PARPTS cross the parameter bounds they are shifted by SHIFT. ! SBK: The sign of BETA(K). ! STP: A small value used to perturb the parameters. stp = 100*epsmac eta = epsmac ! Create points to use in calculating FCN for ETA and NETA do j = -2, 2 if (j == 0) then parpts(0, :) = beta else do k = 1, np if (ifixb(1) < 0) then parpts(j, k) = beta(k) + j*stp*beta(k) elseif (ifixb(k) /= 0) then parpts(j, k) = beta(k) + j*stp*beta(k) else parpts(j, k) = beta(k) end if end do end if end do ! Adjust the points used in calculating FCN to uphold the boundary constraints do k = 1, np sbk = int(sign(one, parpts(2, k) - parpts(-2, k))) if (parpts(sbk*2, k) > upper(k)) then shift = upper(k) - parpts(sbk*2, k) parpts(sbk*2, k) = upper(k) do j = -sbk*2, sbk*1, sbk parpts(j, k) = parpts(j, k) + shift end do if (parpts(-sbk*2, k) < lower(k)) then info = 90010 return end if end if if (parpts(-sbk*2, k) < lower(k)) then shift = lower(k) - parpts(-sbk*2, k) parpts(-sbk*2, k) = lower(k) do j = -sbk*1, sbk*2, sbk parpts(j, k) = parpts(j, k) + shift end do if (parpts(sbk*2, k) > upper(k)) then info = 90010 return end if end if end do ! Evaluate FCN for all points in PARPTS do j = -2, 2 if (all(parpts(j, :) == beta)) then wrk7(j, :) = pv0(nrow, :) else partmp = parpts(j, :) istop = 0 call fcn(partmp, xplusd, ifixb, ifixx, 003, f, fjacb, fjacd, istop) if (istop /= 0) then return else nfev = nfev + 1 end if wrk7(j, :) = f(nrow, :) end if end do ! Calculate ETA and NETA do k = 1, q a = zero b = zero do j = -2, 2 a = a + wrk7(j, k) b = b + j*wrk7(j, k) end do a = p2*a b = p1*b if ((wrk7(0, k) /= zero) .and. (abs(wrk7(1, k) + wrk7(-1, k)) > 100*epsmac)) then fac = 1/abs(wrk7(0, k)) else fac = one end if do j = -2, 2 wrk7(j, k) = abs((wrk7(j, k) - (a + j*b))*fac) eta = max(wrk7(j, k), eta) end do end do neta = int(max(two, p5 - log10(eta))) end subroutine eta_fcn subroutine eval_jac & (fcn, & anajac, cdjac, & n, m, np, q, & betac, beta, stpb, & ifixb, ifixx, ldifx, & x, delta, xplusd, stpd, ldstpd, & ssf, tt, ldtt, neta, fn, & stp, wrk1, wrk2, wrk3, wrk6, tempret, & fjacb, isodr, fjacd, we1, ldwe, ld2we, & njev, nfev, istop, info, & lower, upper) !! Compute the weighted Jacobians wrt `beta` and `delta`. use odrpack_kinds, only: zero procedure(fcn_t) :: fcn !! User-supplied subroutine for evaluating the model. logical, intent(in) :: anajac !! Variable designating whether the Jacobians are computed by finite differences !! (`.false.`) or not (`.true.`). logical, intent(in) :: cdjac !! Variable designating whether the Jacobians are computed by central differences !! (`.true.`) or by forward differences (`.false.`). integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(in) :: betac(np) !! Current estimated values of the unfixed `beta`s. real(wp), intent(out) :: beta(np) !! Function parameters. real(wp), intent(in) :: stpb(np) !! Relative step used for computing finite difference derivatives with respect to `beta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `delta` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: x(n, m) !! Independent variable. real(wp), intent(in) :: delta(n, m) !! Estimated values of `delta`. real(wp), intent(out) :: xplusd(n, m) !! Values of `x + delta`. real(wp), intent(in) :: stpd(ldstpd, m) !! Relative step used for computing finite difference derivatives with respect to `delta`. integer, intent(in) :: ldstpd !! Leading dimension of array `stpd`. real(wp), intent(in) :: ssf(np) !! Scale used for the `beta`s. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values used for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. integer, intent(in) :: neta !! Number of accurate digits in the function results. real(wp), intent(in) :: fn(n, q) !! Predicted values of the function at the current point. real(wp), intent(out) :: stp(n) !! Step used for computing finite difference derivatives with respect to `delta`. real(wp), intent(out) :: wrk1(n, m, q) !! Work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! Work array of `(n, q)` elements. real(wp), intent(out) :: wrk3(np) !! Work array of `(np)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! Work array of `(n, np, q)` elements. real(wp), intent(inout) :: tempret(:, :) !! Temporary work array for holding return values before copying to a lower rank array. real(wp), intent(out) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. real(wp), intent(in) :: we1(ldwe, ld2we, q) !! Square roots of the `epsilon` weights in array `we`. integer, intent(in) :: ldwe !! Leading dimension of arrays `we` and `we1`. integer, intent(in) :: ld2we !! Second dimension of arrays `we` and `we1`. integer, intent(inout) :: njev !! Number of Jacobian evaluations. integer, intent(inout) :: nfev !! Number of function evaluations. integer, intent(out) :: istop !! Variable designating that the user wishes the computations stopped. integer, intent(out) :: info !! Variable designating why the computations were stopped. real(wp), intent(in) :: lower(np) !! Lower bound of `beta`. real(wp), intent(in) :: upper(np) !! Upper bound of `beta`. ! Local scalars integer :: ideval, j, j1 logical :: ferror ! Variable Definitions (alphabetically) ! FERROR: The variable designating whether ODRPACK95 detected nonzero values in array DELTA ! in the OLS case, and thus whether the user may have overwritten important information ! by computing FJACD in the OLS case. ! IDEVAL: The variable designating what computations are to be performed by user-supplied ! subroutine FCN. ! J: An indexing variable. ! J1: An indexing variable. ! Insert current unfixed BETA estimates into BETA call unpack_vec(np, betac, beta, ifixb) ! Compute XPLUSD = X + DELTA xplusd = x + delta ! Compute the Jacobian wrt the estimated BETAS (FJACB) and the Jacobian wrt DELTA (FJACD) istop = 0 if (isodr) then ideval = 110 else ideval = 010 end if if (anajac) then call fcn(beta, xplusd, ifixb, ifixx, ideval, wrk2, fjacb, fjacd, istop) if (istop /= 0) then return else njev = njev + 1 end if ! Make sure fixed elements of FJACD are zero if (isodr) then do j = 1, q call set_ifix(n, m, ifixx, ldifx, fjacd(1, 1, j), n, fjacd(1, 1, j), n) end do end if elseif (cdjac) then call jac_cdiff(fcn, & n, m, np, q, & beta, x, delta, xplusd, ifixb, ifixx, ldifx, & stpb, stpd, ldstpd, & ssf, tt, ldtt, neta, fn, stp, wrk1, wrk2, wrk3, wrk6, & fjacb, isodr, fjacd, nfev, istop, info, & lower, upper) else call jac_fwdiff(fcn, & n, m, np, q, & beta, x, delta, xplusd, ifixb, ifixx, ldifx, & stpb, stpd, ldstpd, & ssf, tt, ldtt, neta, fn, stp, wrk1, wrk2, wrk3, wrk6, & fjacb, isodr, fjacd, nfev, istop, info, & lower, upper) end if if (istop < 0 .or. info >= 10000) then return elseif (.not. isodr) then ! Try to detect whether the user has computed JFACD within FCN in the OLS case ferror = sum(delta**2) /= zero if (ferror) then info = 50300 return end if end if ! Weight the Jacobian wrt the estimated BETAS if (ifixb(1) < 0) then do j = 1, np call scale_mat(n, q, we1, ldwe, ld2we, & reshape(fjacb(:, j, :), [n, q]), & tempret(1:n, 1:q)) fjacb(:, j, :) = tempret(1:n, 1:q) end do else j1 = 0 do j = 1, np if (ifixb(j) >= 1) then j1 = j1 + 1 call scale_mat(n, q, we1, ldwe, ld2we, & reshape(fjacb(:, j, :), [n, q]), & tempret(1:n, 1:q)) fjacb(:, j1, :) = tempret(1:n, 1:q) end if end do end if ! Weight the Jacobian's wrt DELTA as appropriate if (isodr) then do j = 1, m call scale_mat(n, q, we1, ldwe, ld2we, & reshape(fjacd(:, j, :), [n, q]), & tempret(1:n, 1:q)) fjacd(:, j, :) = tempret(1:n, 1:q) end do end if end subroutine eval_jac pure subroutine fctr(oksemi, a, lda, n, info) !! Factor the positive (semi)definite matrix `a` using a modified Cholesky factorization. ! Adapted from LINPACK subroutine DPOFA. use odrpack_kinds, only: zero logical, intent(in) :: oksemi !! Flag indicating whether the factored array can be positive semidefinite !! (`.true.`) or whether it must be found to be positive definite (`.false.`). real(wp), intent(inout) :: a(lda, n) !! Array to be factored. Upon return, `a` contains the upper triangular matrix !! `r` so that `a = trans(r)*r` where the strict lower triangle is set to zero. !! If `info /= 0`, the factorization is not complete. integer, intent(in) :: lda !! Leading dimension of array `a`. integer, intent(in) :: n !! Number of rows and columns of data in array `a`. integer, intent(out) :: info !! Output flag. !! `0`: Factorization was completed. !! `k`: Error condition. The leading minor of order `k` is not positive (semi)definite. ! Local scalars real(wp) :: xi, s, t integer j, k ! Variable Definitions (alphabetically) ! J: An indexing variable. ! K: An indexing variable. ! XI: A value used to test for non positive semidefiniteness. ! Set relative tolerance for detecting non positive semidefiniteness. xi = -10*epsilon(zero) ! Compute factorization, storing in upper triangular portion of A do j = 1, n info = j s = zero do k = 1, j - 1 if (a(k, k) == zero) then t = zero else t = a(k, j) - dot_product(a(1:k - 1, k), a(1:k - 1, j)) t = t/a(k, k) end if a(k, j) = t s = s + t**2 end do s = a(j, j) - s ! Exit if (a(j, j) < zero .or. s < xi*abs(a(j, j))) then return elseif (.not. oksemi .and. s <= zero) then return elseif (s <= zero) then a(j, j) = zero else a(j, j) = sqrt(s) end if end do info = 0 ! Zero out lower portion of A do j = 1, n - 1 a(j + 1:n, j) = zero end do end subroutine fctr pure subroutine fctrw & (n, m, q, npp, & isodr, & we, ldwe, ld2we, wd, ldwd, ld2wd, & wrk0, wrk4, & we1, nnzw, info) !! Check input parameters, indicating errors found using nonzero values of argument `info` !! as described in the ODRPACK95 reference guide. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: npp !! Number of function parameters being estimated. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true`) or by OLS (`.false`). real(wp), intent(in) :: we(ldwe, ld2we, q) !! Squared `epsilon` weights. integer, intent(in) :: ldwe !! Leading dimension of array `we`. integer, intent(in) :: ld2we !! Second dimension of array `we`. real(wp), intent(in) :: wd(ldwd, ld2wd, m) !! Squared `delta` weights. integer, intent(in) :: ldwd !! Leading dimension of array `wd`. integer, intent(in) :: ld2wd !! Second dimension of array `wd`. real(wp), intent(out) :: wrk0(q, q) !! Work array of `(q, q)` elements. real(wp), intent(out) :: wrk4(m, m) !! Work array of `(m, m)` elements. real(wp), intent(out) :: we1(ldwe, ld2we, q) !! Factored `epsilon` weights, such that `trans(we1)*we1 = we`. integer, intent(out) :: nnzw !! Number of nonzero weighted observations. integer, intent(out) :: info !! Variable designating why the computations were stopped. ! Local scalars integer :: i, finfo, j logical :: notzero, exited ! Variable Definitions (alphabetically) ! I: An indexing variable. ! J: An indexing variable. ! L: An indexing variable. ! NOTZERO: The variable designating whether a given component of the weight array WE ! contains a nonzero element (FALSE) or not (TRUE). ! Check EPSILON weights, and store factorization in WE1 exited = .false. if (we(1, 1, 1) < zero) then ! WE contains a scalar we1(1, 1, 1) = -sqrt(abs(we(1, 1, 1))) nnzw = n else nnzw = 0 if (ldwe == 1) then if (ld2we == 1) then ! WE contains a diagonal matrix do j = 1, q if (we(1, 1, j) > zero) then nnzw = n we1(1, 1, j) = sqrt(we(1, 1, j)) elseif (we(1, 1, j) < zero) then info = 30010 exited = .true. exit end if end do else ! WE contains a full Q by Q semidefinite matrix do j = 1, q wrk0(1:j, j) = we(1, 1:j, j) end do call fctr(.true., wrk0, q, q, finfo) if (finfo /= 0) then info = 30010 exited = .true. else do j = 1, q we1(1, :, j) = wrk0(:, j) if (we1(1, j, j) /= zero) then nnzw = n end if end do end if end if else if (ld2we == 1) then ! WE contains an array of diagonal matrix do i = 1, n notzero = .false. do j = 1, q if (we(i, 1, j) > zero) then notzero = .true. we1(i, 1, j) = sqrt(we(i, 1, j)) elseif (we(i, 1, j) < zero) then info = 30010 exited = .true. exit end if end do if (exited) exit if (notzero) then nnzw = nnzw + 1 end if end do else ! WE contains an array of full Q by Q semidefinite matrices do i = 1, n do j = 1, q wrk0(1:j, j) = we(i, 1:j, j) end do call fctr(.true., wrk0, q, q, finfo) if (finfo /= 0) then info = 30010 exited = .true. exit else notzero = .false. do j = 1, q we1(i, :, j) = wrk0(:, j) if (we1(i, j, j) /= zero) then notzero = .true. end if end do end if if (notzero) then nnzw = nnzw + 1 end if end do end if end if end if ! Check for a sufficient number of nonzero EPSILON weights if (.not. exited) then if (nnzw < npp) info = 30020 end if ! Check DELTA weights if (.not. isodr .or. wd(1, 1, 1) < zero) then ! Problem is not ODR, or WD contains a scalar return else if (ldwd == 1) then if (ld2wd == 1) then ! WD contains a diagonal matrix if (any(wd(1, 1, :) <= zero)) then info = max(30001, info + 1) return end if else ! WD contains a full M by M positive definite matrix do j = 1, m wrk4(1:j, j) = wd(1, 1:j, j) end do call fctr(.false., wrk4, m, m, finfo) if (finfo /= 0) then info = max(30001, info + 1) return end if end if else if (ld2wd == 1) then ! WD contains an array of diagonal matrices if (any(wd(:, 1, :) <= zero)) then info = max(30001, info + 1) return end if else ! WD contains an array of full M by M positive definite matrices do i = 1, n do j = 1, m wrk4(1:j, j) = wd(i, 1:j, j) end do call fctr(.false., wrk4, m, m, finfo) if (finfo /= 0) then info = max(30001, info + 1) return end if end do end if end if end if end subroutine fctrw pure subroutine set_flags( & job, restart, initd, dovcv, redoj, anajac, cdjac, chkjac, isodr, implicit) !! Set flags indicating conditions specified by `job`. integer, intent(in) :: job !! Variable controlling problem initialization and computational method. logical, intent(out) :: restart !! Variable designating whether the call is a restart (`.true.`) or not (`.false.`). logical, intent(out) :: initd !! Variable designating whether `delta` is to be initialized to zero (`.true.`) !! or to the first `n` by `m` elements of array `rwork` (`.false.`). logical, intent(out) :: dovcv !! Variable designating whether the covariance matrix is to be computed (`.true.`) !! or not (`.false.`). logical, intent(out) :: redoj !! Variable designating whether the Jacobian matrix is to be recomputed for the !! computation of the covariance matrix (`.true.`) or not (`.false.`). logical, intent(out) :: anajac !! Variable designating whether the Jacobians are computed by finite differences !! (`.false.`) or not (`.true.`). logical, intent(out) :: cdjac !! Variable designating whether the Jacobians are computed by central differences !! (`.true.`) or by forward differences (`.false.`). logical, intent(out) :: chkjac !! Variable designating whether the user-supplied Jacobians are to be checked !! (`.true.`) or not (`.false.`). logical, intent(out) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). logical, intent(out) :: implicit !! Variable designating whether the solution is by implicit ODR (`.true.`) !! or explicit ODR (`.false.`). ! Local scalars integer :: j ! Variable Definitions (alphabetically) ! J: The value of a specific digit of JOB. if (job >= 0) then restart = job >= 10000 initd = mod(job, 10000)/1000 == 0 j = mod(job, 1000)/100 if (j == 0) then dovcv = .true. redoj = .true. elseif (j == 1) then dovcv = .true. redoj = .false. else dovcv = .false. redoj = .false. end if j = mod(job, 100)/10 if (j == 0) then anajac = .false. cdjac = .false. chkjac = .false. elseif (j == 1) then anajac = .false. cdjac = .true. chkjac = .false. elseif (j == 2) then anajac = .true. cdjac = .false. chkjac = .true. else anajac = .true. cdjac = .false. chkjac = .false. end if j = mod(job, 10) if (j == 0) then isodr = .true. implicit = .false. elseif (j == 1) then isodr = .true. implicit = .true. else isodr = .false. implicit = .false. end if else restart = .false. initd = .true. dovcv = .true. redoj = .true. anajac = .false. cdjac = .false. chkjac = .false. isodr = .true. implicit = .false. end if end subroutine set_flags real(wp) pure function hstep(itype, neta, i, j, stp, ldstp) result(res) !! Set relative step size for finite difference derivatives. use odrpack_kinds, only: zero, two, three, ten integer, intent(in) :: itype !! Finite difference method being used. !! `0`: Forward finite differences. !! `1`: Central finite differences. integer, intent(in) :: neta !! Number of good digits in the function results. integer, intent(in) :: i !! Identifier for selecting user-supplied step sizes. integer, intent(in) :: j !! Identifier for selecting user-supplied step sizes. real(wp), intent(in) :: stp(ldstp, j) !! Step size for the finite difference derivative. integer, intent(in) :: ldstp !! Leading dimension of array `stp`. if (stp(1, 1) <= zero) then if (itype == 0) then ! Use default forward finite difference step size res = ten**(-abs(neta)/two - two) else ! Use default central finite difference step size res = ten**(-abs(neta)/three) end if elseif (ldstp == 1) then res = stp(1, j) else res = stp(i, j) end if end function hstep pure subroutine set_ifix(n, m, ifix, ldifix, t, ldt, tfix, ldtfix) !! Set elements of `t` to zero according to `ifix`. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of rows of data in the array. integer, intent(in) :: m !! Number of columns of data in the array. integer, intent(in) :: ifix(ldifix, m) !! Array designating whether an element of `t` is to be set to zero. integer, intent(in) :: ldifix !! Leading dimension of array `ifix`. real(wp), intent(in) :: t(ldt, m) !! Array being set to zero according to the elements of `ifix`. integer, intent(in) :: ldt !! Leading dimension of array `t`. real(wp), intent(out) :: tfix(ldtfix, m) !! Resulting array. integer, intent(in) :: ldtfix !! Leading dimension of array `tfix`. ! Local scalars integer :: j ! Variable Definitions (alphabetically) ! J: an indexing variable. if (n == 0 .or. m == 0) return if (ifix(1, 1) >= zero) then if (ldifix == n) then where (ifix(:, :) == 0) tfix(1:n, :) = zero else where tfix(1:n, :) = t(1:n, :) end where else do j = 1, m if (ifix(1, j) == 0) then tfix(1:n, j) = zero else tfix(1:n, j) = t(1:n, j) end if end do end if end if end subroutine set_ifix pure subroutine loc_iwork & (m, q, np, & msgbi, msgdi, ifix2i, istopi, & nnzwi, nppi, idfi, & jobi, iprinti, lunerri, lunrpti, & nrowi, ntoli, netai, & maxiti, niteri, nfevi, njevi, int2i, iranki, ldtti, & boundi, & liwkmin) !! Get storage locations within integer work space. integer, intent(in) :: m !! Number of columns of data in the independent variable. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: np !! Number of function parameters. integer, intent(out) :: msgbi !! Starting location in array `iwork` of array `msgb`. integer, intent(out) :: msgdi !! Starting location in array `iwork` of array `msgd`. integer, intent(out) :: ifix2i !! Starting location in array `iwork` of array `ifix2`. integer, intent(out) :: istopi !! Location in array `iwork` of variable `istop`. integer, intent(out) :: nnzwi !! Location in array `iwork` of variable `nnzw`. integer, intent(out) :: nppi !! Location in array `iwork` of variable `npp`. integer, intent(out) :: idfi !! Location in array `iwork` of variable `idf`. integer, intent(out) :: jobi !! Location in array `iwork` of variable `job`. integer, intent(out) :: iprinti !! Location in array `iwork` of variable `iprint`. integer, intent(out) :: lunerri !! Location in array `iwork` of variable `lunerr`. integer, intent(out) :: lunrpti !! Location in array `iwork` of variable `lunrpt`. integer, intent(out) :: nrowi !! Location in array `iwork` of variable `nrow`. integer, intent(out) :: ntoli !! Location in array `iwork` of variable `ntol`. integer, intent(out) :: netai !! Location in array `iwork` of variable `neta`. integer, intent(out) :: maxiti !! Location in array `iwork` of variable `maxit`. integer, intent(out) :: niteri !! Location in array `iwork` of variable `niter`. integer, intent(out) :: nfevi !! Location in array `iwork` of variable `nfev`. integer, intent(out) :: njevi !! Location in array `iwork` of variable `njev`. integer, intent(out) :: int2i !! Location in array `iwork` of variable `int2`. integer, intent(out) :: iranki !! Location in array `iwork` of variable `irank`. integer, intent(out) :: ldtti !! Location in array `iwork` of variable `ldtt`. integer, intent(out) :: boundi !! Location in array `iwork` of variable `bound`. integer, intent(out) :: liwkmin !! Minimum acceptable length of array `iwork`. if (np >= 1 .and. m >= 1) then msgbi = 1 msgdi = msgbi + q*np + 1 ifix2i = msgdi + q*m + 1 istopi = ifix2i + np nnzwi = istopi + 1 nppi = nnzwi + 1 idfi = nppi + 1 jobi = idfi + 1 iprinti = jobi + 1 lunerri = iprinti + 1 lunrpti = lunerri + 1 nrowi = lunrpti + 1 ntoli = nrowi + 1 netai = ntoli + 1 maxiti = netai + 1 niteri = maxiti + 1 nfevi = niteri + 1 njevi = nfevi + 1 int2i = njevi + 1 iranki = int2i + 1 ldtti = iranki + 1 boundi = ldtti + 1 liwkmin = boundi + np - 1 else msgbi = 1 msgdi = 1 ifix2i = 1 istopi = 1 nnzwi = 1 nppi = 1 idfi = 1 jobi = 1 iprinti = 1 lunerri = 1 lunrpti = 1 nrowi = 1 ntoli = 1 netai = 1 maxiti = 1 niteri = 1 nfevi = 1 njevi = 1 int2i = 1 iranki = 1 ldtti = 1 boundi = 1 liwkmin = 1 end if end subroutine loc_iwork pure subroutine loc_rwork & (n, m, q, np, ldwe, ld2we, isodr, & deltai, epsi, xplusdi, fni, sdi, vcvi, & rvari, wssi, wssdeli, wssepsi, rcondi, etai, & olmavgi, taui, alphai, actrsi, pnormi, rnormsi, prersi, & partoli, sstoli, taufaci, epsmaci, & beta0i, betaci, betasi, betani, si, ssi, ssfi, qrauxi, ui, & fsi, fjacbi, we1i, diffi, & deltasi, deltani, ti, tti, omegai, fjacdi, & wrk1i, wrk2i, wrk3i, wrk4i, wrk5i, wrk6i, wrk7i, & loweri, upperi, & lrwkmin) !! Get storage locations within real work space. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: ldwe !! Leading dimension of array `we`. integer, intent(in) :: ld2we !! Second dimension of array `we`. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). integer, intent(out) :: deltai !! Starting location in array `rwork` of array `delta(n, m)`. integer, intent(out) :: epsi !! Starting location in array `rwork` of array `eps(n, q)`. integer, intent(out) :: xplusdi !! Starting location in array `rwork` of array `xplusd(n, m)`. integer, intent(out) :: fni !! Starting location in array `rwork` of array `fn(n, q)`. integer, intent(out) :: sdi !! Starting location in array `rwork` of array `sd(np)`. integer, intent(out) :: vcvi !! Starting location in array `rwork` of array `vcv(np**2)`. integer, intent(out) :: rvari !! Location in array `rwork` of variable `rvar`. integer, intent(out) :: wssi !! Location in array `rwork` of variable `wss`. integer, intent(out) :: wssdeli !! Location in array `rwork` of variable `wssdel`. integer, intent(out) :: wssepsi !! Location in array `rwork` of variable `wsseps`. integer, intent(out) :: rcondi !! Location in array `rwork` of variable `rcondi`. integer, intent(out) :: etai !! Location in array `rwork` of variable `eta`. integer, intent(out) :: olmavgi !! Location in array `rwork` of variable `olmavg`. integer, intent(out) :: taui !! Location in array `rwork` of variable `tau`. integer, intent(out) :: alphai !! Location in array `rwork` of variable `alpha`. integer, intent(out) :: actrsi !! Location in array `rwork` of variable `actrs`. integer, intent(out) :: pnormi !! Location in array `rwork` of variable `pnorm`. integer, intent(out) :: rnormsi !! Location in array `rwork` of variable `rnorms`. integer, intent(out) :: prersi !! Location in array `rwork` of variable `prers`. integer, intent(out) :: partoli !! Location in array `rwork` of variable `partol`. integer, intent(out) :: sstoli !! Location in array `rwork` of variable `sstol`. integer, intent(out) :: taufaci !! Location in array `rwork` of variable `taufac`. integer, intent(out) :: epsmaci !! Location in array `rwork` of variable `epsmac`. integer, intent(out) :: beta0i !! Starting location in array `rwork` of array `beta0(np)`. integer, intent(out) :: betaci !! Starting location in array `rwork` of array `betac(np)`. integer, intent(out) :: betasi !! Starting location in array `rwork` of array `betas(np)`. integer, intent(out) :: betani !! Starting location in array `rwork` of array `betan(np)`. integer, intent(out) :: si !! Starting location in array `rwork` of array `s(np)`. integer, intent(out) :: ssi !! Starting location in array `rwork` of array `ss(np)`. integer, intent(out) :: ssfi !! Starting location in array `rwork` of array `ssf(np)`. integer, intent(out) :: qrauxi !! Starting location in array `rwork` of array `qraux(np)`. integer, intent(out) :: ui !! Starting location in array `rwork` of array `u(np)`. integer, intent(out) :: fsi !! Starting location in array `rwork` of array `fs(n, q)`. integer, intent(out) :: fjacbi !! Starting location in array `rwork` of array `fjacb(n, np, q)`. integer, intent(out) :: we1i !! Starting location in array `rwork` of array `we1(ldwe, ldwe2, q)`. integer, intent(out) :: diffi !! Starting location in array `rwork` of array `diff(q*(np + m))`. integer, intent(out) :: deltasi !! Starting location in array `rwork` of array `deltas(n, m)`. integer, intent(out) :: deltani !! Starting location in array `rwork` of array `deltan(n, m)`. integer, intent(out) :: ti !! Starting location in array `rwork` of array `t(n, m)`. integer, intent(out) :: tti !! Starting location in array `rwork` of array `tt(n, m)`. integer, intent(out) :: omegai !! Starting location in array `rwork` of array `omega(q**2)`. integer, intent(out) :: fjacdi !! Starting location in array `rwork` of array `fjacd(n, m, q)`. integer, intent(out) :: wrk1i !! Starting location in array `rwork` of array `wrk1(n, m, q)`. integer, intent(out) :: wrk2i !! Starting location in array `rwork` of array `wrk2(n, q)`. integer, intent(out) :: wrk3i !! Starting location in array `rwork` of array `wrk3(np)`. integer, intent(out) :: wrk4i !! Starting location in array `rwork` of array `wrk4(m, m)`. integer, intent(out) :: wrk5i !! Starting location in array `rwork` of array `wrk5(m)`. integer, intent(out) :: wrk6i !! Starting location in array `rwork` of array `wrk6(n, np, q)`. integer, intent(out) :: wrk7i !! Starting location in array `rwork` of array `wrk7(5*q)`. integer, intent(out) :: loweri !! Starting location in array `rwork` of array `lower(np)`. integer, intent(out) :: upperi !! Starting location in array `rwork` of array `upper(np)`. integer, intent(out) :: lrwkmin !! Minimum acceptable length of vector `rwork`. ! Local scalars integer :: next ! Variable Definitions (alphabetically) ! NEXT: The next available location with RWORK. if (n >= 1 .and. m >= 1 .and. np >= 1 .and. q >= 1 .and. ldwe >= 1 & .and. ld2we >= 1) then deltai = 1 epsi = deltai + n*m xplusdi = epsi + n*q fni = xplusdi + n*m sdi = fni + n*q vcvi = sdi + np rvari = vcvi + np*np wssi = rvari + 1 wssdeli = wssi + 1 wssepsi = wssdeli + 1 rcondi = wssepsi + 1 etai = rcondi + 1 olmavgi = etai + 1 taui = olmavgi + 1 alphai = taui + 1 actrsi = alphai + 1 pnormi = actrsi + 1 rnormsi = pnormi + 1 prersi = rnormsi + 1 partoli = prersi + 1 sstoli = partoli + 1 taufaci = sstoli + 1 epsmaci = taufaci + 1 beta0i = epsmaci + 1 betaci = beta0i + np betasi = betaci + np betani = betasi + np si = betani + np ssi = si + np ssfi = ssi + np qrauxi = ssfi + np ui = qrauxi + np fsi = ui + np fjacbi = fsi + n*q we1i = fjacbi + n*np*q diffi = we1i + ldwe*ld2we*q next = diffi + q*(np + m) if (isodr) then deltasi = next deltani = deltasi + n*m ti = deltani + n*m tti = ti + n*m omegai = tti + n*m fjacdi = omegai + q**2 wrk1i = fjacdi + n*m*q next = wrk1i + n*m*q else deltasi = deltai deltani = deltai ti = deltai tti = deltai omegai = deltai fjacdi = deltai wrk1i = deltai end if wrk2i = next wrk3i = wrk2i + n*q wrk4i = wrk3i + np wrk5i = wrk4i + m*m wrk6i = wrk5i + m wrk7i = wrk6i + n*q*np loweri = wrk7i + 5*q upperi = loweri + np next = upperi + np lrwkmin = next else deltai = 1 epsi = 1 xplusdi = 1 fni = 1 sdi = 1 vcvi = 1 rvari = 1 wssi = 1 wssdeli = 1 wssepsi = 1 rcondi = 1 etai = 1 olmavgi = 1 taui = 1 alphai = 1 actrsi = 1 pnormi = 1 rnormsi = 1 prersi = 1 partoli = 1 sstoli = 1 taufaci = 1 epsmaci = 1 beta0i = 1 betaci = 1 betasi = 1 betani = 1 si = 1 ssi = 1 ssfi = 1 qrauxi = 1 fsi = 1 ui = 1 fjacbi = 1 we1i = 1 diffi = 1 deltasi = 1 deltani = 1 ti = 1 tti = 1 fjacdi = 1 omegai = 1 wrk1i = 1 wrk2i = 1 wrk3i = 1 wrk4i = 1 wrk5i = 1 wrk6i = 1 wrk7i = 1 loweri = 1 upperi = 1 lrwkmin = 1 end if end subroutine loc_rwork pure subroutine init_work & (n, m, np, rwork, lrwork, iwork, liwork, & x, ifixx, ldifx, scld, ldscld, & beta, sclb, & sstol, partol, maxit, taufac, & job, iprint, lunerr, lunrpt, & lower, upper, & epsmai, sstoli, partli, maxiti, taufci, & jobi, iprini, luneri, lunrpi, & ssfi, tti, ldtti, deltai, & loweri, upperi, boundi) !! Initialize work vectors as necessary. use odrpack_kinds, only: zero, one, two, three integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. integer, intent(in) :: np !! Number of function parameters. real(wp), intent(out) :: rwork(lrwork) !! Real work space. integer, intent(in) :: lrwork !! Length of vector `rwork`. integer, intent(out) :: iwork(liwork) !! Integer work space. integer, intent(in) :: liwork !! Length of vector `iwork`. real(wp), intent(in) :: x(n, m) !! Independent variable. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: scld(ldscld, m) !! Scaling values for `delta`. integer, intent(in) :: ldscld !! Leading dimension of array `scld`. real(wp), intent(in) :: beta(np) !! Function parameters. real(wp), intent(in) :: sclb(np) !! Scaling values for `beta`. real(wp), intent(in) :: sstol !! Sum-of-squares convergence stopping criteria. real(wp), intent(in) :: partol !! Parameter convergence stopping criteria. integer, intent(in) :: maxit !! Maximum number of iterations allowed. real(wp), intent(in) :: taufac !! Factor used to compute the initial trust region diameter. integer, intent(in) :: job !! Variable controlling problem initialization and computational method. integer, intent(in) :: iprint !! Print control variable. integer, intent(in) :: lunerr !! Logical unit number used for error messages. integer, intent(in) :: lunrpt !! Logical unit number used for computation reports. real(wp), intent(in) :: lower(np) !! Lower bounds for the function parameters. real(wp), intent(in) :: upper(np) !! Upper bounds for the function parameters. integer, intent(in) :: epsmai !! Location in array `rwork` of variable `epsmac`. integer, intent(in) :: sstoli !! Location in array `rwork` of variable `sstol`. integer, intent(in) :: partli !! Location in array `rwork` of variable `partol`. integer, intent(in) :: maxiti !! Location in array `iwork` of variable `maxit`. integer, intent(in) :: taufci !! Location in array `rwork` of variable `taufac`. integer, intent(in) :: jobi !! Location in array `iwork` of variable `job`. integer, intent(in) :: iprini !! Location in array `iwork` of variable `iprint`. integer, intent(in) :: luneri !! Location in array `iwork` of variable `lunerr`. integer, intent(in) :: lunrpi !! Location in array `iwork` of variable `lunrpt`. integer, intent(in) :: ssfi !! Starting location in array `rwork` of array `ssf`. integer, intent(in) :: tti !! Starting location in array `rwork` of the array `tt`. integer, intent(in) :: ldtti !! Leading dimension of array `tt`. integer, intent(in) :: deltai !! Starting location in array `rwork` of array `delta`. integer, intent(in) :: loweri !! Starting location in array `iwork` of array `lower`. integer, intent(in) :: upperi !! Starting location in array `iwork` of array `upper`. integer, intent(in) :: boundi !! Location in array `iwork` of variable `bound`. ! Local scalars integer :: i, j, istart logical :: anajac, cdjac, chkjac, dovcv, implicit, initd, isodr, redoj, restart ! Variable Definitions (alphabetically) ! ANAJAC: The variable designating whether the Jacobians are computed by finite differences ! (FALSE) or not (TRUE). ! CDJAC: The variable designating whether the Jacobians are computed by central differences ! (TRUE) or by forward differences (FALSE). ! CHKJAC: The variable designating whether the user-supplied Jacobians are to be checked ! (TRUE) or not (FALSE). ! DOVCV: The variable designating whether the covariance matrix is to be computed (TRUE) ! or not (FALSE). ! I: An indexing variable. ! IMPLICIT: The variable designating whether the solution is by implicit ODR (TRUE) or ! explicit ODR (FALSE). ! INITD: The variable designating whether DELTA is to be initialized to zero (TRUE) or ! to the values in the first N by M elements of array RWORK (FALSE). ! ISODR: The variable designating whether the solution is by ODR (TRUE) or by OLS ! (FALSE). ! J: An indexing variable. ! REDOJ: The variable designating whether the Jacobian matrix is to be recomputed for the ! computation of the covariance matrix (TRUE) or not (FALSE). ! RESTART: The variable designating whether the call is a restart (TRUE) or not ! (FALSE). call set_flags(job, restart, initd, dovcv, redoj, anajac, cdjac, chkjac, isodr, implicit) ! Store value of machine precision in work vector rwork(epsmai) = epsilon(zero) ! Set tolerance for stopping criteria based on the change in the parameters (see also ! subprogram DODCNT) if (partol < zero) then rwork(partli) = rwork(epsmai)**(two/three) else rwork(partli) = min(partol, one) end if ! Set tolerance for stopping criteria based on the change in the sum of squares of the ! weighted observational errors if (sstol < zero) then rwork(sstoli) = sqrt(rwork(epsmai)) else rwork(sstoli) = min(sstol, one) end if ! Set factor for computing trust region diameter at first iteration if (taufac <= zero) then rwork(taufci) = one else rwork(taufci) = min(taufac, one) end if ! Set maximum number of iterations if (maxit < 0) then iwork(maxiti) = 50 else iwork(maxiti) = maxit end if ! Store problem initialization and computational method control variable if (job <= 0) then iwork(jobi) = 0 else iwork(jobi) = job end if ! Set print control if (iprint < 0) then iwork(iprini) = 2001 else iwork(iprini) = iprint end if ! Set logical unit number for error messages iwork(luneri) = lunerr ! Set logical unit number for computation reports iwork(lunrpi) = lunrpt ! Compute scaling for BETA's and DELTA's if (sclb(1) <= zero) then call scale_beta(np, beta, rwork(ssfi)) else rwork(ssfi:ssfi + (np - 1)) = sclb end if if (isodr) then if (scld(1, 1) <= zero) then iwork(ldtti) = n call scale_delta(n, m, x, rwork(tti), iwork(ldtti)) else if (ldscld == 1) then iwork(ldtti) = 1 rwork(tti:tti + (m - 1)) = scld(1:m, 1) else iwork(ldtti) = n do j = 1, m istart = tti + (j - 1)*iwork(ldtti) rwork(istart:istart + (n - 1)) = scld(1:n, j) end do end if end if end if ! Initialize DELTA's as necessary if (isodr) then if (initd) then rwork(deltai:deltai + (n*m - 1)) = zero else if (ifixx(1, 1) >= 0) then if (ldifx == 1) then do j = 1, m if (ifixx(1, j) == 0) then istart = deltai + (j - 1)*n rwork(istart:istart + (n - 1)) = zero end if end do else do j = 1, m do i = 1, n if (ifixx(i, j) == 0) then rwork(deltai - 1 + i + (j - 1)*n) = zero end if end do end do end if end if end if else rwork(deltai:deltai + (n*m - 1)) = zero end if ! Copy bounds into RWORK rwork(loweri:loweri + np - 1) = lower rwork(upperi:upperi + np - 1) = upper ! Initialize parameters on bounds in IWORK iwork(boundi:boundi + np - 1) = 0 end subroutine init_work subroutine jac_cdiff & (fcn, & n, m, np, q, & beta, x, delta, xplusd, ifixb, ifixx, ldifx, & stpb, stpd, ldstpd, & ssf, tt, ldtt, neta, fn, stp, wrk1, wrk2, wrk3, wrk6, & fjacb, isodr, fjacd, nfev, istop, info, & lower, upper) !! Compute central difference approximations to the Jacobian wrt the estimated `beta`s and !! wrt the `delta`s. use odrpack_kinds, only: zero, one procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(in) :: x(n, m) !! Explanatory variable. real(wp), intent(in) :: delta(n, m) !! Estimated errors in the explanatory variables. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: stpb(np) !! Relative step used for computing finite difference derivatives with respect to each `beta`. real(wp), intent(in) :: stpd(ldstpd, m) !! Relative step used for computing finite difference derivatives with respect to each `delta`. integer, intent(in) :: ldstpd !! Leading dimension of array `stpd`. real(wp), intent(in) :: ssf(np) !! Scaling values used for `beta`. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values used for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. integer, intent(in) :: neta !! Number of good digits in the function results. real(wp), intent(in) :: fn(n, q) !! New predicted values from the function. Used when parameter is on a boundary. real(wp), intent(out) :: stp(n) !! Step used for computing finite difference derivatives with respect to each `delta`. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk3(np) !! A work array of `(np)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. real(wp), intent(out) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(out) :: info !! Variable designating why the computations were stopped. real(wp), intent(in) :: lower(np) !! Lower bound on `beta`. real(wp), intent(in) :: upper(np) !! Upper bound on `beta`. ! Local scalars real(wp) :: betak, typj integer :: i, j, k, l logical :: doit, setzero ! Variable Definitions (alphabetically) ! BETAK: The K-th function parameter. ! DOIT: The variable designating whether the derivative wrt a given BETA or DELTA ! needs to be computed (TRUE) or not (FALSE). ! I: An indexing variable. ! J: An indexing variable. ! K: An indexing variable. ! L: An indexing variable. ! SETZERO: The variable designating whether the derivative wrt some DELTA needs to be ! set to zero (TRUE) or not (FALSE). ! TYPJ: The typical size of the J-th unknown BETA or DELTA. ! Compute the Jacobian wrt the estimated BETAS do k = 1, np if (ifixb(1) >= 0) then if (ifixb(k) == 0) then doit = .false. else doit = .true. end if else doit = .true. end if if (.not. doit) then fjacb(:, k, :) = zero else betak = beta(k) wrk3(k) = betak + derstep(1, k, betak, ssf, stpb, neta) wrk3(k) = wrk3(k) - betak beta(k) = betak + wrk3(k) if (beta(k) > upper(k)) then beta(k) = upper(k) elseif (beta(k) < lower(k)) then beta(k) = lower(k) end if if (beta(k) - 2*wrk3(k) < lower(k)) then beta(k) = lower(k) + 2*wrk3(k) elseif (beta(k) - 2*wrk3(k) > upper(k)) then beta(k) = upper(k) + 2*wrk3(k) end if if (beta(k) > upper(k) .or. beta(k) < lower(k)) then info = 60001 return end if istop = 0 if (beta(k) == betak) then wrk2 = fn else call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop) if (istop /= 0) then return else nfev = nfev + 1 end if end if fjacb(:, k, :) = wrk2 beta(k) = beta(k) - 2*wrk3(k) if (beta(k) > upper(k)) then info = 60001 return end if if (beta(k) < lower(k)) then info = 60001 return end if istop = 0 if (beta(k) == betak) then wrk2 = fn else call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop) if (istop /= 0) then return else nfev = nfev + 1 end if end if fjacb(:, k, :) = (fjacb(:, k, :) - wrk2)/(2*wrk3(k)) beta(k) = betak end if end do ! Compute the Jacobian wrt the X's if (isodr) then do j = 1, m if (ifixx(1, 1) < 0) then doit = .true. setzero = .false. elseif (ldifx == 1) then if (ifixx(1, j) == 0) then doit = .false. else doit = .true. end if setzero = .false. else doit = any(ifixx(:, j) /= 0) setzero = any(ifixx(:, j) == 0) end if if (.not. doit) then fjacd(:, j, :) = zero else do i = 1, n if (xplusd(i, j) == zero) then if (tt(1, 1) < zero) then typj = 1/abs(tt(1, 1)) elseif (ldtt == 1) then typj = 1/tt(1, j) else typj = 1/tt(i, j) end if else typj = abs(xplusd(i, j)) end if stp(i) = xplusd(i, j) & + sign(one, xplusd(i, j))*typj*hstep(1, neta, i, j, stpd, ldstpd) stp(i) = stp(i) - xplusd(i, j) xplusd(i, j) = xplusd(i, j) + stp(i) end do istop = 0 call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop) if (istop /= 0) then return else nfev = nfev + 1 fjacd(:, j, :) = wrk2 end if xplusd(:, j) = x(:, j) + delta(:, j) - stp istop = 0 call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop) if (istop /= 0) then return else nfev = nfev + 1 end if if (setzero) then do i = 1, n if (ifixx(i, j) == 0) then fjacd(i, j, :) = zero else fjacd(i, j, :) = (fjacd(i, j, :) - wrk2(i, :))/(2*stp(i)) end if end do else do l = 1, q fjacd(:, j, l) = (fjacd(:, j, l) - wrk2(:, l))/(2*stp(:)) end do end if xplusd(:, j) = x(:, j) + delta(:, j) end if end do end if end subroutine jac_cdiff subroutine jac_fwdiff & (fcn, & n, m, np, q, & beta, x, delta, xplusd, ifixb, ifixx, ldifx, & stpb, stpd, ldstpd, & ssf, tt, ldtt, neta, fn, stp, wrk1, wrk2, wrk3, wrk6, & fjacb, isodr, fjacd, nfev, istop, info, & lower, upper) !! Compute forward difference approximations to the Jacobian wrt the estimated `beta`s and !! wrt the `delta`s. use odrpack_kinds, only: zero, one procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(in) :: x(n, m) !! Explanatory variable. real(wp), intent(in) :: delta(n, m) !! Estimated errors in the explanatory variables. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: stpb(np) !! Relative step used for computing finite difference derivatives with respect to each `beta`. real(wp), intent(in) :: stpd(ldstpd, m) !! Relative step used for computing finite difference derivatives with respect to each `delta`. integer, intent(in) :: ldstpd !! Leading dimension of array `stpd`. real(wp), intent(in) :: ssf(np) !! Scaling values used for `beta`. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values used for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. integer, intent(in) :: neta !! Number of good digits in the function results. real(wp), intent(in) :: fn(n, q) !! New predicted values from the function. Used when parameter is on a boundary. real(wp), intent(out) :: stp(n) !! Step used for computing finite difference derivatives with respect to each `delta`. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk3(np) !! A work array of `(np)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. real(wp), intent(out) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(out) :: info !! Variable designating why the computations were stopped. real(wp), intent(in) :: lower(np) !! Lower bound on `beta`. real(wp), intent(in) :: upper(np) !! Upper bound on `beta`. ! Local scalars real(wp) :: betak, step, typj integer :: i, j, k, l logical :: doit, setzero ! Variable Definitions (alphabetically) ! BETAK: The K-th function parameter. ! DOIT: The variable designating whether the derivative wrt a given BETA or DELTA needs ! to be computed (TRUE) or not (FALSE). ! I: An indexing variable. ! J: An indexing variable. ! K: An indexing variable. ! L: An indexing variable. ! SETZRO: The variable designating whether the derivative wrt some DELTA needs to be set ! to zero (TRUE) or not (FALSE). ! TYPJ: The typical size of the J-th unknown BETA or DELTA. ! Compute the Jacobian wrt the estimated BETAS do k = 1, np if (ifixb(1) >= 0) then if (ifixb(k) == 0) then doit = .false. else doit = .true. end if else doit = .true. end if if (.not. doit) then fjacb(:, k, :) = zero else betak = beta(k) step = derstep(0, k, betak, ssf, stpb, neta) wrk3(k) = betak + step wrk3(k) = wrk3(k) - betak beta(k) = betak + wrk3(k) if (beta(k) > upper(k)) then step = -step wrk3(k) = betak + step wrk3(k) = wrk3(k) - betak beta(k) = betak + wrk3(k) end if if (beta(k) < lower(k)) then step = -step wrk3(k) = betak + step wrk3(k) = wrk3(k) - betak beta(k) = betak + wrk3(k) if (beta(k) > upper(k)) then info = 60001 return end if end if istop = 0 call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop) if (istop /= 0) then return else nfev = nfev + 1 end if fjacb(:, k, :) = (wrk2 - fn)/wrk3(k) beta(k) = betak end if end do ! Compute the Jacobian wrt the X'S if (isodr) then do j = 1, m if (ifixx(1, 1) < 0) then doit = .true. setzero = .false. elseif (ldifx == 1) then if (ifixx(1, j) == 0) then doit = .false. else doit = .true. end if setzero = .false. else doit = any(ifixx(:, j) /= 0) setzero = any(ifixx(:, j) == 0) end if if (.not. doit) then fjacd(:, j, :) = zero else do i = 1, n if (xplusd(i, j) == zero) then if (tt(1, 1) < zero) then typj = 1/abs(tt(1, 1)) elseif (ldtt == 1) then typj = 1/tt(1, j) else typj = 1/tt(i, j) end if else typj = abs(xplusd(i, j)) end if stp(i) = xplusd(i, j) & + sign(one, xplusd(i, j))*typj*hstep(0, neta, i, j, stpd, ldstpd) stp(i) = stp(i) - xplusd(i, j) xplusd(i, j) = xplusd(i, j) + stp(i) end do istop = 0 call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop) if (istop /= 0) then return else nfev = nfev + 1 fjacd(:, j, :) = wrk2 end if if (setzero) then do i = 1, n if (ifixx(i, j) == 0) then fjacd(i, j, :) = zero else fjacd(i, j, :) = (fjacd(i, j, :) - fn(i, :))/stp(i) end if end do else do l = 1, q fjacd(:, j, l) = (fjacd(:, j, l) - fn(:, l))/stp(:) end do end if xplusd(:, j) = x(:, j) + delta(:, j) end if end do end if end subroutine jac_fwdiff subroutine check_jac & (fcn, & n, m, np, q, & beta, betaj, xplusd, & ifixb, ifixx, ldifx, stpb, stpd, ldstpd, & ssf, tt, ldtt, & eta, neta, ntol, nrow, isodr, epsmac, & pv0i, fjacb, fjacd, & msgb, msgd, diff, istop, nfev, njev, & wrk1, wrk2, wrk6, & interval) !! Driver routine for the derivative checking process. ! Adapted from STARPAC subroutine DCKCNT. use odrpack_kinds, only: zero, one, p5 => half procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(inout) :: betaj(np) !! Function parameters offset such that steps don't cross bounds. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: stpb(np) !! Step size for finite difference derivatives wrt `beta`. real(wp), intent(in) :: stpd(ldstpd, m) !! Step size for finite difference derivatives wrt `delta`. integer, intent(in) :: ldstpd !! Leading dimension of array `stpd`. real(wp), intent(in) :: ssf(np) !! Scaling values used for `beta`. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values used for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. real(wp), intent(in) :: eta !! Relative noise in the function results. integer, intent(in) :: neta !! Number of reliable digits in the model results. integer, intent(out) :: ntol !! Number of digits of agreement required between the numerical derivatives and the !! user supplied derivatives. integer, intent(in) :: nrow !! Row number of the explanatory variable array at which the derivative is checked. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(in) :: epsmac !! Value of machine precision. real(wp), intent(in) :: pv0i(n, q) !! Predicted values using the user supplied parameter estimates. real(wp), intent(out) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. integer, intent(out) :: msgb(1 + q*np) !! Error checking results for the Jacobian wrt `beta`. integer, intent(out) :: msgd(1 + q*m) !! Error checking results for the Jacobian wrt `delta`. real(wp), intent(out) :: diff(q, np + m) !! Relative differences between the user supplied and finite difference derivatives !! for each derivative checked. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta` integer, intent(inout) :: nfev !! Number of function evaluations. integer, intent(inout) :: njev !! Number of Jacobian evaluations. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. integer, intent(in) :: interval(np) !! Specifies which checks can be performed when checking derivatives based on the !! interval of the bound constraints. ! Local scalars real(wp) :: diffj, h0, hc0, pv, tol, typj integer :: ideval, j, lq, msgb1, msgd1 logical :: isfixd, iswrtb ! Local arrays real(wp) :: pv0(n, q) ! Variable Definitions (alphabetically) ! DIFFJ: The relative differences between the user supplied and finite difference ! derivatives for the derivative being checked. ! H0: The initial relative step size for forward differences. ! HC0: The initial relative step size for central differences. ! IDEVAL: The variable designating what computations are to be performed by user supplied ! subroutine FCN. ! ISFIXD: The variable designating whether the parameter is fixed (TRUE) or not (FALSE). ! ISWRTB: The variable designating whether the derivatives wrt BETA (TRUE) or DELTA ! (FALSE) are being checked. ! J: An index variable. ! LQ: The response currently being examined. ! MSGB1: The error checking results for the Jacobian wrt BETA. ! MSGD1: The error checking results for the Jacobian wrt DELTA. ! PV: The scalar in which the predicted value from the model for row NROW is stored. ! PV0: The predicted values using the current parameter estimates (possibly offset from ! the user supplied estimates to create distance between parameters and the bounds ! on the parameters). ! TOL: The agreement tolerance. ! TYPJ: The typical size of the J-th unknown BETA or DELTA. ! Set tolerance for checking derivatives tol = eta**(0.25E0_wp) ntol = int(max(one, p5 - log10(tol))) ! Compute, if necessary, PV0 pv0 = pv0i if (any(beta /= betaj)) then istop = 0 ideval = 001 call fcn(betaj, xplusd, ifixb, ifixx, ideval, pv0, fjacb, fjacd, istop) if (istop /= 0) then return else njev = njev + 1 end if end if ! Compute user-supplied derivative values istop = 0 if (isodr) then ideval = 110 else ideval = 010 end if call fcn(betaj, xplusd, ifixb, ifixx, ideval, wrk2, fjacb, fjacd, istop) if (istop /= 0) then return else njev = njev + 1 end if ! Check derivatives wrt BETA for each response of observation NROW msgb1 = 0 msgd1 = 0 do lq = 1, q ! Set predicted value of model at current parameter estimates pv = pv0(nrow, lq) iswrtb = .true. do j = 1, np if (ifixb(1) < 0) then isfixd = .false. elseif (ifixb(j) == 0) then isfixd = .true. else isfixd = .false. end if if (isfixd) then msgb(1 + lq + (j - 1)*q) = -1 else if (beta(j) == zero) then if (ssf(1) < zero) then typj = one/abs(ssf(1)) else typj = one/ssf(j) end if else typj = abs(beta(j)) end if h0 = hstep(0, neta, 1, j, stpb, 1) hc0 = h0 ! Check derivative wrt the J-th parameter at the NROW-th row if (interval(j) >= 1) then call check_jac_value(fcn, & n, m, np, q, & betaj, xplusd, & ifixb, ifixx, ldifx, & eta, tol, nrow, epsmac, j, lq, typj, h0, hc0, & iswrtb, pv, fjacb(nrow, j, lq), & diffj, msgb1, msgb(2), istop, nfev, & wrk1, wrk2, wrk6, interval) if (istop /= 0) then msgb(1) = -1 return else diff(lq, j) = diffj end if else msgb(1 + j) = 9 end if end if end do ! Check derivatives wrt X for each response of observation NROW if (isodr) then iswrtb = .false. do j = 1, m if (ifixx(1, 1) < 0) then isfixd = .false. elseif (ldifx == 1) then if (ifixx(1, j) == 0) then isfixd = .true. else isfixd = .false. end if else isfixd = .false. end if if (isfixd) then msgd(1 + lq + (j - 1)*q) = -1 else if (xplusd(nrow, j) == zero) then if (tt(1, 1) < zero) then typj = one/abs(tt(1, 1)) elseif (ldtt == 1) then typj = one/tt(1, j) else typj = one/tt(nrow, j) end if else typj = abs(xplusd(nrow, j)) end if h0 = hstep(0, neta, nrow, j, stpd, ldstpd) hc0 = hstep(1, neta, nrow, j, stpd, ldstpd) ! Check derivative wrt the J-th column of DELTA at row NROW call check_jac_value(fcn, & n, m, np, q, & betaj, xplusd, & ifixb, ifixx, ldifx, & eta, tol, nrow, epsmac, j, lq, typj, h0, hc0, & iswrtb, pv, fjacd(nrow, j, lq), & diffj, msgd1, msgd(2), istop, nfev, & wrk1, wrk2, wrk6, interval) if (istop /= 0) then msgd(1) = -1 return else diff(lq, np + j) = diffj end if end if end do end if end do msgb(1) = msgb1 msgd(1) = msgd1 end subroutine check_jac subroutine check_jac_curv & (fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & eta, tol, nrow, epsmac, j, lq, hc, iswrtb, & fd, typj, pvpstp, stp0, & pv, d, & diffj, msg, istop, nfev, & wrk1, wrk2, wrk6) !! Check whether high curvature could be the cause of the disagreement between the numerical !! and analytic derviatives. ! Adapted from STARPAC subroutine DCKCRV. use odrpack_kinds, only: one, two, ten procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x` + `delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: eta !! Relative noise in the model. real(wp), intent(in) :: tol !! Agreement tolerance. integer, intent(in) :: nrow !! Row number of the explanatory variable array at which the derivative is to be checked. real(wp), intent(in) :: epsmac !! Value of machine precision. integer, intent(in) :: j !! Index of the partial derivative being examined. integer, intent(in) :: lq !! Response currently being examined. real(wp), intent(in) :: hc !! Relative step size for central finite differences. logical, intent(in) :: iswrtb !! Variable designating whether the derivatives wrt `beta` (`iswrtb = .true.`) or !! `delta` (`iswrtb = .false.`) are being checked. real(wp), intent(out) :: fd !! Forward difference derivative wrt the `j`-th parameter. real(wp), intent(in) :: typj !! Typical size of the `j`-th unknown `beta` or `delta`. real(wp), intent(out) :: pvpstp !! Predicted value for row `nrow` of the model based on the current parameter estimates !! for all but the `j`-th parameter value, which is `beta(j) + stp0`. real(wp), intent(in) :: stp0 !! Initial step size for the finite difference derivative. real(wp), intent(in) :: pv !! Predicted value of the model for row `nrow`. real(wp), intent(in) :: d !! Derivative with respect to the `j`-th unknown parameter. real(wp), intent(out) :: diffj !! Relative differences between the user supplied and finite difference derivatives !! for the derivative being checked. integer, intent(out) :: msg(q, j) !! Error checking results. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. ! Local scalars real(wp), parameter :: p01 = 0.01_wp real(wp) :: curve, pvmcrv, pvpcrv, stp, stpcrv ! Variable Definitions (alphabetically) ! CURVE: A measure of the curvature in the model. ! PVMCRV: The predicted value for row NROW of the model based on the current parameter ! estimates for all but the Jth parameter value, which is BETA(J)-STPCRV. ! PVPCRV: The predicted value for row NROW of the model based on the current parameter ! estimates for all but the Jth parameter value, which is BETA(J)+STPCRV. ! STP: A step size for the finite difference derivative. ! STPCRV: The step size selected to check for curvature in the model. if (iswrtb) then ! Perform central difference computations for derivatives wrt BETA stpcrv = (hc*typj*sign(one, beta(j)) + beta(j)) - beta(j) call fpvb(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stpcrv, & istop, nfev, pvpcrv, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if call fpvb(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, -stpcrv, & istop, nfev, pvmcrv, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if else ! Perform central difference computations for derivatives wrt DELTA stpcrv = (hc*typj*sign(one, xplusd(nrow, j)) + xplusd(nrow, j)) - xplusd(nrow, j) call fpvd(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stpcrv, & istop, nfev, pvpcrv, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if call fpvd(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, -stpcrv, & istop, nfev, pvmcrv, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if end if ! Estimate curvature by second derivative of model curve = abs((pvpcrv - pv) + (pvmcrv - pv))/(stpcrv*stpcrv) curve = curve + eta*(abs(pvpcrv) + abs(pvmcrv) + two*abs(pv))/(stpcrv**2) ! Check if finite precision arithmetic could be the culprit. call check_jac_fp(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & eta, tol, nrow, j, lq, iswrtb, & fd, typj, pvpstp, stp0, curve, pv, d, & diffj, msg, istop, nfev, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if if (msg(lq, j) == 0) then return end if ! Check if high curvature could be the problem. stp = two*max(tol*abs(d)/curve, epsmac) if (stp < abs(ten*stp0)) then stp = min(stp, p01*abs(stp0)) end if if (iswrtb) then ! Perform computations for derivatives wrt BETA stp = (stp*sign(one, beta(j)) + beta(j)) - beta(j) call fpvb(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp, & istop, nfev, pvpstp, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if else ! Perform computations for derivatives wrt DELTA stp = (stp*sign(one, xplusd(nrow, j)) + xplusd(nrow, j)) - xplusd(nrow, j) call fpvd(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp, & istop, nfev, pvpstp, & wrk1, wrk2, wrk6) if (istop /= 0) then return end if end if ! Compute the new numerical derivative fd = (pvpstp - pv)/stp diffj = min(diffj, abs(fd - d)/abs(d)) ! Check whether the new numerical derivative is ok if (abs(fd - d) <= tol*abs(d)) then msg(lq, j) = 0 ! Check if finite precision may be the culprit (fudge factor = 2) elseif (abs(stp*(fd - d)) < two*eta*(abs(pv) + abs(pvpstp)) + & curve*(epsmac*typj)**2) then msg(lq, j) = 5 end if end subroutine check_jac_curv subroutine check_jac_fp & (fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & eta, tol, nrow, j, lq, iswrtb, & fd, typj, pvpstp, stp0, curve, pv, d, & diffj, msg, istop, nfev, & wrk1, wrk2, wrk6) !! Check whether finite precision arithmetic could be the cause of the disagreement between !! the derivatives. ! Adapted from STARPAC subroutine DCKFPA. use odrpack_kinds, only: one, two, hundred procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: eta !! Relative noise in the model. real(wp), intent(in) :: tol !! Agreement tolerance. integer, intent(in) :: nrow !! Row number of the explanatory variable array at which the derivative is to be checked. integer, intent(in) :: j !! Index of the partial derivative being examined. integer, intent(in) :: lq !! Response currently being examined. logical, intent(in) :: iswrtb !! Variable designating whether the derivatives wrt `beta` (`iswrtb = .true.`) !! or `delta` (`iswrtb = .false.`) are being checked. real(wp), intent(out) :: fd !! Forward difference derivative wrt the `j`-th parameter. real(wp), intent(in) :: typj !! Typical size of the `j`-th unknown `beta` or `delta`. real(wp), intent(out) :: pvpstp !! Predicted value for row `nrow` of the model based on the current parameter !! estimates for all but the `j`-th parameter value, which is `beta(j) + stp0`. real(wp), intent(in) :: stp0 !! Step size for the finite difference derivative. real(wp), intent(inout) :: curve !! A measure of the curvature in the model. real(wp), intent(in) :: pv !! Predicted value for row `nrow`. real(wp), intent(in) :: d !! Derivative with respect to the `j`-th unknown parameter. real(wp), intent(out) :: diffj !! Relative differences between the user supplied and finite difference derivatives !! for the derivative being checked. integer, intent(out) :: msg(q, j) !! Error checking results. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. ! Local scalars real(wp), parameter :: p1 = 0.1_wp real(wp) :: stp logical :: large ! Variable Definitions (alphabetically) ! LARGE: The value designating whether the recommended increase in the step size would ! be greater than TYPJ. ! STP: A step size for the finite difference derivative. ! Finite precision arithmetic could be the problem. ! Try a larger step size based on estimate of condition error. stp = eta*(abs(pv) + abs(pvpstp))/(tol*abs(d)) if (stp > abs(p1*stp0)) then stp = max(stp, hundred*abs(stp0)) end if if (stp > typj) then stp = typj large = .true. else large = .false. end if if (iswrtb) then ! Perform computations for derivatives wrt BETA stp = (stp*sign(one, beta(j)) + beta(j)) - beta(j) call fpvb(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp, & istop, nfev, pvpstp, & wrk1, wrk2, wrk6) else ! Perform computations for derivatives wrt DELTA stp = (stp*sign(one, xplusd(nrow, j)) + xplusd(nrow, j)) - xplusd(nrow, j) call fpvd(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp, & istop, nfev, pvpstp, & wrk1, wrk2, wrk6) end if if (istop /= 0) then return end if fd = (pvpstp - pv)/stp diffj = min(diffj, abs(fd - d)/abs(d)) ! Check for agreement if ((abs(fd - d)) <= tol*abs(d)) then ! Forward difference quotient and analytic derivatives agree. msg(lq, j) = 0 elseif ((abs(fd - d) <= abs(two*curve*stp)) .or. large) then ! Curvature may be the culprit (fudge factor = 2) if (large) then msg(lq, j) = 4 else msg(lq, j) = 5 end if end if end subroutine check_jac_fp subroutine check_jac_value & (fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & eta, tol, nrow, epsmac, j, lq, typj, h0, hc0, & iswrtb, pv, d, & diffj, msg1, msg, istop, nfev, & wrk1, wrk2, wrk6, interval) !! Check user supplied analytic derivatives against numerical derivatives. ! Adapted from STARPAC subroutine DCKMN. use odrpack_kinds, only: zero, one, two, three, ten, hundred procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. real(wp), intent(in) :: eta !! Relative noise in the model. real(wp), intent(in) :: tol !! Agreement tolerance. integer, intent(in) :: nrow !! Row number of the explanatory variable array at which the derivative is to be checked. real(wp), intent(in) :: epsmac !! Value of machine precision. integer, intent(in) :: j !! Index of the partial derivative being examined. integer, intent(in) :: lq !! Response currently being examined. real(wp), intent(in) :: typj !! Typical size of the `j`-th unknown `beta` or `delta`. real(wp), intent(in) :: h0 !! Initial step size for the finite difference derivative. real(wp), intent(in) :: hc0 !! Relative step size for central finite differences. logical, intent(in) :: iswrtb !! Variable designating whether the derivatives wrt `beta` (`iswrtb = .true.`) !! or `delta` (`iswrtb = .false.`) are being checked. real(wp), intent(in) :: pv !! Predicted value for row `nrow`. real(wp), intent(in) :: d !! Derivative with respect to the `j`-th unknown parameter. real(wp), intent(out) :: diffj !! Relative differences between the user supplied and finite difference derivatives !! for the derivative being checked. integer, intent(out) :: msg1 !! First set of error checking results. integer, intent(out) :: msg(q, j) !! Error checking results. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. integer, intent(in) :: interval(np) !! Specifies which checks can be performed when checking derivatives based on the !! interval of the bound constraints. ! Local scalars real(wp), parameter :: p01 = 0.01_wp, p1 = 0.1_wp real(wp), parameter :: big = 1.0E19_wp, tol2 = 5.0E-2_wp real(wp) :: fd, h, hc, h1, hc1, pvpstp, stp0 integer :: i ! Variable Definitions (alphabetically) ! BIG: A big value, used to initialize DIFFJ. ! FD: The forward difference derivative wrt the Jth parameter. ! H: The relative step size for forward differences. ! H1: The default relative step size for forward differences. ! HC: The relative step size for central differences. ! HC1: The default relative step size for central differences. ! PVPSTP: The predicted value for row NROW of the model using the current parameter ! estimates for all but the Jth parameter value, which is BETA(J) + STP0. ! STP0: The initial step size for the finite difference derivative. ! TOL2: A minimum agreement tolerance. ! Calculate the Jth partial derivative using forward difference quotients and decide if ! it agrees with user supplied values h1 = sqrt(eta) hc1 = eta**(1/three) msg(lq, j) = 7 diffj = big do i = 1, 3 if (i == 1) then ! Try initial relative step size h = h0 hc = hc0 elseif (i == 2) then ! Try larger relative step size h = max(ten*h1, min(hundred*h0, one)) hc = max(ten*hc1, min(hundred*hc0, one)) elseif (i == 3) then ! Try smaller relative step size h = min(p1*h1, max(p01*h, two*epsmac)) hc = min(p1*hc1, max(p01*hc, two*epsmac)) end if if (iswrtb) then ! Perform computations for derivatives wrt BETA stp0 = (h*typj*sign(one, beta(j)) + beta(j)) - beta(j) call fpvb(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp0, & istop, nfev, pvpstp, & wrk1, wrk2, wrk6) else ! Perform computations for derivatives wrt DELTA stp0 = (h*typj*sign(one, xplusd(nrow, j)) + xplusd(nrow, j)) - xplusd(nrow, j) call fpvd(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp0, & istop, nfev, pvpstp, & wrk1, wrk2, wrk6) end if if (istop /= 0) then return end if fd = (pvpstp - pv)/stp0 ! Check for agreement ! Numerical and analytic derivatives agree if (abs(fd - d) <= tol*abs(d)) then ! Set relative difference for derivative checking report if ((d == zero) .or. (fd == zero)) then diffj = abs(fd - d) else diffj = abs(fd - d)/abs(d) end if ! Set message flag if (d == zero) then ! JTH analytic and numerical derivatives are both zero. msg(lq, j) = 1 else ! JTH analytic and numerical derivatives are both nonzero. msg(lq, j) = 0 end if else ! Numerical and analytic derivatives disagree. Check why if ((d == zero) .or. (fd == zero)) then if (interval(j) >= 10 .or. .not. iswrtb) then call check_jac_zero(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, epsmac, j, lq, iswrtb, & tol, d, fd, typj, pvpstp, stp0, pv, & diffj, msg, istop, nfev, & wrk1, wrk2, wrk6) else msg(lq, j) = 8 end if else if (interval(j) >= 100 .or. .not. iswrtb) then call check_jac_curv(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & eta, tol, nrow, epsmac, j, lq, hc, iswrtb, & fd, typj, pvpstp, stp0, pv, d, & diffj, msg, istop, nfev, & wrk1, wrk2, wrk6) else msg(lq, j) = 8 end if end if if (msg(lq, j) <= 2) then exit end if end if end do ! Set summary flag to indicate questionable results if ((msg(lq, j) >= 7) .and. (diffj <= tol2)) then msg(lq, j) = 6 end if if ((msg(lq, j) >= 1) .and. (msg(lq, j) <= 6)) then msg1 = max(msg1, 1) elseif (msg(lq, j) >= 7) then msg1 = 2 end if end subroutine check_jac_value subroutine check_jac_zero & (fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, epsmac, j, lq, iswrtb, & tol, d, fd, typj, pvpstp, stp0, pv, & diffj, msg, istop, nfev, & wrk1, wrk2, wrk6) !! Recheck the derivatives in the case where the finite difference derivative disagrees with !! the analytic derivative and the analytic derivative is zero. ! Adapted from STARPAC subroutine DCKZRO. use odrpack_kinds, only: zero, three procedure(fcn_t) :: fcn !! User supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. integer, intent(in) :: nrow !! Row number of the explanatory variable array at which the derivative is to be checked. real(wp), intent(in) :: epsmac !! Value of machine precision. integer, intent(in) :: j !! Index of the partial derivative being examined. integer, intent(in) :: lq !! Response currently being examined. logical, intent(in) :: iswrtb !! Variable designating whether the derivatives wrt `beta` (`iswrtb = .true.`) !! or `delta` (`iswrtb = .false.`) are being checked. real(wp), intent(in) :: tol !! Agreement tolerance. real(wp), intent(in) :: d !! Derivative with respect to the `j`-th unknown parameter. real(wp), intent(in) :: fd !! Forward difference derivative wrt the `j`-th parameter. real(wp), intent(in) :: typj !! Typical size of the `j`-th unknown `beta` or `delta`. real(wp), intent(in) :: pvpstp !! Predicted value for row `nrow` of the model using the current parameter estimates !! for all but the `j`-th parameter value, which is `beta(j) + stp0`. real(wp), intent(in) :: stp0 !! Initial step size for the finite difference derivative. real(wp), intent(in) :: pv !! Predicted value from the model for row `nrow`. real(wp), intent(out) :: diffj !! Relative differences between the user supplied and finite difference derivatives !! for the derivative being checked. integer, intent(out) :: msg(q, j) !! Error checking results. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: wrk1(n, m, q) !! A work array of `(n, m, q)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk6(n, np, q) !! A work array of `(n, np, q)` elements. ! Local scalars real(wp) :: cd, pvmstp ! Variable Definitions (alphabetically) ! CD: The central difference derivative wrt the Jth parameter. ! PVMSTP: The predicted value for row NROW of the model using the current parameter ! estimates for all but the Jth parameter value, which is BETA(J) - STP0. ! Recalculate numerical derivative using central difference and step size of 2*STP0 if (iswrtb) then ! Perform computations for derivatives wrt BETA call fpvb(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, -stp0, & istop, nfev, pvmstp, & wrk1, wrk2, wrk6) else ! Perform computations for derivatives wrt DELTA call fpvd(fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, -stp0, & istop, nfev, pvmstp, & wrk1, wrk2, wrk6) end if if (istop /= 0) then return end if cd = (pvpstp - pvmstp)/(2*stp0) diffj = min(abs(cd - d), abs(fd - d)) ! Check for agreement if (diffj <= tol*abs(d)) then ! Finite difference and analytic derivatives now agree if (d == zero) then msg(lq, j) = 1 else msg(lq, j) = 0 end if elseif (diffj*typj <= abs(pv*epsmac**(1/three))) then ! Derivatives are both close to zero msg(lq, j) = 2 else ! Derivatives are not both close to zero msg(lq, j) = 3 end if end subroutine check_jac_zero pure subroutine check_inputs & (n, m, np, q, & isodr, anajac, & beta, ifixb, & ldifx, ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, & lrwork, lrwkmin, liwork, liwkmin, & sclb, scld, stpb, stpd, & info, & lower, upper) !! Check input parameters, indicating errors found using nonzero values of argument `info`. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). logical, intent(in) :: anajac !! Variable designating whether the Jacobians are computed by finite differences !! (`.false.`) or not (`.true.`). real(wp), intent(in) :: beta(np) !! Function parameters. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. integer, intent(in) :: ldscld !! Leading dimension of array `scld`. integer, intent(in) :: ldstpd !! Leading dimension of array `stpd`. integer, intent(in) :: ldwe !! Leading dimension of array `we`. integer, intent(in) :: ld2we !! Second dimension of array `we`. integer, intent(in) :: ldwd !! Leading dimension of array `wd`. integer, intent(in) :: ld2wd !! Second dimension of array `wd`. integer, intent(in) :: lrwork !! Length of vector `rwork`. integer, intent(in) :: lrwkmin !! Minimum acceptable length of array `rwork`. integer, intent(in) :: liwork !! Length of vector `iwork`. integer, intent(in) :: liwkmin !! Minimum acceptable length of array `iwork`. real(wp), intent(in) :: sclb(np) !! Scaling values for `beta`. real(wp), intent(in) :: scld(ldscld, m) !! Scaling value for `delta`. real(wp), intent(in) :: stpb(np) !! Step for the finite difference derivative wrt `beta`. real(wp), intent(in) :: stpd(ldstpd, m) !! Step for the finite difference derivative wrt `delta`. integer, intent(out) :: info !! Variable designating why the computations were stopped. real(wp), intent(in) :: lower(np) !! Lower bound on `beta`. real(wp), intent(in) :: upper(np) !! Upper bound on `beta`. ! Local scalars integer :: last, npp ! Variable Definitions (alphabetically) ! LAST: The last row of the array to be accessed. ! NPP: The number of function parameters being estimated. ! Find actual number of parameters being estimated if ((np <= 0) .or. (ifixb(1) < 0)) then npp = np else npp = count(ifixb /= 0) end if ! Check problem specification parameters if ((n <= 0) .or. (m <= 0) .or. (npp <= 0 .or. npp > n) .or. (q <= 0)) then info = 10000 if (n <= 0) then info = info + 1000 end if if (m <= 0) then info = info + 100 end if if (npp <= 0 .or. npp > n) then info = info + 10 end if if (q <= 0) then info = info + 1 end if return end if ! Check dimension specification parameters if (((ldwe /= 1) .and. (ldwe < n)) .or. & ((ld2we /= 1) .and. (ld2we < q)) .or. & (isodr .and. ((ldwd /= 1) .and. (ldwd < n))) .or. & (isodr .and. ((ld2wd /= 1) .and. (ld2wd < m))) .or. & (isodr .and. ((ldifx /= 1) .and. (ldifx < n))) .or. & (isodr .and. ((ldstpd /= 1) .and. (ldstpd < n))) .or. & (isodr .and. ((ldscld /= 1) .and. (ldscld < n))) .or. & (lrwork < lrwkmin) .or. & (liwork < liwkmin)) then info = 20000 if ((ldwe /= 1 .and. ldwe < n) .or. (ld2we /= 1 .and. ld2we < q)) then info = info + 100 end if if (isodr .and. & ((ldwd /= 1 .and. ldwd < n) .or. (ld2wd /= 1 .and. ld2wd < m))) then info = info + 200 end if if (isodr .and. (ldifx /= 1 .and. ldifx < n)) then info = info + 10 end if if (isodr .and. (ldstpd /= 1 .and. ldstpd < n)) then info = info + 20 end if if (isodr .and. & (ldscld /= 1 .and. ldscld < n)) then info = info + 40 end if if (lrwork < lrwkmin) then info = info + 1 end if if (liwork < liwkmin) then info = info + 2 end if end if ! Check DELTA scaling if (isodr .and. scld(1, 1) > zero) then if (ldscld >= n) then last = n else last = 1 end if if (any(scld(1:last, :) <= zero)) then info = 30200 end if end if ! Check BETA scaling if (sclb(1) > zero) then if (any(sclb <= zero)) then if (info == 0) then info = 30100 else info = info + 100 end if end if end if ! Check DELTA finite difference step sizes if ((.not. anajac) .and. isodr .and. stpd(1, 1) > zero) then if (ldstpd >= n) then last = n else last = 1 end if if (any(stpd(1:last, :) <= zero)) then if (info == 0) then info = 32000 else info = info + 2000 end if end if end if ! Check BETA finite difference step sizes if ((.not. anajac) .and. stpb(1) > zero) then if (any(stpb <= zero)) then if (info == 0) then info = 31000 else info = info + 1000 end if end if end if ! Check bounds if (any(upper < lower)) then if (info == 0) then info = 91000 end if end if if (any( & ((upper < beta) .or. (lower > beta)) & .and. .not. (upper < lower))) then if (info >= 90000) then info = info + 100 else info = 90100 end if end if end subroutine check_inputs subroutine lcstep & (n, m, np, q, npp, & f, fjacb, fjacd, & wd, ldwd, ld2wd, ss, tt, ldtt, delta, & alpha, epsfcn, isodr, & tfjacb, omega, u, qraux, kpvt, & s, t, phi, irank, rcond, forvcv, & wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc) !! Compute locally constrained steps `s` and `t`, and `phi(alpha)`. ! @note: This is one of the most time-consuming subroutines in ODRPACK (~25% of total). use odrpack_kinds, only: zero, one use linpack, only: dchex, dqrdc, dqrsl, dtrco, dtrsl use blas_interfaces, only: drot, drotg integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: npp !! Number of function parameters being estimated. real(wp), intent(in) :: f(n, q) !! Weighted estimated values of `epsilon`. real(wp), intent(in) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. real(wp), intent(in) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. real(wp), intent(in) :: wd(ldwd, ld2wd, m) !! Squared `delta` weights. integer, intent(in) :: ldwd !! Leading dimension of array `wd`. integer, intent(in) :: ld2wd !! Second dimension of array `wd`. real(wp), intent(in) :: ss(np) !! Scaling values for the unfixed `beta`s. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. real(wp), intent(in) :: delta(n, m) !! Estimated errors in the explanatory variables. real(wp), intent(in) :: alpha !! Levenberg-Marquardt parameter. real(wp), intent(in) :: epsfcn !! Function's precision. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(out) :: tfjacb(n, q, np) !! Array `omega*fjacb`. real(wp), intent(out) :: omega(q, q) !! Array defined such that: !! `omega*trans(omega) = inv(I + fjacd*inv(e)*trans(fjacd)) !! = (I - fjacd*inv(p)*trans(fjacd))` !! where `e = d**2 + alpha*tt**2` and !! `p = trans(fjacd)*fjacd + d**2 + alpha*tt**2`. real(wp), intent(out) :: u(np) !! Approximate null vector for `tfjacb`. real(wp), intent(out) :: qraux(np) !! Array required to recover the orthogonal part of the Q-R decomposition. integer, intent(out) :: kpvt(np) !! Pivot vector. real(wp), intent(out) :: s(np) !! Step for `beta`. real(wp), intent(out) :: t(n, m) !! Step for `delta`. real(wp), intent(out) :: phi !! Difference between the norm of the scaled step and the trust region diameter. integer, intent(out) :: irank !! Rank deficiency of the Jacobian wrt `beta`. real(wp), intent(out) :: rcond !! Approximate reciprocal condition number of `tfjacb`. logical, intent(in) :: forvcv !! Variable designating whether this subroutine was called to set up for the !! covariance matrix computations (`forvcv = .true.`) or not (`forvcv = .false.`). real(wp), intent(out) :: wrk1(n, q, m) !! A work array of `(n, q, m)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk3(np) !! A work array of `(np)` elements. real(wp), intent(out) :: wrk4(m, m) !! A work array of `(m, m)` elements. real(wp), intent(out) :: wrk5(m) !! A work array of `(m)` elements. real(wp), intent(out) :: wrk(lwrk) !! A work array of `(lwrk)` elements, _equivalenced_ to `wrk1` and `wrk2`. integer, intent(in) :: lwrk !! Length of vector `wrk`. integer, intent(inout) :: istopc !! Variable designating whether the computations were stopped due to a numerical !! error within subroutine `dodstp`. ! Local scalars real(wp) :: co, si, temp integer :: i, imax, inf, ipvt, j, k, k1, k2, kp, l logical :: elim ! Local arrays real(wp) :: dum(2) ! Variable definitions (alphabetically) ! CO: The cosine from the plane rotation. ! DUM: A dummy array. ! ELIM: The variable designating whether columns of the Jacobian ! I: An indexing variable. ! IMAX: The index of the element of U having the largest absolute value. ! INF: The return code from LINPACK routines. ! IPVT: The variable designating whether pivoting is to be done. ! J: An indexing variable. ! K: An indexing variable. ! K1: An indexing variable. ! K2: An indexing variable. ! KP: The rank of the Jacobian wrt BETA. ! L: An indexing variable. ! SI: The sine from the plane rotation. ! TEMP: A temporary storage LOCATION. ! Compute loop parameters which depend on weight structure ! Set up KPVT if ALPHA = 0 if (alpha == zero) then kp = npp do k = 1, np kpvt(k) = k end do else if (npp >= 1) then kp = npp - irank else kp = npp end if end if if (isodr) then ! T = WD * DELTA = D*G2 call scale_mat(n, m, wd, ldwd, ld2wd, delta, t) do i = 1, n ! Compute WRK4, such that TRANS(WRK4)*WRK4 = E = (D**2 + ALPHA*TT**2) call esubi(n, m, wd, ldwd, ld2wd, alpha, tt, ldtt, i, wrk4) call fctr(.false., wrk4, m, m, inf) if (inf /= 0) then istopc = 60000 return end if ! Compute OMEGA, such that ! trans(OMEGA)*OMEGA = I+FJACD*inv(E)*trans(FJACD) ! inv(trans(OMEGA)*OMEGA) = I-FJACD*inv(P)*trans(FJACD) call vevtr(m, q, i, fjacd, n, m, wrk4, m, wrk1, n, q, omega, q, wrk5) do l = 1, q omega(l, l) = one + omega(l, l) end do call fctr(.false., omega, q, q, inf) if (inf /= 0) then istopc = 60000 return end if ! Compute WRK1 = trans(FJACD)*(I-FJACD*inv(P)*trans(JFACD)) ! = trans(FJACD)*inv(trans(OMEGA)*OMEGA) do j = 1, m wrk1(i, :, j) = fjacd(i, j, :) call solve_trl(q, omega, q, wrk1(i, 1:q, j), 4) call solve_trl(q, omega, q, wrk1(i, 1:q, j), 2) end do ! Compute WRK5 = inv(E)*D*G2 wrk5 = t(i, :) call solve_trl(m, wrk4, m, wrk5, 4) call solve_trl(m, wrk4, m, wrk5, 2) ! Compute TFJACB = inv(trans(OMEGA))*FJACB do k = 1, kp tfjacb(i, :, k) = fjacb(i, kpvt(k), :) call solve_trl(q, omega, q, tfjacb(i, 1:q, k), 4) if (ss(1) > zero) then tfjacb(i, :, k) = tfjacb(i, :, k)/ss(kpvt(k)) else tfjacb(i, :, k) = tfjacb(i, :, k)/abs(ss(1)) end if end do ! Compute WRK2 = (V*inv(E)*D**2*G2 - G1) do l = 1, q wrk2(i, l) = dot_product(fjacd(i, :, l), wrk5) - f(i, l) end do ! Compute WRK2 = inv(trans(OMEGA))*(V*inv(E)*D**2*G2 - G1) call solve_trl(q, omega, q, wrk2(i, 1:q), 4) end do else do l = 1, q do i = 1, n do k = 1, kp tfjacb(i, l, k) = fjacb(i, kpvt(k), l) if (ss(1) > zero) then tfjacb(i, l, k) = tfjacb(i, l, k)/ss(kpvt(k)) else tfjacb(i, l, k) = tfjacb(i, l, k)/abs(ss(1)) end if end do wrk2(i, l) = -f(i, l) end do end do end if ! Compute S ! Do QR factorization (with column pivoting of TFJACB if ALPHA = 0) if (alpha == zero) then ipvt = 1 kpvt = 0 else ipvt = 0 end if call dqrdc(tfjacb, n*q, n*q, kp, qraux, kpvt, wrk3, ipvt) call dqrsl(tfjacb, n*q, n*q, kp, qraux, wrk2, dum, wrk2, dum, dum, dum, 1000, inf) if (inf /= 0) then istopc = 60000 return end if ! Eliminate alpha part using givens rotations if (alpha /= zero) then s(1:npp) = zero do k1 = 1, kp wrk3(1:kp) = zero wrk3(k1) = sqrt(alpha) do k2 = k1, kp call drotg(tfjacb(k2, 1, k2), wrk3(k2), co, si) if (kp - k2 >= 1) then call drot(kp - k2, tfjacb(k2, 1, k2 + 1), n*q, wrk3(k2 + 1), 1, co, si) end if temp = co*wrk2(k2, 1) + si*s(kpvt(k1)) s(kpvt(k1)) = -si*wrk2(k2, 1) + co*s(kpvt(k1)) wrk2(k2, 1) = temp end do end do end if ! Compute solution - eliminate variables if necessary if (npp >= 1) then if (alpha == zero) then kp = npp elim = .true. do while (elim .and. kp >= 1) ! Estimate RCOND - U will contain approx null vector call dtrco(tfjacb, n*q, kp, rcond, u, 1) if (rcond <= epsfcn) then elim = .true. imax = maxloc(u(1:kp), dim=1) ! IMAX is the column to remove - use DCHEX and fix KPVT if (imax /= kp) then call dchex(tfjacb, n*q, kp, imax, kp, wrk2, n*q, 1, qraux, wrk3, 2) k = kpvt(imax) do i = imax, kp - 1 kpvt(i) = kpvt(i + 1) end do kpvt(kp) = k end if kp = kp - 1 else elim = .false. end if end do irank = npp - kp end if end if if (forvcv) return ! Backsolve and unscramble if (npp >= 1) then do i = kp + 1, npp wrk2(i, 1) = zero end do if (kp >= 1) then call dtrsl(tfjacb, n*q, kp, wrk2, 01, inf) if (inf /= 0) then istopc = 60000 return end if end if do i = 1, npp if (ss(1) > zero) then s(kpvt(i)) = wrk2(i, 1)/ss(kpvt(i)) else s(kpvt(i)) = wrk2(i, 1)/abs(ss(1)) end if end do end if if (isodr) then ! NOTE: T and WRK1 have been initialized above, ! where T = WD * DELTA = D*G2 ! WRK1 = trans(FJACD)*(I-FJACD*inv(P)*trans(JFACD)) do i = 1, n ! Compute WRK4, such that trans(WRK4)*WRK4 = E = (D**2 + ALPHA*TT**2) call esubi(n, m, wd, ldwd, ld2wd, alpha, tt, ldtt, i, wrk4) call fctr(.false., wrk4, m, m, inf) if (inf /= 0) then istopc = 60000 return end if ! Compute WRK5 = inv(E)*D*G2 wrk5 = t(i, :) call solve_trl(m, wrk4, m, wrk5, 4) call solve_trl(m, wrk4, m, wrk5, 2) do l = 1, q wrk2(i, l) = f(i, l) & + dot_product(fjacb(i, 1:npp, l), s(1:npp)) & - dot_product(fjacd(i, :, l), wrk5) end do do j = 1, m wrk5(j) = dot_product(wrk1(i, :, j), wrk2(i, :)) t(i, j) = -(wrk5(j) + t(i, j)) end do call solve_trl(m, wrk4, m, t(i, 1:m), 4) call solve_trl(m, wrk4, m, t(i, 1:m), 2) end do end if ! Compute PHI(ALPHA) from scaled S and T call scale_mat(npp, 1, ss, npp, 1, s, wrk(1:npp)) if (isodr) then call scale_mat(n, m, tt, ldtt, 1, t, wrk(npp + 1:npp + 1 + n*m - 1)) phi = norm2(wrk(1:npp + n*m)) else phi = norm2(wrk(1:npp)) end if end subroutine lcstep subroutine vcv_beta & (n, m, np, q, npp, & f, fjacb, fjacd, & wd, ldwd, ld2wd, ssf, ss, tt, ldtt, delta, & epsfcn, isodr, & vcv, sd, & wrk6, omega, u, qraux, jpvt, & s, t, irank, rcond, rss, idf, rvar, ifixb, & wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc) !! Compute covariance matrix of estimated parameters. use odrpack_kinds, only: zero use linpack, only: dpodi integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the explanatory variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: npp !! Number of function parameters being estimated. real(wp), intent(in) :: f(n, q) !! Weighted estimated values of `epsilon`. real(wp), intent(in) :: fjacb(n, np, q) !! Jacobian with respect to `beta`. real(wp), intent(in) :: fjacd(n, m, q) !! Jacobian with respect to `delta`. real(wp), intent(in) :: wd(ldwd, ld2wd, m) !! `delta` weights. integer, intent(in) :: ldwd !! Leading dimension of array `wd`. integer, intent(in) :: ld2wd !! Second dimension of array `wd`. real(wp), intent(in) :: ssf(np) !! Scaling values used for `beta`. real(wp), intent(in) :: ss(np) !! Scaling values for the unfixed `beta`s. real(wp), intent(in) :: tt(ldtt, m) !! Scaling values for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. real(wp), intent(in) :: delta(n, m) !! Estimated errors in the explanatory variables. real(wp), intent(in) :: epsfcn !! Function's precision. logical, intent(in) :: isodr !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`). real(wp), intent(out) :: vcv(np, np) !! Covariance matrix of the estimated `beta`s. real(wp), intent(out) :: sd(np) !! Standard deviations of the estimated `beta`s. real(wp), intent(out) :: wrk6(n*q, np) !! A work array of `(n*q, np)` elements. real(wp), intent(out) :: omega(q, q) !! Array defined such that `omega*trans(omega) = inv(I + fjacd*inv(e)*trans(fjacd)) !! = (I - fjacd*inv(p)*trans(fjacd))`. real(wp), intent(out) :: u(np) !! Approximate null vector for `fjacb`. real(wp), intent(out) :: qraux(np) !! Array required to recover the orthogonal part of the Q-R decomposition. integer, intent(out) :: jpvt(np) !! Pivot vector. real(wp), intent(out) :: s(np) !! Step for `beta`. real(wp), intent(out) :: t(n, m) !! Step for `delta`. integer, intent(out) :: irank !! Rank deficiency of the Jacobian wrt `beta`. real(wp), intent(out) :: rcond !! Approximate reciprocal condition of `fjacb`. real(wp), intent(inout) :: rss !! Residual sum of squares. integer, intent(out) :: idf !! Degrees of freedom of the fit, equal to the number of observations with nonzero !! weighted derivatives minus the number of parameters being estimated. real(wp), intent(out) :: rvar !! Residual variance. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. real(wp), intent(out) :: wrk1(n, q, m) !! A work array of `(n, q, m)` elements. real(wp), intent(out) :: wrk2(n, q) !! A work array of `(n, q)` elements. real(wp), intent(out) :: wrk3(np) !! A work array of `(np)` elements. real(wp), intent(out) :: wrk4(m, m) !! A work array of `(m, m)` elements. real(wp), intent(out) :: wrk5(m) !! A work array of `(m)` elements. real(wp), intent(out) :: wrk(lwrk) !! A work array of `(lwrk)` elements, _equivalenced_ to `wrk1` and `wrk2`. integer, intent(in) :: lwrk !! Length of vector `lwrk`. integer, intent(out) :: istopc !! Variable designating whether the computations were stoped due to a numerical !! error within subroutine `dodstp`. ! Local scalars real(wp) :: temp integer :: i, iunfix, j, junfix, kp logical :: forvcv ! Variable definitions (alphabetically) ! FORVCV: The variable designating whether subroutine DODSTP is called to set up for ! the covariance matrix computations (TRUE) or not (FALSE). ! I: An indexing variable. ! IUNFIX: The index of the next unfixed parameter. ! J: An indexing variable. ! JUNFIX: The index of the next unfixed parameter. ! KP: The rank of the Jacobian wrt BETA. ! TEMP: A temporary storage location forvcv = .true. istopc = 0 call lcstep(n, m, np, q, npp, & f, fjacb, fjacd, & wd, ldwd, ld2wd, ss, tt, ldtt, delta, & zero, epsfcn, isodr, & wrk6, omega, u, qraux, jpvt, & s, t, temp, irank, rcond, forvcv, & wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc) if (istopc /= 0) then return end if kp = npp - irank call dpodi(wrk6, n*q, kp, wrk3, 1) idf = 0 do i = 1, n if (any(fjacb(i, :, :) /= zero)) then idf = idf + 1 cycle end if if (isodr) then if (any(fjacd(i, :, :) /= zero)) then idf = idf + 1 cycle end if end if end do if (idf > kp) then idf = idf - kp rvar = rss/idf else idf = 0 rvar = rss end if ! Store variances in SD, restoring original order sd = zero do i = 1, kp sd(jpvt(i)) = wrk6(i, i) end do if (np > npp) then junfix = npp do j = np, 1, -1 if (ifixb(j) == 0) then sd(j) = zero else sd(j) = sd(junfix) junfix = junfix - 1 end if end do end if ! Store covariance matrix in VCV, restoring original order vcv = zero do i = 1, kp do j = i + 1, kp if (jpvt(i) > jpvt(j)) then vcv(jpvt(i), jpvt(j)) = wrk6(i, j) else vcv(jpvt(j), jpvt(i)) = wrk6(i, j) end if end do end do if (np > npp) then iunfix = npp do i = np, 1, -1 if (ifixb(i) == 0) then do j = i, 1, -1 vcv(i, j) = zero end do else junfix = npp do j = np, 1, -1 if (ifixb(j) == 0) then vcv(i, j) = zero else vcv(i, j) = vcv(iunfix, junfix) junfix = junfix - 1 end if end do iunfix = iunfix - 1 end if end do end if do i = 1, np vcv(i, i) = sd(i) sd(i) = sqrt(rvar*sd(i)) do j = 1, i vcv(j, i) = vcv(i, j) end do end do ! Unscale standard errors and covariance matrix do i = 1, np if (ssf(1) > zero) then sd(i) = sd(i)/ssf(i) else sd(i) = sd(i)/abs(ssf(1)) end if do j = 1, np if (ssf(1) > zero) then vcv(i, j) = vcv(i, j)/(ssf(i)*ssf(j)) else vcv(i, j) = vcv(i, j)/(ssf(1)*ssf(1)) end if end do end do end subroutine vcv_beta pure subroutine pack_vec(n2, n1, v1, v2, ifix) !! Select the unfixed elements of `v2` and return them in `v1`. integer, intent(in) :: n2 !! Number of items in `v2`. integer, intent(out) :: n1 !! Number of items in `v1`. real(wp), intent(out) :: v1(n2) !! Vector of the unfixed items from `v2`. real(wp), intent(in) :: v2(n2) !! Vector of the fixed and unfixed items from which the unfixed elements are to be extracted. integer, intent(in) :: ifix(n2) !! Values designating whether the elements of `v2` are fixed at their input values or not. ! Local scalars integer :: i ! Variable definitions (alphabetically) ! I: An indexing variable. n1 = 0 if (ifix(1) >= 0) then do i = 1, n2 if (ifix(i) /= 0) then n1 = n1 + 1 v1(n1) = v2(i) end if end do else n1 = n2 v1 = v2 end if end subroutine pack_vec pure subroutine unpack_vec(n2, v1, v2, ifix) !! Copy the elements of `v1` into the locations of `v2` which are unfixed. integer, intent(in) :: n2 !! Number of items in `v2`. real(wp), intent(in) :: v1(n2) !! Vector of the unfixed items. real(wp), intent(out) :: v2(n2) !! Vector of the fixed and unfixed items into which the elements of `v1` are to !! be inserted. integer, intent(in) :: ifix(n2) !! Values designating whether the elements of `v2` are fixed at their input values !! or not. ! Local scalars integer :: i, n1 ! Variable Definitions (alphabetically) ! I: An indexing variable. ! N1: The number of items in V1. n1 = 0 if (ifix(1) >= 0) then do i = 1, n2 if (ifix(i) /= 0) then n1 = n1 + 1 v2(i) = v1(n1) end if end do else n1 = n2 v2 = v1 end if end subroutine unpack_vec real(wp) pure function ppf_normal(p) result(res) !! Compute the percent point function value for the normal (Gaussian) distribution with !! mean 0 and standard deviation 1, and with probability density function: !! !! `f(x) = (1/sqrt(2*pi))*exp(-x^2/2)` !! ! Adapted from DATAPAC subroutine TPPF, with modifications to facilitate conversion to ! real(wp) automatically. !***Author Filliben, James J., ! Statistical Engineering Division ! National Bureau of Standards ! Washington, D. C. 20234 ! (Original Version--June 1972. ! (Updated --September 1975, ! November 1975, AND ! October 1976. !***Description ! --The coding as presented below is essentially ! identical to that presented by Odeh and Evans ! as Algortihm 70 of Applied Statistics. ! --As pointed out by Odeh and Evans in Applied ! Statistics, their algorithm representes a ! substantial improvement over the previously employed ! Hastings approximation for the normal percent point ! function, with accuracy improving from 4.5*(10**-4) ! to 1.5*(10**-8). !***References Odeh and Evans, the Percentage Points of the Normal ! Distribution, Algortihm 70, Applied Statistics, 1974, ! Pages 96-97. ! Evans, Algorithms for Minimal Degree Polynomial and ! Rational Approximation, M. Sc. Thesis, 1972, ! University of Victoria, B. C., Canada. ! Hastings, Approximations for Digital Computers, 1955, ! Pages 113, 191, 192. ! National Bureau of Standards Applied Mathematics ! Series 55, 1964, Page 933, Formula 26.2.23. ! Filliben, Simple and Robust Linear Estimation of the ! Location Parameter of a Symmetric Distribution ! (Unpublished Ph.D. Dissertation, Princeton ! University), 1969, Pages 21-44, 229-231. ! Filliben, "The Percent Point Function", ! (Unpublished Manuscript), 1970, Pages 28-31. ! Johnson and Kotz, Continuous Univariate Distributions, ! Volume 1, 1970, Pages 40-111. ! Kelley Statistical Tables, 1948. ! Owen, Handbook of Statistical Tables, 1962, Pages 3-16. ! Pearson and Hartley, Biometrika Tables for ! Statisticians, Volume 1, 1954, Pages 104-113. use odrpack_kinds, only: zero, half, one real(wp), intent(in) :: p !! Probability at which the percent point is to be evaluated; `0 < p < 1`. ! Local scalars real(wp), parameter :: p0 = -0.322232431088E0_wp, & p1 = -1.0E0_wp, & p2 = -0.342242088547E0_wp, & p3 = -0.204231210245E-1_wp, & p4 = -0.453642210148E-4_wp, & q0 = 0.993484626060E-1_wp, & q1 = 0.588581570495E0_wp, & q2 = 0.531103462366E0_wp, & q3 = 0.103537752850E0_wp, & q4 = 0.38560700634E-2_wp real(wp) :: aden, anum, r, t ! Variable Definitions (alphabetically) ! ADEN: A value used in the approximation. ! ANUM: A value used in the approximation. ! P0: A parameter used in the approximation. ! P1: A parameter used in the approximation. ! P2: A parameter used in the approximation. ! P3: A parameter used in the approximation. ! P4: A parameter used in the approximation. ! Q0: A parameter used in the approximation. ! Q1: A parameter used in the approximation. ! Q2: A parameter used in the approximation. ! Q3: A parameter used in the approximation. ! Q4: A parameter used in the approximation. ! R: The probability at which the percent point is evaluated. ! T: A value used in the approximation. if (p == half) then res = zero else r = p if (p > half) r = one - r t = sqrt(-2*log(r)) anum = ((((t*p4 + p3)*t + p2)*t + p1)*t + p0) aden = ((((t*q4 + q3)*t + q2)*t + q1)*t + q0) res = t + (anum/aden) if (p < half) res = -res end if end function ppf_normal real(wp) pure function ppf_tstudent(p, idf) result(res) !! Compute the percent point function value for the student's T distribution with `idf` !! degrees of freedom. ! Adapted from DATAPAC subroutine TPPF, with modifications to facilitate conversion to ! real(wp) automatically. !***Author Filliben, James J., ! Statistical Engineering Division ! National Bureau of Standards ! Washington, D. C. 20234 ! (Original Version--October 1975.) ! (Updated --November 1975.) !***Description ! --For IDF = 1 AND IDF = 2, the percent point function ! for the T distribution exists in simple closed form ! and so the computed percent points are exact. ! --For IDF between 3 and 6, inclusively, the approximation ! is augmented by 3 iterations of Newton's method to ! improve the accuracy, especially for P near 0 or 1. !***References National Bureau of Standards Applied Mathmatics ! Series 55, 1964, Page 949, Formula 26.7.5. ! Johnson and Kotz, Continuous Univariate Distributions, ! Volume 2, 1970, Page 102, Formula 11. ! Federighi, "Extended Tables of the Percentage Points ! of Student"S T Distribution, Journal of the American ! Statistical Association, 1969, Pages 683-688. ! Hastings and Peacock, Statistical Distributions, A ! Handbook for Students and Practitioners, 1975, ! Pages 120-123. use odrpack_kinds, only: pi, zero, half, one, two, three, eight, fiftn real(wp), intent(in) :: p !! Probability at which the percent point is to be evaluated; `0 < p < 1`. integer, intent(in) :: idf !! Degrees of freedom. ! Local scalars real(wp), parameter :: b21 = 4.0E0_wp, & b31 = 96.0E0_wp, & b32 = 5.0E0_wp, & b33 = 16.0E0_wp, & b34 = 3.0E0_wp, & b41 = 384.0E0_wp, & b42 = 3.0E0_wp, & b43 = 19.0E0_wp, & b44 = 17.0E0_wp, & b45 = -15.0E0_wp, & b51 = 9216.0E0_wp, & b52 = 79.0E0_wp, & b53 = 776.0E0_wp, & b54 = 1482.0E0_wp, & b55 = -1920.0E0_wp, & b56 = -945.0E0_wp real(wp) :: arg, c, con, d1, d3, d5, d7, d9, df, ppfn, s, term1, term2, term3, & term4, term5, z integer :: ipass, maxit ! Variable definitions (alphabetically) ! ARG: A value used in the approximation. ! B21: A parameter used in the approximation. ! B31: A parameter used in the approximation. ! B32: A parameter used in the approximation. ! B33: A parameter used in the approximation. ! B34: A parameter used in the approximation. ! B41: A parameter used in the approximation. ! B42: A parameter used in the approximation. ! B43: A parameter used in the approximation. ! B44: A parameter used in the approximation. ! B45: A parameter used in the approximation. ! B51: A parameter used in the approximation. ! B52: A parameter used in the approximation. ! B53: A parameter used in the approximation. ! B54: A parameter used in the approximation. ! B55: A parameter used in the approximation. ! B56: A parameter used in the approximation. ! C: A value used in the approximation. ! CON: A value used in the approximation. ! DF: The degrees of freedom. ! D1: A value used in the approximation. ! D3: A value used in the approximation. ! D5: A value used in the approximation. ! D7: A value used in the approximation. ! D9: A value used in the approximation. ! IPASS: A value used in the approximation. ! MAXIT: The maximum number of iterations allowed for the approx. ! PPFN: The normal percent point value. ! S: A value used in the approximation. ! TERM1: A value used in the approximation. ! TERM2: A value used in the approximation. ! TERM3: A value used in the approximation. ! TERM4: A value used in the approximation. ! TERM5: A value used in the approximation. ! Z: A value used in the approximation. df = idf maxit = 5 if (idf <= 0) then !Treat the IDF < 1 case res = zero elseif (idf == 1) then !Treat the IDF = 1 (Cauchy) case arg = pi*p res = -1/tan(arg) elseif (idf == 2) then ! Treat the IDF = 2 case term1 = sqrt(two)/2 term2 = 2*p - one term3 = sqrt(p*(one - p)) res = term1*term2/term3 elseif (idf >= 3) then ! Treat the IDF greater than or equal to 3 case ppfn = ppf_normal(p) d1 = ppfn d3 = ppfn**3 d5 = ppfn**5 d7 = ppfn**7 d9 = ppfn**9 term1 = d1 term2 = (one/b21)*(d3 + d1)/df term3 = (one/b31)*(b32*d5 + b33*d3 + b34*d1)/(df**2) term4 = (one/b41)*(b42*d7 + b43*d5 + b44*d3 + b45*d1)/(df**3) term5 = (one/b51)*(b52*d9 + b53*d7 + b54*d5 + b55*d3 + b56*d1)/(df**4) res = term1 + term2 + term3 + term4 + term5 if (idf == 3) then ! Augment the results for the IDF = 3 case con = pi*(p - half) arg = res/sqrt(df) z = atan(arg) do ipass = 1, maxit s = sin(z) c = cos(z) z = z - (z + s*c - con)/(2*c**2) end do res = sqrt(df)*s/c elseif (idf == 4) then ! Augment the results for the IDF = 4 case con = 2*(p - half) arg = res/sqrt(df) z = atan(arg) do ipass = 1, maxit s = sin(z) c = cos(z) z = z - ((one + half*c**2)*s - con)/((one + half)*c**3) end do res = sqrt(df)*s/c elseif (idf == 5) then ! Augment the results for the IDF = 5 case con = pi*(p - half) arg = res/sqrt(df) z = atan(arg) do ipass = 1, maxit s = sin(z) c = cos(z) z = z - (z + (c + (two/three)*c**3)*s - con)/((eight/three)*c**4) end do res = sqrt(df)*s/c elseif (idf == 6) then ! Augment the results for the IDF = 6 case con = 2*(p - half) arg = res/sqrt(df) z = atan(arg) do ipass = 1, maxit s = sin(z) c = cos(z) z = z - ((one + half*c**2 + (three/eight)*c**4)*s - con)/((fiftn/eight)*c**5) end do res = sqrt(df)*s/c end if end if end function ppf_tstudent subroutine fpvb & (fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp, & istop, nfev, pvb, & fjacd, f, fjacb) !! Compute the `nrow`-th function value using `beta(j) + stp`. procedure(fcn_t) :: fcn !! User-supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(in) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. integer, intent(in) :: nrow !! Row number of the independent variable array at which the derivative is to be checked. integer, intent(in) :: j !! Index of the partial derivative being examined. integer, intent(in) :: lq !! Response currently being examined. real(wp), intent(in) :: stp !! Step size for the finite difference derivative. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: pvb !! Function value for the selected observation & response. real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian wrt delta. real(wp), intent(out) :: f(n, q) !! Predicted function values. real(wp), intent(out) :: fjacb(n, np, q) !! Jocobian wrt beta. ! Local scalars real(wp) :: betaj ! Variable Definitions (alphabetically) ! BETAJ: The current estimate of the jth parameter. betaj = beta(j) beta(j) = beta(j) + stp istop = 0 call fcn(beta, xplusd, ifixb, ifixx, 003, f, fjacb, fjacd, istop) if (istop == 0) then nfev = nfev + 1 else return end if beta(j) = betaj pvb = f(nrow, lq) end subroutine fpvb subroutine fpvd & (fcn, & n, m, np, q, & beta, xplusd, ifixb, ifixx, ldifx, & nrow, j, lq, stp, & istop, nfev, pvd, & fjacd, f, fjacb) !! Compute `nrow`-th function value using `x(nrow, j) + delta(nrow, j) + stp`. procedure(fcn_t) :: fcn !! User-supplied subroutine for evaluating the model. integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. integer, intent(in) :: np !! Number of function parameters. integer, intent(in) :: q !! Number of responses per observation. real(wp), intent(in) :: beta(np) !! Function parameters. real(wp), intent(inout) :: xplusd(n, m) !! Values of `x + delta`. integer, intent(in) :: ifixb(np) !! Values designating whether the elements of `beta` are fixed at their input values or not. integer, intent(in) :: ifixx(ldifx, m) !! Values designating whether the elements of `x` are fixed at their input values or not. integer, intent(in) :: ldifx !! Leading dimension of array `ifixx`. integer, intent(in) :: nrow !! Row number of the independent variable array at which the derivative is to be checked. integer, intent(in) :: j !! Index of the partial derivative being examined. integer, intent(in) :: lq !! Response currently being examined. real(wp), intent(in) :: stp !! Step size for the finite difference derivative. integer, intent(out) :: istop !! Variable designating whether there are problems computing the function at the !! current `beta` and `delta`. integer, intent(inout) :: nfev !! Number of function evaluations. real(wp), intent(out) :: pvd !! Function value for the selected observation & response. real(wp), intent(out) :: fjacd(n, m, q) !! Jacobian wrt delta. real(wp), intent(out) :: f(n, q) !! Predicted function values. real(wp), intent(out) :: fjacb(n, np, q) !! Jocobian wrt beta. ! Local scalars real(wp) :: xpdj ! Variable Definitions (alphabetically) ! XPDJ: The (NROW,J)th element of XPLUSD. xpdj = xplusd(nrow, j) xplusd(nrow, j) = xplusd(nrow, j) + stp istop = 0 call fcn(beta, xplusd, ifixb, ifixx, 003, f, fjacb, fjacd, istop) if (istop == 0) then nfev = nfev + 1 else return end if xplusd(nrow, j) = xpdj pvd = f(nrow, lq) end subroutine fpvd pure subroutine scale_vec(n, m, scl, ldscl, t, ldt, sclt, ldsclt) !! Scale `t` by the inverse of `scl`, i.e., compute `t/scl`. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of rows of data in `t`. integer, intent(in) :: m !! Number of columns of data in `t`. real(wp), intent(in) :: scl(ldscl, m) !! Scale values. integer, intent(in) :: ldscl !! Leading dimension of array `scl`. real(wp), intent(in) :: t(ldt, m) !! Array to be inversely scaled by `scl`. integer, intent(in) :: ldt !! Leading dimension of array `t`. real(wp), intent(out) :: sclt(ldsclt, m) !! Inversely scaled matrix. integer, intent(in) :: ldsclt !! Leading dimension of array `sclt`. ! Local scalars integer :: j ! Variable Definitions (alphabetically) ! J: An indexing variable. if (n == 0 .or. m == 0) return if (scl(1, 1) >= zero) then if (ldscl >= n) then sclt(1:n, 1:m) = t(1:n, 1:m)/scl(1:n, 1:m) else do j = 1, m sclt(1:n, j) = t(1:n, j)/scl(1, j) end do end if else sclt(1:n, 1:m) = t(1:n, 1:m)/abs(scl(1, 1)) end if end subroutine scale_vec pure subroutine scale_beta(np, beta, ssf) !! Select scaling values for `beta` according to the algorithm given in the ODRPACK95 !! reference guide. use odrpack_kinds, only: zero, one integer, intent(in) :: np !! Number of function parameters. real(wp), intent(in) :: beta(np) !! Function parameters. real(wp), intent(out) :: ssf(np) !! Scaling values for `beta`. ! Local scalars real(wp) :: bmax, bmin integer :: k logical ::bigdif ! Variable Definitions (alphabetically) ! BIGDIF: The variable designating whether there is a significant difference in the ! magnitudes of the nonzero elements of BETA (TRUE) or not (FALSE). ! BMAX: The largest nonzero magnitude. ! BMIN: The smallest nonzero magnitude. ! K: An indexing variable. bmax = abs(beta(1)) do k = 2, np bmax = max(bmax, abs(beta(k))) end do if (bmax == zero) then ! All input values of BETA are zero ssf(1:np) = one else ! Some of the input values are nonzero bmin = bmax do k = 1, np if (beta(k) /= zero) then bmin = min(bmin, abs(beta(k))) end if end do bigdif = log10(bmax) - log10(bmin) >= one do k = 1, np if (beta(k) == zero) then ssf(k) = 10/bmin else if (bigdif) then ssf(k) = 1/abs(beta(k)) else ssf(k) = 1/bmax end if end if end do end if end subroutine scale_beta pure subroutine scale_delta(n, m, x, tt, ldtt) !! Select scaling values for `delta` according to the algorithm given in the ODRPACK95 !! reference guide. use odrpack_kinds, only: zero, one integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. real(wp), intent(in) :: x(n, m) !! Independent variable. real(wp), intent(out) :: tt(ldtt, m) !! Scaling values for `delta`. integer, intent(in) :: ldtt !! Leading dimension of array `tt`. ! Local scalars real(wp) :: xmax, xmin integer :: i, j logical :: bigdif ! Variable Definitions (alphabetically) ! BIGDIF: The variable designating whether there is a significant difference in the ! magnitudes of the nonzero elements of X (TRUE) or not (FALSE). ! I: An indexing variable. ! J: An indexing variable. ! XMAX: The largest nonzero magnitude. ! XMIN: THE SMALLEST NONZERO MAGNITUDE. do j = 1, m xmax = abs(x(1, j)) do i = 2, n xmax = max(xmax, abs(x(i, j))) end do if (xmax == zero) then ! All input values of X(I,J), I=1,...,N, are zero tt(1:n, j) = one else ! Some of the input values are nonzero xmin = xmax do i = 1, n if (x(i, j) /= zero) then xmin = min(xmin, abs(x(i, j))) end if end do bigdif = log10(xmax) - log10(xmin) >= one do i = 1, n if (x(i, j) /= zero) then if (bigdif) then tt(i, j) = 1/abs(x(i, j)) else tt(i, j) = 1/xmax end if else tt(i, j) = 10/xmin end if end do end if end do end subroutine scale_delta pure subroutine select_row(n, m, x, nrow) !! Select the row at which the derivative will be checked. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of observations. integer, intent(in) :: m !! Number of columns of data in the independent variable. real(wp), intent(in) :: x(n, m) !! Independent variable. integer, intent(inout) :: nrow !! Selected row number of the independent variable. ! Local scalars integer :: i ! Variable Definitions (alphabetically) ! I: An index variable. ! J: An index variable. ! M: The number of columns of data in the independent variable. ! N: The number of observations. ! NROW: The selected row number of the independent variable. ! X: The independent variable. if ((nrow >= 1) .and. (nrow <= n)) return ! Select first row of independent variables which contains no zeros ! if there is one, otherwise first row is used. nrow = 1 do i = 1, n if (all(x(i, :) /= zero)) then nrow = i return end if end do end subroutine select_row pure subroutine solve_trl(n, t, ldt, b, job) !! Solve systems of the form: !! !! `t * x = b or trans(t) * x = b` !! !! where `t` is an upper or lower triangular matrix of order `n`, and the solution `x` !! overwrites the RHS `b`. ! Adapted from LINPACK subroutine DTRSL. ! @note: This is one of the most time-consuming subroutines in ODRPACK (~25% of total). use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of rows and columns of data in array `t`. real(wp), intent(in) :: t(ldt, n) !! Upper or lower tridiagonal system. integer, intent(in) :: ldt !! Leading dimension of array `t`. real(wp), intent(inout) :: b(:) !! On input: the right hand side; On exit: the solution. integer, intent(in) :: job !! What kind of system is to be solved: !! 1: Solve `t * x = b`, where `t` is lower triangular, !! 2: Solve `t * x = b`, where `t` is upper triangular, !! 3: Solve `trans(t) * x = b`, where `t` is lower triangular, !! 4: Solve `trans(t) * x = b`, where `t` is upper triangular. ! Local scalars real(wp) :: temp integer :: j1, j, jn ! Variable Definitions (alphabetically) ! B: On input: the right hand side; On exit: the solution ! J1: The first nonzero entry in T. ! J: An indexing variable. ! JN: The last nonzero entry in T. ! JOB: What kind of system is to be solved, where if JOB is ! 1 Solve T*X=B, T lower triangular, ! 2 Solve T*X=B, T upper triangular, ! 3 Solve trans(T)*X=B, T lower triangular, ! 4 Solve trans(T)*X=B, T upper triangular. ! LDT: The leading dimension of array T. ! N: The number of rows and columns of data in array T. ! T: The upper or lower tridiagonal system. ! Find first nonzero diagonal entry in T j1 = 0 do j = 1, n if (j1 == 0 .and. t(j, j) /= zero) then j1 = j elseif (t(j, j) == zero) then b(j) = zero end if end do if (j1 == 0) return ! Find last nonzero diagonal entry in T jn = 0 do j = n, j1, -1 if (jn == 0 .and. t(j, j) /= zero) then jn = j elseif (t(j, j) == zero) then b(j) = zero end if end do select case (job) case (1) ! Solve T*X=B for T lower triangular b(j1) = b(j1)/t(j1, j1) do j = j1 + 1, jn temp = -b(j - 1) b(j:jn) = b(j:jn) + temp*t(j:jn, j - 1) if (t(j, j) /= zero) then b(j) = b(j)/t(j, j) else b(j) = zero end if end do case (2) ! Solve T*X=B for T upper triangular. b(jn) = b(jn)/t(jn, jn) do j = jn - 1, j1, -1 temp = -b(j + 1) b(1:j) = b(1:j) + temp*t(1:j, j + 1) if (t(j, j) /= zero) then b(j) = b(j)/t(j, j) else b(j) = zero end if end do case (3) ! Solve trans(T)*X=B for T lower triangular. b(jn) = b(jn)/t(jn, jn) do j = jn - 1, j1, -1 b(j) = b(j) - dot_product(t(j + 1:jn + 1, j), b(j + 1:jn + 1)) if (t(j, j) /= zero) then b(j) = b(j)/t(j, j) else b(j) = zero end if end do case (4) ! Solve trans(T)*X=B for T upper triangular b(j1) = b(j1)/t(j1, j1) do j = j1 + 1, jn b(j) = b(j) - dot_product(t(1:j - 1, j), b(1:j - 1)) if (t(j, j) /= zero) then b(j) = b(j)/t(j, j) else b(j) = zero end if end do case default error stop "Invalid value of JOB." end select end subroutine solve_trl pure subroutine vevtr & (m, q, indx, v, ldv, ld2v, e, lde, ve, ldve, ld2ve, vev, ldvev, wrk5) !! Compute `v*e*trans(v)` for the (`indx`)th `m` by `q` array in `v`. integer, intent(in) :: m !! Number of columns of data in the independent variable. integer, intent(in) :: q !! Number of responses per observation. integer, intent(in) :: indx !! Row in `v` in which the `m` by `q` array is stored. integer, intent(in) :: ldv !! Leading dimension of array `v`. integer, intent(in) :: ld2v !! Second dimension of array `v`. integer, intent(in) :: lde !! Leading dimension of array `e`. integer, intent(in) :: ldve !! Leading dimension of array `ve`. integer, intent(in) :: ldvev !! Leading dimension of array `vev`. integer, intent(in) :: ld2ve !! Second dimension of array `ve`. real(wp), intent(in) :: v(ldv, ld2v, q) !! An array of `q` by `m` matrices. real(wp), intent(in) :: e(lde, m) !! Matrix of the factors, so `ete = (d**2 + alpha*t**2)`. real(wp), intent(out) :: ve(ldve, ld2ve, m) !! Array `ve = v * inv(e)`. real(wp), intent(out) :: vev(ldvev, q) !! Array `vev = v * inv(ete) * trans(v)`. real(wp), intent(out) :: wrk5(m) !! Work vector. ! Local scalars integer :: l1, l2 ! Variable Definitions (alphabetically) ! J: An indexing variable. ! L1: An indexing variable. ! L2: An indexing variable. if (q == 0 .or. m == 0) return do l1 = 1, q wrk5 = v(indx, 1:m, l1) call solve_trl(m, e, lde, wrk5, 4) ve(indx, l1, :) = wrk5 end do do l1 = 1, q do l2 = 1, l1 vev(l1, l2) = dot_product(ve(indx, l1, :), ve(indx, l2, :)) vev(l2, l1) = vev(l1, l2) end do end do end subroutine vevtr pure subroutine scale_mat(n, m, wt, ldwt, ld2wt, t, wtt) !! Scale matrix `t` using `wt`, i.e., compute `wtt = wt*t`. use odrpack_kinds, only: zero integer, intent(in) :: n !! Number of rows of data in `t`. integer, intent(in) :: m !! Number of columns of data in `t`. real(wp), intent(in), target :: wt(..) !! Array of shape conformable to `(ldwt,ld2wt,m)` holding the weights. integer, intent(in) :: ldwt !! Leading dimension of array `wt`. integer, intent(in) :: ld2wt !! Second dimension of array `wt`. real(wp), intent(in), target :: t(..) !! Array of shape conformable to `(n,m)` being scaled by `wt`. real(wp), intent(out), target :: wtt(..) !! Array of shape conformable to `(n,m)` holding the result of weighting array `t` by !! array `wt`. Array `wtt` can be the same as `t` only if the arrays in `wt` are upper !! triangular with zeros below the diagonal. ! Local scalars integer :: i, j real(wp), pointer :: wt_(:, :, :), t_(:, :), wtt_(:, :) ! Variable Definitions (alphabetically) ! I: An indexing variable. ! J: An indexing variable. if (n == 0 .or. m == 0) return select rank (wt) rank (1); wt_(1:ldwt, 1:ld2wt, 1:m) => wt rank (2); wt_(1:ldwt, 1:ld2wt, 1:m) => wt rank (3); wt_(1:ldwt, 1:ld2wt, 1:m) => wt rank default; error stop "Invalid rank of `wt`." end select select rank (t) rank (1); t_(1:n, 1:m) => t rank (2); t_(1:n, 1:m) => t rank default; error stop "Invalid rank of `t`." end select select rank (wtt) rank (1); wtt_(1:n, 1:m) => wtt rank (2); wtt_(1:n, 1:m) => wtt rank default; error stop "Invalid rank of `wtt`." end select if (wt_(1, 1, 1) >= zero) then if (ldwt >= n) then if (ld2wt >= m) then ! WT is an N-array of M by M matrices do j = 1, m do i = 1, n wtt_(i, j) = dot_product(wt_(i, j, :), t_(i, :)) end do end do else ! WT is an N-array of diagonal matrices do j = 1, m wtt_(:, j) = wt_(:, 1, j)*t_(:, j) end do end if else if (ld2wt >= m) then ! WT is an M by M matrix do j = 1, m do i = 1, n wtt_(i, j) = dot_product(wt_(1, j, :), t_(i, :)) end do end do else ! WT is a diagonal matrix do j = 1, m wtt_(:, j) = wt_(1, 1, j)*t_(:, j) end do end if end if else ! WT is a scalar wtt_ = abs(wt_(1, 1, 1))*t_ end if end subroutine scale_mat pure subroutine move_beta( & np, beta, lower, upper, ssf, stpb, neta, eta, interval) !! Ensure range of bounds is large enough for derivative checking. !! Move beta away from bounds so that derivatives can be calculated. use odrpack_kinds, only: zero, one, three, ten, hundred integer, intent(in) :: np !! Number of parameters `np`. real(wp), intent(inout) :: beta(np) !! Function parameters. real(wp), intent(in) :: lower(np) !! !! Lower bound on `beta`. real(wp), intent(in) :: upper(np) !! Upper bound on `beta`. real(wp), intent(in) :: ssf(np) !! Scale used for the `beta`s. real(wp), intent(in) :: stpb(np) !! Relative step used for computing finite difference derivatives with respect to `beta`. integer, intent(in) :: neta !! Number of good digits in the function results. real(wp), intent(in) :: eta !! Relative noise in the function results. integer, intent(out) :: interval(np) !! Specifies which difference methods and step sizes are supported by the current !! interval `upper-lower`. ! Local scalars integer :: k real(wp) :: h, h0, h1, hc, hc0, hc1, stpl, stpr, typj ! VARIABLE DEFINITIONS (ALPHABETICALLY) ! H: Relative step size for forward differences. ! H0: Initial relative step size for forward differences. ! H1: Default relative step size for forward differences. ! HC: Relative step size for center differences. ! HC0: Initial relative step size for center differences. ! HC1: Default relative step size for center differences. ! K: Index variable for BETA. ! STPL: Maximum step to the left of BETA (-) the derivative checker will use. ! STPR: Maximum step to the right of BETA (+) the derivative checker will use. ! TYPJ: The typical size of the J-th unkonwn BETA. interval = 111 do k = 1, np h0 = hstep(0, neta, 1, k, stpb, 1) hc0 = h0 h1 = sqrt(eta) hc1 = eta**(one/three) h = max(ten*h1, min(hundred*h0, one)) hc = max(ten*hc1, min(hundred*hc0, one)) if (beta(k) == zero) then if (ssf(1) < zero) then typj = one/abs(ssf(1)) else typj = one/ssf(k) end if else typj = abs(beta(k)) end if stpr = (h*typj*sign(one, beta(k)) + beta(k)) - beta(k) stpl = (hc*typj*sign(one, beta(k)) + beta(k)) - beta(k) ! Check outer interval if (lower(k) + 2*abs(stpl) > upper(k)) then if (interval(k) >= 100) then interval(k) = interval(k) - 100 end if elseif (beta(k) + stpl > upper(k) .or. beta(k) - stpl > upper(k)) then beta(k) = upper(k) - abs(stpl) elseif (beta(k) + stpl < lower(k) .or. beta(k) - stpl < lower(k)) then beta(k) = lower(k) + abs(stpl) end if ! Check middle interval if (lower(k) + 2*abs(stpr) > upper(k)) then if (mod(interval(k), 100) >= 10) then interval(k) = interval(k) - 10 end if elseif (beta(k) + stpr > upper(k) .or. beta(k) - stpr > upper(k)) then beta(k) = upper(k) - abs(stpr) elseif (beta(k) + stpr < lower(k) .or. beta(k) - stpr < lower(k)) then beta(k) = lower(k) + abs(stpr) end if ! Check inner interval if (lower(k) + abs(stpr) > upper(k)) then interval(k) = 0 elseif (beta(k) + stpr > upper(k)) then beta(k) = upper(k) - stpr elseif (beta(k) + stpr < lower(k)) then beta(k) = lower(k) - stpr end if end do end subroutine move_beta end module odrpack_core