!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File:  ct11prnt.f90
!   Subroutines for printing the control solution.
!
! 24 Dec 2008: Updated for v4 of snctrl.
! 09 Feb 2010: v5.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

module ct11prnt
  use precision,    only : ip, rp
  use ctUsrModule,  only : ctProb
  use ctWorkModule, only : ctWork, ctGrid
  use ctGlobal,     only : ncOde, ncAlg
  implicit none

  private
  public  :: c1prnt
  private :: c1psol

contains

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine c1prnt ( prob, grid, x, work, iw, leniw )
    integer(ip), intent(in)    :: leniw
    real(rp),    intent(in)    :: x(:)
    integer(ip), intent(inout) :: iw(leniw)
    type(ctProb) :: prob
    type(ctGrid) :: grid
    type(ctWork) :: work

    !===========================================================================
    ! Control option to print control related output to Print file,
    ! using snPRNT.
    !
    ! 21 Aug 2008: Current version of sncPrnt.
    ! 19 Nov 2008: Added refinement limit option print.
    ! 24 Dec 2008: Updated for v4 of snctrl.
    ! 28 Jan 2009: Updated for new phase/sparse struct implementation.
    ! 09 Feb 2010: v5.
    !===========================================================================
    integer(ip) :: iIntf, iDisc, iRefn, iRefL, nY, nU, nP, nC, nPhs, nNodes
    real(rp)    :: refTol
    character   :: str*120, out*3

    iIntf  = prob%iIntf
    iDisc  = prob%iDisc
    iRefn  = prob%iRefn
    iRefL  = prob%iRefL
    refTol = prob%refTol

    nY   = prob%nY
    nU   = prob%nU
    nP   = prob%nP
    nC   = prob%nC
    nPhs = prob%nPhs

    call snPRNT (13, ' -----------------------------', iw, leniw)
    call snPRNT ( 3, '     SNCTRL Summary Output    ', iw, leniw)
    call snPRNT ( 3, ' -----------------------------', iw, leniw)

    ! Discretization Opt
    if ( iDisc == 0 ) then
       out = ' TR'
    else
       out = ' HS'
    end if
    write (str, 1000) out
    call snPRNT (3, str, iw, leniw)

    ! Problem Data
    write (str, 1010) nY, nU
    call snPRNT (3, str, iw, leniw)
    write (str, 1020) nP, nC
    call snPRNT (3, str, iw, leniw)
    write (str, 1040) nPhs
    call snPRNT (3, str, iw, leniw)

    nNodes = grid%nNodes - (nPhs-1)
    write (str, 1030) grid%nInt, nNodes
    call snPRNT (3, str, iw, leniw)

    ! Refinement Opts
    if (iRefn == 0) then
       out = ' No'
    else
       out = 'Yes'
    end if
    call snPRNT (3, '', iw, leniw)
    write (str, 1050) out, refTol
    call snPRNT (3, str, iw, leniw)
    write (str, 1060) iRefL
    call snPRNT (3, str, iw, leniw)
    call snPRNT (3, '', iw, leniw)

    if (iIntf == 0) then
       out = '  S'
    else if ( iIntf == 1 ) then
       out = '  D'
    else if ( iIntf == 2 ) then
       out = '  A'
    end if

    ! Internal stats
    write (str, 1100) out
    call snPRNT (3, str, iw, leniw)
    write (str, 1110) ncOde, ncAlg
    call snPRNT (3, str, iw, leniw)

    call c1psol ( iDisc, nY, nU, nP, nPhs, grid, x, work%n, iw, leniw )

    return

1000 format ('  Discretization               ', a6)
1010 format ('  Number of states             ', i6, 2x, &
             '  Number of controls           ', i10 )
1020 format ('  Number of parameters         ', i6, 2x, &
             '  Number of add. constraints   ', i10 )
1030 format ('  Total number of intervals    ', i6, 2x, &
             '  Total number of nodes        ', i10)
1040 format ('  Number of phases             ', i6 )

1050 format ('  Refinement                   ', a6, 2x, &
             '  Refinement Tolerance       '  , e12.4)
1060 format ('  Refinement Limit             ', i6 )

