odr Subroutine

   impure subroutine odr &
      (fcn, &
       n, m, np, nq, &
       beta, &
       y, x, &
       delta, &
       we, wd, &
       ifixb, ifixx, &
       job, ndigit, taufac, &
       sstol, partol, maxit, &
       iprint, lunerr, lunrpt, &
       stpb, stpd, &
       sclb, scld, &
       work, iwork, &
       info, &
       lower, upper)
   !! Driver routine for finding the weighted explicit or implicit orthogonal distance
   !! regression (ODR) or ordinary linear or nonlinear least squares (OLS) solution.
      ! Date Written:   860529   (YYMMDD)
      ! Revision Date:  20040301 (YYYYMMDD)
      ! Category No.:  G2E,I1B1
      ! Keywords:  Orthogonal distance regression,
      !   Nonlinear least squares,
      !   Measurement error models,
      !   Errors in variables
      ! Authors:
      !   Boggs, Paul T.
      !     Applied and Computational Mathematics Division
      !     National Institute of Standards and Technology
      !     Gaithersburg, MD 20899
      !   Byrd, Richard H.
      !     Department of Computer Science
      !     University of Colorado, Boulder, CO 80309
      !   Rogers, Janet E.
      !     Applied and Computational Mathematics Division
      !     National Institute of Standards and Technology
      !     Boulder, CO 80303-3328
      !   Schnabel, Robert B.
      !     Department of Computer Science
      !     University of Colorado, Boulder, CO 80309
      !     and
      !     Applied and Computational Mathematics Division
      !     National Institute of Standards and Technology
      !     Boulder, CO 80303-3328
      ! Purpose:  REAL(wp) driver routine for finding the weighted explicit or implicit
      !   orthogonal distance regression (ODR) or ordinary linear or nonlinear least squares (OLS)
      !   solution (long call statement)
      ! Description:
      !   For details, see ODRPACK95 User's Reference Guide.
      ! References:
      !   Boggs, P. T., R. H. Byrd, J. R. Donaldson, and R. B. Schnabel (1989),
      !     "Algorithm 676 --- ODRPACK: Software for Weighted Orthogonal Distance Regression,"
      !     ACM Trans. Math. Software., 15(4):348-364.
      !       Boggs, P. T., R. H. Byrd, J. E. Rogers, and
      !   R. B. Schnabel (1992),
      !     "User's Reference Guide for ODRPACK Version 2.01,
      !     Software for Weighted Orthogonal Distance Regression,"
      !     National Institute of Standards and Technology
      !     Internal Report Number 92-4834.
      !  Boggs, P. T., R. H. Byrd, and R. B. Schnabel (1987),
      !    "A Stable and Efficient Algorithm for Nonlinear Orthogonal Distance Regression,"
      !    SIAM J. Sci. Stat. Comput., 8(6):1052-1078.

      use odrpack_kinds, only: negone, zero
      use odrpack_core, only: fcn_t
      use odrpack_reports, only: dodphd, dodpe1

      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) :: nq
         !! Number of responses per observation.
      real(wp), intent(inout) :: beta(:)
         !! Function parameters. `Shape: (np)`.
      real(wp), intent(in) :: y(:, :)
         !! Dependent variable. `Shape: (n, nq)`. Unused when the model is implicit.
      real(wp), intent(in) :: x(:, :)
         !! Explanatory variable. `Shape: (n, m)`.
      real(wp), intent(inout), optional :: delta(:, :)
         !! Error in the `x` data. `Shape: (n, m)`. Initial guess on input and estimated value
         !! on output.
      real(wp), intent(in), optional :: we(:, :, :)
         !! `epsilon` weights. `Shape: (1<=ldwe<=n, 1<=ld2we<=nq, nq)`. See p. 25.
      real(wp), intent(in), optional :: wd(:, :, :)
         !! `delta` weights. `Shape: (1<=ldwd<=n, 1<=ld2wd<=m, m)`. See p. 26.
      integer, intent(in), optional :: ifixb(:)
         !! Values designating whether the elements of `beta` are fixed at their input values
         !! or not. `Shape: (np)`.
      integer, intent(in), optional :: ifixx(:, :)
         !! Values designating whether the elements of `x` are fixed at their input values
         !! or not. `Shape: (1<=ldifx<=n, m)`. See p. 27.
      integer, intent(in), optional :: job
         !! Variable controlling problem initialization and computational method.
      integer, intent(in), optional :: ndigit
         !! Number of accurate digits in the function results, as supplied by the user.
      real(wp), intent(in), optional :: taufac
         !! Factor used to compute the initial trust region diameter.
      real(wp), intent(in), optional :: sstol
         !! Sum-of-squares convergence stopping tolerance.
      real(wp), intent(in), optional :: partol
      !! Parameter convergence stopping tolerance.
      integer, intent(in), optional :: maxit
         !! Maximum number of iterations allowed.
      integer, intent(in), optional :: iprint
         !! Print control variable.
      integer, intent(in), optional :: lunerr
         !! Logical unit number for error messages. Available options are:
         !!   0 => no output. 
         !!   6 => output to standard error.
         !!   other => output to logical unit number `lunerr`.  
      integer, intent(in), optional :: lunrpt
         !! Logical unit number for computation reports. Available options are:
         !!   0 => no output. 
         !!   6 => output to standard error.
         !!   other => output to logical unit number `lunrpt`.  
      real(wp), intent(in), optional :: stpb(:)
         !! Relative step for computing finite difference derivatives with respect to `beta`.
         !! `Shape: (np)`.
      real(wp), intent(in), optional :: stpd(:, :)
         !! Relative step for computing finite difference derivatives with respect to `delta`.
         !! `Shape: (1<=ldstpd<=n, m)`. See p. 31.
      real(wp), intent(in), optional :: sclb(:)
      !! Scaling values for `beta`. `Shape: (np)`.
      real(wp), intent(in), optional :: scld(:, :)
         !! Scaling values for `delta`. `Shape: (1<=ldscld<=n, m)`. See p. 32.
      real(wp), intent(inout), optional, target :: work(:)
         !! Real work space.
      integer, intent(inout), optional, target :: iwork(:)
         !! Integer work space.
      integer, intent(out), optional :: info
         !! Variable designating why the computations were stopped.
      real(wp), intent(in), optional :: lower(:)
         !! Lower bound on `beta`. `Shape: (np)`.
      real(wp), intent(in), optional :: upper(:)
         !! Upper bound on `beta`. `Shape: (np)`.

      ! Local variables
      integer :: ldwe, ld2we, ldwd, ld2wd, ldifx, ldscld, ldstpd, job_, ndigit_, maxit_, &
                 iprint_, lunerr_, lunrpt_, info_, lwork, liwork, info1_, info2_, info3_, &
                 info4_, info5_
      integer :: ifixb_(np), ifixx_(n, m)
      real(wp) :: taufac_, sstol_, partol_
      real(wp) :: lower_(np), we_(n, nq, nq), wd_(n, m, m), stpb_(np), stpd_(n, m), &
                  sclb_(np), scld_(n, m), upper_(np), wd1(1, 1, 1)
      real(wp), allocatable :: tempret(:, :)

      real(wp), allocatable, target :: work_local(:)
      real(wp), pointer :: work_(:)
      integer, allocatable, target :: iwork_local(:)
      integer, pointer :: iwork_(:)

      logical :: head, isodr, restart

      ! Set LINFO to zero indicating no errors have been found thus far
      info_ = 0
      info1_ = 0
      info2_ = 0
      info3_ = 0
      info4_ = 0
      info5_ = 0

      ! Set all scalar variable defaults except JOB
      ldwe = 1
      ld2we = 1
      ldwd = 1
      ld2wd = 1
      ldifx = 1
      ldscld = 1
      ldstpd = 1
      iprint_ = -1
      lunerr_ = 6
      lunrpt_ = 6
      maxit_ = -1
      ndigit_ = -1
      partol_ = negone
      sstol_ = negone
      taufac_ = negone
      head = .true.

      !  Check for the option arguments for printing (so error messages can be
      !  printed appropriately from here on out
      if (present(iprint)) then
         iprint_ = iprint
      end if

      if (present(lunrpt)) then
         lunrpt_ = lunrpt
      end if
      if (lunrpt_ == 6) then
         lunrpt_ = output_unit
      end if

      if (present(lunerr)) then
         lunerr_ = lunerr
      end if
      if (lunerr_ == 6) then
         lunerr_ = error_unit
      end if

      ! Ensure the problem size is valid
      if (n <= 0) then
         info5_ = 1
         info4_ = 1
      end if

      if (m <= 0) then
         info5_ = 1
         info3_ = 1
      end if

      if (np <= 0) then
         info5_ = 1
         info2_ = 1
      end if

      if (nq <= 0) then
         info5_ = 1
         info1_ = 1
      end if

      if (info5_ /= 0) then
         info_ = 10000*info5_ + 1000*info4_ + 100*info3_ + 10*info2_ + info1_
         if (lunerr_ /= 0 .and. iprint_ /= 0) then
            call dodphd(head, lunrpt_)
            call dodpe1( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, nq, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lwork, liwork)
         end if
         if (present(info)) then
            info = info_
         end if
         return
      end if

      ! Define LJOB and check that necessary arguments are passed for JOB
      if (present(job)) then
         job_ = job
         if (mod(job, 10000)/1000 >= 1) then
            if (.not. present(delta)) then
               info5_ = 7
               info4_ = 1
            end if
         end if
         if (job >= 10000) then
            if (.not. present(iwork)) then
               info5_ = 7
               info2_ = 1
            end if
            if (.not. present(work)) then
               info5_ = 7
               info3_ = 1
            end if
         end if
      else
         job_ = -1
      end if

      if (info5_ /= 0) then
         info_ = 10000*info5_ + 1000*info4_ + 100*info3_ + 10*info2_ + info1_
         if (lunerr_ /= 0 .and. iprint_ /= 0) then
            call dodphd(head, lunrpt_)
            call dodpe1( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, nq, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lwork, liwork)
         end if
         if (present(info)) then
            info = info_
         end if
         return
      end if

      ! Determine the size of the work arrays
      isodr = (job_ < 0 .or. mod(job_, 10) <= 1)
      call workspace_dimensions(n, m, np, nq, isodr, lwork, liwork)
      
      ! Allocate the local work arrays, if not provided by user
      ! The complementary case is treated below, because of the way the info flags are designed 
      if (.not. present(iwork)) then
         allocate (iwork_local(liwork), stat=info2_)
         iwork_ => iwork_local
      end if

      if (.not. present(work)) then
         allocate (work_local(lwork), stat=info3_)
         work_ => work_local
      end if
     
      allocate (tempret(max(n, np), max(nq, m)), stat=info4_)

      if (info4_ /= 0 .or. info3_ /= 0 .or. info2_ /= 0) then
         info5_ = 8
      end if

      if (info5_ /= 0) then
         info_ = 10000*mod(info5_, 10) + 1000*mod(info4_, 10) + &
                 100*mod(info3_, 10) + 10*mod(info2_, 10) + mod(info1_, 10)
         if (lunerr_ /= 0 .and. iprint_ /= 0) then
            call dodphd(head, lunrpt_)
            call dodpe1( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, nq, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lwork, liwork)
         end if
         if (present(info)) then
            info = info_
         end if
         return
      end if

      ! Set internal array variable defaults except IWORK and WORK
      ifixb_(1) = -1
      ifixx_(1, 1) = -1
      lower_(1:np) = -huge(zero)
      sclb_(1) = negone
      scld_(1, 1) = negone
      stpb_(1) = negone
      stpd_(1, 1) = negone
      upper_(1:np) = huge(zero)
      we_(1, 1, 1) = negone
      wd_(1, 1, 1) = negone

      ! Check the size of required arguments and return errors if they are too small
      if (size(beta) < np) then
         info1_ = info1_ + 1
      end if

      if (any(size(y) < [n, nq])) then
         info1_ = info1_ + 2
      end if

      if (any(size(x) < [n, m])) then
         info1_ = info1_ + 4
      end if

      ! Check the presence of optional arguments and copy their values internally or
      ! report errors as necessary
      if (present(ifixb)) then
         if (size(ifixb) < np) then
            info1_ = info1_ + 64
         end if
         if (ifixb(1) < 0) then
            ifixb_(1) = ifixb(1)
         else
            ifixb_(1:np) = ifixb(1:np)
         end if
      end if

      if (present(ifixx)) then
         ldifx = size(ifixx, 1)
         if (any(size(ifixx) <= [0, 0])) then
            info1_ = info1_ + 128
         end if
         if (.not. (ifixx(1, 1) < 0 .or. ldifx == 1 .or. ldifx >= n) &
             .or. size(ifixx, 2) < m) then
            info1_ = info1_ + 128
         end if
         if (ldifx > n) then
            ldifx = n
         end if
         if (ifixx(1, 1) < 0) then
            ifixx_(1, 1) = ifixx(1, 1)
         else
            ifixx_(1:ldifx, 1:m) = ifixx(1:ldifx, 1:m)
         end if
      end if
    
      if (present(iwork)) then
         if (size(iwork) >= liwork) then
            iwork_ => iwork(1:liwork)
         else
             info1_ = info1_ + 8192           
         end if
      end if

      if (present(maxit)) then
         maxit_ = maxit
      end if

      if (present(ndigit)) then
         ndigit_ = ndigit
      end if

      if (present(partol)) then
         partol_ = partol
      end if

      if (present(sclb)) then
         if (size(sclb) < np) then
            info1_ = info1_ + 1024
         end if
         if (sclb(1) <= zero) then
            sclb_(1) = sclb(1)
         else
            sclb_(1:np) = sclb(1:np)
         end if
      end if

      if (present(scld)) then
         ldscld = size(scld, 1)
         if (any(size(scld) <= [0, 0])) then
            info1_ = info1_ + 2048
         end if
         if (.not. (scld(1, 1) <= zero .or. ldscld == 1 .or. ldscld >= n) &
             .or. size(scld, 2) < m) then
            info1_ = info1_ + 2048
         end if
         if (ldscld > n) then
            ldscld = n
         end if
         if (scld(1, 1) <= zero) then
            scld_(1, 1) = scld(1, 1)
         else
            scld_(1:ldscld, 1:m) = scld(1:ldscld, 1:m)
         end if
      end if

      if (present(sstol)) then
         sstol_ = sstol
      end if

      if (present(stpb)) then
         if (size(stpb) < np) then
            info1_ = info1_ + 256
         end if
         if (stpb(1) <= zero) then
            stpb_(1) = stpb(1)
         else
            stpb_(1:np) = stpb(1:np)
         end if
      end if

      if (present(stpd)) then
         ldstpd = size(stpd, 1)
         if (any(size(stpd) <= [0, 0])) then
            info1_ = info1_ + 512
         end if
         if (.not. (stpd(1, 1) <= zero .or. ldstpd == 1 .or. ldstpd >= n) &
             .or. size(stpd, 2) < m) then
            info1_ = info1_ + 512
         end if
         if (ldstpd > n) then
            ldstpd = n
         end if
         if (stpd(1, 1) <= zero) then
            stpd_(1, 1) = stpd(1, 1)
         else
            stpd_(1:ldstpd, 1:m) = stpd(1:ldstpd, 1:m)
         end if
      end if

      if (present(taufac)) then
         taufac_ = taufac
      end if

      if (present(we)) then
         ldwe = size(we, 1)
         ld2we = size(we, 2)
         if (any(size(we) <= [0, 0, 0])) then
            info1_ = info1_ + 16
         end if
         if (.not. (we(1, 1, 1) < zero .or. &
                    ((ldwe == 1 .or. ldwe >= n) .and. &
                     (ld2we == 1 .or. ld2we >= nq))) .or. size(we, 3) < nq) then
            info1_ = info1_ + 16
         end if
         if (ldwe > n) then
            ldwe = n
         end if
         if (ld2we > nq) then
            ld2we = nq
         end if
         if (we(1, 1, 1) < zero) then
            we_(1, 1, 1) = we(1, 1, 1)
         else
            we_(1:ldwe, 1:ld2we, 1:nq) = we(1:ldwe, 1:ld2we, 1:nq)
         end if
      end if

      if (present(wd)) then
         ldwd = size(wd, 1)
         ld2wd = size(wd, 2)
         if (any(size(wd) <= [0, 0, 0])) then
            info1_ = info1_ + 32
         end if
         if (.not. (wd(1, 1, 1) < zero .or. &
                    ((ldwd == 1 .or. ldwd >= n) .and. &
                     (ld2wd == 1 .or. ld2wd >= m))) .or. size(wd, 3) < m) then
            info1_ = info1_ + 32
         end if
         if (ldwd > n) then
            ldwd = n
         end if
         if (ld2wd > m) then
            ld2wd = m
         end if
         if (wd(1, 1, 1) <= zero) then
            wd_(1, 1, 1) = wd(1, 1, 1)
         else
            wd_(1:ldwd, 1:ld2wd, 1:m) = wd(1:ldwd, 1:ld2wd, 1:m)
         end if
      end if

      if (present(work)) then
         if (size(work) >= lwork) then
            work_ => work(1:lwork)
         else
            info1_ = info1_ + 4096
         end if
      end if

      restart = (mod(job_/10000, 10) >= 1)
      if (.not. restart) then
         work_ = zero
         iwork_ = 0
      end if

      if (present(delta) .and. (.not. restart)) then
         if (any(shape(delta) < [n, m])) then
            info1_ = info1_ + 8
         end if
         work_(1:n*m) = reshape(delta(1:n, 1:m), [n*m])
      end if

      if (present(lower)) then
         if (size(lower) < np) then
            info1_ = info1_ + 32768
         end if
         lower_(1:np) = lower(1:np)
      end if

      if (present(upper)) then
         if (size(upper) < np) then
            info1_ = info1_ + 16384
         end if
         upper_(1:np) = upper(1:np)
      end if

      ! Report an error if any of the array sizes didn't match.
      if (info1_ /= 0) then
         info_ = 100000 + info1_
         info1_ = 0
         if (lunerr_ /= 0 .and. iprint_ /= 0) then
            call dodphd(head, lunrpt_)
            call dodpe1( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, nq, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lwork, liwork)
         end if
         if (present(info)) then
            info = info_
         end if
         return
      end if

      if (wd_(1, 1, 1) /= 0) then
         call dodcnt &
            (fcn, &
             n, m, np, nq, &
             beta(1:np), &
             y(1:n, 1:nq), n, x(1:n, 1:m), n, &
             we_(1:ldwe, 1:ld2we, 1:nq), ldwe, ld2we, &
             wd_(1:ldwd, 1:ld2wd, 1:m), ldwd, ld2wd, &
             ifixb_, ifixx_(1:ldifx, 1:m), ldifx, &
             job_, ndigit_, taufac_, &
             sstol_, partol_, maxit_, &
             iprint_, lunerr_, lunrpt_, &
             stpb_, stpd_(1:ldstpd, 1:m), ldstpd, &
             sclb_, scld_(1:ldscld, 1:m), ldscld, &
             work_, lwork, tempret, iwork_, liwork, &
             info_, &
             lower_, upper_)
      else
         wd1(1, 1, 1) = negone
         call dodcnt &
            (fcn, &
             n, m, np, nq, &
             beta(1:np), &
             y(1:n, 1:nq), n, x(1:n, 1:m), n, &
             we_(1:ldwe, 1:ld2we, 1:nq), ldwe, ld2we, &
             wd1, 1, 1, &
             ifixb_, ifixx_(1:ldifx, 1:m), ldifx, &
             job_, ndigit_, taufac_, &
             sstol_, partol_, maxit_, &
             iprint_, lunerr_, lunrpt_, &
             stpb_, stpd_(1:ldstpd, 1:m), ldstpd, &
             sclb_, scld_(1:ldscld, 1:m), ldscld, &
             work_, lwork, tempret, iwork_, liwork, &
             info_, &
             lower_, upper_)
      end if

      if (present(delta)) then
         delta(1:n, 1:m) = reshape(work_(1:n*m), [n, m])
      end if

      if (present(info)) then
         info = info_
      end if

   end subroutine odr