odr Subroutine

   impure subroutine odr &
      (fcn, &
       n, m, q, np, &
       beta, &
       y, x, &
       delta, &
       we, wd, &
       ifixb, ifixx, &
       job, ndigit, taufac, &
       sstol, partol, maxit, &
       iprint, lunerr, lunrpt, &
       stpb, stpd, &
       sclb, scld, &
       rwork, 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.
      ! 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: print_header, print_error_inputs

      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) :: q
         !! Number of responses per observation.
      integer, intent(in) :: np
         !! Number of function parameters.
      real(wp), intent(inout) :: beta(:)
         !! Function parameters. `Shape: (np)`.
      real(wp), intent(in) :: y(:, :)
         !! Dependent variable. `Shape: (n, q)`. 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, target :: we(:, :, :)
         !! `epsilon` weights. `Shape: ({1,n}, {1,q}, q)`. See p. 25.
      real(wp), intent(in), optional, target :: wd(:, :, :)
         !! `delta` weights. `Shape: ({1,n}, {1,m}, m)`. See p. 26.
      integer, intent(in), optional, target :: ifixb(:)
         !! Values designating whether the elements of `beta` are fixed at their input values
         !! or not. `Shape: (np)`.
      integer, intent(in), optional, target :: ifixx(:, :)
         !! Values designating whether the elements of `x` are fixed at their input values
         !! or not. `Shape: ({1,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 output.
         !!  `k /= 0,6`: Output to logical unit `k`.
      integer, intent(in), optional :: lunrpt
         !! Logical unit number for computation reports. Available options are:
         !!  `0`: No output.
         !!  `6`: Output to standard output.
         !!  `k /= 0,6`: Output to logical unit `k`.
      real(wp), intent(in), optional, target :: stpb(:)
         !! Relative step for computing finite difference derivatives with respect to `beta`.
         !! `Shape: (np)`.
      real(wp), intent(in), optional, target :: stpd(:, :)
         !! Relative step for computing finite difference derivatives with respect to `delta`.
         !! `Shape: ({1,n}, m)`. See p. 31.
      real(wp), intent(in), optional, target :: sclb(:)
      !! Scaling values for `beta`. `Shape: (np)`.
      real(wp), intent(in), optional, target :: scld(:, :)
         !! Scaling values for `delta`. `Shape: ({1,n}, m)`. See p. 32.
      real(wp), intent(inout), optional, target :: rwork(:)
         !! 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, target :: lower(:)
         !! Lower bound on `beta`. `Shape: (np)`.
      real(wp), intent(in), optional, target :: upper(:)
         !! Upper bound on `beta`. `Shape: (np)`.

      ! Local variables
      integer :: ldwe, ld2we, ldwd, ld2wd, ldifx, ldscld, ldstpd, job_, ndigit_, maxit_, &
                 iprint_, lunerr_, lunrpt_, info_, lrwork, liwork, info1_, info2_, info3_, &
                 info4_, info5_
      integer, allocatable, target :: ifixb_local(:), ifixx_local(:, :), iwork_local(:)
      integer, pointer :: ifixb_(:), ifixx_(:, :), iwork_(:)

      real(wp) :: taufac_, sstol_, partol_
      real(wp), allocatable :: tempret(:, :)
      real(wp), allocatable, target :: rwork_local(:), lower_local(:), upper_local(:), &
                                       sclb_local(:), scld_local(:, :), stpb_local(:), &
                                       stpd_local(:, :), we_local(:, :, :), wd_local(:, :, :)
      real(wp), pointer :: rwork_(:), lower_(:), upper_(:), sclb_(:), scld_(:, :), stpb_(:), &
                           stpd_(:, :), we_(:, :, :), wd_(:, :, :)

      logical :: head, isodr, restart

      ! Nullify pointers
      nullify (iwork_, rwork_, lower_, upper_, ifixb_, ifixx_, sclb_, scld_, stpb_, stpd_, &
               we_, wd_)

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

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

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

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

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

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

      if (np < 1 .or. np > n) then
         info5_ = 1
         info2_ = 1
      end if

      if (q < 1) 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 print_header(head, lunrpt_)
            call print_error_inputs( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, q, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lrwork, 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(rwork)) 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 print_header(head, lunrpt_)
            call print_error_inputs( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, q, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lrwork, 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, q, np, isodr, lrwork, 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(rwork)) then
         allocate (rwork_local(lrwork), stat=info3_)
         rwork_ => rwork_local
      end if

      allocate (tempret(max(n, np), max(q, 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 print_header(head, lunrpt_)
            call print_error_inputs( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, q, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
               lrwork, liwork)
         end if
         if (present(info)) then
            info = info_
         end if
         return
      end if

      if (present(iwork)) then
         if (size(iwork) == liwork) then
            iwork_ => iwork
         else
            info1_ = info1_ + 8192
         end if
      end if

      if (present(rwork)) then
         if (size(rwork) == lrwork) then
            rwork_ => rwork
         else
            info1_ = info1_ + 4096
         end if
      end if

      ! Check the size of required arguments and return errors if they are too small
      if (any(shape(x) /= [n, m])) then
         info1_ = info1_ + 4
      end if

      if (any(shape(y) /= [n, q])) then
         info1_ = info1_ + 2
      end if

      if (size(beta) /= np) then
         info1_ = info1_ + 1
      end if

      ! Check the presence of optional array arguments
      ! Arrays are pointed to or allocated if necessary

      if (present(delta)) then
         if (any(shape(delta) /= [n, m])) then
            info1_ = info1_ + 8
         end if
      end if

      if (present(we)) then
         ldwe = size(we, 1)
         ld2we = size(we, 2)
         if ((ldwe == 1 .or. ldwe == n) .and. &
             (ld2we == 1 .or. ld2we == q) .and. &
             (size(we, 3) == q)) then
            we_ => we
         else
            info1_ = info1_ + 16
         end if
      else
         ldwe = 1
         ld2we = 1
         allocate (we_local(ldwe, ld2we, q))
         we_local = negone
         we_ => we_local
      end if

      if (present(wd)) then
         ldwd = size(wd, 1)
         ld2wd = size(wd, 2)
         if ((ldwd == 1 .or. ldwd == n) .and. &
             (ld2wd == 1 .or. ld2wd == m) .and. &
             (size(wd, 3) == m)) then
            wd_ => wd
         else
            info1_ = info1_ + 32
         end if
      else
         ldwd = 1
         ld2wd = 1
         allocate (wd_local(ldwd, ld2wd, m))
         wd_local = negone
         wd_ => wd_local
      end if

      if (present(ifixb)) then
         if (size(ifixb) == np) then
            ifixb_ => ifixb
         else
            info1_ = info1_ + 64
         end if
      else
         allocate (ifixb_local(np))
         ifixb_local = 1
         ifixb_local(1) = -1
         ifixb_ => ifixb_local
      end if

      if (present(ifixx)) then
         ldifx = size(ifixx, 1)
         if ((ldifx == 1 .or. ldifx == n) .and. (size(ifixx, 2) == m)) then
            ifixx_ => ifixx
         else
            info1_ = info1_ + 128
         end if
      else
         ldifx = 1
         allocate (ifixx_local(ldifx, m))
         ifixx_local = 1
         ifixx_local(1, 1) = -1
         ifixx_ => ifixx_local
      end if

      if (present(stpb)) then
         if (size(stpb) == np) then
            stpb_ => stpb
         else
            info1_ = info1_ + 256
         end if
      else
         allocate (stpb_local(np))
         stpb_local = negone
         stpb_ => stpb_local
      end if

      if (present(stpd)) then
         ldstpd = size(stpd, 1)
         if ((ldstpd == 1 .or. ldstpd == n) .and. (size(stpd, 2) == m)) then
            stpd_ => stpd
         else
            info1_ = info1_ + 512
         end if
      else
         ldstpd = 1
         allocate (stpd_local(ldstpd, m))
         stpd_local = negone
         stpd_ => stpd_local
      end if

      if (present(sclb)) then
         if (size(sclb) == np) then
            sclb_ => sclb
         else
            info1_ = info1_ + 1024
         end if
      else
         allocate (sclb_local(np))
         sclb_local = negone
         sclb_ => sclb_local
      end if

      if (present(scld)) then
         ldscld = size(scld, 1)
         if ((ldscld == 1 .or. ldscld == n) .and. (size(scld, 2) == m)) then
            scld_ => scld
         else
            info1_ = info1_ + 2048
         end if
      else
         ldscld = 1
         allocate (scld_local(ldscld, m))
         scld_local = negone
         scld_ => scld_local
      end if

      if (present(upper)) then
         if (size(upper) == np) then
            upper_ => upper
         else
            info1_ = info1_ + 16384
         end if
      else
         allocate (upper_local(np))
         upper_local = huge(zero)
         upper_ => upper_local
      end if

      if (present(lower)) then
         if (size(lower) == np) then
            lower_ => lower
         else
            info1_ = info1_ + 32768
         end if
      else
         allocate (lower_local(np))
         lower_local = -huge(zero)
         lower_ => lower_local
      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 print_header(head, lunrpt_)
            call print_error_inputs( &
               lunerr_, info_, info5_, info4_, info3_, info2_, info1_, &
               n, m, q, &
               ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, lrwork, liwork)
         end if
         if (present(info)) then
            info = info_
         end if
         return
      end if

      ! Now it is safe to use the array pointers
      restart = (mod(job_/10000, 10) >= 1)
      if (.not. restart) then
         rwork_ = zero
         iwork_ = 0
      end if

      if (present(delta) .and. (.not. restart)) then
         rwork_(1:n*m) = reshape(delta, [n*m])
      end if

      ! Set default values for optional scalar arguments
      maxit_ = -1
      if (present(maxit)) then
         maxit_ = maxit
      end if

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

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

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

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

      ! Call the driver routine
      call odrdrive( &
         fcn, &
         n, m, np, q, &
         beta, y, x, &
         we_, ldwe, ld2we, &
         wd_, ldwd, ld2wd, &
         ifixb_, ifixx_, ldifx, &
         job_, ndigit_, taufac_, &
         sstol_, partol_, maxit_, &
         iprint_, lunerr_, lunrpt_, &
         stpb_, stpd_, ldstpd, &
         sclb_, scld_, ldscld, &
         rwork_, lrwork, tempret, iwork_, liwork, &
         info_, &
         lower_, upper_)

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

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

   end subroutine odr