1100 format ('  Control Interface            ', a6)
1110 format ('  Number of calls to odecon    ', i6, 2x, &
             '  Number of calls to algcon    ', i10 )

  end subroutine c1prnt

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine c1psol ( iDisc, nY, nU, nP, nPhs, grid, x, n, iw, leniw )
    integer(ip), intent(in)    :: iDisc, nY, nU, nP, nPhs, n, leniw
    real(rp),    intent(in)    :: x(n)
    integer(ip), intent(inout) :: iw(leniw)
    type(ctGrid) :: grid

    !===========================================================================
    ! Prints control interface output to the SNOPT output file.
    !
    ! Prints a column with the time nodes and then 5 columns of the
    ! state/control values.
    !
    ! 21 Aug 2008: Current version of c1psol.
    ! 24 Dec 2008: Updated for v4 of snctrl.
    ! 09 Feb 2010: v5.
    !===========================================================================
    integer(ip) :: nNodes, nVar, pind, tind, ind, vS, vE, which, &
                   i, j, p, c1, c2, ct
    integer(ip),  pointer :: intPtr(:), ndPtr(:)
    real(rp),     pointer :: step(:)
    real(rp), allocatable :: time(:)
    character(120)        :: str


    nNodes =  grid%nNodes
    intPtr => grid%intPtr
    ndPtr  => grid%ndPtr
    step   => grid%step

    allocate ( time(nNodes+1) )
    nVar = nY+nU

    ! Set up time array
    tind = 1
    time(tind) = grid%t0

    if ( iDisc == 0 ) then
       do p = 1, nPhs
          do j = intPtr(p), intPtr(p+1)-1
             tind = tind + 1
             time(tind) = time(tind-1) + step(j)
          end do
          tind = tind + 1
          time(tind) = time(tind-1)
       end do
       time(nNodes+1) = grid%tf

    else
       do p = 1, nPhs
          do j = intPtr(p), intPtr(p+1)-1
             tind = tind + 1
             time(tind) = time(tind-1) + step(j) / 2.0

             tind = tind + 1
             time(tind) = time(tind-1) + step(j) / 2.0
          end do
          tind = tind + 1
          time(tind) = time(tind-1)
       end do
       time(nNodes+1) = grid%tf

    end if


    ! Start with states from the beginning...
    which = 0
    call snPRNT(11, ' State Variables', iw, leniw)
    call snPRNT( 1, ' ---------------', iw, leniw)

    vS  = 1
    vE  = min (5, nY)

100 tind = 0
    ind   = 0

    do p = 1, nPhs
       str = 's'

       do j = ndPtr(p), ndPtr(p+1)
          c1 = 2
          c2 = 16
          tind = tind + 1
          write (str(c1:c2), '(f12.5)') time(tind)

          ct = 0
          do i = vS, vE
             ct = ct + 1
             c1 = 16*ct
             c2 = 16*(ct+1)
             write (str(c1:c2), '(f12.5)') x(i + ind)
          end do
          ind = ind + nVar
          call snPRNT (1, str, iw, leniw)
          str = ' '
       end do

       str = ''
       call snPRNT (1, str, iw, leniw)
    end do


    ! Check where we're at...
    if ( which == 0 ) then
       if ( vE < nY ) then
          call snPRNT(11, ' State Variables', iw, leniw)
          call snPRNT( 1, ' ---------------', iw, leniw)
          vS = vE + 1
          vE = min (5, nY-vE) + vS - 1
          go to 100
       else
          str = ' '
          call snPRNT(11, ' Control Variables', iw, leniw)
          call snPRNT( 1, ' -----------------', iw, leniw)
          which = 1
          vS = 1 + nY
          vE = min (5, nU) + vS - 1
          go to 100
       end if
    end if

    if ( which == 1 ) then
       if ( vE < nU+nY ) then
          call snPRNT(11, ' Control Variables', iw, leniw)
          call snPRNT( 1, ' -----------------', iw, leniw)
          vS = vE + 1
          vE = min (5, nY+nU-vE) + vS - 1
          go to 100
       end if
    end if

   ! Print out parameters if any
   if ( nP > 0 ) then
      pind = n - nP
      call snPRNT(11, ' Parameters     ', iw, leniw)
      call snPRNT( 1, ' ---------------', iw, leniw)
      do i = 1, nP
          write(str, '(f12.5)') x(pind + i)
          call snPRNT(1, str, iw, leniw)
       end do
    end if

    deallocate ( time )

  end subroutine c1psol

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

end module ct11prnt
