text:
!##############################################################################
!# ****************************************************************************
!# <name> ccpostprocessing </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains postprocessing routines for the CC2D solver.
!#
!# The following routines can be found here:
!#
!# 1.) cc_initpostprocessing
!#     -> Initialise the postprocessing.
!#
!# 2.) cc_donepostprocessing
!#     -> Clean up the postprocessing.
!#
!# 3.) cc_postprocessingStationary
!#     -> Postprocessing of a solution in a stationary solution process.
!#        Evaluate the solution of the stationary solver, write GMV-file.
!#
!# 4.) cc_postprocessingNonstat
!#     -> Postprocessing of a solution in a nonstationary solution process.
!#        Evaluate the solution of the nonstationary solver, write GMV-file.
!#
!# Auxiliary routines:
!#
!# 1.) cc_errorAnalysis
!#     -> Perform error analysis (comparison to analytic solution;
!#        as configured in the DAT files).
!#
!# 2.) cc_calculateBodyForces
!#     -> Calculate the body forces on an object.
!#
!# 3.) cc_calculateDivergence
!#     -> Calculate the divergence.
!#
!# 4.) cc_writeUCD
!#     -> Write UCD (AVS, GMV,...) output.
!#
!# 5.) cc_writeFilm
!#     -> Write film output (=raw solution vectors).
!#
!# </purpose>
!##############################################################################

module ccpostprocessing

  use fsystem
  use storage
  use linearsolver
  use boundary
  use linearalgebra
  use cubature
  use matrixfilters
  use vectorfilters
  use discretebc
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use coarsegridcorrection
  use spdiscprojection
  use nonlinearsolver
  use paramlist
  use bilinearformevaluation
  use linearformevaluation
  use statistics
  use element
  use multilevelprojection
  use vectorio

  use collection
  use convection

  use ucd

  use pprocnavierstokes
  use pprocerror

  use ccboundaryconditionparser
  use ccmatvecassembly
  use ccbasic
  use cccallback

  implicit none

!<types>

!<typeblock>

  ! Postprocessing structure. Whenever the postprocessing calculation routines are
  ! called, they fill this structure using the current solution vector (and other
  ! information if necessary). The information in this structure can then be used
  ! for GMV output e.g.
  type t_c2d2postprocessing

    ! A discretisation structure that describes a piecewise constant discretisation
    ! (usually P0 or Q0).
    type(t_spatialDiscretisation) :: rdiscrConstant

    ! A discretisation structure that describes a piecewise linear discretisation
    ! (usually P1 or Q1).
    type(t_spatialDiscretisation) :: rdiscrLinear

    ! A discretisation structure that describes a piecewise quadratic discretisation
    ! (usually P2 or Q2).
    type(t_spatialDiscretisation) :: rdiscrQuadratic

    ! Whether nonstationary postprocessing should be used or not.
    logical              :: bnonstationaryPostprocessing

    ! Point in time when the next UCD file is to be written out
    real(DP)             :: dnextTimeUCD = 0.0_DP

    ! Next file extension for UCD output file.
    integer              :: inextFileSuffixUCD = 0

    ! Point in time when the next Film file is to be written out
    real(DP)             :: dnextTimeFilm = 0.0_DP

    ! Next file extension for Film output file.
    integer              :: inextFileSuffixFilm = 0

    ! A vector that describes the X-velocity field in the vertices
    type(t_vectorScalar) :: rvectorVelX

    ! A vector that describes the Y-velocity field in the vertices
    type(t_vectorScalar) :: rvectorVelY

    ! A vector that describes the pressure field in the vertices
    type(t_vectorScalar) :: rvectorPressure

    ! A vector that describes the pressure field in the cells
    type(t_vectorScalar) :: rvectorPressureCells

    ! A vector that describes the streamfunction
    type(t_vectorScalar) :: rvectorStreamfunction

    ! A vector that describes the H1-error of the velocity field in the vertices
    type(t_vectorScalar) :: rvectorH1err

    ! A vector that describes the H1-error of the pressure in the cells
    type(t_vectorScalar) :: rvectorH1errCells

  end type

!</typeblock>

!</types>

contains

  ! ***************************************************************************

!<subroutine>

  subroutine cc_postprocessingStationary (rproblem,rvector,rpostprocessing)

!<description>
  ! Postprocessing of solutions of stationary simulations.
  ! Writes the solution into a GMV file, calculates forces,...
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem

  ! Postprocessing structure.
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!<input>
  ! The solution vector which is to be evaluated by the postprocessing routines.
  type(t_vectorBlock), intent(IN) :: rvector
!</input>

!</subroutine>

    ! local variables
    type(t_timer) :: rtimer

    call stat_clearTimer(rtimer)
    call stat_startTimer(rtimer)

    ! Calculate body forces.
    call cc_calculateBodyForces (rvector,rproblem)

    ! Calculate the divergence
    call cc_calculateDivergence (rvector,rproblem)

    ! Error analysis, comparison to reference function.
    call cc_errorAnalysis (rvector,rproblem)

    ! Write the UCD export file (GMV, AVS,...) as configured in the DAT file.
    call cc_writeUCD (rpostprocessing, rvector, rproblem)

    ! Gather statistics
    call stat_stopTimer(rtimer)
    rproblem%rstatistics%dtimePostprocessing = &
      rproblem%rstatistics%dtimePostprocessing + rtimer%delapsedReal

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine cc_postprocessingNonstat (rproblem,rvector,rpostprocessing)

!<description>
  ! Postprocessing of solutions of stationary simulations.
  ! Writes the solution into a GMV file, calculates forces,...
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem

  ! Postprocessing structure. Defines what to do with solution vectors.
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!<input>
  ! The solution vector which is to be evaluated by the postprocessing routines.
  type(t_vectorBlock), intent(IN) :: rvector
  !TYPE(t_vectorBlock), INTENT(IN) :: rconcentration
!</input>

!</subroutine>

    ! local variables
    type(t_timer) :: rtimer

    call stat_clearTimer(rtimer)
    call stat_startTimer(rtimer)

    ! Calculate body forces.
    call cc_calculateBodyForces (rvector,rproblem)
    !
    ! Calculate the divergence
    call cc_calculateDivergence (rvector,rproblem)

    ! Error analysis, comparison to reference function.
    call cc_errorAnalysis (rvector,rproblem)

    ! Write the UCD export file (GMV, AVS,...) as configured in the DAT file.
    call cc_writeUCD (rpostprocessing, rvector, rproblem, &
        rproblem%rtimedependence%dtime)
    ! Write film output (raw data vectors)
    call cc_writeFilm (rpostprocessing, rvector, rproblem, &
        rproblem%rtimedependence%dtime)

    ! Gather statistics
    call stat_stopTimer(rtimer)
    rproblem%rstatistics%dtimePostprocessing = &
      rproblem%rstatistics%dtimePostprocessing + rtimer%delapsedReal

  end subroutine

!******************************************************************************

!<subroutine>

  subroutine cc_errorAnalysis (rsolution,rproblem)

!<description>
  ! Performs error analysis on a given solution rsolution as specified
  ! in the .DAT file.
  ! The result of the error analysis is written to the standard output.
!</description>

!<input>
  ! Solution vector to compute the norm/error from.
  type(t_vectorBlock), intent(IN) :: rsolution
!</input>

!<inputoutput>
  ! Problem structure.
  type(t_problem), intent(INOUT) :: rproblem
!</inputoutput>

!</subroutine>

    ! local variables
    real(DP),dimension(4) :: Derr
    real(DP) :: derrorVel, derrorP, derrorConc
    integer :: icalcL2,icalcH1

    call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                     'IERRORANALYSISL2', icalcL2, 0)
    call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                     'IERRORANALYSISH1', icalcH1, 0)

    if ((icalcL2 .ne. 0) .or. (icalcH1 .ne. 0)) then
      call output_lbrk()
      call output_line ('Error Analysis')
      call output_line ('--------------')
    end if

    if (icalcL2 .ne. 0) then

      call cc_initCollectForAssembly (rproblem,rproblem%rcollection)

      ! Perform error analysis to calculate and add 1/2||u-z||_{L^2}.
      call pperr_scalar (rsolution%RvectorBlock(1),PPERR_L2ERROR,Derr(1),&
                         ffunction_TargetX,rproblem%rcollection)

      call pperr_scalar (rsolution%RvectorBlock(2),PPERR_L2ERROR,Derr(2),&
                         ffunction_TargetY,rproblem%rcollection)

      derrorVel = (0.5_DP*(Derr(1)**2+Derr(2)**2))

      call pperr_scalar (rsolution%RvectorBlock(3),PPERR_L2ERROR,Derr(3),&
                         ffunction_TargetP,rproblem%rcollection)

      derrorP = Derr(3)

      call pperr_scalar (rsolution%RvectorBlock(4),PPERR_L2ERROR,Derr(4),&
                         ffunction_TargetC,rproblem%rcollection)

      derrorConc = Derr(4)

      call output_line ('||u-reference||_L2 = '//trim(sys_sdEP(derrorVel,15,6)) )
      call output_line ('||p-reference||_L2 = '//trim(sys_sdEP(derrorP,15,6)) )
      call output_line ('||c-reference||_L2 = '//trim(sys_sdEP(derrorConc,15,6)) )

      call cc_doneCollectForAssembly (rproblem,rproblem%rcollection)

    end if

    if (icalcH1 .ne. 0) then

      call cc_initCollectForAssembly (rproblem,rproblem%rcollection)

      ! Perform error analysis to calculate and add 1/2||u-z||_{L^2}.
      call pperr_scalar (rsolution%RvectorBlock(1),PPERR_H1ERROR,Derr(1),&
                         ffunction_TargetX,rproblem%rcollection)

      call pperr_scalar (rsolution%RvectorBlock(2),PPERR_H1ERROR,Derr(2),&
                         ffunction_TargetY,rproblem%rcollection)

      derrorVel = (0.5_DP*(Derr(1)**2+Derr(2)**2))

      call output_line ('||u-reference||_H1 = '//trim(sys_sdEP(derrorVel,15,6)) )

      call cc_doneCollectForAssembly (rproblem,rproblem%rcollection)

    end if

  end subroutine

!******************************************************************************

!<subroutine>

  subroutine cc_calculateBodyForces (rsolution,rproblem)

!<description>
  ! Calculates body forces as configured in the .DAT file.
  ! The result is written to the standard output.
!</description>

!<input>
  ! Solution vector to compute the norm/error from.
  type(t_vectorBlock), intent(IN) :: rsolution
!</input>

!<inputoutput>
  ! Problem structure.
  type(t_problem), intent(INOUT) :: rproblem
!</inputoutput>

!</subroutine>

    ! local variables

    ! Forces on the object
    real(DP), dimension(NDIM2D) :: Dforces
    real(DP) :: df1,df2
    type(t_boundaryRegion) :: rregion

    ! If we have a uniform discretisation, calculate the body forces on the
    ! 2nd boundary component - if it exists.
    if ((rsolution%p_rblockDiscr%RspatialDiscr(1)% &
         ccomplexity .eq. SPDISC_UNIFORM) .and. &
        (boundary_igetNBoundComp(rproblem%rboundary) .ge. 2)) then

      ! Calculate drag-/lift coefficients on the 2nd boundary component.
      ! This is for the benchmark channel!
      call boundary_createRegion (rproblem%rboundary, &
          2, 0, rregion)
      rregion%iproperties = BDR_PROP_WITHSTART+BDR_PROP_WITHEND
      df1 = 1.0_DP/1000.0_DP
      df2 = 0.1_DP * 0.2_DP**2
      call ppns2D_bdforces_uniform (rsolution,rregion,Dforces,CUB_G1_1D,df1,df2)

      call output_lbrk()
      call output_line ('Body forces')
      call output_line ('-----------')
      call output_line ('Body forces real bd., bdc/horiz/vert')
      call output_line (' 2 / ' &
          //trim(sys_sdEP(Dforces(1),15,6)) // ' / '&
          //trim(sys_sdEP(Dforces(2),15,6)) )

    endif

  end subroutine

!******************************************************************************

!<subroutine>

  subroutine cc_calculateDivergence (rsolution,rproblem)

!<description>
  ! Calculates the divergence of a solution.
  ! The result is written to the standard output.
!</description>

!<input>
  ! Solution vector to compute the norm/error from.
  type(t_vectorBlock), intent(IN) :: rsolution
!</input>

!<inputoutput>
  ! Problem structure.
  type(t_problem), intent(INOUT) :: rproblem
!</inputoutput>

!</subroutine>

    ! local variables
    integer(I32) :: ieltype
    type(t_matrixScalar) :: rBmatrix
    type(t_vectorScalar), target :: rtempVector

    if (rsolution%p_rblockDiscr%RspatialDiscr(1)% &
        ccomplexity .eq. SPDISC_UNIFORM) then

      ieltype = rsolution%p_rblockDiscr%RspatialDiscr(1)% &
                RelementDistr(1)%celement

      if (elem_getPrimaryElement(ieltype) .eq. EL_Q1T) then

        ! Create a temporary vector
        call lsyssc_createVecByDiscr (rsolution%RvectorBlock(3)%p_rspatialDiscr,&
            rtempVector,.true.)

        ! Calculate divergence = B1^T u1 + B2^T u2
        call lsyssc_transposeMatrix (rproblem%RlevelInfo(rproblem%nlmax)%rmatrixB1,&
            rBmatrix,LSYSSC_TR_VIRTUAL)
        call lsyssc_scalarMatVec (&
            rBmatrix, rsolution%RvectorBlock(1), &
            rtempVector, 1.0_DP, 0.0_DP)
        call lsyssc_transposeMatrix (rproblem%RlevelInfo(rproblem%nlmax)%rmatrixB2,&
            rBmatrix,LSYSSC_TR_VIRTUAL)
        call lsyssc_scalarMatVec (&
            rBmatrix, rsolution%RvectorBlock(2), &
            rtempVector, 1.0_DP, 1.0_DP)

        call output_lbrk()
        call output_line ('Divergence')
        call output_line ('----------')
        call output_line ('Divergence = ' &
            //trim(sys_sdEP(lsyssc_vectorNorm(rtempVector,LINALG_NORML2),15,6)) )

        call lsyssc_releaseVector (rtempVector)

      end if

    end if

  end subroutine

!******************************************************************************

!<subroutine>

  subroutine cc_writeUCD (rpostprocessing,rvector,rproblem,dtime)

!<description>
  ! Writes an UCD postprocessing file as configured in the DAT file.
  ! (-> GMV, AVS, Paraview,...)
!</description>

!<input>
  ! Solution vector.
  type(t_vectorBlock), intent(IN) :: rvector
  !TYPE(t_vectorBlock), INTENT(IN),OPTIONAL :: rconcentration
  ! OPTIONAL: Simulation time.
  ! Must be ommitted in stationary simulations.
  real(DP), intent(IN), optional :: dtime
!</input>

!<inputoutput>
  ! Problem structure.
  type(t_problem), intent(INOUT), target :: rproblem

  ! Postprocessing structure. Must have been initialised prior
  ! to calling this routine.
  ! The time stamp of the last written out GMV is updated.
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!</subroutine>

    ! local variables

    ! We need some more variables for postprocessing - i.e. writing
    ! a GMV file.
    real(DP), dimension(:), pointer :: p_Ddata,p_Ddata2

    ! A pointer to the triangulation.
    type(t_triangulation), pointer :: p_rtriangulation

    ! A vector accepting Q1 data
    type(t_vectorBlock) :: rprjVector

    ! A discretisation structure for Q1
    type(t_blockDiscretisation) :: rprjDiscretisation

    ! Discrete boundary conditions for the output vector
    type(t_discreteBC), target :: rdiscreteBC
    type(t_discreteFBC), target :: rdiscreteFBC

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport

    real(DP) :: dminTime, dmaxTime, dtimeDifferenceUCD
    integer :: ioutputUCD,ilevelUCD
    integer(i32) :: ieltype

    character(SYS_STRLEN) :: sfile,sfilename

    if (present(dtime)) then
      ! In a nonstationary simulation, first check if we are allowed
      ! to write something.

      call parlst_getvalue_double (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                      'DMINTIMEUCD', dminTime, -1.E100_DP)
      call parlst_getvalue_double (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                      'DMAXTIMEUCD', dmaxTime, 1.E100_DP)
      call parlst_getvalue_double (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                      'DTIMEDIFFERENCEUCD', dtimeDifferenceUCD, 0.0_DP)

      if ((dtime .lt. dminTime) .or. (dtime .gt. dmaxTime)) return

      if (dtimeDifferenceUCD .gt. 0.0_DP) then
        if (rpostprocessing%dnextTimeUCD .gt.  dtime) return
      else if (dtimeDifferenceUCD .lt. 0.0_DP) then
        if (rpostprocessing%dnextTimeUCD .lt. dtime) return
      end if
      ! Otherwise: Always write!

    end if

    ! Type of output:
    call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
                              'IOUTPUTUCD', ioutputUCD, 0)
    if (ioutputUCD .eq. 0) return

    ! Level of output:
    call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
                              'ILEVELUCD', ilevelUCD, 0)
    if (ilevelUCD .le. 0) then
      ilevelUCD = rproblem%NLMAX+ilevelUCD
    end if

    ilevelUCD = min(rproblem%NLMAX,max(rproblem%NLMIN,ilevelUCD))

    ! The solution vector is probably not in the way, GMV likes it!
    ! GMV for example does not understand Q1~ vectors!
    ! Therefore, we first have to convert the vector to a form that
    ! GMV understands.
    ! GMV understands only Q1 solutions! So the task is now to create
    ! a Q1 solution from rvector and write that out.
    !
    ! For this purpose, first create a 'derived' simple discretisation
    ! structure based on Q1 by copying the main guiding block discretisation
    ! structure and modifying the discretisation structures of the
    ! two velocity subvectors:

    call spdiscr_duplicateBlockDiscr(rvector%p_rblockDiscr,rprjDiscretisation)

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(1), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(1))

    call spdiscr_deriveSimpleDiscrSc (&
                 rvector%p_rblockDiscr%RspatialDiscr(2), &
                 EL_Q1, CUB_G2X2, &
                 rprjDiscretisation%RspatialDiscr(2))

    ! The pressure discretisation substructure stays the old.
    !
    ! Now set up a new solution vector based on this discretisation,
    ! allocate memory.
    call lsysbl_createVecBlockByDiscr (rprjDiscretisation,rprjVector,.false.)

    ! Then take our original solution vector and convert it according to the
    ! new discretisation:
    call spdp_projectSolution (rvector,rprjVector)

    ! Initialise the collection for the assembly process with callback routines.
    ! Basically, this stores the simulation time in the collection if the
    ! simulation is nonstationary.
    call cc_initCollectForAssembly (rproblem,rproblem%rcollection)

    ! Initialise the discrete BC structure
    call bcasm_initDiscreteBC(rdiscreteBC)

    ! Discretise the boundary conditions according to the Q1/Q1/Q0
    ! discretisation for implementing them into a solution vector.
    call cc_assembleBDconditions (rproblem,rprjDiscretisation,&
        rdiscreteBC,rproblem%rcollection,.true.)

    ! Connect the vector to the BC's
    rprjVector%p_rdiscreteBC => rdiscreteBC

    ! The same way, discretise boundary conditions of fictitious boundary components.
    call bcasm_initDiscreteFBC(rdiscreteFBC)
    call cc_assembleFBDconditions (rproblem,rprjDiscretisation,&
        rdiscreteFBC,rproblem%rcollection)
    rprjVector%p_rdiscreteBCfict => rdiscreteFBC

    ! Filter the solution vector to implement discrete BC's.
    call vecfil_discreteBCsol (rprjVector)

    ! Filter the solution vector to implement discrete BC's for fictitious
    ! boundary components.
    call vecfil_discreteFBCsol (rprjVector)

    ! Clean up the collection (as we are done with the assembly, that's it.
    call cc_doneCollectForAssembly (rproblem,rproblem%rcollection)

    ! Basic filename
    call parlst_getvalue_string (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                 'SFILENAMEUCD', sfile, '')

    ! Remove possible ''-characters
    read(sfile,*) sfilename

    ! Create the actual filename
    sfile = trim(adjustl(sfilename))//'.'//sys_si0(rpostprocessing%inextFileSuffixUCD,5)

    ! Now we have a Q1/Q1/Q0 solution in rprjVector -- on the level NLMAX.
    ! The next step is to project it down to level ilevelUCD.
    ! Due to the fact that solutions are usually 2-level-ordered,
    ! this can be shortened by taking only the first NVT vertices
    ! of the solution vector!

    ! From the attached discretisation, get the underlying triangulation
    ! of that level
    p_rtriangulation => rproblem%RlevelInfo(ilevelUCD)%rtriangulation

    ! Start UCD export to GMV file:
    call output_lbrk ()
    call output_line ('Writing GMV file: '//sfile)

    select case (ioutputUCD)
    case (1)
      call ucd_startGMV (rexport,UCD_FLAG_STANDARD,p_rtriangulation,sfile)

    case (2)
      call ucd_startAVS (rexport,UCD_FLAG_STANDARD,p_rtriangulation,sfile)

    case (3)
      call ucd_startVTK (rexport,UCD_FLAG_STANDARD,p_rtriangulation,sfile)

    case DEFAULT
      call output_line ('Invalid UCD ooutput type.', &
                        OU_CLASS_ERROR,OU_MODE_STD,'mysubroutine')
      stop
    end select

    ! Is there a simulation time?
    if (present(dtime)) &
      call ucd_setSimulationTime (rexport,dtime)

    ! Write the configuration of the application as comment block
    ! to the output file.
    call ucd_addCommentLine (rexport,'Configuration:')
    call ucd_addCommentLine (rexport,'---------------')
    call ucd_addParameterList (rexport,rproblem%rparamList)
    call ucd_addCommentLine (rexport,'---------------')

    ! Write velocity field
    call lsyssc_getbase_double (rprjVector%RvectorBlock(1),p_Ddata)
    call lsyssc_getbase_double (rprjVector%RvectorBlock(2),p_Ddata2)

    ! CALL ucd_addVariableVertexBased (rexport,'X-vel',UCD_VAR_XVELOCITY, &
    !     p_Ddata(1:p_rtriangulation%NVT))
    ! CALL ucd_addVariableVertexBased (rexport,'Y-vel',UCD_VAR_YVELOCITY, &
    !     p_Ddata2(1:p_rtriangulation%NVT))
    call ucd_addVarVertBasedVec (rexport,'velocity',&
        p_Ddata(1:p_rtriangulation%NVT),p_Ddata2(1:p_rtriangulation%NVT))

    ! Write pressure
    call lsyssc_getbase_double (rprjVector%RvectorBlock(3),p_Ddata)
    call ucd_addVariableElementBased (rexport,'pressure',UCD_VAR_STANDARD, &
        p_Ddata(1:p_rtriangulation%NEL))

    ! Write concentration
    !IF (PRESENT(rconcentration)) THEN
    CALL lsyssc_getbase_double (rprjVector%RvectorBlock(4),p_Ddata)
    CALL ucd_addVariableVertexBased (rexport,'concentration',UCD_VAR_STANDARD, &
        p_Ddata)
    !END IF
    ! If we have a simple Q1~ discretisation, calculate the streamfunction.
    if (rvector%p_rblockDiscr%RspatialDiscr(1)% &
        ccomplexity .eq. SPDISC_UNIFORM) then

      ieltype = rvector%p_rblockDiscr%RspatialDiscr(1)% &
                RelementDistr(1)%celement

      if (elem_getPrimaryElement(ieltype) .eq. EL_Q1T) then

        call ppns2D_streamfct_uniform (rvector,rprjVector%RvectorBlock(1))

        call lsyssc_getbase_double (rprjVector%RvectorBlock(1),p_Ddata)
        call ucd_addVariableVertexBased (rexport,'streamfunction',&
            UCD_VAR_STANDARD, p_Ddata(1:p_rtriangulation%NVT))

      end if

    end if

    ! Write the file to disc, that's it.
    call ucd_write (rexport)
    call ucd_release (rexport)

    ! Release the auxiliary vector
    call lsysbl_releaseVector (rprjVector)

    ! Release the discretisation structure.
    call spdiscr_releaseBlockDiscr (rprjDiscretisation)

    ! Throw away the discrete BC's - not used anymore.
    call bcasm_releaseDiscreteBC (rdiscreteBC)
    call bcasm_releaseDiscreteFBC (rdiscreteFBC)

    if (present(dtime)) then
      ! Update time stamp of last written out GMV.
      rpostprocessing%dnextTimeUCD = rpostprocessing%dnextTimeUCD+dtimeDifferenceUCD
      rpostprocessing%inextFileSuffixUCD = rpostprocessing%inextFileSuffixUCD + 1
    end if

  end subroutine

!******************************************************************************

!<subroutine>

  subroutine cc_writeFilm (rpostprocessing,rvector,rproblem,dtime)

!<description>
  ! Writes Film output (raw data vectors) to a file as configured in the
  ! DAT file.
  !
  ! Note: This file is usually only used in a nonstationary simulation.
  ! In a stationary simulation, Film output makes no sense!
!</description>

!<input>
  ! Solution vector.
  type(t_vectorBlock), intent(IN) :: rvector

  ! Simulation time.
  real(DP), intent(IN) :: dtime
!</input>

!<inputoutput>
  ! Problem structure.
  type(t_problem), intent(INOUT), target :: rproblem

  ! Postprocessing structure. Must have been initialised prior
  ! to calling this routine.
  ! The time stamp of the last written out Film file is updated.
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!</subroutine>

    ! local variables
    real(DP) :: dminTime, dmaxTime, dtimeDifferenceFilm
    integer :: ioutputFilm,ilevelFilm

    type(t_vectorBlock) :: rvector1,rvector2
    type(t_vectorScalar) :: rvectorTemp
    character(LEN=SYS_STRLEN) :: sfile,sfilename
    integer :: ilev
    integer :: NEQ
    type(t_interlevelProjectionBlock) :: rprojection
    logical :: bformatted

    ! Type of output:
    call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
                              'IOUTPUTFILM', ioutputFilm, 0)
    if (ioutputFilm .eq. 0) return

    ! First check if we are allowed to write something.
    call parlst_getvalue_double (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                    'DMINTIMEFILM', dminTime, -1.E100_DP)
    call parlst_getvalue_double (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                    'DMAXTIMEFILM', dmaxTime, 1.E100_DP)
    call parlst_getvalue_double (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                    'DTIMEDIFFERENCEFILM', dtimeDifferenceFilm, 0.0_DP)

    if ((dtime .lt. dminTime) .or. (dtime .gt. dmaxTime)) return

    if (dtimeDifferenceFilm .gt. 0.0_DP) then
      if (rpostprocessing%dnextTimeFilm .gt.  dtime) return
    else if (dtimeDifferenceFilm .lt. 0.0_DP) then
      if (rpostprocessing%dnextTimeFilm .lt. dtime) return
    end if

    ! Basic filename
    call parlst_getvalue_string (rproblem%rparamList, 'CC-POSTPROCESSING', &
                                 'SFILENAMEFILM', sfile, '')

    ! Remove possible ''-characters
    read(sfile,*) sfilename

    ! Create the actual filename
    sfile = trim(adjustl(sfilename))//'.'//sys_si0(rpostprocessing%inextFileSuffixFilm,5)

    ! Level of output:
    call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
                              'ILEVELFILM', ilevelFilm, 0)
    if (ilevelFilm .le. 0) then
      ilevelFilm = rproblem%NLMAX+ilevelFilm
    end if

    ilevelFilm = min(rproblem%NLMAX,max(rproblem%NLMIN,ilevelFilm))

    if (ilevelFilm .lt. rproblem%NLMIN) then
      call output_line ('Warning: Level for solution vector is < NLMIN! &
          &Writing out at level NLMIN!', &
          OU_CLASS_WARNING,OU_MODE_STD,'cc_releasePreconditioner')
      call sys_halt()
      ilevelFilm = rproblem%NLMIN
    end if

    ! Write formatted output?
    bformatted = ioutputFilm .ne. 2

    ! Interpolate the solution down to level istart.
    call lsysbl_copyVector (rvector,rvector1)   ! creates new rvector1!

    do ilev = rproblem%NLMAX,ilevelFilm+1,-1

      ! Initialise a vector for the lower level and a prolongation structure.
      call lsysbl_createVectorBlock (&
          rproblem%RlevelInfo(ilev-1)%rdiscretisation,rvector2,.false.)

      call mlprj_initProjectionVec (rprojection,rvector2)

      ! Interpolate to the next higher level.
      ! (Don't 'restrict'! Restriction would be for the dual space = RHS vectors!)

      NEQ = mlprj_getTempMemoryVec (rprojection,rvector2,rvector1)
      if (NEQ .ne. 0) call lsyssc_createVector (rvectorTemp,NEQ,.false.)
      call mlprj_performInterpolation (rprojection,rvector2,rvector1, &
                                       rvectorTemp)
      if (NEQ .ne. 0) call lsyssc_releaseVector (rvectorTemp)

      ! Swap rvector1 and rvector2. Release the fine grid vector.
      call lsysbl_swapVectors (rvector1,rvector2)
      call lsysbl_releaseVector (rvector2)

      call mlprj_doneProjection (rprojection)

    end do

    call output_lbrk ()
    call output_line ('Writing Film file: '//sfile)

    ! Write out the solution.
    if (bformatted) then
      call vecio_writeBlockVectorHR (rvector1, 'SOLUTION', .true.,&
         0, sfile, '(E22.15)')
    else
      call vecio_writeBlockVectorHR (rvector1, 'SOLUTION', .true.,0, sfile)
    end if

    ! Release temp memory.
    call lsysbl_releaseVector (rvector1)

    ! Update time stamp of last written out Film file.
    rpostprocessing%dnextTimeFilm = rpostprocessing%dnextTimeFilm+dtimeDifferenceFilm
    rpostprocessing%inextFileSuffixFilm = rpostprocessing%inextFileSuffixFilm + 1

  end subroutine

  !****************************************************************************

!<subroutine>

  subroutine cc_initpostprocessing (rproblem,rpostprocessing)

!<description>
  ! Initialises the given postprocessing structure rpostprocessing
  ! according to the main problem rproblem. The structure can then be used
  ! to generate postprocessing data.
!</description>

!<input>
  ! A problem structure that describes the main problem to solve.
  type(t_problem), intent(IN),target :: rproblem
!</input>

!<output>
  ! Postprocessing structure.
  type(t_c2d2postprocessing), intent(OUT) :: rpostprocessing
!</output>

!</subroutine>

  ! local variables
  type(t_blockDiscretisation), pointer :: p_rdiscr

    ! For postprocessing, we need discretisation structures in the Q0 and Q1 space,
    ! later perhaps in the Q2 space. For this purpose, derive the corresponding
    ! discretisation structure using the 'main' discretisation structure on the
    ! maximum level.
    !
    ! For simplicity, we use only the discretisation structure of the X-velocity
    ! to derive everything.

    p_rdiscr => rproblem%RlevelInfo(rproblem%NLMAX)%rdiscretisation

    ! Piecewise constant space:
    call spdiscr_deriveSimpleDiscrSc (&
                 p_rdiscr%RspatialDiscr(1), &
                 EL_Q0, CUB_G1X1, &
                 rpostprocessing%rdiscrConstant)

    ! Piecewise linear space:
    call spdiscr_deriveSimpleDiscrSc (&
                 p_rdiscr%RspatialDiscr(1), &
                 EL_Q1, CUB_G2X2, &
                 rpostprocessing%rdiscrLinear)

    ! Piecewise quadratic space:
    call spdiscr_deriveSimpleDiscrSc (&
                 p_rdiscr%RspatialDiscr(1), &
                 EL_Q2, CUB_G3X3, &
                 rpostprocessing%rdiscrQuadratic)

    ! Initialise the time/file suffix when the first UCD file is to be written out.
    rpostprocessing%bnonstationaryPostprocessing = (rproblem%itimedependence .ne. 0)
    if (rproblem%itimedependence .ne. 0) then
      rpostprocessing%dnextTimeUCD = rproblem%rtimedependence%dtimeInit
      call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
         'ISTARTSUFFIXUCD', rpostprocessing%inextFileSuffixUCD, 1)

      rpostprocessing%dnextTimeFilm = rproblem%rtimedependence%dtimeInit
      call parlst_getvalue_int (rproblem%rparamList, 'CC-POSTPROCESSING', &
         'ISTARTSUFFIXFILM', rpostprocessing%inextFileSuffixFilm, 1)
    end if

  end subroutine

  !****************************************************************************

!<subroutine>

  subroutine cc_copyPostprocessing (rpostprocessingSrc,rpostprocessingDst)

!<description>
  ! Copies the state of the postprocessing from rpostprocessingSrc
  ! to rpostprocessingDst.
  ! This is typically called to back up or to restore the current state
  ! of a running postprocessing structure.
!</description>

!<input>
  ! Source Postprocessing structure.
  type(t_c2d2postprocessing), intent(IN) :: rpostprocessingSrc
!</input>

!<inputoutput>
  ! Destination Postprocessing structure.
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessingDst
!</inputoutput>

!</subroutine>

    ! Initialise the time/file suffix when the first UCD file is to be written out.
    if (rpostprocessingSrc%bnonstationaryPostprocessing) then
      rpostprocessingDst%bnonstationaryPostprocessing = &
          rpostprocessingSrc%bnonstationaryPostprocessing

      rpostprocessingDst%dnextTimeUCD = rpostprocessingSrc%dnextTimeUCD
      rpostprocessingDst%inextFileSuffixUCD = rpostprocessingSrc%inextFileSuffixUCD

      rpostprocessingDst%dnextTimeFilm = rpostprocessingSrc%dnextTimeFilm
      rpostprocessingDst%inextFileSuffixFilm = rpostprocessingSrc%inextFileSuffixFilm
    end if

  end subroutine

  !****************************************************************************

!<subroutine>

  subroutine cc_clearpostprocessing (rpostprocessing)

!<description>
  ! Releases all calculated data in the given postprocessing structure
  ! so that it can be allocated again in the calculation routines.
  ! This routine must be called at the end of the postprocessing routines
  ! to release the temporary memory that was allocated for the vectors
  ! in the postprocessing structure.
!</description>

!<inputoutput>
  ! Postprocessing structure.
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!</subroutine>

    ! Release all vectors which might be allocated.
    if (rpostprocessing%rvectorVelX%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorVelX)
    if (rpostprocessing%rvectorVelY%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorVelY)
    if (rpostprocessing%rvectorPressure%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorPressure)
    if (rpostprocessing%rvectorPressureCells%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorPressureCells)
    if (rpostprocessing%rvectorStreamfunction%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorStreamfunction)
    if (rpostprocessing%rvectorH1err%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorH1err)
    if (rpostprocessing%rvectorH1errCells%NEQ .ne. 0) &
      call lsyssc_releaseVector (rpostprocessing%rvectorH1errCells)

  end subroutine

  !****************************************************************************

!<subroutine>

  subroutine cc_donepostprocessing (rpostprocessing)

!<description>
  ! Releases a given problem structure. All allocated memory of this structure
  ! is released.
!</description>

!<inputoutput>
  type(t_c2d2postprocessing), intent(INOUT) :: rpostprocessing
!</inputoutput>

!</subroutine>

    ! Release all vectors -- if there are still some allocated
    call cc_clearpostprocessing (rpostprocessing)

    ! Release the discretisation structures allocated above.
    call spdiscr_releaseDiscr(rpostprocessing%rdiscrQuadratic)
    call spdiscr_releaseDiscr(rpostprocessing%rdiscrLinear)
    call spdiscr_releaseDiscr(rpostprocessing%rdiscrConstant)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
Let Nat = forall A. (A -> A) -> A -> A in
let Z: Nat = Lambda A. lambda f: A -> A. lambda x: A. x in
let S: Nat -> Nat = lambda n: Nat. Lambda A. lambda f: A -> A. lambda x: A. f (n [A] f x) in
let add: Nat -> Nat -> Nat = lambda n: Nat. lambda m: Nat. n [Nat] S m in
let mult: Nat -> Nat -> Nat = lambda n: Nat. lambda m: Nat. n [Nat] (add m) Z in
let two: Nat = mult (S Z) (S (S Z)) in
two

ext: f

max_stars_repo_name: eperdew/System-F

max_issues_repo_name: eperdew/System-F

max_forks_repo_name: eperdew/System-F


text:
      subroutine dgefa(a,lda,n,ipvt,info)
      integer lda,n,ipvt(n),info
      double precision a(lda,n)
c
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c
c     internal variables
c
      double precision t
      integer idamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0d0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end
      subroutine dgesl(a,lda,n,ipvt,b,job)
      integer lda,n,ipvt(n),job
      double precision a(lda,n),b(n)
c
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output from dgeco or dgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from dgeco or dgefa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c
c     internal variables
c
      double precision ddot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call daxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            t = ddot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end
      subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
      integer lda,n,ml,mu,ipvt(n),info
      double precision abd(lda,n)
c
c     dgbfa factors a double precision band matrix by elimination.
c
c     dgbfa is usually called by dgbco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c
c     on entry
c
c        abd     double precision(lda, n)
c                contains the matrix in band storage.  the columns
c                of the matrix are stored in the columns of  abd  and
c                the diagonals of the matrix are stored in rows
c                ml+1 through 2*ml+mu+1 of  abd .
c                see the comments below for details.
c
c        lda     integer
c                the leading dimension of the array  abd .
c                lda must be .ge. 2*ml + mu + 1 .
c
c        n       integer
c                the order of the original matrix.
c
c        ml      integer
c                number of diagonals below the main diagonal.
c                0 .le. ml .lt. n .
c
c        mu      integer
c                number of diagonals above the main diagonal.
c                0 .le. mu .lt. n .
c                more efficient if  ml .le. mu .
c     on return
c
c        abd     an upper triangular matrix in band storage and
c                the multipliers which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgbsl will divide by zero if
c                     called.  use  rcond  in dgbco for a reliable
c                     indication of singularity.
c
c     band storage
c
c           if  a  is a band matrix, the following program segment
c           will set up the input.
c
c                   ml = (band width below the diagonal)
c                   mu = (band width above the diagonal)
c                   m = ml + mu + 1
c                   do 20 j = 1, n
c                      i1 = max0(1, j-mu)
c                      i2 = min0(n, j+ml)
c                      do 10 i = i1, i2
c                         k = i - j + m
c                         abd(k,j) = a(i,j)
c                10    continue
c                20 continue
c
c           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
c           in addition, the first  ml  rows in  abd  are used for
c           elements generated during the triangularization.
c           the total number of rows needed in  abd  is  2*ml+mu+1 .
c           the  ml+mu by ml+mu  upper left triangle and the
c           ml by ml  lower right triangle are not referenced.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c     fortran max0,min0
c
c     internal variables
c
      double precision t
      integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
c
c
      m = ml + mu + 1
      info = 0
c
c     zero initial fill-in columns
c
      j0 = mu + 2
      j1 = min0(n,m) - 1
      if (j1 .lt. j0) go to 30
      do 20 jz = j0, j1
         i0 = m + 1 - jz
         do 10 i = i0, ml
            abd(i,jz) = 0.0d0
   10    continue
   20 continue
   30 continue
      jz = j1
      ju = 0
c
c     gaussian elimination with partial pivoting
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 130
      do 120 k = 1, nm1
         kp1 = k + 1
c
c        zero next fill-in column
c
         jz = jz + 1
         if (jz .gt. n) go to 50
         if (ml .lt. 1) go to 50
            do 40 i = 1, ml
               abd(i,jz) = 0.0d0
   40       continue
   50    continue
c
c        find l = pivot index
c
         lm = min0(ml,n-k)
         l = idamax(lm+1,abd(m,k),1) + m - 1
         ipvt(k) = l + k - m
c
c        zero pivot implies this column already triangularized
c
         if (abd(l,k) .eq. 0.0d0) go to 100
c
c           interchange if necessary
c
            if (l .eq. m) go to 60
               t = abd(l,k)
               abd(l,k) = abd(m,k)
               abd(m,k) = t
   60       continue
c
c           compute multipliers
c
            t = -1.0d0/abd(m,k)
            call dscal(lm,t,abd(m+1,k),1)
c
c           row elimination with column indexing
c
            ju = min0(max0(ju,mu+ipvt(k)),n)
            mm = m
            if (ju .lt. kp1) go to 90
            do 80 j = kp1, ju
               l = l - 1
               mm = mm - 1
               t = abd(l,j)
               if (l .eq. mm) go to 70
                  abd(l,j) = abd(mm,j)
                  abd(mm,j) = t
   70          continue
               call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
   80       continue
   90       continue
         go to 110
  100    continue
            info = k
  110    continue
  120 continue
  130 continue
      ipvt(n) = n
      if (abd(m,n) .eq. 0.0d0) info = n
      return
      end
      subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
      integer lda,n,ml,mu,ipvt(n),job
      double precision abd(lda,n),b(n)
c
c     dgbsl solves the double precision band system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgbco or dgbfa.
c
c     on entry
c
c        abd     double precision(lda, n)
c                the output from dgbco or dgbfa.
c
c        lda     integer
c                the leading dimension of the array  abd .
c
c        n       integer
c                the order of the original matrix.
c
c        ml      integer
c                number of diagonals below the main diagonal.
c
c        mu      integer
c                number of diagonals above the main diagonal.
c
c        ipvt    integer(n)
c                the pivot vector from dgbco or dgbfa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b , where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgbco has set rcond .gt. 0.0
c        or dgbfa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c     fortran min0
c
c     internal variables
c
      double precision ddot,t
      integer k,kb,l,la,lb,lm,m,nm1
c
      m = mu + ml + 1
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve l*y = b
c
         if (ml .eq. 0) go to 30
         if (nm1 .lt. 1) go to 30
            do 20 k = 1, nm1
               lm = min0(ml,n-k)
               l = ipvt(k)
               t = b(l)
               if (l .eq. k) go to 10
                  b(l) = b(k)
                  b(k) = t
   10          continue
               call daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
   20       continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/abd(m,k)
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = -b(k)
            call daxpy(lm,t,abd(la,k),1,b(lb),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = ddot(lm,abd(la,k),1,b(lb),1)
            b(k) = (b(k) - t)/abd(m,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (ml .eq. 0) go to 90
         if (nm1 .lt. 1) go to 90
            do 80 kb = 1, nm1
               k = n - kb
               lm = min0(ml,n-k)
               b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
               l = ipvt(k)
               if (l .eq. k) go to 70
                  t = b(l)
                  b(l) = b(k)
                  b(k) = t
   70          continue
   80       continue
   90    continue
  100 continue
      return
      end
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      integer i,incx,incy,ix,iy,m,mp1,n
      double precision dx(*),dy(*),da
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
      subroutine dcopy(n,sx,incx,sy,incy)
c
c     copies a vector, x, to a vector, y.
c     uses unrolled loops for increments equal to 1.
c     jack dongarra, linpack, 3/11/78.
c
      double precision sx(*),sy(*)
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,7)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sx(i)
   30 continue
      if( n .lt. 7 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,7
        sy(i) = sx(i)
        sy(i + 1) = sx(i + 1)
        sy(i + 2) = sx(i + 2)
        sy(i + 3) = sx(i + 3)
        sy(i + 4) = sx(i + 4)
        sy(i + 5) = sx(i + 5)
        sy(i + 6) = sx(i + 6)
   50 continue
      return
      end
      subroutine dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loop for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision da,dx(*)
      integer i,incx,m,mp1,n,nincx
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(*),dy(*),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
      integer function idamax(n,dx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(*),dmax
      integer i,incx,ix,n
c
      idamax = 0
      if(n.lt.1) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      dmax = dabs(dx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(dabs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = dabs(dx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 dmax = dabs(dx(1))
      do 30 i = 2,n
         if(dabs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = dabs(dx(i))
   30 continue
      return
      end
      DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX )
*     .. Scalar Arguments ..
      INTEGER                           INCX, N
*     .. Array Arguments ..
      DOUBLE PRECISION                  X( * )
*     ..
*
*  DNRM2 returns the euclidean norm of a vector via the function
*  name, so that
*
*     DNRM2 := sqrt( x'*x )
*
*
*
*  -- This version written on 25-October-1982.
*     Modified on 14-October-1993 to inline the call to DLASSQ.
*     Sven Hammarling, Nag Ltd.
*
*
*     .. Parameters ..
      DOUBLE PRECISION      ONE         , ZERO
      PARAMETER           ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     .. Local Scalars ..
      INTEGER               IX
      DOUBLE PRECISION      ABSXI, NORM, SCALE, SSQ
*     .. Intrinsic Functions ..
      INTRINSIC             ABS, SQRT
*     ..
*     .. Executable Statements ..
      IF( N.LT.1 .OR. INCX.LT.1 )THEN
         NORM  = ZERO
      ELSE IF( N.EQ.1 )THEN
         NORM  = ABS( X( 1 ) )
      ELSE
         SCALE = ZERO
         SSQ   = ONE
*        The following loop is equivalent to this call to the LAPACK
*        auxiliary routine:
*        CALL DLASSQ( N, X, INCX, SCALE, SSQ )
*
         DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX
            IF( X( IX ).NE.ZERO )THEN
               ABSXI = ABS( X( IX ) )
               IF( SCALE.LT.ABSXI )THEN
                  SSQ   = ONE   + SSQ*( SCALE/ABSXI )**2
                  SCALE = ABSXI
               ELSE
                  SSQ   = SSQ   +     ( ABSXI/SCALE )**2
               END IF
            END IF
   10    CONTINUE
         NORM  = SCALE * SQRT( SSQ )
      END IF
*
      DNRM2 = NORM
      RETURN
*
*     End of DNRM2.
*
      END

max_stars_repo_name: davidbrochart/daetk

max_issues_repo_name: davidbrochart/daetk

max_forks_repo_name: davidbrochart/daetk


module constrained
    use       ls, only: lsquad
    use   trisys, only: backcol, forwcol
    use cholesky, only: cholcol
implicit none
    integer :: n           ! Dimensão do problema
    integer :: m           ! Número de restrições do problema
    double precision :: mu ! Peso da penalidade

    interface
        function f_t(x)
            double precision, intent(in) :: x(:)
            double precision             :: f_t
        end function f_t

        subroutine c_t(cx, x)
            double precision, intent(out) :: cx(:)
            double precision, intent(in)  :: x(:)
        end subroutine c_t
    end interface

    procedure (f_t), pointer :: f_ptr => null()
    procedure (c_t), pointer :: c_ptr => null()
contains

    function Q(x)
        double precision, intent(in) :: x(:)
        double precision             :: Q

        double precision :: cx(size(x))

        call c_ptr(cx, x)
        Q = f_ptr(x) + 1/(2*mu)*sum(cx**2)
    end function Q

    subroutine dQ(dQx, gx, cx, dcx)
        double precision, intent(out) :: dQx(:)
        double precision, intent(in)  :: gx(:)
        double precision, intent(in)  :: cx(:)
        double precision, intent(in)  :: dcx(:, :)
        integer :: k, i

        dQx(:) = gx(:)

        do k = 1, m
            do i = 1, n
                dQx(i) = dQx(i) + cx(k)*dcx(k, i)/mu
            end do
        end do
    end subroutine dQ

    subroutine S(Sx, hx, cx, dcx, d2cx)
        double precision, intent(out) :: Sx(:, :)
        double precision, intent(in)  :: hx(:, :)
        double precision, intent(in)  :: cx(:)
        double precision, intent(in)  :: dcx(:, :)
        double precision, intent(in)  :: d2cx(:, :, :)
        integer :: i, j, k

        do j = 1, n
            do i = 1, n
                Sx(i, j) = hx(i, j)

                do k = 1, m
                    Sx(i, j) = Sx(i, j) + cx(k)*d2cx(k, i, j)/mu
                end do
            end do

            do i = n+1, n+m
                Sx(i, j) = dcx(i-n, j)
            end do
        end do

        do j = n+1, n+m
            do i = 1, n
                Sx(i, j) = dcx(j-n, i)
            end do

            do i = n+1, n+m
                if (i == j) then
                    Sx(i, j) = -mu
                else
                    Sx(i, j) = 0.d0
                end if
            end do
        end do
    end subroutine S

    subroutine v(vx, dQx)
        double precision, intent(out) :: vx(:)
        double precision, intent(in) :: dQx(:)
        integer i

        do i = 1, n
            vx(i) = -dQx(i)
        end do

        do i = n+1, n+m
            vx(i) = 0.d0
        end do
    end subroutine v

    subroutine newt(x, f, g, h, c, dc, d2c, gamma, sigma, theta, eps, iteration, armijo, norm, angle)
        ! Saída
        double precision, intent(inout) :: x(:)

        ! Entrada
        double precision, intent(in) :: eps     ! Tolerância
        double precision, intent(in) :: gamma   ! γ > 0, Constante para condição de Armijo
        double precision, intent(in) :: sigma   ! σ > 0, constante para condição da norma
        double precision, intent(in) :: theta   ! Θ > 0, constante para condição do ângulo

        ! Função, gradiente, hessiana e busca linear
        interface
            function f(x)
                double precision, intent(in) :: x(:)
                double precision             :: f
            end function f

            subroutine g(gx, x)
                double precision, intent(out) :: gx(:)
                double precision, intent(in)  :: x(:)
            end subroutine g

            subroutine h(hx, x)
                double precision, intent(out) :: hx(:, :)
                double precision, intent(in)  :: x(:)
            end subroutine h

            subroutine c(cx, x)
                double precision, intent(out) :: cx(:)
                double precision, intent(in)  :: x(:)
            end subroutine c

            subroutine dc(dcx, x)
                double precision, intent(out) :: dcx(:, :)
                double precision, intent(in)  :: x(:)
            end subroutine dc

            subroutine d2c(d2cx, x)
                double precision, intent(out) :: d2cx(:, :, :)
                double precision, intent(in)  :: x(:)
            end subroutine d2c

            subroutine iteration()
            end subroutine iteration

            subroutine armijo()
            end subroutine armijo

            subroutine norm()
            end subroutine norm

            subroutine angle()
            end subroutine angle
        end interface

        ! Variáveis locais
        double precision :: d(n)    ! Direção
        double precision :: dz(n+m) ! Direção aumentada
        integer          :: p(n)    ! Vetor de permutação
        double precision :: alpha   ! Passo
        double precision :: rho     ! Constante para shift nos autovalores
                                    ! da hessiana

        integer          :: i, status, tmp

        ! Variáveis para evitar chamadas de função desnecessárias
        double precision :: xd(n)
        double precision :: Sr(n+m, n+m) ! S(x)
        double precision :: dQTd         ! ∇ᵀf(x)d
        double precision :: ndQ          ! ‖∇f(x)‖
        double precision :: nd           ! ‖d‖

        double precision :: fx            ! f(x)
        double precision :: gx(n)         ! ∇f(x)
        double precision :: hx(n, n)      ! ∇²f(x)

        double precision :: cx(m)         ! [c₁(x) ..   cₘ(x)]ᵀ
        double precision :: dcx(m, n)     ! [∇c₁(x) ..  ∇cₘ(x)]ᵀ = A(x)
        double precision :: d2cx(m, n, n) ! [∇²c₁(x) .. ∇²cₘ(x)]ᵀ

        double precision :: Qx            !   Q(x, μ) =     f(x) + 1/2μ*∑cᵢ(x)²
        double precision :: dQx(n)        !  ∇Q(x, μ) =   ∇f(x) + 1/μ*∑cᵢ(x)∇cᵢ(x)
        double precision :: d2Qx(n, n)    ! ∇²Q(x, μ) =  ∇²f(x) + 1/μ*∑cᵢ(x)∇²cᵢ(x) + 1/μ*A(x)ᵀA(x)

        double precision :: vx(n+m), Sx(n+m, n+m) ! S(x)*[d z]ᵀ = v(x)

        ! Valores iniciais
        call g(gx, x)
        call c(cx, x)
        call dc(dcx, x)
        call dQ(dQx, gx, cx, dcx)
        call v(vx, dQx)

        do while (norm2(dQx) > eps)
            call iteration()

            ! Calcula fator de Cholesky de uma aproximação
            ! definida positiva de S(x)
            !     GGᵀ = S(x) + ρI
            call h(hx, x)
            call d2c(d2cx, x)
            call S(Sx, hx, cx, dcx, d2cx)

            ! Verifica se a própria S(x) não é definida positiva
            Sr  = Sx
            status = cholcol(n+m, Sr)
            if (status == 0) then
                ! É definida positiva

                dz = vx
                ! Resolve Gᵀdz solução da equação
                !    GGᵀdz = v(x)
                tmp = forwcol(n+m, Sr, dz, unit = .false.)

                ! Resolve dz solução da equação
                !     Gᵀdz = G⁻¹v(x)
                tmp = backcol(n+m, Sr, dz, trans = .true.)
            end if

            ! Extrai a direção de descida do vetor aumentado
            do i = 1, n
                d(i) = dz(i)
            end do

            ndQ  = norm2(dQx)
            nd  = norm2(d)
            dQTd = dot_product(dQx, d)
            rho = 1e-3

            ! Condição do ângulo
            do while (status == -1 .or. dQTd > -theta*ndQ*nd)
                call angle()
                ! Faz shift nos autovalores e tenta novamente
                ! Sr = S(x) + ρI
                Sr = Sx
                do i = 1, n+m
                    Sr(i, i) = Sr(i, i) + rho
                end do

                status = cholcol(n+m, Sr)
                if (status == 0) then
                    ! É definida positiva

                    dz = vx
                    ! Resolve Gᵀd solução da equação
                    !    GGᵀdz = v(x)
                    tmp = forwcol(n+m, Sr, dz, unit = .false.)

                    ! Resolve d solução da equação
                    !     Gᵀdz = G⁻¹v(x)
                    tmp = backcol(n+m, Sr, dz, trans = .true.)
                end if

                ! Extrai a direção de descida do vetor aumentado
                do i = 1, n
                    d(i) = dz(i)
                end do

                ndQ  = norm2(dQx)
                nd   = norm2(d)
                dQTd = dot_product(dQx, d)
                rho  = 10*rho
            end do

            ! Condição da norma
            if (nd < sigma*ndQ) then
                call norm()
                d = sigma*ndQ/nd*d
            end if

            dQTd = dot_product(dQx, d)

            ! Busca linear
            alpha = 1.d0
            call lsquad(alpha, x, d, Q, gamma, dQTd, xd, armijo)

            ! Atualiza x com o novo valor
            x  = x + alpha*d
            call g(gx, x)
            call c(cx, x)
            call dc(dcx, x)
            call dQ(dQx, gx, cx, dcx)
            call v(vx, dQx)
        end do
    end subroutine newt

    subroutine penalty(x, x0, n_, m_, mu_, f, g, h, c, dc, d2c, &
                        gamma, sigma, theta, eps, subproblem, iteration, armijo, norm, angle)
        ! Entrada
        double precision, intent(in)    :: x0(:)   ! Ponto inicial
        integer,          intent(in)    :: n_, m_
        double precision, intent(inout) :: mu_     ! μ > 0, valor inicial para o peso da penalidade
        double precision, intent(in)    :: gamma   ! γ > 0, Constante para condição de Armijo
        double precision, intent(in)    :: sigma   ! σ > 0, constante para condição da norma
        double precision, intent(in)    :: theta   ! Θ > 0, constante para condição do ângulo
        double precision, intent(in)    :: eps     ! Tolerância

        ! Função, gradiente, hessiana e busca linear
        interface
            function f(x)
                double precision, intent(in) :: x(:)
                double precision             :: f
            end function f

            subroutine g(gx, x)
                double precision, intent(out) :: gx(:)
                double precision, intent(in)  :: x(:)
            end subroutine g

            subroutine h(hx, x)
                double precision, intent(out) :: hx(:, :)
                double precision, intent(in)  :: x(:)
            end subroutine h

            subroutine c(cx, x)
                double precision, intent(out) :: cx(:)
                double precision, intent(in)  :: x(:)
            end subroutine c

            subroutine dc(dcx, x)
                double precision, intent(out) :: dcx(:, :)
                double precision, intent(in)  :: x(:)
            end subroutine dc

            subroutine d2c(d2cx, x)
                double precision, intent(out) :: d2cx(:, :, :)
                double precision, intent(in)  :: x(:)
            end subroutine d2c

            subroutine subproblem()
            end subroutine subproblem

            subroutine iteration()
            end subroutine iteration

            subroutine armijo()
            end subroutine armijo

            subroutine norm()
            end subroutine norm

            subroutine angle()
            end subroutine angle
        end interface

        ! Saída
        double precision, intent(out) :: x(:)

        ! Variáveis locais
        double precision :: cx(m_)
        double precision :: dcx(m_, n_)
        double precision :: gx(n_)
        double precision :: dL(n_) ! ∇L(x, -c(x)/μ) = ∇f(x) + 1/μ*∑cₖ(x)∇cₖ(x)
        integer          :: i, j, k
        logical          :: kkt

        n = n_
        m = m_
        mu = mu_
        f_ptr => f
        c_ptr => c

        ! Valores iniciais
        x(:) = x0(:)
        call c(cx, x)
        call dc(dcx, x)
        call g(gx, x)

        ! Lagrangeana
        dL(:) = gx(:)
        do k = 1, m
            do i = 1, n
                dL(i) = dL(i) + cx(k)*dcx(k, i)/mu
            end do
        end do

        write (19, *), x
        if (norm2(cx) < eps .and. norm2(dL) < eps) then
            kkt = .true.
        else
            kkt = .false.
        end if

        do while(.not. kkt)
            ! Minimiza Q(x, μ)
            call subproblem()
            call newt(x, f, g, h, c, dc, d2c, gamma, sigma, theta, eps, iteration, armijo, norm, angle)

            ! Verifica KKT
            call c(cx, x)
            call dc(dcx, x)
            call g(gx, x)

            dL(:) = gx(:)
            do k = 1, m
                do i = 1, n
                    dL(i) = dL(i) + cx(k)*dcx(k, i)/mu
                end do
            end do

            write (19, *), x
            write (38, '(e10.2, a, e10.2, a)'), norm2(cx), " & ", norm2(dL), " \\"
            if (norm2(cx) < eps .and. norm2(dL) < eps) then
                kkt = .true.
            else
                ! Atualiza μ
                mu = mu/2
            end if
        end do

        mu_ = mu
    end subroutine penalty

end module constrained


text:
! mimic NWChem tgt_sd_t_s1_1 kernel

#define simd schedule(static,1)

  implicit integer (a-z)
  l1 = 1; l2 = 1; l3 = 1; l4 = 1; l5 = 1; l6 = 1;
  u1 = 24; u2 = 24; u3 = 24; u4 = 24; u5 = 24; u6 = 24;
  call tgt_sd_t_s1_1(l1,l2,l3,l4,l5,l6, u1,u2,u3,u4,u5,u6)
  end

subroutine tgt_sd_t_s1_1(l1,l2,l3,l4,l5,l6, u1,u2,u3,u4,u5,u6)
  implicit integer (a-z)
  real a(24,24,24,24,24,24)
  real b(24,24,24,24,24,24)

  a=3.0
  b=0.0

  n1 = 24
  n2 = 24
  n3 = 24
  n4 = 24
  n5 = 24
  n6 = 24

  nn = n1*n2*n3*n4*n5*n6

!$omp target teams distribute parallel do simd
  do ii = 1, nn
     t = ii-1
     q = t / n6
     r = t - q*n6
     i6 = l6 + r

     t = q
     q = t / n5
     r = t - q*n5
     i5 = l5 + r

     t = q
     q = t / n4
     r = t - q*n4
     i4 = l4 + r

     t = q
     q = t / n3
     r = t - q*n3
     i3 = l3 + r

     t = q
     q = t / n2
     r = t - q*n2
     i2 = l2 + r

     t = q
     q = t / n1
     r = t - q*n1
     i1 = l1 + r

     b(i6,i5,i4,i3,i2,i1) = a(i6,i5,i4,i3,i2,i1) + i3
enddo
!$omp end target teams distribute parallel do

!  write(6,*) b(1,1,1,1,1,1)
!  write(6,*) a(1,1,1,1,1,1)

  write(6,*) ((b(k,j,1,1,1,1),j=1,4),k=1,4)
  return
end

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
In addition to Locally Owned Websites local websites there are a whole bunch of websites designed to facilitate social interaction.

Some specifically targeted people in Davis. They are all gone now, but they once included the following:
  DavisCoeds  dead
  UCDavisParties.com  dead
  UDavisLife.com  dead

Others have users from all over the world:
  Facebook
  Livejournal
  Meetup
  MySpace
  Nextdoor
  http://www.orkut.com Orkut  Had a fairly large UCD community at one point.
  Xanga
  Blogger
  Flickr
  Twitter
  Instagram Snap photos and share them

These succeed to various degrees in their objective.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module int8_m
  !!
  !! Provide functionality for comparing scalar integers of `kind=int8`
  !!
  use, intrinsic :: iso_fortran_env,  only : WP => int8
  use test_base_m,                    only : testline
  use test_base_m,                    only : tests
  implicit none
  private
  public :: is

  interface is
    module procedure is_int
  end interface

  ! integer, parameter :: WP = int8

contains

  include 'is_int.inc'

end module int8_m

max_stars_repo_name: jbdv-no/fortapt

max_issues_repo_name: jbdv-no/fortapt

max_forks_repo_name: jbdv-no/fortapt


module dyn_init
! initialize dynamics (including reading in initial conditions).
! also does IO of spectral data.
! public subroutines:
! init_dyn:  allocate spectral arrays (init_specdata),
!  read in initial conditions, initalize pressure variables (init_pressdata,
!  get ak,bk from IC file header), compute gradient of orography,
!  set up linear damping operators (disspec,diff_prof,damp_prof), initialize arrays
!  for semi-implicit time stepping.
! wrtout_sig: write out spectral data.
! readin_sig: read in spectral data.
 use kinds, only: r_kind, r_single, r_double
! JOY removing sigio dependency
! use sigio_module, only: sigio_sclose,sigio_swohdc,&
!  sigio_srohdc,sigio_aldata,sigio_data,sigio_head,sigio_sropen,sigio_srdata,sigio_axdata
 use params, only: &
 nlons,nlats,nlevs,ndimspec,ntrunc,initfile,dry,ndiss,efold,dcmip,polar_opt,&
 heldsuarez,explicit,tstart,idate_start,hdif_fac,hdif_fac2,fshk,ntrac,taustratdamp,&
 ntoz,ntclw,pdryini,massfix,toa_pressure
 !sighead
 !JOY removed sighead to remove sigio dependency
 use shtns, only: shtns_init, spectogrd, grdtospec, getgrad, getvrtdivspec, lap, lats, lons
 use spectral_data, only: vrtspec,divspec,virtempspec,tracerspec,topospec,lnpsspec,&
                          disspec,diff_prof,dmp_prof,init_specdata
 use pressure_data, only: ak,bk,ck,dbk,bkl,sl,si,psg,prs,pk,dpk,calc_pressdata
 use grid_data, only: init_griddata, dphisdx, dphisdy, phis, ug, vg, virtempg, &
 tracerg,lnpsg,dlnpdtg
! JOY removing stochastic routines
! use stoch_data, only: init_stochdata
 use physcons, only: rerth => con_rerth, rd => con_rd, cp => con_cp, &
                     omega => con_omega, grav => con_g, pi => con_pi, &
                     fv => con_fvirt
 use semimp_data, only: init_semimpdata
 use iso_c_binding

 implicit none
!namelist/dynamics_list/ak,bk
 private
 !public :: init_dyn, wrtout_sig, readin_sig, wrtout_gfsgrb
 !! removing all output calls
 ! JOY remove routines which read initial conditions from file
 !public :: init_dyn, wrtout_sig, readin_sig
 public :: init_dyn

 real(r_kind) :: slrd0=0
 contains
 !JOY subroutine to set topography from arbitrary array
 subroutine set_topography(surf_geop) bind(c, name='gfs_set_topography')

    real(r_kind), intent(in), dimension(nlons,nlats) :: surf_geop

    call grdtospec(surf_geop/grav, topospec)
    call spectogrd(grav*topospec, phis)
    call getgrad(grav*topospec, dphisdx, dphisdy, rerth)

 end subroutine set_topography

 subroutine init_dyn(py_damping_level, py_taustratdamp) bind(c, name='gfs_init_dynamics')
    integer k
    ! allocate arrays
    ! JOY data array inits will be done from python

    !call init_specdata()
    !call init_griddata()
    integer(c_int), intent(in):: py_damping_level
    real(c_double), intent(in):: py_taustratdamp

    ! initialize spherical harmonic lib
    !print *,'Initialising shtns'
    call shtns_init(nlons,nlats,ntrunc,nthreads=1,polar_opt=polar_opt)
        ! convert to ln(ps) in Pa.
    !print *,'calling spectogrd'
    !print *, 'Grid size', nlats, nlons, ntrunc
    !print *, shape(lnpsspec)
    !call spectogrd(lnpsspec, psg)
    !psg = 1000.*exp(psg) ! convert to Pa
    !call grdtospec(log(psg), lnpsspec) ! back to spectral.

    do k=1,nlevs
       dbk(k) = bk(k+1)-bk(k)
       bkl(k) = 0.5*(bk(k+1)+bk(k))
       ck(k)  = ak(k+1)*bk(k)-ak(k)*bk(k+1)
    enddo
         ! compute sigma coordinate quantities (bottom to top).
    do k=1,nlevs+1
        si(nlevs+2-k)= (ak(k) - toa_pressure)/(pdryini - toa_pressure)+bk(k) ! si are now sigmas
    enddo
    do k=1,nlevs
        sl(k) = 0.5*(si(k)+si(k+1))
         !sl(k) = ((si(k)**(rk+1.) - si(k+1)**(rk+1.))/&
         !        ((rk+1.)*(si(k)-si(k+1))))**(1./rk)
    enddo

    slrd0 = si(nlevs+1 - py_damping_level)
    taustratdamp = py_taustratdamp

    !call spectogrd(grav*topospec, phis)
    !print *,'min/max surface geopotential',minval(phis),maxval(phis)
    !call spectogrd(lnpsspec, psg)
    !psg = exp(psg)
    !print *,'min/max sfc pressure (hPa)',minval(psg/100.),maxval(psg/100.)
    ! initialize model interface and level pressures, related variables.
    !call calc_pressdata(lnpsg)
    ! compute gradient of surface orography
    !call getgrad(grav*topospec, dphisdx, dphisdy, rerth)
    ! hyper-diffusion operator (plus rayleigh damping in upper stratosphere)

    call setdampspec(ndiss,efold,hdif_fac,hdif_fac2,fshk,disspec,diff_prof,dmp_prof)

    ! initialize arrays for semi-implicit adjustments.
    if (.not. explicit) call init_semimpdata()
    ! initialize stochastic data.
    ! JOY removing stochastic routines
    !call init_stochdata()
 end subroutine init_dyn

 subroutine copyspecin(rspecdata,cspecdata)
    real(r_single), intent(in) :: rspecdata(2*ndimspec)
    complex(r_kind), intent(out) :: cspecdata(ndimspec)
    integer n,nn
    nn = 1
    ! factor of sqrt(2.*pi) accounts for difference in normalization
    ! between ncep libs and shtns (which uses orthonormalized norm)
    do n=1,ndimspec
       cspecdata(n) = sqrt(2.*pi)*cmplx(rspecdata(nn),rspecdata(nn+1))
       nn = nn + 2
    enddo
 end subroutine copyspecin

 subroutine copyspecout(cspecdata, rspecdata)
    real(r_single), intent(out) :: rspecdata(2*ndimspec)
    complex(r_kind), intent(in) :: cspecdata(ndimspec)
    integer n,nn
    nn = 1
    ! factor of sqrt(2.*pi) accounts for difference in normalization
    ! between ncep libs and shtns (which uses orthonormalized norm)
    do n=1,ndimspec
       rspecdata(nn) = real(cspecdata(n))/sqrt(2.*pi)
       rspecdata(nn+1) = imag(cspecdata(n))/sqrt(2.*pi)
       nn = nn + 2
    enddo
 end subroutine copyspecout

 subroutine twodtooned(data2,data1)
   real(r_kind), intent(in) :: data2(nlons,nlats)
   real(r_kind), intent(out) :: data1(nlons*nlats)
   integer i,j,n
   do n=1,nlons*nlats
      j = 1+(n-1)/nlons
      i = n-(j-1)*nlons
      data1(n) = data2(i,j)
   enddo
 end subroutine twodtooned

! subroutine dcmip_ics(testcase,dry)
!   ! DCMIP 2012 test case initial conditions (test cases 4 and 5)
!   ! (idealized baroclinic instability or trop cyclone test case initial conditions)
!   use dcmip_initial_conditions_test_4, only: test4_baroclinic_wave
!   use dcmip_initial_conditions_test_5, only: test5_tropical_cyclone
!   integer, intent(in) :: testcase
!   logical, intent(in) :: dry
!   real(8) lon,lat,z,X,p,u,v,w,t,phi_s,ps,rho,q,q1,q2
!   integer moist,zcoords,i,j,k,nt
!
!   if (testcase/10 .eq. 5) then
!      if (dry) then
!         print *,'cannot run dcmip trop cyclone test case with dry=.true.'
!         stop
!      endif
!      print *,'replacing initial conds with dcmip tropical cyclone case..'
!   else if (testcase/10 .eq. 4) then
!      print *,'replacing initial conds with dcmip baroclinic wave test case..'
!   else
!      print *,'dcmip must be 40 something or 50 something'
!      stop
!   endif
!   zcoords = 0 ! p, not z, is specified.
!   if (dry) then
!     moist = 0
!   else
!     moist = 1
!   endif
!   if (testcase/10 .eq. 4) psg = 1.e5
!   if (testcase/10 .eq. 5) then
!      psg = 1.015e5
!      lnpsg = log(psg)
!      call grdtospec(lnpsg,lnpsspec)
!      call calc_pressdata(lnpsg)
!      ! run once, just to calculate ps
!      ! ps then used to compute pressure on model levels for subsequent call
!      do k=1,nlevs
!      do j=1,nlats
!      do i=1,nlons
!         lon = lons(i,j); lat = lats(i,j); p = prs(i,j,k)
!         call test5_tropical_cyclone (lon,lat,p,z,zcoords,u,v,w,t,phi_s,ps,rho,q)
!         psg(i,j) = ps
!      enddo
!      enddo
!      enddo
!   endif
!   ! calculate pressure on model levels given ps.
!   lnpsg = log(psg)
!   call grdtospec(lnpsg,lnpsspec)
!   call calc_pressdata(lnpsg)
!   ! calculate u,v,t,q,phis and other tracers.
!   do k=1,nlevs
!      do j=1,nlats
!      do i=1,nlons
!         lon = lons(i,j); lat = lats(i,j); p = prs(i,j,k)
!         if (testcase/10 .eq. 4) then
!            call test4_baroclinic_wave (moist,X,lon,lat,p,z,zcoords,u,v,w,t,phi_s,ps,rho,q,q1,q2)
!         else if (testcase/10 .eq. 5) then
!            call test5_tropical_cyclone (lon,lat,p,z,zcoords,u,v,w,t,phi_s,ps,rho,q)
!         endif
!         ug(i,j,k) = u; vg(i,j,k) = v; phis(i,j) = phi_s
!         if (dry) then
!            virtempg(i,j,k) = t
!            tracerg(i,j,k,1) = 0.
!            tracerg(i,j,k,2) = q1
!            tracerg(i,j,k,3) = q2
!         else
!            virtempg(i,j,k) = t*(1. + fv*q)
!            tracerg(i,j,k,1) = q
!            if (testcase/10 .eq. 4) then
!               tracerg(i,j,k,2) = q1
!               tracerg(i,j,k,3) = q2
!            else if (testcase .eq. 52) then
!               tracerg(i,j,k,2:ntrac) = 0.
!            endif
!         endif
!       enddo
!       enddo
!       call getvrtdivspec(ug(:,:,k),vg(:,:,k),vrtspec(:,k),divspec(:,k),rerth)
!       call grdtospec(virtempg(:,:,k),virtempspec(:,k))
!       do nt=1,ntrac
!          call grdtospec(tracerg(:,:,k,nt),tracerspec(:,k,nt))
!       enddo
!   enddo
!   call grdtospec(phis/grav, topospec)
!
! end subroutine dcmip_ics

 subroutine heldsuarez_ics(toa_pressure)
   ! jablonowski and williamson (2006, QJR, p. 2943, doi: 10.1256/qj.06.12)
   ! initial perturbation on an isothermal state.
   real(r_kind), dimension(nlons,nlats) :: rnh,xnh,rsh,xsh, rand_vel
   real(r_kind) :: lonc,latc,up,pertrad
   real(c_double), intent(in):: toa_pressure
   integer k
   print *,'replacing initial conds with held and suarez test case..'
   lonc = pi/9.
   up = 1.
   pertrad = rerth/10.
   !JOY changing ICs
   latc = 0.
   !latc = 2.*pi/9.
   !xnh = sin(latc)*sin(lats) + cos(latc)*cos(lats)*cos(lons-lonc)
   xnh = sin(latc)*sin(lats) + cos(latc)*cos(lats)
   latc = -2.*pi/9.
   !xsh = sin(latc)*sin(lats) + cos(latc)*cos(lats)*cos(lons-lonc)
   xsh = sin(latc)*sin(lats) + cos(latc)*cos(lats)
   rnh = rerth*acos(xnh)
   rsh = 0.*rerth*acos(xsh)
   virtempg = 300.  ! isothermal state.
   vg = 0.
   psg = 1.e5
   lnpsg = log(psg)
   call grdtospec(lnpsg,lnpsspec)
   call calc_pressdata(lnpsg)
   do k=1,nlevs
      ! add a barotropic zonal wind perturbation (opp sign in each hemisphere)
      !JOY changed to single tropical jet
      call random_number(rand_vel)
      ug(:,:,k) = 0.8*rand_vel
      !ug(:,:,k) = up*(exp(-(rnh/pertrad)**2)+exp(-(rsh/pertrad)**2))
      call getvrtdivspec(ug(:,:,k),vg(:,:,k),vrtspec(:,k),divspec(:,k),rerth)
      call grdtospec(virtempg(:,:,k),virtempspec(:,k))
   enddo
   topospec = 0.
 end subroutine heldsuarez_ics

 subroutine setdampspec(ndiss,efold,hdif_fac,hdif_fac2,fshk,disspec,diff_prof,dmp_prof)
   ! set hyper-diffusion and upper level rayleigh damping parameters/structure.
   ! if efold and/or ndiss are zero, default GFS parameters are used.
   integer, intent(inout) :: ndiss
   real(r_kind), intent(inout) :: efold,fshk
   real(r_kind), intent(in) :: hdif_fac,hdif_fac2
   real(r_kind), intent(out),dimension(ndimspec) :: disspec
   real(r_kind), dimension(nlevs),intent(out) :: diff_prof,dmp_prof
   integer k
   real(r_kind) dmp_prof1
   ! if ndiss=0, use default value of 8
   if (ndiss == 0) ndiss = 8
   ! if efold <= 0, use GFS defaults.
   if (efold .le. 0) then
      if (ntrunc > 170) then
         efold = 3600./(hdif_fac2*(ntrunc/170.)**4*1.1)
      else if (ntrunc == 126 .or. ntrunc == 170) then
         efold = 1./(hdif_fac2*12.E15/(RERTH**4)*(80.*81.)**2)
      else
         efold = 1./(hdif_fac2*3.E15/(RERTH**4)*(80.*81.)**2)
      end if
      efold = 2.*efold ! mysterious factor 2 in deldifs.f
      !print *,'ndiss,efold =',ndiss,efold
   end if
   if (fshk .le. 0) then
      ! factor to multiply diffusion per scale height.
      fshk = 1.0*hdif_fac
      if (ntrunc > 170) fshk = 2.2*hdif_fac
      if (ntrunc == 126) fshk = 1.5*hdif_fac
   end if
   !slrd0=0.00        ! SIGMA LEVEL AT WHICH TO BEGIN RAYLEIGH MOMTUM DAMPING
   dmp_prof1=1./taustratdamp ! RECIPROCAL OF TIME SCALE PER SCALE HEIGHT
                      ! ABOVE BEGINNING SIGMA LEVEL FOR RAYLEIGH DAMPING
   dmp_prof = 0.
   diff_prof = sl**(log(1./fshk))
   !print *,'profiles for diffusion and linear damping:'
   !print *,'(level, sigma, diffusion enhancment, linear drag coeff)'
   do k=1,nlevs
      if (sl(k) .lt. slrd0) dmp_prof(k)=dmp_prof1*log(slrd0/sl(k))
   !   print *,k,sl(k),diff_prof(k),dmp_prof(k)
   enddo
   disspec(:) = -(1./efold)*(lap(:)/minval(lap))**(ndiss/2)
   return
 end subroutine setdampspec

end module dyn_init


text:
      subroutine getdef (token, toksiz, defn, defsiz)
      integer token (100), defn (2048)
      integer toksiz, defsiz
      common /cdefio/ bp, buf (4096)
      integer bp
      integer buf
      common /cfname/ fcname (30)
      integer fcname
      common /cfor/ fordep, forstk (200)
      integer fordep
      integer forstk
      common /cgoto/ xfer
      integer xfer
      common /clabel/ label, retlab, memflg, col, logic0
      integer label
      integer retlab
      integer memflg
      integer col
      integer logic0
      common /cline/ dbgout, dbglev, level, linect (5), infile (5), fnam
     *p, fnames ( 150)
      integer dbgout
      integer dbglev
      integer level
      integer linect
      integer infile
      integer fnamp
      integer fnames
      common /cmacro/ cp, ep, evalst (500), deftbl
      integer cp
      integer ep
      integer evalst
      integer deftbl
      common /coutln/ outp, outbuf (74)
      integer outp
      integer outbuf
      common /csbuf/ sbp, sbuf(2048), smem(240)
      integer sbp
      integer sbuf
      integer smem
      common /cswtch/ swtop, swlast, swstak(1000), swvnum, swvlev, swvst
     *k(10), swinrg
      integer swtop
      integer swlast
      integer swstak
      integer swvnum
      integer swvlev
      integer swvstk
      integer swinrg
      common /ckword/ rkwtbl
      integer rkwtbl
      common /clname/ fkwtbl, namtbl, gentbl, errtbl, xpptbl
      integer fkwtbl
      integer namtbl
      integer gentbl
      integer errtbl
      integer xpptbl
      common /erchek/ ername, body, esp, errstk(30)
      integer ername
      integer body
      integer esp
      integer errstk
      integer mem( 60000)
      common/cdsmem/mem
      integer c, t, ptoken (100)
      integer gtok, ngetch
      integer i, nlpar
      call skpblk
      c = gtok (ptoken, 100)
      if (.not.(c .eq. 40))goto 23000
      t = 40
      goto 23001
23000 continue
      t = 32
      call pbstr (ptoken)
23001 continue
      call skpblk
      if (.not.(gtok (token, toksiz) .ne. -9))goto 23002
      call baderr (22Hnon-alphanumeric name.)
23002 continue
      call skpblk
      c = gtok (ptoken, 100)
      if (.not.(t .eq. 32))goto 23004
      call pbstr (ptoken)
      i = 1
23006 continue
      c = ngetch (c)
      if (.not.(i .gt. defsiz))goto 23009
      call baderr (20Hdefinition too long.)
23009 continue
      defn (i) = c
      i = i + 1
23007 if (.not.(c .eq. 35 .or. c .eq. 10 .or. c .eq. -1))goto 23006
23008 continue
      if (.not.(c .eq. 35))goto 23011
      call putbak (c)
23011 continue
      goto 23005
23004 continue
      if (.not.(t .eq. 40))goto 23013
      if (.not.(c .ne. 44))goto 23015
      call baderr (24Hmissing comma in define.)
23015 continue
      nlpar = 0
      i = 1
23017 if (.not.(nlpar .ge. 0))goto 23019
      if (.not.(i .gt. defsiz))goto 23020
      call baderr (20Hdefinition too long.)
      goto 23021
23020 continue
      if (.not.(ngetch (defn (i)) .eq. -1))goto 23022
      call baderr (20Hmissing right paren.)
      goto 23023
23022 continue
      if (.not.(defn (i) .eq. 40))goto 23024
      nlpar = nlpar + 1
      goto 23025
23024 continue
      if (.not.(defn (i) .eq. 41))goto 23026
      nlpar = nlpar - 1
23026 continue
23025 continue
23023 continue
23021 continue
23018 i = i + 1
      goto 23017
23019 continue
      goto 23014
23013 continue
      call baderr (19Hgetdef is confused.)
23014 continue
23005 continue
      defn (i - 1) = -2
      return
      end
c     logic0  logical_column

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


!---------------------------------------------------!
! Copyright (c) 2016 Shunsuke A. Sato               !
! Released under the MIT license                    !
! https://opensource.org/licenses/mit-license.php   !
!---------------------------------------------------!
!-------10--------20--------30--------40--------50--------60--------70--------80--------90
! Two band model (3D)
module global_variables_2d
  use global_variables
  use mpi_module

  real(8) :: dir_pol_1(3),dir_pol_2(3)

end module global_variables_2d


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine lhsx(c)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c This function computes the left hand side for the three x-factors
c---------------------------------------------------------------------

       include 'header.h'

       double precision ru1
       integer          i, j, k, c

c---------------------------------------------------------------------
c      treat only cell c
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c      first fill the lhs for the u-eigenvalue
c---------------------------------------------------------------------
       do  k = start(3,c), cell_size(3,c)-end(3,c)-1
          do  j = start(2,c), cell_size(2,c)-end(2,c)-1
             do  i = start(1,c)-1, cell_size(1,c)-end(1,c)
                ru1 = c3c4*rho_i(i,j,k,c)
                cv(i) = us(i,j,k,c)
                rhon(i) = dmax1(dx2+con43*ru1,
     >                          dx5+c1c5*ru1,
     >                          dxmax+ru1,
     >                          dx1)
             end do

             do  i = start(1,c), cell_size(1,c)-end(1,c)-1
                lhs(i,j,k,1,c) =   0.0d0
                lhs(i,j,k,2,c) = - dttx2 * cv(i-1) - dttx1 * rhon(i-1)
                lhs(i,j,k,3,c) =   1.0d0 + c2dttx1 * rhon(i)
                lhs(i,j,k,4,c) =   dttx2 * cv(i+1) - dttx1 * rhon(i+1)
                lhs(i,j,k,5,c) =   0.0d0
             end do
          end do
       end do

c---------------------------------------------------------------------
c      add fourth order dissipation
c---------------------------------------------------------------------
       if (start(1,c) .gt. 0) then
          i = 1
          do   k = start(3,c), cell_size(3,c)-end(3,c)-1
             do   j = start(2,c), cell_size(2,c)-end(2,c)-1
                lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz5
                lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
                lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1

                lhs(i+1,j,k,2,c) = lhs(i+1,j,k,2,c) - comz4
                lhs(i+1,j,k,3,c) = lhs(i+1,j,k,3,c) + comz6
                lhs(i+1,j,k,4,c) = lhs(i+1,j,k,4,c) - comz4
                lhs(i+1,j,k,5,c) = lhs(i+1,j,k,5,c) + comz1
             end do
          end do
       endif

       do   k = start(3,c), cell_size(3,c)-end(3,c)-1
          do   j = start(2,c), cell_size(2,c)-end(2,c)-1
             do   i=3*start(1,c), cell_size(1,c)-3*end(1,c)-1
                lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
                lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
                lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
                lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
                lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
             end do
          end do
       end do

       if (end(1,c) .gt. 0) then
          i = cell_size(1,c)-3
          do   k = start(3,c), cell_size(3,c)-end(3,c)-1
             do   j = start(2,c), cell_size(2,c)-end(2,c)-1
                lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
                lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
                lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
                lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4

                lhs(i+1,j,k,1,c) = lhs(i+1,j,k,1,c) + comz1
                lhs(i+1,j,k,2,c) = lhs(i+1,j,k,2,c) - comz4
                lhs(i+1,j,k,3,c) = lhs(i+1,j,k,3,c) + comz5
             end do
          end do
       endif

c---------------------------------------------------------------------
c      subsequently, fill the other factors (u+c), (u-c) by a4ing to
c      the first
c---------------------------------------------------------------------
       do   k = start(3,c), cell_size(3,c)-end(3,c)-1
          do   j = start(2,c), cell_size(2,c)-end(2,c)-1
             do   i = start(1,c), cell_size(1,c)-end(1,c)-1
                lhs(i,j,k,1+5,c)  = lhs(i,j,k,1,c)
                lhs(i,j,k,2+5,c)  = lhs(i,j,k,2,c) -
     >                            dttx2 * speed(i-1,j,k,c)
                lhs(i,j,k,3+5,c)  = lhs(i,j,k,3,c)
                lhs(i,j,k,4+5,c)  = lhs(i,j,k,4,c) +
     >                            dttx2 * speed(i+1,j,k,c)
                lhs(i,j,k,5+5,c) = lhs(i,j,k,5,c)
                lhs(i,j,k,1+10,c) = lhs(i,j,k,1,c)
                lhs(i,j,k,2+10,c) = lhs(i,j,k,2,c) +
     >                            dttx2 * speed(i-1,j,k,c)
                lhs(i,j,k,3+10,c) = lhs(i,j,k,3,c)
                lhs(i,j,k,4+10,c) = lhs(i,j,k,4,c) -
     >                            dttx2 * speed(i+1,j,k,c)
                lhs(i,j,k,5+10,c) = lhs(i,j,k,5,c)
             end do
          end do
       end do

       return
       end


!> \file read_write.f90
!! \brief Using dictionary_t to store a configuration file

!> Read a "equal sign" separated input file
!!
!! Usage:
!!    ./read_write name_of_input_file [dict_size]
!!
!! dict_size is an optional integer argument that sets the size of the hash table. The
!! default value is 1024.
!!
!! Output: the program returns one line per key, with the key and the value separated by an
!! equal sign.

program use_ssdm
  use dictionary_m
  implicit none

  type(dictionary_t) :: b
  character(len=128) :: line
  character(len=:), allocatable :: k, v, fname
  integer :: iostat, eq_location, funit, dict_size

  if (command_argument_count() < 1) then
     stop 'missing argument for parameter file'
  end if

  call get_command_argument(1, line)
  fname = trim(adjustl(line))

  if (command_argument_count() < 2) then
     dict_size = 1024
  else
     call get_command_argument(2, line)
     read(line, *) dict_size
  end if

  call b%init(dict_size)

  open(file=fname, newunit=funit)
  read_loop: do
     read(funit, '(a)', iostat=iostat) line
     if (iostat < 0) exit read_loop
     if (iostat > 0) cycle read_loop
     eq_location = index(line, '=')
     if (eq_location == 0) cycle read_loop

     k = trim(adjustl(line(1:eq_location-1)))
     v = trim(adjustl(line(eq_location+1:)))

     if (len(k) ==0) then
        write(*,*) 'empty key while reading data'
        stop
     end if
     if (len(v) ==0) then
        write(*,*) 'empty value while reading data'
        stop
     end if

     call b%set(k, v)

  end do read_loop

  close(funit)

  call fileshow(b)

contains

  subroutine fileshow(d)
    class(dictionary_t), intent(in) :: d

    integer :: i, j, s

    do i = 1, d%dict_size
       s = d%buckets(i)%current_idx
       if (s > 0) then
          do j = 1, s
             write(*,'(a)', advance='no') d%buckets(i)%entries(j)%key
             write(*,'(a)', advance='no') ' = '
             write(*,'(a)', advance='no') d%buckets(i)%entries(j)%value
             write(*,*)
          end do
       end if
    end do

  end subroutine fileshow

end program use_ssdm