module putils

use, intrinsic :: iso_fortran_env

implicit none

! Fortran 2008
integer, parameter, private :: sp = REAL32
integer, parameter, private :: dp = REAL64
integer, parameter, private :: qp = REAL128

contains

      subroutine pwrite(ptcls,n1in,nrayspin,nunit,idmin,idmax,ncols,nprecision)
! routine to perform parallel write of ptcls array
      use mpi
      implicit none
!     real*8, dimension(n1,*) :: ptcls
      real*8, dimension(:,:) :: ptcls
      real*8, allocatable, dimension(:,:) :: tblock
      integer :: n1in,nrayspin,nunit,idmin,idmax,ncols,nprecision
      integer :: n1,nraysp
      character*17 :: myformat='(99(1x,1pe16.9 ))'
      character*2 :: a2
      integer :: nrayspmax,nlength,mycount,nraysw,i,l,mprocs,myrank,ierr
      integer ::  mpistat(MPI_STATUS_SIZE)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,mprocs,ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)
      n1=ubound(ptcls,1)
      nraysp=ubound(ptcls,2)
!     if(myrank.eq.0)write(6,*)'(pwrite) n1,nraysp=',n1,nraysp
!
!---format:
      nlength=nprecision+7
      call num2stringit(nlength,a2,2)
      myformat(11:12)=a2
      call num2stringit(nprecision,a2,2)
      myformat(14:15)=a2
!---
!
      if(idmax.le.nraysp)then
        do i=1,idmax
          if(myrank.eq.0)write(nunit,myformat)ptcls(1:ncols,i)
        enddo
        return
      endif
!
!------------------------------------------------------------------
      call MPI_REDUCE(nraysp,nrayspmax,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ierr)
      if(myrank.eq.0)then
        allocate(tblock(n1,nrayspmax))
        mycount=0
        do i=1,nraysp
         mycount=mycount+1
         if(mycount.lt.idmin)cycle
         if(mycount.gt.idmax)exit
         write(nunit,myformat)ptcls(1:ncols,i)
        enddo
        do l=1,mprocs-1
          call MPI_RECV(tblock,n1*nrayspmax,MPI_DOUBLE,l,98,MPI_COMM_WORLD,mpistat,ierr)
          call MPI_GET_COUNT(mpistat,MPI_DOUBLE,nraysw,ierr)
          nraysw=nraysw/n1
          do i=1,nraysw
           mycount=mycount+1
           if(mycount.lt.idmin)cycle
           if(mycount.gt.idmax)exit
           write(nunit,myformat)tblock(1:ncols,i)
          enddo
        enddo
        deallocate(tblock)
      else
        call MPI_SEND(ptcls,n1*nraysp,MPI_DOUBLE,0,98,MPI_COMM_WORLD,ierr)
      endif
      return
      end
!
      subroutine num2stringit(num,a,ndigits)
! converts an integer "num" with at most "ndigits" digits into a character string "a"
! if ndigits exceeds the actual number of digits in num, then leading zeroes are inserted in the string
      implicit none
      integer num,ndigits
      character a*(*)
      integer n,m,k
      m=num
      do n=1,ndigits
      k=m/10**(ndigits-n)
      if(k.eq.0)a(n:n)='0'
      if(k.eq.1)a(n:n)='1'
      if(k.eq.2)a(n:n)='2'
      if(k.eq.3)a(n:n)='3'
      if(k.eq.4)a(n:n)='4'
      if(k.eq.5)a(n:n)='5'
      if(k.eq.6)a(n:n)='6'
      if(k.eq.7)a(n:n)='7'
      if(k.eq.8)a(n:n)='8'
      if(k.eq.9)a(n:n)='9'
      m=m-k*10**(ndigits-n)
      enddo
      return
      end

end module putils


text:
// This is for simulation (unified hdl/hvl)

// +define+PRINTF_COND=0

// Sadly, this seems important
+define+RANDOMIZE_REG_INIT
+define+RANDOMIZE_MEM_INIT
+define+RANDOMIZE_GARBAGE_ASSIGN
+define+RANDOMIZE_INVALID_ASSIGN

-f ${ROCKET_SOC}/ve/rocket_soc_uvm/sim/scripts/rocket_soc_uvm_hdl.f
-f ${ROCKET_SOC}/ve/rocket_soc_uvm/sim/scripts/rocket_soc_uvm_hvl.f

ext: f


text: The http://wric.ucdavis.edu Weed Research & Information Center is a cooperative extension program that fosters research in weed management. They host Weed Day each year.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! ----------------------------------------- !
!      Set up filtro periodico              !
!             F. Califano, Marzo 2001       !
! ----------------------------------------- !

!**************************************************
!  MPI PARALLEL VERSION: VALENTINI-FAGANELLO 2009
!   3D PARALLEL VERSION: FAGANELLO 2010
!**************************************************

           subroutine init_filtro_per

use box_mod
use parameter_mod
use filtro_per_mod

IMPLICIT NONE

  integer  :: kk
  real(dp) :: tr1, tr2, tr3, tr4, tr5, w0, w11, w22

!   Set up coefficienti filtro periodico in x !!!!!!!!!!!!!!!!refuso!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!   filter parameters

       tr1 = (2.0 + 3.0 * omega_x) / 4.0
       tr2 = (6.0 + 7.0 * omega_x) / 8.0
       tr3 = (6.0 + omega_x)       / 20.0
       tr4 = (2.0 - 3.0 * omega_x) / 40.0
       tr5 = (3.0 - 2.0 * omega_x) / 10.0

 do kk = 1, n2x
  w0 = 2.d0 * pgreco * dfloat(kk) / dfloat(nx)
  w11 = tr1 + tr2*cos(w0) + tr3*cos(2.d0*w0) + tr4*cos(3.d0*w0)
  w22 = 1.0 +2.0 * omega_x*cos(w0) + 2.0 * tr5*cos(2.d0*w0)
  work4x(kk) = w11 / w22
 enddo

!   Set up coefficienti filtro periodico in y

!   filter parameters

       tr1 = (2.0 + 3.0 * omega_y) / 4.0
       tr2 = (6.0 + 7.0 * omega_y) / 8.0
       tr3 = (6.0 + omega_y)       / 20.0
       tr4 = (2.0 - 3.0 * omega_y) / 40.0
       tr5 = (3.0 - 2.0 * omega_y) / 10.0

 do kk = 1, n2y
  w0 = 2.d0 * pgreco * dfloat(kk) / dfloat(ny)
  w11 = tr1 + tr2*cos(w0) + tr3*cos(2.d0*w0) + tr4*cos(3.d0*w0)
  w22 = 1.0 +2.0 * omega_y*cos(w0) + 2.0 * tr5*cos(2.d0*w0)
  work4y(kk) = w11 / w22
 enddo

!   Set up coefficienti filtro periodico in z

!   filter parameters

       tr1 = (2.0 + 3.0 * omega_z) / 4.0
       tr2 = (6.0 + 7.0 * omega_z) / 8.0
       tr3 = (6.0 + omega_z)       / 20.0
       tr4 = (2.0 - 3.0 * omega_z) / 40.0
       tr5 = (3.0 - 2.0 * omega_z) / 10.0

 do kk = 1, n2z
  w0 = 2.d0 * pgreco * dfloat(kk) / dfloat(nz)
  w11 = tr1 + tr2*cos(w0) + tr3*cos(2.d0*w0) + tr4*cos(3.d0*w0)
  w22 = 1.0 +2.0 * omega_z*cos(w0) + 2.0 * tr5*cos(2.d0*w0)
  work4z(kk) = w11 / w22
 enddo

 end subroutine

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
Jump to Timeline #Navigation Navigation

  Future barber, farmer and member of City Council, Warren R. Pugh is born in Davis in 1877.
  State Senator  J. H. Harlan
  Assemblyman  D. N. Hershey
  Sheriff and tax collector  F. M. Rahm


text:
subroutine case_OldDomain(a,itask)
   !-----------------------------------------------------------------------
   !    This routine loads an old mesh
   !-----------------------------------------------------------------------
   use typre
   use Mod_GeneralCase
   use Mod_CaseVariables
   use Mod_DomainVariables
   implicit none
   class(GeneralCase), target :: a
   type(masterVariables), pointer :: m => NULL()
   type(domainVariables), pointer :: d => NULL()

   integer(ip) :: itask

   real(rp)    :: cputim2,cputim1

   m => a%caseVars%masterVars
   d => a%caseVars%domainVars

   !Initialize OldMesh
   if (itask == 1) then

      call m%cpu_start(2)%Tic
      !OldMesh Input data
      call d%OldMesh%SetMPI(m%MPIcomm,m%MPIsize,m%MPIroot,m%MPIrank)
      call d%OldMesh%SetReadType(m%ReadTypeString)
      call d%OldMesh%SetInputFile(m%oldnamda)
      call d%OldMesh%SetInputFolder(m%OldDataFolder)
      call d%OldMesh%SetOutputFolder(m%ResultsFolder)
      call d%OldMesh%SetOutputFiles(m%lun_memor,m%lun_outpu)
      call d%OldMesh%SetParallelLibrary(m%ParallelLibrary)
      call d%OldMesh%SetFlush(m%kfl_flush)

      !Read data
      call d%OldMesh%Readat

      !Initial operations (Graph and Scatter)
      call d%OldMesh%Initialize

      call m%cpu_start(2)%Toc

   !Deallocate Read Globals
   elseif (itask == 2) then
      call m%cpu_start(2)%Tic

      call d%OldMesh%DeallocateReadGlobals

      call m%cpu_start(2)%Toc

   !Deallocate Locals
   elseif (itask == 3) then
      call d%OldMesh%Turnof

   endif

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


module histogram_module

  use multifab_module
  use ml_layout_module
  use probin_common_module, only: n_cells
  use probin_reactdiff_module, only: nspecies

  implicit none

  private
  public :: init_histogram, destroy_histogram, write_histogram

  type(multifab), save ::  n_serial
  type(layout), save :: la_serial

contains

  subroutine init_histogram(mla)

    type(ml_layout), intent(in   ) :: mla

    type(box)  :: bx_serial
    type(boxarray)  :: bxa_serial
    integer :: lo(mla%dim), hi(mla%dim)
    integer :: dm

    dm = mla%dim

    ! first set lo and hi to the problem domain
    lo = 0
    hi(1:dm) = n_cells(1:dm) - 1

    ! build n_serial
    ! set bx_serial to reduced dimension problem domain
    call box_build_2(bx_serial,lo,hi)
    ! build a boxarray containing only one box
    call boxarray_build_bx(bxa_serial,bx_serial)
    call layout_build_ba(la_serial, bxa_serial, bx_serial, mla%pmask, &
                         explicit_mapping=(/parallel_IOProcessorNode()/) )
    call destroy(bxa_serial)
    call multifab_build(n_serial, la_serial, nspecies, 0)

  end subroutine init_histogram

  subroutine destroy_histogram()

    call multifab_destroy(n_serial)
    call destroy(la_serial)

  end subroutine destroy_histogram

  subroutine write_histogram(n_in,time,step)

    type(multifab) , intent(in) :: n_in(:)
    real(kind=dp_t), intent(in) :: time
    integer        , intent(in) :: step

    real(kind=dp_t), pointer :: np(:,:,:,:)
    integer :: i,j,k

    character(len=12) :: hist_name

    write(unit=hist_name,fmt='("hist",i8.8)') step

    OPEN(100,FILE=hist_name)

    call copy(n_serial,1,n_in(1),1,nspecies)

    if (parallel_IOProcessor()) then

       np => dataptr(n_serial,1)

       do k=lbound(np,3), ubound(np,3)
          do j=lbound(np,2), ubound(np,2)
             do i=lbound(np,1), ubound(np,1)
                write(100,*) real(time), real(np(i,j,k,:))
             end do
          end do
       end do

    end if

  end subroutine write_histogram

end module histogram_module


text:
program main

!*****************************************************************************80
!
!! MAIN is the main program for HEATED_PLATE_OPENMP.
!
!  Discussion:
!
!    This code solves the steady state heat equation on a rectangular region.
!
!    The sequential version of this program needs approximately
!    18/eps iterations to complete.
!
!
!    The physical region, and the boundary conditions, are suggested
!    by this diagram;
!
!                   W = 0
!             +------------------+
!             |                  |
!    W = 100  |                  | W = 100
!             |                  |
!             +------------------+
!                   W = 100
!
!    The region is covered with a grid of M by N nodes, and an N by N
!    array W is used to record the temperature.  The correspondence between
!    array indices and locations in the region is suggested by giving the
!    indices of the four corners:
!
!                  I = 0
!          [0][0]-------------[0][N-1]
!             |                  |
!      J = 0  |                  |  J = N-1
!             |                  |
!        [M-1][0]-----------[M-1][N-1]
!                  I = M-1
!
!    The steady state solution to the discrete heat equation satisfies the
!    following condition at an interior grid point:
!
!      W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
!
!    where "Central" is the index of the grid point, "North" is the index
!    of its immediate neighbor to the "north", and so on.
!
!    Given an approximate solution of the steady state heat equation, a
!    "better" solution is given by replacing each interior point by the
!    average of its 4 neighbors - in other words, by using the condition
!    as an ASSIGNMENT statement:
!
!      W[Central]  <=  (1/4) * ( W[North] + W[South] + W[East] + W[West] )
!
!    If this process is repeated often enough, the difference between
!    successive estimates of the solution will go to zero.
!
!    This program carries out such an iteration, using a tolerance specified by
!    the user, and writes the final estimate of the solution to a file that can
!    be used for graphic processing.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    18 October 2011
!
!  Author:
!
!    Original FORTRAN90 version by Michael Quinn.
!    This version by John Burkardt.
!
!  Reference:
!
!    Michael Quinn,
!    Parallel Programming in C with MPI and OpenMP,
!    McGraw-Hill, 2004,
!    ISBN13: 978-0071232654,
!    LC: QA76.73.C15.Q55.
!
!  Local parameters:
!
!    Local, real ( kind = 8 ) DIFF, the norm of the change in the solution from
!    one iteration to the next.
!
!    Local, real ( kind = 8 ) MEAN, the average of the boundary values, used
!    to initialize the values of the solution in the interior.
!
!    Local, real ( kind = 8 ) U(M,N), the solution at the previous iteration.
!
!    Local, real ( kind = 8 ) W(M,N), the solution computed at the latest
!    iteration.
!
  use omp_lib

  implicit none

  integer ( kind = 4 ), parameter :: m = 500
  integer ( kind = 4 ), parameter :: n = 500

  real ( kind = 8 ) diff
  real ( kind = 8 ) :: eps = 0.001D+00
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iterations
  integer ( kind = 4 ) iterations_print
  integer ( kind = 4 ) j
  real ( kind = 8 ) mean
  real ( kind = 8 ) u(m,n)
  real ( kind = 8 ) w(m,n)
  real ( kind = 8 ) wtime

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'HEATED_PLATE_OPENMP'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) &
    '  A program to solve for the steady state temperature distribution'
  write ( *, '(a)' ) '  over a rectangular plate.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8,a,i8,a)' ) '  Spatial grid of ', m, ' by ', n, ' points.'
  write ( *, '(a,g14.6)' ) &
    '  The iteration will repeat until the change is <= ', eps
  write ( *, '(a,i8)' ) &
    '  The number of processors available = ', omp_get_num_procs ( )
  write ( *, '(a,i8)' ) &
    '  The number of threads available    = ', omp_get_max_threads ( )
!
!  Set the boundary values, which don't change.
!
!  OpenMP Note:
!  You CANNOT set MEAN to zero inside the parallel region.
!
  mean = 0.0D+00
!
!$omp parallel shared ( w ) private ( i, j )

  !$omp do
  do i = 2, m - 1
    w(i,1) = 100.0D+00
    w(i,n) = 100.0D+00
  end do
  !$omp end do

  !$omp do
  do j = 1, n
    w(m,j) = 100.0D+00
    w(1,j) =   0.0D+00
  end do
  !$omp end do
!
!  Average the boundary values, to come up with a reasonable
!  initial value for the interior.
!
  !$omp do reduction ( + : mean )
  do i = 2, m - 1
    mean = mean + w(i,1) + w(i,n)
  end do
  !$omp end do

  !$omp do reduction ( + : mean )
  do j = 1, n
    mean = mean + w(1,j) + w(m,j)
  end do
  !$omp end do

!$omp end parallel
!
!  OpenMP note:
!  You cannot normalize MEAN inside the parallel region.  It
!  only gets its correct value once you leave the parallel region.
!  So we interrupt the parallel region, set MEAN, and go back in.
!
  mean = mean / dble ( 2 * m + 2 * n - 4 )
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  MEAN = ', mean
!
!  Initialize the interior solution to the mean value.
!
!$omp parallel shared ( mean, w ) private ( i, j )

  !$omp do
  do j = 2, n - 1
    do i = 2, m - 1
      w(i,j) = mean
    end do
  end do
  !$omp end do

!$omp end parallel
!
!  Iterate until the  new solution W differs from the old solution U
!  by no more than EPS.
!
  iterations = 0
  iterations_print = 1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' Iteration  Change'
  write ( *, '(a)' ) ' '

  wtime = omp_get_wtime ( )

  diff = eps

  do while ( eps <= diff )
!
!  OpenMP node: You CANNOT set DIFF to 0.0 inside the parallel region.
!
    diff = 0.0D+00

!$omp parallel shared ( u, w ) private ( i, j )

    !$omp do
    do j = 1, n
      do i = 1, m
        u(i,j) = w(i,j)
      end do
    end do
    !$omp end do

    !$omp do
    do j = 2, n - 1
      do i = 2, m - 1
        w(i,j) = 0.25D+00 * ( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) )
      end do
    end do
    !$omp end do

    !$omp do reduction ( max : diff )
    do j = 1, n
      do i = 1, m
        diff = max ( diff, abs ( u(i,j) - w(i,j) ) )
      end do
    end do
    !$omp end do

!$omp end parallel

    iterations = iterations + 1

    if ( iterations == iterations_print ) then
      write ( *, '(2x,i8,2x,g14.6)' ) iterations, diff
      iterations_print = 2 * iterations_print
    end if

  end do

  wtime = omp_get_wtime ( ) - wtime

  write ( *, '(a)' ) ' '
  write ( *, '(2x,i8,2x,g14.6)' ) iterations, diff
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Error tolerance achieved.'
  write ( *, '(a,g14.6)' ) '  Wall clock time = ', wtime
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'HEATED_PLATE_OPENMP:'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end

max_forks_repo_name: vlkale/lw-sched


!***********************************/
!	Name:ローカル格子同士の接触判定プログラム
!	Alias:OContactDetermination
!	Description:ローカル格子の接触を2段階で判定する
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2017.12.20
!	Update:2017.12.20
!	Other:
!***********************************/
subroutine OContactDetermination(UG1,UG2,MG1,MG2,iTimeCombination)
    use StructVar_Mod
!    use LoopVar_Mod
    implicit none
    type(UnstructuredGrid), intent(inout) :: UG1,UG2 !Target,Donor
    type(MoveGrid), intent(in) :: MG1,MG2 !Target,Donor
    integer, intent(in) :: iTimeCombination !N+1段階とN+1段階のとき0，N段階が入っていると1
    integer :: iToF, BoundarySphereMethod

!    write(6,*) "Determination ""CONTACT"" between TargetGrid and DonorGrid"
!	1.Lp格子と，n+1段階のLx格子(x<p)またはn段階のLy格子(p=<y)格子について接触判定を行う
!   L1格子は2~N，L2格子は3~Nのように接触判定を行う
    iToF = BoundarySphereMethod(UG1%RoughRadius,UG2%RoughRadius,MG1%RC,MG2%RC,iTimeCombination)

    !if(iToF == 1) then
        !N+1段階とN+1段階の接触判定ならiKind=0
        !N+1段階とN段階の接触判定ならiKind=1
        !call GJKAlgorithm(iTimeCombination,MG1,MG2,UG1,UG2,iToF)
        !どうも間違いがあるようなのであとで修正する(自格子同士のときは問題なく接触判定できてるので，あからさまに接触していないケースが怪しいのだが...)

    !end if

    if(iToF == 1) then !格子が重なっているとき
        if(iTimeCombination == 0) then !格子番号的にUG2 < UG1
            UG1%Overlap(iTimeCombination+1,UG2%Number) = 1 !n+1段階の上位格子と被っているとき1
            UG2%Overlap(iTimeCombination+1,UG1%Number) = 2 !一括判定を行う内部境界を持ち，下側から参照する側として

        else !格子番号的にUG1 <= UG2 N+1とNの接触判定
            UG1%Overlap(iTimeCombination+1,UG2%Number) = 3 !n段階の下位格子と被っているとき2
        end if
    else
        UG1%Overlap(iTimeCombination+1,UG2%Number) = 0 !格子が重なっていない(or考慮しなくてよい)とき0
        if(iTimeCombination == 0) then !格子番号的にUG2 < UG1
            UG2%Overlap(iTimeCombination+1,UG1%Number) = 0 !一括判定を行う内部境界を持ち，下側から参照する側として
        end if
    end if

return
end subroutine OContactDetermination


text:
      subroutine filterSD(n,m,x,al,f,fmin,cstype,bl,bu,ws,lws,v,nv,
     *  maxa,maxla,maxu,maxiu,kmax,maxg,rho,htol,rgtol,maxit,iprint,
     *  nout,ifail)
      implicit double precision (a-h,o-z)
      dimension x(*),al(*),bl(*),bu(*),ws(*),lws(*),v(*)
      character cstype(*)

c  Copyright (C) 2010 Roger Fletcher

c  Current version dated 5 October 2011

c  THE ACCOMPANYING PROGRAM IS PROVIDED UNDER THE TERMS OF THE ECLIPSE PUBLIC
c  LICENSE ("AGREEMENT"). ANY USE, REPRODUCTION OR DISTRIBUTION OF THE PROGRAM
c  CONSTITUTES RECIPIENT'S ACCEPTANCE OF THIS AGREEMENT

c  Solves an NLP problem of the form: find a local solution x to
c
c          minimize    f(x)
c                            [  x   ]
c          subject to  bl <= [      ] <= bu
c                            [ c(x) ]
c
c  where f(x) is a given function of n variables x and c(x) is a vector of m
c  given constraint functions. f(x) is minimized subject to lower and upper
c  bounds bl and bu on x and c(x). f(x) and c(x) are defined by a user supplied
c  subroutine 'functions'. The user is also required to supply a subroutine
c  'gradients' which calculates gradients of f(x) and c(x) with respect to x.
c
c  If linearized constraints and the trust region are incompatible, the code
c  enters 'phase 1' in which an 'l1 feasibility problem' is solved. If this
c  is unsuccessful in resolving the situation then the code exits with
c  ifail=3 and returns a 'locally infeasible point' in x.
c
c  Parameter List  (variables in a line starting with C must be set on entry)
c  ==============
C  n     number of variables
C  m     number of general constraints
C  x(n+m)  x(1:n) stores the vector of variables. Initially an estimate of the
c            solution must be set, replaced by the solution (if it exists)
c            on exit. The rest of x is workspace
c  al(n+m) stores Lagrange multipliers at the solution on exit. A positive
c            multiplier indicates that the lower bound is active, and a
c            negative multiplier indicates that the upper bound is active.
c            Inactive constraints have a zero multiplier.
c  f     returns the value of f(x) when x is a feasible solution
C  fmin  set a strict lower bound on f(x) for feasible x (used to identify
c          an unbounded NLP)
c  cstype(m) character workspace: if ifail=3, cstype indicates constraints
c              that are infeasible in the L1 solution. cstype(i)='A' if the
c              lower bound on c/s i is infeasible, 'Z' if the upper bound is
c              infeasible, else 'N' if feasible.
C  bl(n+m) lower bounds on x and c(x) (use numbers no less than -ainfty
c            (see below) and where possible supply realistic bounds on x)
C  bu(n+m) upper bounds on x and c(x) (use numbers no greater than ainfty)
c  ws(*)   double precision workspace
c  lws(*)  integer workspace
C  v(maxg) stores nv Ritz values (estimates of eigenvalues of reduced Hessian)
C            supply the setting from a previous run of filterSD, or set
C            nv=1 and v(1)=1.D0 in absence of other information
C  nv      number of values set in v
C  maxa    maximum number of entries in the Jacobian a(*) set by gradients
C  maxla   number of entries required for sparse matrix indices and pointers
C            la(0:*) to be set up in lws(*) (maxla>=maxa+m+3).
C            Set maxla=1 if using dense matrix format
C  maxu    length of workspace user(*) passed through to user subroutines
C            'functions' and 'gradients'
C  maxiu   length of workspace iuser(*) passed through to user subroutines
C  kmax    maximum dimension of null space allowed for (kmax<=n)
C  maxg    maximum number of reduced gradient vectors stored by the
C            limited memory method (typically 6 or 7)
C  rho     initial trust region radius (typically 1.D1)
C  htol    tolerance allowed in sum h of constraint feasibilities (e.g. 1.D-6)
C  rgtol   tolerance allowed in reduced gradient l2 norm (typically 1.D-4)
C  maxit   maximum number of major iterations allowed
C  iprint  verbosity of printing (0=none, 1=one line per iteration,
C            2=additional text information given)
C  nout    output channel for printing
c  ifail   returns failure indication as follows
c              0 = successful run
c              1 = unbounded NLP (f <= fmin at an htol-feasible point)
c              2 = bounds on x are inconsistent
c              3 = local minimum of feasibility problem and h > htol
c                  (nonlinear constraints are locally inconsistent)
c              4 = initial point x has h > ubd (reset ubd or x and re-enter)
c              5 = maxit major iterations have been carried out
c              6 = termination with rho <= htol
c              7 = not enough workspace in ws or lws (see message)
c              8 = insufficient space for filter (increase mxf and re-enter)
c             >9 = unexpected fail in LCP solver (10 has been added to ifail)

c  User Routines
c  =============
c  The user must supply two subroutines to calculate f(x), c(x) and their
c  first derivatives as follows
c
c     subroutine functions(n,m,x,f,c,user,iuser)
c     implicit double precision (a-h,o-z)
c     dimension x(*),c(*),user(*),iuser(*)
c     ...
c     Statements to calculate f(x) and the m-vector c(x). The user is
c     responsible for ensuring that any failures such as IEEE errors
c     (overflow, NaN's etc.) are trapped and not returned to filterSD.
c     The same holds for gradients.
c     ...
c     return
c     end
c
c     subroutine gradients(n,m,x,a,user,iuser)
c     implicit double precision(a-h,o-z)
c     dimension x(*),a(*),user(*),iuser(*)
c     ...
c     Statements to calculate gradients of f(x) and c(x) and set in a(*).
c     The column vector grad(f) must be followed by the column vectors
c     grad(c_i), i=1,2,...,m, in the one dimensional array a(*). Either a
c     dense or sparse data structure may be used. If using the sparse data
c     structure, only stucturally non-zero entries are set. Pointers etc. for
c     the data structure are set once and for all in lws as described below.
c     The user may assume that a call of 'gradients' immediately follows one
c      of 'functions' with the same vector x.)
c     ...
c     return
c     end
c
c  The user must also supply a driver routine which calls filterSD. This must
c  set parameters and common blocks of filterSD as appropriate.
c  Space for x,al,bl,bu,ws,lws,v and cstype must be assigned.
c  If using the sparse data structure for setting gradients, indices and
c  pointers la(0:maxla-1) must be set in the driver in lws, immediately
c  following any user workspace in lws(1:maxiu). No changes in this
c  data structure are allowed during the operation of filterSD.
c  More details of the format of a(*) and la(*) are given in the file sparseA.f
c  For dense format just set maxla=1 and set lws(maxiu+1) to the
c  'stride' (>=n) used in setting the columns of grad(f) and grad(c_i).
c  For efficiency, constant entries in the gradients may be set in the driver.
c  However two copies of the gradients are kept by filterSD . These reside in
c  ws(maxu+1:maxu+maxa) and ws(maxu+maxa+1:maxu+2*maxa). Any constant entries
c  must be set in both copies.

C  Common blocks
c  =============
c      common/wsc/kk,ll,kkk,lll,mxws,mxlws
c  The user must specify the length of the workspace arrays ws(*) and lws(*)
c  in mxws and mxlws respectively. It may not be easy to specify a-priori how
c  large these arrays should be. Set a suitable large estimate, and filterSD
c  will prompt if larger values are required. As a guide, ws(*) contains
c  first user workspace, then workspace for filterSD, then workspace for glcpd,
c  and finally workspace for denseL.f or schurQR.f. lws(*) contains user
c  workspace, then maxla+n+m+mlp locations for filterSD and additional locationsc  for denseL.f or schurQR.f.
c     common/defaultc/ainfty,ubd,mlp,mxf
c  Default values of some control parameters are set here. ainfty is used to
c  represent infinity. ubd provides an upper bound on the allowed constraint
c  violation. mlp is the maximum length of arrays used in degeneracy control.
c  mxf is the maximum length of filter arrays. Default values are 1.D20, 1.D4,
c  50,  50 respectively.
c     common/ngrc/mxgr
c  The user can limit the time spent in each call of the LCP solver by setting
c  an upper limit on the number of gradient calls in mxgr (default=1000000)
c     common/mxm1c/mxm1
c  When using denseL.f, mxm1 must be set to the maximum number of general
c  constraints allowed in the active set. mxm1=min(m+1,n) is always sufficient
c     common/epsc/eps,tol,emin
c     common/repc/sgnf,nrep,npiv,nres
c  These common blocks provide default parameters that control glcpd
c     common/statsc/dnorm,h,hJt,hJ,ipeq,k,itn,nft,ngt
c  This common block returns information about the outcome of filterSD.
c  dnorm=final step length, h=final c/s violation, hJt=ditto for 'N' c/s,
c  hJ=ditto for 'A' and 'Z' c/s, ipeq=number of active equations, k=number of
c  free variables, itn=number of iterations, nft=total number of function
c  calls, ngt=total number of gradient calls

      common/wsc/kk,ll,kkk,lll,mxws,mxlws
      common/defaultc/ainfty,ubd,mlp,mxf
      common/functc/fxd,alc,m_,iph,last1,next1,nx,nx1,
     *  nal,nal1,naal,naal1,nxd,nxd1,ncx,ncx1,ncxd,ncxd1,nla1

      m_=m
      nm=n+m
c  set real storage map for ws
c  first maxu locations are user storage for functions and gradients
c  vectors required by funct: two slots of length maxa for a(*)
      last1=maxu+1
      next1=last1+maxa
c  slot of length n for x
      nx1=next1+maxa
      nx=nx1-1
c  slot of length m for lambda
      nal=nx+n
      nal1=nal+1
c  slot of length n for Ak.al, or for storing x
      naal=nal+m
      naal1=naal+1
c  slot of length n for x at x+d
      nxd=naal+n
      nxd1=nxd+1
c  slot of length m for c at x
      ncx=nxd+n
      ncx1=ncx+1
c  slot of length m for c at x+d
      ncxd=ncx+m
      ncxd1=ncxd+1
c  local storage for filter_SD
c  slot of length n for d
      id1=ncxd1+m
c  slot of length n+m for dl
      idl1=id1+n
c  slot of length n+m for du
      idu1=idl1+nm
c  slot of length n for g
      ig1=idu1+nm
c  slot of length n+m for e
      ie1=ig1+n
c  slot of length mlp for alp
      ialp1=ie1+nm
c  slot of length mxf for filh
      ifilh1=ialp1+mlp
c  slot of length mxf for filf
      ifilf1=ifilh1+mxf
c  total length of ws so far
      kk=ifilf1+mxf-1

c  set integer storage map for lws
c  first maxiu locations are user storage for functions and gradients
c  storage of length maxla for la(0:*)
      nla1=maxiu+1
c  local storage for filter_SD
c  slot of length n+m for ls
      ils1=nla1+maxla
      ils=ils1-1
c  slot of length mlp for lp
      ilp1=ils1+nm
c  total length of lws so far
      ll=ilp1+mlp-1

      do i=1,n
        if(bl(i).gt.bu(i))then
          if(iprint.gt.1)write(nout,*)'simple bounds infeasible'
          ifail=2
          return
        endif
        ws(nx+i)=min(max(bl(i),x(i)),bu(i))
      enddo

c  note x and al are just used as workspace: the true values are those in ws

      call filter_SD(n,f,fmin,cstype,bl,bu,ws,lws,v,nv,
     *  maxa,kmax,maxg,
     *  ws(id1),ws(idl1),ws(idu1),ws(ig1),x,al,ws(ie1),lws(ils1),
     *  ws(ialp1),lws(ilp1),ws(ifilh1),ws(ifilf1),rho,htol,rgtol,
     *  maxit,iprint,nout,ifail)

      if(ifail.ge.7)return
c  scatter ws(nx.. and ws(nal.. and bound multipliers into x and al
      do i=1,n
        al(i)=0.D0
      enddo
      do i=1,m
        al(n+i)=ws(nal+i)
      enddo
      do j=1,n
        i=abs(lws(ils+j))
        if(i.le.n)then
          if(ws(nx+i).eq.bl(i))then
            al(i)=x(i)
          elseif(ws(nx+i).eq.bu(i))then
            al(i)=-x(i)
          endif
        endif
      enddo
      do i=1,n
        x(i)=ws(nx+i)
      enddo

      return
    1 format(A,15I5)
      end

      subroutine filter_SD(n,f,fmin,cstype,bl,bu,ws,lws,v,nv,
     *  maxa,kmax,maxg,
     *  d,dl,du,g,r,w,e,ls,alp,lp,filh,filf,rho,htol,rgtol,
     *  maxit,iprint,nout,ifail)
      implicit double precision (a-h,o-z)
      dimension bl(*),bu(*),ws(*),lws(*),v(*),
     *  d(*),dl(*),du(*),g(*),r(*),w(*),e(*),ls(*),alp(*),lp(*),
     *  filh(*),filf(*)
      character cstype(*)

      parameter (sigma=1.D-1,infty=100000000)

      common/defaultc/ainfty,ubd,mlp,mxf
      common/epsc/eps,tol,emin
      common/repc/sgnf,nrep,npiv,nres
      common/functc/fxd,alc,m,iph,last1,next1,nx,nx1,
     *  nal,nal1,naal,naal1,nxd,nxd1,ncx,ncx1,ncxd,ncxd1,nla1
      common/infoc/rgnorm,vstep,iter,npv,nfn,ngr
      common/ngrc/mxgr
      common/statsc/dnorm,h,hJt,hJ,ipeq,k,itn,nft,ngt

    1 format(A,15I5)
    2 format(A,6E15.7)
    3 format(A/(15I5))
c   4 format(A/(6E13.5))
    4 format(A/(5E15.7))
    5 format((6E15.7))
    6 format(A,2E15.7,I2)
 1000 format(I4,1X,E14.6,E16.8,'  < reset J          ',11X,E12.4)
 1001 format(I4,1X,E14.6,E16.8,'  < LCP',E13.5,2E12.4)
 1002 format(I4,1X,E14.6,E16.8,'  < project')
 2000 format(I4,E14.6,E16.8,'  << feasible LP    ',13X,E12.4)
 2001 format(I4,E14.6,E16.8,'  << LCP',E13.5,2E12.4)
 2002 format(I4,E14.6,E16.8,'  << project')

      n1=n+1
      nm=n+m
      mode=0
      nrep=0
      iph=2
      nfil=0
      itn=0
      nft=0
      ngt=0

c  evaluate f,c and a
      call functions(n,m,ws(nx1),f,ws(ncx1),ws,lws)
      call gradients(n,m,ws(nx1),ws(last1),ws,lws)
c  evaluate h
      h=0.D0
      do i=1,m
        ci=ws(ncx+i)
        h=h+max(0.D0,bl(n+i)-ci,ci-bu(n+i))
      enddo
      if(h.gt.ubd)then
        if(iprint.gt.1)write(nout,2)'h.gt.ubd: h =',h
        ifail=4
        return
      endif
      if(iprint.ge.1)
     *  write(nout,*)' itn     h/hJt          f/hJ      ',
     *  '           rgnorm       dnorm        rho'
   10 continue
      if(rho.lt.htol)then
        if(iprint.gt.1)write(nout,2)'rho less than htol: rho =',rho
        ifail=6
        return
      endif
c     print 4,'x =',(ws(nx+i),i=1,n)
c     print 4,'c =',(ws(ncx+i),i=1,m)
c  set up LP subproblem
      do i=1,n
        dl(i)=max(-rho,bl(i)-ws(nx+i))
        du(i)=min(rho,bu(i)-ws(nx+i))
        d(i)=0.D0
      enddo
      if(abs(iph).eq.1)then
        do j=1,n
          i=abs(ls(j))-n
          if(i.gt.0)then
            if(cstype(i).eq.'A')then
              ls(j)=n+i
            elseif(cstype(i).eq.'Z')then
              ls(j)=-n-i
            endif
          endif
        enddo
      endif
      do i=1,m
        ci=ws(ncx+i)
        dl(n+i)=bl(n+i)-ci
        du(n+i)=bu(n+i)-ci
        ws(nal+i)=0.D0
        cstype(i)='N'
      enddo
      iph_=iph
      iph=0
      k=0
c     print *,'solve LP subproblem',itn
      iii=0
c     if(itn.eq.14)iii=1
      call glcpd(n,m,k,kmax,maxg,ws(last1),lws(nla1),d,dl,du,phi,
     *  -ainfty,g,r,w,e,ls,alp,lp,mlp,ipeq,ws,lws,cstype,v,nv,rgtol,
     *  mode,ifail,infty,iii,0)
c     print 4,'d =',(d(i),i=1,n)
c     print 4,'r =',(r(i),i=1,nm)
c     print 3,'ls =',(ls(i),i=1,nm)
c     print 1,'ipeq,k,ifail',ipeq,k,ifail
      mode=2
      iph=iph_
      if(ifail.eq.0.or.ifail.eq.4)then
        if(iprint.ge.1)write(nout,2000)itn,h,f,rho
        if(abs(iph).eq.1)nfil=nfil1-1
        iph=2
        goto50
      elseif(ifail.ne.3)then
c       print 1,'itn =',itn
        if(iprint.gt.1)write(nout,*)'unexpected fail in LP subproblem'
        goto99
      endif
   15 continue
      if(h.le.htol)then
        if(iprint.gt.1)write(nout,*)'htol-feasible but LP is infeasible'
        ifail=0
        return
      endif
c  infeasibility: enter feasibility restoration
c     print *,'LP is infeasible: solve l1 subproblem',itn
      iii=0
c     if(itn.eq.9)iii=2
      call l1sold(n,m,k,kmax,maxg,ws(last1),lws(nla1),d,dl,du,phi,
     *  g,r,w,e,ls,alp,lp,mlp,ipeq,ws,lws,cstype,v,nv,rgtol,ifail,iii,0)
c     print 4,'d1 =',(d(i),i=1,n)
c     print 4,'r =',(r(i),i=1,nm)
c     print 3,'ls =',(ls(i),i=1,nm)
      if(ifail.ne.0)then
        if(print.gt.1)print *,'unexpected fail in l1 subproblem'
        goto99
      endif
      if(abs(iph).eq.2)then
        call addfil(h,f,filh,filf,1,nfil,mxf,ifail)
        if(ifail.gt.0)return
        nfil1=nfil+1
      endif
c  relax infeasible c/s
      hJ=0.D0
      hJt=0.D0
      do j=1,n
        i=abs(ls(j))
        if(i.gt.n)hJt=hJt+max(0.D0,dl(i),-du(i))
      enddo
      do j=n1,nm
        i=abs(ls(j))
        if(i.gt.n)then
          if(r(i).lt.0.D0)then
            hJ=hJ+max(0.D0,dl(i),-du(i))
            if(ls(j).ge.0)then
              du(i)=dl(i)
              dl(i)=-ainfty
              cstype(i-n)='A'
            else
              dl(i)=du(i)
              du(i)=ainfty
              cstype(i-n)='Z'
            endif
          else
            hJt=hJt+max(0.D0,dl(i),-du(i))
          endif
        endif
      enddo
c     print *,'phase 1 filter entries followed by (hJt,hJ)'
c     do i=nfil1,nfil
c       print 5,filh(i),filf(i)
c     enddo
c     print 5,hJt,hJ
c     print *,'cstype = ',(cstype(i),i=1,m)
c     print 2,'hJt,hJ',hJt,hJ
      if(iprint.ge.1)write(nout,1000)itn,hJt,hJ,rho
c     if(hJt.gt.tol)then
c       call addfil(hJt,hJ,filh,filf,nfil1,nfil,mxf,ifail)
c       if(ifail.gt.0)return
        call testfil(hJt,hJ,filh,filf,nfil1,nfil,ifail)
        if(ifail.eq.1)then
          if(iprint.gt.1)write(nout,*)'l1 solution not acceptable'
          dnorm=0.D0
          do i=1,n
            dnorm=max(dnorm,abs(d(i)))
          enddo
          rho=5.D-1*dnorm
          goto10
        endif
c     endif
c  collect multipliers from l1 subproblem
      do j=1,n
        i=abs(ls(j))
        if(i.gt.n)then
          if(ls(j).gt.0)then
            ws(nal+i-n)=r(i)
          else
            ws(nal+i-n)=-r(i)
          endif
        endif
        ws(naal+j)=0.D0
      enddo
c     print 4,'al =',(ws(i),i=nal1,nal+m)
      do i=1,m
        call saipy(ws(nal+i),ws(last1),lws(nla1),i,ws(naal1),n)
      enddo
   20 continue
      if(itn.eq.maxit)then
        if(iprint.gt.1)write(nout,*)'itn.ge.maxit'
        ifail=5
        return
      endif
      iph=1
c     mode=2
      k=0
c  solve LCP subproblem
c     print 4,'x =',(ws(i),i=nx1,nx+n)
c     print 4,'c =',(ws(i),i=ncx1,ncx+m)
c     print 4,'al =',(ws(i),i=nal1,nal+m)
      do i=1,n
        dl(i)=max(-rho,bl(i)-ws(nx+i))
        du(i)=min(rho,bu(i)-ws(nx+i))
        d(i)=0.D0
      enddo
      alc=scpr(0.D0,ws(nal1),ws(ncx1),m)
c     print *,'solve phase 1 LCP subproblem',itn
      iii=0
c     if(itn.eq.11)iii=2
      call glcpd(n,m,k,kmax,maxg,ws(last1),lws(nla1),d,dl,du,phi,
     *  fmin,g,r,w,e,ls,alp,lp,mlp,ipeq,ws,lws,cstype,v,nv,rgtol,
     *  mode,ifail,mxgr,iii,0)
      nft=nft+nfn
      ngt=ngt+ngr
c     print 1,'nfn,ngr',nfn,ngr
c     print 4,'d =',(d(i),i=1,n)
c     print 4,'r =',(r(i),i=1,nm)
c     print 3,'ls =',(ls(i),i=1,nm)
c     print 1,'ipeq,k,ifail =',ipeq,k,ifail
      itn=itn+1
      dnorm=0.D0
      do i=1,n
        dnorm=max(dnorm,abs(d(i)))
      enddo
c     print 2,'dnorm,rho',dnorm,rho
      if(ifail.eq.3)then
        if(iprint.gt.1)
     *    write(nout,*)'phase 1 LCP problem is infeasible'
        if(dnorm.le.htol)goto10
        rho=5.D-1*dnorm
        goto10
      elseif(ifail.eq.1)then
        if(iprint.gt.1)
     *    write(nout,*)'phase 1 LCP subproblem is unbounded'
        goto99
      elseif(ifail.gt.5)then
        if(iprint.gt.1)
     *    write(nout,*)'malfunction in phase 1 LCP subproblem'
        goto99
      endif
      hxdJt=0.D0
      hxdJ=0.D0
      do i=1,m
        ci=ws(ncxd+i)
        if(cstype(i).eq.'N')then
          hxdJt=hxdJt+max(0.D0,bl(n+i)-ci,ci-bu(n+i))
        else
          hxdJ=hxdJ+max(0.D0,bl(n+i)-ci,ci-bu(n+i))
        endif
      enddo
      if(iprint.ge.1)write(nout,1001)itn,hxdJt,hxdJ,rgnorm,dnorm,rho
c     print 4,'x+d =',(ws(nxd+i),i=1,n)
c     print 4,'c at x+d =',(ws(ncxd+i),i=1,m)
      if(hxdJt.le.htol.and.dnorm.le.htol)goto40
c     print *,'phase 1 filter entries followed by (hJt,hJ)'
c     do i=nfil1,nfil
c       print 5,filh(i),filf(i)
c     enddo
c     print 5,hJt,hJ
      hxd=hxdJt+hxdJ
      if(hxd.ge.ubd)then
        if(iprint.gt.1)write(nout,*)'upper bound on h exceeded (1)'
        rho=max(1.D-1,5.D-1*h/hxd)*dnorm
        goto10
      endif
      dq=hJ-phi
      df=hJ-hxdJ
c     print 2,'dq,df',dq,df
c  filter test for LCP solution
      call testfil(hxdJt,hxdJ,filh,filf,nfil1,nfil,ifail)
      if(ifail.eq.0)call testfil(hxdJt,hxdJ,hJt,hJ,1,1,ifail)
c     print 6,'hxdJt,hxdJ,ifail',hxdJt,hxdJ,ifail
      if(ifail.eq.1.or.(dq.ge.tol.and.df.lt.sigma*dq))then
        if(hxdJt.eq.0.D0.or.dq.lt.tol)then
          if(iprint.gt.1)write(nout,*)'hxdJt.eq.0.D0.or.dq.lt.tol'
          rho=max(1.D-1,min(5.D-1*h/hxd,5.D-1))*dnorm
          goto10
        endif
c  projection step
        nv=1
        v(1)=1.D0
        iph=-1
        do i=1,n
          ws(naal+i)=ws(nx+i)
        enddo
   30   continue
        hxJt=hxdJt
        hxJ=hxdJ
        do i=1,n
          ws(nx+i)=ws(nxd+i)
          dl(i)=max(-rho,bl(i)-ws(nx+i))
          du(i)=min(rho,bu(i)-ws(nx+i))
c         dl(i)=bl(i)-ws(nx+i)
c         du(i)=bu(i)-ws(nx+i)
          d(i)=0.D0
        enddo
        do i=1,m
          ci=ws(ncxd+i)
          if(cstype(i).eq.'A')then
            dl(n+i)=-ainfty
            du(n+i)=bl(n+i)-ci
          elseif(cstype(i).eq.'Z')then
            dl(n+i)=bu(n+i)-ci
            du(n+i)=ainfty
          else
            dl(n+i)=bl(n+i)-ci
            du(n+i)=bu(n+i)-ci
          endif
        enddo
c       mode=2
        k=0
c  solve projection subproblem
c       print *,'solve phase 1 projection subproblem',itn
        iii=0
c       if(itn.eq.68)iii=1
        call glcpd(n,m,k,kmax,maxg,ws(next1),lws(nla1),d,dl,du,phi,0.D0,
     *    g,r,w,e,ls,alp,lp,mlp,ipeq,ws,lws,cstype,v,nv,rgtol,
     *    mode,ifail,mxgr,iii,0)
c       print 4,'d =',(d(i),i=1,n)
c       print 4,'r =',(r(i),i=1,nm)
c       print 3,'ls =',(ls(i),i=1,nm)
c       print 1,'ipeq,k,ifail =',ipeq,k,ifail
        if(ifail.eq.3)then
          if(iprint.gt.1)
     *      write(nout,*)'phase 1 projection problem is infeasible'
          do i=1,n
            ws(nx+i)=ws(naal+i)
          enddo
          rho=max(1.D-1,min(5.D-1*h/hxd,5.D-1))*dnorm
          goto10
        elseif(ifail.gt.5)then
          if(iprint.gt.1)
     *      write(nout,*)'malfunction in phase 1 projection subproblem'
          goto99
        endif
        do i=1,n
          ws(nxd+i)=ws(nx+i)+d(i)
        enddo
        call functions(n,m,ws(nxd1),fxd,ws(ncxd1),ws,lws)
        call gradients(n,m,ws(nxd1),ws(next1),ws,lws)
        hxdJt=0.D0
        hxdJ=0.D0
        do i=1,m
          if(cstype(i).eq.'N')then
            hxdJt=hxdJt+max(0.D0,bl(n+i)-ws(ncxd+i),ws(ncxd+i)-bu(n+i))
          else
            hxdJ=hxdJ+max(0.D0,bl(n+i)-ws(ncxd+i),ws(ncxd+i)-bu(n+i))
          endif
        enddo
        if(iprint.ge.1)write(nout,1002)itn,hxdJt,hxdJ
c       print 4,'c at x+d =',(ws(ncxd+i),i=1,m)
c  filter test for projection solution
        hxd=hxdJt+hxdJ
        if(hxd.ge.ubd)then
          if(iprint.gt.1)write(nout,*)'upper bound on h exceeded (2)'
          rho=max(1.D-1,5.D-1*h/hxd)*dnorm
          goto10
        endif
        df=hJ-hxdJ
c       print 2,'dq,df',dq,df
        call testfil(hxdJt,hxdJ,filh,filf,nfil1,nfil,ifail)
        if(ifail.eq.0)call testfil(hxdJt,hxdJ,hJt,hJ,1,1,ifail)
c       if(ifail.eq.1)print 2,'hxdJt/hxJt =',hxdJt/hxJt
c       print 6,'project: hxdJt,hxdJ,ifail',hxdJt,hxdJ,ifail
        if(ifail.eq.1.or.df.lt.sigma*dq)then
          if(hxdJt.le.8.D-1*hxJt)then
            df=hJ-(hxJt*hxdJ-hxdJt*hxJ)/(hxJt-hxdJt)
            if(df.ge.sigma*dq)goto30
          endif
          do i=1,n
            ws(nx+i)=ws(naal+i)
          enddo
          rho=max(1.D-1,min(5.D-1*h/hxd,5.D-1))*dnorm
          if(iprint.gt.1)write(nout,*)'phase 1 projection step fails'
          goto10
        endif
c       print *,'accept projection step (1)'
        do i=1,n
          ws(nx+i)=ws(naal+i)
        enddo
      endif
   40 continue
c  accept LCP (iph=1) or projection (iph=-1) solution
      if(dq.lt.tol)then
        call addfil(hJt,hJ,filh,filf,nfil1,nfil,mxf,ifail)
        if(ifail.gt.0)return
      endif
      do i=1,n
        ws(nx+i)=ws(nxd+i)
        ws(naal+i)=0.D0
      enddo
      do i=1,m
        ws(ncx+i)=ws(ncxd+i)
      enddo
      call iexch(last1,next1)
      h=hxdJt+hxdJ
      f=fxd
      if(dnorm.eq.rho)rho=2.D0*rho
      if(h.le.htol)goto10
c  check for situations where the l1 partition needs recalculating ...
c  if there are any active relaxed c/s
      do j=1,n
        i=abs(ls(j))-n
c       if(i.gt.0.and.cstype(i).ne.'N')print 1,'active relaxed c/s',i
        if(i.gt.0.and.cstype(i).ne.'N')goto10
      enddo
c  or any infeasible relaxed c/s
      do i=1,m
c       if((cstype(i).eq.'A'.and.ws(ncx+i).ge.bl(n+i)).or.
c    *    (cstype(i).eq.'Z'.and.ws(ncx+i).le.bu(n+i)))
c    *    print 1,'infeasible relaxed c/s',i
        if((cstype(i).eq.'A'.and.ws(ncx+i).ge.bl(n+i)).or.
     *    (cstype(i).eq.'Z'.and.ws(ncx+i).le.bu(n+i)))goto10
      enddo
      hJt=hxdJt
      hJ=hxdJ
      if(iph.eq.-1)goto10
      if(hxdJt.le.htol.and.dnorm.le.htol)then
        if(iprint.gt.1)write(nout,*)'locally infeasible problem'
c       print 2,'hxdJt.le.htol.and.dnorm.le.htol'
        ifail=3
        return
      endif
c  collect LCP multipliers
      do i=1,m
        ws(nal+i)=0.D0
      enddo
      do j=1,n
        i=abs(ls(j))
        if(i.gt.n)then
          if(ls(j).gt.0)then
            ws(nal+i-n)=r(i)
          else
            ws(nal+i-n)=-r(i)
          endif
        endif
      enddo
c     print 4,'al =',(ws(i),i=nal1,nal+m)
      do i=1,m
        ci=ws(ncx+i)
        if(cstype(i).eq.'A')then
          dl(n+i)=-ainfty
          du(n+i)=bl(n+i)-ci
        elseif(cstype(i).eq.'Z')then
          dl(n+i)=bu(n+i)-ci
          du(n+i)=ainfty
        else
          dl(n+i)=bl(n+i)-ci
          du(n+i)=bu(n+i)-ci
        endif
        call saipy(ws(nal+i),ws(last1),lws(nla1),i,ws(naal1),n)
      enddo
      goto20

   50 continue
c  Phase 2 code
c  collect multipliers from LP subproblem
      do j=1,n
        i=abs(ls(j))
        if(i.gt.n)then
          if(ls(j).gt.0)then
            ws(nal+i-n)=r(i)
          else
            ws(nal+i-n)=-r(i)
          endif
        endif
        ws(naal+j)=0.D0
      enddo
      do i=1,m
        call saipy(ws(nal+i),ws(last1),lws(nla1),i,ws(naal1),n)
      enddo
   60 continue
      if(itn.eq.maxit)then
        if(iprint.gt.1)write(nout,*)'itn.ge.maxit'
        ifail=5
        return
      endif
c     print 2,'h,f =',h,f
      iph=2
      k=0
c     mode=2
c  solve LCP subproblem
c     print 4,'x =',(ws(i),i=nx1,nx+n)
c     print 4,'al =',(ws(i),i=nal1,nal+m)
c     print 4,'c =',(ws(i),i=ncx1,ncx+m)
      do i=1,n
        dl(i)=max(-rho,bl(i)-ws(nx+i))
        du(i)=min(rho,bu(i)-ws(nx+i))
        d(i)=0.D0
      enddo
      alc=scpr(0.D0,ws(nal1),ws(ncx1),m)
c     print *,'solve phase 2 LCP subproblem',itn
      iii=0
c     if(itn.eq.164)iii=1
      call glcpd(n,m,k,kmax,maxg,ws(last1),lws(nla1),d,dl,du,phi,fmin,
     *  g,r,w,e,ls,alp,lp,mlp,ipeq,ws,lws,cstype,v,nv,rgtol,
     *  mode,ifail,mxgr,iii,0)
      nft=nft+nfn
      ngt=ngt+ngr
c     print 1,'nfn,ngr',nfn,ngr
c     print 4,'d =',(d(i),i=1,n)
c     print 4,'r =',(r(i),i=1,nm)
c     print 3,'ls =',(ls(i),i=1,nm)
c     print 1,'ipeq,k,ifail =',ipeq,k,ifail
      itn=itn+1
      dnorm=0.D0
      do i=1,n
        dnorm=max(dnorm,abs(d(i)))
      enddo
c     print 2,'dnorm,rho',dnorm,rho
      if(ifail.eq.3)then
        if(iprint.gt.1)
     *    write(nout,*)'phase 2 LCP problem is infeasible'
c       mode=2
        goto15
c       if(dnorm.le.htol)goto10
c       rho=5.D-1*dnorm
c       goto10
      elseif(ifail.eq.1)then
        if(iprint.gt.1)
     *    write(nout,*)'phase 2 LCP subproblem is unbounded'
        goto99
      elseif(ifail.gt.5)then
        if(iprint.gt.1)
     *    write(nout,*)'malfunction in phase 2 LCP subproblem'
        goto99
      endif
      hxd=0.D0
      do i=1,m
        ci=ws(ncxd+i)
        hxd=hxd+max(0.D0,bl(n+i)-ci,ci-bu(n+i))
      enddo
      if(iprint.ge.1)write(nout,2001)itn,hxd,fxd,rgnorm,dnorm,rho
c     print 4,'c at x+d =',(ws(ncxd+i),i=1,m)
      if(hxd.le.htol.and.(fxd.le.fmin.or.dnorm.le.htol))goto80
c     print *,'phase 2 filter entries followed by (h,f)'
c     do i=1,nfil
c       print 5,filh(i),filf(i)
c     enddo
c     print 5,h,f
c  filter test for LCP solution
      if(hxd.ge.ubd)then
        if(iprint.gt.1)write(nout,*)'upper bound on h exceeded (3)'
        rho=max(1.D-1,5.D-1*h/hxd)*dnorm
        goto10
      endif
      dq=f-phi
      df=f-fxd
c     print 2,'dq,df',dq,df
      call testfil(hxd,fxd,filh,filf,1,nfil,ifail)
      if(ifail.eq.0)call testfil(hxd,fxd,h,f,1,1,ifail)
c     print 6,'hxd,fxd,ifail',hxd,fxd,ifail
      if(ifail.eq.1.or.(dq.ge.tol.and.df.lt.sigma*dq))then
        if(hxd.eq.0.D0.or.dq.lt.tol)then
          rho=5.D-1*dnorm
          if(iprint.gt.1)write(nout,*)'hxd.eq.0.D0.or.dq.lt.tol'
          rho=max(1.D-1,min(5.D-1*h/hxd,5.D-1))*dnorm
          goto10
        endif
c  projection step
        nv=1
        v(1)=1.D0
        iph=-2
        do i=1,n
          ws(naal+i)=ws(nx+i)
        enddo
   70   continue
        hx=hxd
        fx=fxd
        do i=1,n
          ws(nx+i)=ws(nxd+i)
          dl(i)=max(-rho,bl(i)-ws(nx+i))
          du(i)=min(rho,bu(i)-ws(nx+i))
c         dl(i)=bl(i)-ws(nx+i)
c         du(i)=bu(i)-ws(nx+i)
          d(i)=0.D0
        enddo
        do i=1,m
          ci=ws(ncxd+i)
          dl(n+i)=bl(n+i)-ci
          du(n+i)=bu(n+i)-ci
        enddo
c       mode=2
        k=0
c       print 4,'x =',(ws(nx+i),i=1,n)
c       print 4,'v =',(v(i),i=1,nv)
c  solve projection subproblem
c       print *,'solve phase 2 projection subproblem'
        iii=0
c       if(itn.eq.9)iii=3
        call glcpd(n,m,k,kmax,maxg,ws(next1),lws(nla1),d,dl,du,phi,0.D0,
     *    g,r,w,e,ls,alp,lp,mlp,ipeq,ws,lws,cstype,v,nv,rgtol,
     *    mode,ifail,mxgr,iii,0)
c       print 4,'d =',(d(i),i=1,n)
c       print 4,'r =',(r(i),i=1,nm)
c       print 3,'ls =',(ls(i),i=1,nm)
c       print 1,'ipeq,k,ifail =',ipeq,k,ifail
        if(ifail.eq.3)then
          if(iprint.gt.1)
     *      write(nout,*)'phase 2 projection problem is infeasible'
          do i=1,n
            ws(nx+i)=ws(naal+i)
          enddo
          rho=5.D-1*dnorm
          goto10
        elseif(ifail.gt.5)then
          if(iprint.gt.1)
     *     write(nout,*),'malfunction in phase 2 projection subproblem'
          goto99
        endif
        do i=1,n
          ws(nxd+i)=ws(nx+i)+d(i)
        enddo
c       print 4,'xd =',(ws(nxd+i),i=1,n)
        call functions(n,m,ws(nxd1),fxd,ws(ncxd1),ws,lws)
        call gradients(n,m,ws(nxd1),ws(next1),ws,lws)
        hxd=0.D0
        do i=1,m
          ci=ws(ncxd+i)
          hxd=hxd+max(0.D0,bl(n+i)-ci,ci-bu(n+i))
        enddo
        if(iprint.ge.1)write(nout,2002)itn,hxd,fxd
c       print 4,'x+d =',(ws(nxd+i),i=1,n)
c       print 4,'c at x+d =',(ws(ncxd+i),i=1,m)
c  filter test for projection solution
        if(hxd.ge.ubd)then
          if(iprint.gt.1)write(nout,*)'upper bound on h exceeded (4)'
          rho=max(1.D-1,5.D-1*h/hxd)*dnorm
          goto10
        endif
        df=f-fxd
        call testfil(hxd,fxd,filh,filf,1,nfil,ifail)
        if(ifail.eq.0)call testfil(hxd,fxd,h,f,1,1,ifail)
c       if(ifail.eq.1)print 2,'hxd/hx =',hxd/hx
c       print 6,'hxd,fxd,ifail',hxd,fxd,ifail
        if(ifail.eq.1.or.df.lt.sigma*dq)then
          if(hxd.le.8.D-1*hx)then
            df=f-(hx*fxd-hxd*fx)/(hx-hxd)
            if(df.ge.sigma*dq)goto70
          endif
          do i=1,n
            ws(nx+i)=ws(naal+i)
          enddo
          rho=5.D-1*dnorm
          if(iprint.gt.1)write(nout,*)'phase 2 projection step fails'
          goto10
        endif
c       print *,'accept phase 2 projection step'
        do i=1,n
          ws(nx+i)=ws(naal+i)
        enddo
      endif
   80 continue
c  accept LCP (iph=2) or projection (iph=-2) solution
      if(dq.lt.tol)then
        call addfil(h,f,filh,filf,1,nfil,mxf,ifail)
        if(ifail.gt.0)return
      endif
      h=hxd
      f=fxd
      if(dnorm.eq.rho)rho=2.D0*rho
      do i=1,n
        ws(nx+i)=ws(nxd+i)
        ws(naal+i)=0.D0
      enddo
      do i=1,m
        ws(ncx+i)=ws(ncxd+i)
      enddo
      call iexch(last1,next1)
      if(iph.eq.-2)goto10
      if(h.le.htol)then
        if(f.le.fmin)then
          if(iprint.gt.1)
     *      write(nout,*)'phase 2 LCP unbounded'
          ifail=1
          return
        elseif(dnorm.le.htol)then
          if(iprint.gt.1)write(nout,*)'local NLP solution found'
c         print 2,'h.le.htol.and.dnorm.le.htol'
          ifail=0
          return
        endif
      endif
c  collect LCP multipliers
      do i=1,m
        ws(nal+i)=0.D0
      enddo
      do j=1,n
        i=abs(ls(j))
        if(i.gt.n)then
          if(ls(j).gt.0)then
            ws(nal+i-n)=r(i)
          else
            ws(nal+i-n)=-r(i)
          endif
        endif
      enddo
c     print 4,'al =',(ws(i),i=nal1,nal+m)
      do i=1,m
        ci=ws(ncx+i)
        dl(n+i)=bl(n+i)-ci
        du(n+i)=bu(n+i)-ci
        call saipy(ws(nal+i),ws(last1),lws(nla1),i,ws(naal1),n)
      enddo
c     iph=2
      goto60
   99 continue
      if(ifail.eq.7)return
      ifail=ifail+10
      return
      end

      block data nlp_defaults
      implicit double precision (a-h,o-z)
      common/defaultc/ainfty,ubd,mlp,mxf
      common/ngrc/mxgr
      data  ainfty, ubd, mlp, mxf,  mxgr
     *    / 1.D20,  1.D4, 50,  50, 1000000/
      end

      subroutine testfil(h,f,filh,filf,nfil1,nfil,ifail)
      implicit double precision (a-h,o-z)
      dimension filh(*),filf(*)
      parameter (beta=99999.D-5,gamma=1.D-5)
      ifail=0
c     if(h.eq.0.D0)return
      hd=h/beta
      fp=f+gamma*h
      do i=nfil1,nfil
        if(hd.ge.filh(i).and.fp.gt.filf(i))then
c       if(hd.gt.filh(i).and.fp.gt.filf(i))then
          ifail=1
          return
        endif
      enddo
      return
      end

      subroutine addfil(h,f,filh,filf,nfil1,nfil,mxf,ifail)
      implicit double precision (a-h,o-z)
      dimension filh(*),filf(*)
      do i=nfil,nfil1,-1
        if(h.le.filh(i).and.f.le.filf(i))then
          filh(i)=filh(nfil)
          filf(i)=filf(nfil)
          nfil=nfil-1
        endif
      enddo
      if(nfil.ge.mxf)then
        ifail=8
        return
      endif
      nfil=nfil+1
      filh(nfil)=h
      filf(nfil)=f
      ifail=0
      return
      end

      subroutine funct(n,d,phi,ws,lws,cstype)
      implicit double precision (a-h,o-z)
      dimension d(*),ws(*),lws(*)
      character cstype(*)
      common/defaultc/ainfty,ubd,mlp,mxf
      common/functc/fxd,alc,m,iph,last1,next1,nx,nx1,
     *  nal,nal1,naal,naal1,nxd,nxd1,ncx,ncx1,ncxd,ncxd1,nla1
    2 format(A,6E15.7)
    4 format(A/(5E15.7))
      if(iph.lt.0)then
c  projection subproblem
        phi=5.D-1*scpr(0.D0,d,d,n)
        return
      elseif(iph.eq.0)then
c  LP subproblem
        phi=aiscpr(n,ws(last1),lws(nla1),0,d,0.D0)
        return
      endif
      do i=1,n
        ws(nxd+i)=ws(nx+i)+d(i)
      enddo
c     print 4,'$x =',(ws(nx+i),i=1,n)
c     print 4,'$d =',(d(i),i=1,n)
c     print 4,'$x+d =',(ws(nxd+i),i=1,n)
      call functions(n,m,ws(nxd1),fxd,ws(ncxd1),ws,lws)
c     print 2,'$fxd =',fxd
c     print 4,'$al =',(ws(nal+i),i=1,m)
c     print 4,'$A.al =',(ws(naal+i),i=1,n)
c     print 4,'$cxd =',(ws(ncxd+i),i=1,m)
c     print *,'$cstype =',(cstype(i),i=1,m)
      if(iph.eq.2)then
        phi=scpr(-scpr(-fxd-alc,ws(nal1),ws(ncxd1),m),ws(naal1),d,n)
      else
        phi=scpr(-scpr(-alc,ws(nal1),ws(ncxd1),m),ws(naal1),d,n)
        do i=1,m
          if(cstype(i).eq.'A')then
            phi=phi-ws(ncxd+i)
          elseif(cstype(i).eq.'Z')then
            phi=phi+ws(ncxd+i)
          endif
        enddo
      endif
c     print 2,'$fxd,cxd,phi =',fxd,ws(ncxd1),phi
c     print 2,'$phi =',phi
      return
      end

      subroutine grad(n,d,g,ws,lws,cstype)
      implicit double precision (a-h,o-z)
      dimension d(*),g(*),ws(*),lws(*)
      character cstype(*)
      common/defaultc/ainfty,ubd,mlp,mxf
      common/functc/fxd,alc,m,iph,last1,next1,nx,nx1,
     *  nal,nal1,naal,naal1,nxd,nxd1,ncx,ncx1,ncxd,ncxd1,nla1
      common/maxac/maxa
      if(iph.lt.0)then
        do i=1,n
          g(i)=d(i)
        enddo
        return
      elseif(iph.eq.0)then
        do i=1,n
          g(i)=0.D0
        enddo
        call saipy(1.D0,ws(last1),lws(nla1),0,g,n)
        return
      endif
c     print 4,'£x+d =',(ws(nxd+i),i=1,n)
      call gradients(n,m,ws(nxd1),ws(next1),ws,lws)
c     print 4,'£a =',(ws(next1+i),i=0,maxa-1)
c     print 4,'£al =',(ws(nal+i),i=1,m)
      do i=1,n
        g(i)=ws(naal+i)
      enddo
c     print 4,'£Ak.al =',(g(i),i=1,n)
      do i=1,m
        call saipy(-ws(nal+i),ws(next1),lws(nla1),i,g,n)
      enddo
      if(iph.eq.2)then
        call saipy(1.D0,ws(next1),lws(nla1),0,g,n)
      else
        do i=1,m
          if(cstype(i).eq.'A')then
            call saipy(-1.D0,ws(next1),lws(nla1),i,g,n)
          elseif(cstype(i).eq.'Z')then
            call saipy(1.D0,ws(next1),lws(nla1),i,g,n)
          endif
        enddo
      endif
c     print 4,'£g =',(g(i),i=1,n)
    1 format(A,15I5)
    3 format(A/(15I5))
    4 format(A/(5E15.7))
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ############################################################
c     ##                                                        ##
c     ##  subroutine extra3  --  user defined extra potentials  ##
c     ##                                                        ##
c     ############################################################
c
c
c     "extra3" calculates any additional user defined potential
c     contribution and also partitions the energy among the atoms
c
c
      subroutine extra3
      implicit none
      include 'sizes.i'
      include 'action.i'
      include 'analyz.i'
      include 'atoms.i'
      include 'energi.i'
      integer i
c
c
c     zero out the energy due to extra potential terms
c
      nex = 0
      ex = 0.0d0
      do i = 1, n
         aex(i) = 0.0d0
      end do
c
c     add any user-defined extra potentials and partitioning
c
c     e = ......
c     nex = nex + 1
c     ex = ex + e
c     aex(i) = ......
c
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
!!# FUNCTION MODULE: <OldFile>
MODULE FUN_OldFile

!!## PURPOSE
!! Return a unit corresponding to a file.

!!## USAGE
!
!  Unit = OldFile( File [, IOSTAT] )
!
!! where <File> is a filename and <IOSTAT> is an optional
!! status variable.

!!## EXTERNAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_S !!((01-A-KND_IntrinsicTypes.f90))
USE FUN_NewFile                     !!((05-B-FUN_NewFile.f90))

!!## HISTORY
! 1.[waw] Author   = William A. Wieselquist
!         Modified = 99.2007
!         Contact  = william.wieselquist AT gmail.com

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## DEFAULT ACCESS
PRIVATE

!!## PUBLIC ACCESS LIST
PUBLIC :: OldFile

!!## MODULE PROCEDURES
CONTAINS

!!### FUNCTION: <OldFile>
FUNCTION OldFile(File,IOSTAT) RESULT(Unit)

!!#### REQUIRED INPUT
CHARACTER(LEN=*,KIND=KIND_S),INTENT(IN) :: File

!!#### REQUIRED OUTPUT
INTEGER :: Unit

!!#### OPTIONAL OUTPUT
INTEGER,INTENT(OUT),OPTIONAL :: IOSTAT

!!#### LOCAL VARIABLES
LOGICAL :: Opened
INTEGER :: IOSTAT_

!!--begin--
INQUIRE( File=File , Opened=Opened , IOSTAT=IOSTAT_ )

IF( .NOT.Opened )THEN
 Unit = NewFile(File)
ELSE
 INQUIRE( File=File , Number=Unit , IOSTAT=IOSTAT_ )
END IF

IF( PRESENT(IOSTAT) )THEN
 IOSTAT = IOSTAT_
END IF

!!--end--
END FUNCTION

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
module feval_module
  use iso_c_binding
  use pf_mod_dtype
  use encap_module
  use pf_mod_utils
  use pf_mod_misdcQ
  use pf_mod_rkstepper
  implicit none

  ! Define the derived sweeper type
  type, extends(pf_misdcQ_t) :: sweet_sweeper_t
     type(c_ptr)    :: ctx = c_null_ptr ! c pointer to PlaneDataCtx/SphereDataCtx
     integer        :: nnodes           ! number of nodes
     integer        :: sweep_niter      ! number of the current sweep
     integer        :: sweep_niter_max  ! max number of sweeps
     integer        :: level            ! level on which the sweeper is acting
     real(c_double) :: dt               ! full timestep size
   contains
     procedure :: f_eval                => sweet_f_eval_sweeper
     procedure :: f_comp                => sweet_f_comp_sweeper
     procedure :: destroy               => sweet_sweeper_destroy
     procedure :: compute_dt            => sweet_sweeper_compute_dt
  end type sweet_sweeper_t

  ! Define the derived stepper type
  type, extends(pf_ark_stepper_t) :: sweet_stepper_t
     type(c_ptr)    :: ctx = c_null_ptr ! c pointer to PlaneDataCtx/SphereDataCtx
     integer        :: nnodes           ! number of nodes
     integer        :: sweep_niter      ! number of the current sweep
     integer        :: sweep_niter_max  ! max number of sweeps
     integer        :: level            ! level on which the sweeper is acting
     real(c_double) :: dt               ! full timestep size
   contains
     procedure :: f_eval                => sweet_f_eval_stepper
     procedure :: f_comp                => sweet_f_comp_stepper
     procedure :: f_finalize            => sweet_f_finalize_stepper
     procedure :: destroy               => sweet_stepper_destroy
  end type sweet_stepper_t

  ! prototypes of the C functions

  interface
     subroutine cinitial(i_ctx, i_t, i_dt, o_Y) bind(c, name="cinitial")
       use iso_c_binding
       type(c_ptr),    value :: i_ctx, o_Y
       real(c_double), value :: i_t, i_dt
     end subroutine cinitial

     subroutine cfinal(i_ctx, i_Y, i_nnodes, i_niter) bind(c, name="cfinal")
       use iso_c_binding
       type(c_ptr), value :: i_ctx, i_Y
       integer,     value :: i_nnodes
       integer,     value :: i_niter
     end subroutine cfinal

     subroutine ceval_f1(i_Y, i_t, i_ctx, o_F1) bind(c, name="ceval_f1")
       use iso_c_binding
       type(c_ptr),    value :: i_Y, i_ctx, o_F1
       real(c_double), value :: i_t
     end subroutine ceval_f1

     subroutine ceval_f2(i_Y, i_t, i_ctx, o_F2) bind(c, name="ceval_f2")
       use iso_c_binding
       type(c_ptr),    value :: i_Y, i_ctx, o_F2
       real(c_double), value :: i_t
     end subroutine ceval_f2

     subroutine ccomp_f2(io_Y, i_t, i_dtq, i_Rhs, i_ctx, o_F2) bind(c, name="ccomp_f2")
       use iso_c_binding
       type(c_ptr),    value :: io_Y, i_Rhs, i_ctx, o_F2
       real(c_double), value :: i_t, i_dtq
     end subroutine ccomp_f2

     subroutine ceval_f3(i_Y, i_t, i_level, i_ctx, o_F3) bind(c, name="ceval_f3")
       use iso_c_binding
       type(c_ptr),    value :: i_Y, i_ctx, o_F3
       integer,        value :: i_level
       real(c_double), value :: i_t
     end subroutine ceval_f3

     subroutine ccomp_f3(i_Y, i_t, i_dtq, i_level, i_Rhs, i_ctx, o_F3) bind(c, name="ccomp_f3")
       use iso_c_binding
       type(c_ptr),    value :: i_Y, i_Rhs, i_ctx, o_F3
       real(c_double), value :: i_t, i_dtq
       integer,        value :: i_level
     end subroutine ccomp_f3

     subroutine cfinalize(i_Y, i_t, i_dt, i_ctx) bind(c, name="cfinalize")
       use iso_c_binding
       type(c_ptr),    value :: i_Y, i_ctx
       real(c_double), value :: i_t, i_dt
     end subroutine cfinalize

  end interface

contains

  ! function to cast the pf_sweeper_t object into a sweet_sweeper_t object
  function as_sweet_sweeper(sweeper) result(r)
    class(pf_sweeper_t),    intent(in), target  :: sweeper
    class(sweet_sweeper_t),             pointer :: r

    select type(sweeper)
    type is (sweet_sweeper_t)
       r => sweeper
    class default
       stop "TYPE ERROR"
    end select

  end function as_sweet_sweeper

  ! function to cast the pf_rkstepper_t object into a sweet_rkstepper_t object
  function as_sweet_stepper(stepper) result(r)
    class(pf_stepper_t),    intent(in), target  :: stepper
    class(sweet_stepper_t),             pointer :: r

    select type(stepper)
    type is (sweet_stepper_t)
       r => stepper
    class default
       stop "TYPE ERROR"
    end select

  end function as_sweet_stepper

  ! initialization of the sweet data vector using z

  subroutine finitial(sweeper, sd, z, t, dt)
    class(pf_sweeper_t),       intent(inout) :: sweeper
    class(pf_encap_t),         intent(inout) :: sd
    class(pf_encap_t),         intent(inout) :: z
    real(c_double),            intent(in)    :: t, dt

    class(sweet_sweeper_t),    pointer       :: sweet_sweeper_ptr
    class(sweet_data_encap_t), pointer       :: sd_ptr

    sweet_sweeper_ptr => as_sweet_sweeper(sweeper)
    sd_ptr            => as_sweet_data_encap(sd)

    call cinitial(sweet_sweeper_ptr%ctx, &
                  t,                     &
                  dt,                    &
                  sd_ptr%c_sweet_data_ptr)

    call z%copy(sd)

  end subroutine finitial

  subroutine ffinal(sweeper, sd, nnodes, niter)
    use mpi

    class(pf_sweeper_t),       intent(in)    :: sweeper
    class(pf_encap_t),         intent(inout) :: sd
    integer,                   intent(in)    :: nnodes, niter

    class(sweet_sweeper_t),    pointer       :: sweet_sweeper_ptr
    class(sweet_data_encap_t), pointer       :: sd_ptr

    sweet_sweeper_ptr => as_sweet_sweeper(sweeper)
    sd_ptr            => as_sweet_data_encap(sd)

    call cfinal(sweet_sweeper_ptr%ctx,   &
                sd_ptr%c_sweet_data_ptr, &
                nnodes,                  &
                niter)

  end subroutine ffinal

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!             SWEEPER             !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! evaluate the right-hand side

  subroutine sweet_f_eval_sweeper(this, y, t, level_index, f, piece)
    class(sweet_sweeper_t),    intent(inout) :: this
    class(pf_encap_t),         intent(in)    :: y
    real(pfdp),                intent(in)    :: t
    integer,                   intent(in)    :: level_index
    class(pf_encap_t),         intent(inout) :: f
    integer,                   intent(in)    :: piece

    class(sweet_data_encap_t), pointer       :: y_sd_ptr
    class(sweet_data_encap_t), pointer       :: f_sd_ptr

    y_sd_ptr  => as_sweet_data_encap(y)
    f_sd_ptr  => as_sweet_data_encap(f)

    select case (piece)

       case (1) ! misdc rhs
          call ceval_f1(y_sd_ptr%c_sweet_data_ptr, &
                        t,                         &
                        this%ctx,                  &
                        f_sd_ptr%c_sweet_data_ptr)

       case (2) ! first implicit rhs
          call ceval_f2(y_sd_ptr%c_sweet_data_ptr, &
                        t,                         &
                        this%ctx,                  &
                        f_sd_ptr%c_sweet_data_ptr)

       case (3) ! second implicit rhs
          call ceval_f3(y_sd_ptr%c_sweet_data_ptr, &
                        t,                         &
                        this%level-1,              &
                        this%ctx,                  &
                        f_sd_ptr%c_sweet_data_ptr)

       case DEFAULT
          print *, 'Piece argument in f_eval can only be 1, 2'!, or 3'
          call exit(0)
    end select

  end subroutine sweet_f_eval_sweeper

  ! solve for y in the implicit system + update f

  subroutine sweet_f_comp_sweeper(this, y, t, dtq, rhs, level_index, f, piece)
    class(sweet_sweeper_t),   intent(inout) :: this
    class(pf_encap_t),        intent(inout) :: y, f
    real(pfdp),               intent(in)    :: t, dtq
    class(pf_encap_t),        intent(in)    :: rhs
    integer,                  intent(in)    :: level_index
    integer,                  intent(in)    :: piece

    class(sweet_data_encap_t), pointer      :: y_sd_ptr
    class(sweet_data_encap_t), pointer      :: f_sd_ptr
    class(sweet_data_encap_t), pointer      :: rhs_sd_ptr

    y_sd_ptr   => as_sweet_data_encap(y)
    f_sd_ptr   => as_sweet_data_encap(f)
    rhs_sd_ptr => as_sweet_data_encap(rhs)

    select case (piece)

         case (2) ! first implicit solve
          call ccomp_f2(y_sd_ptr%c_sweet_data_ptr,   &
                        t,                           &
                        dtq,                         &
                        rhs_sd_ptr%c_sweet_data_ptr, &
                        this%ctx,                    &
                        f_sd_ptr%c_sweet_data_ptr)

        case (3) ! second implicit solve
           call ccomp_f3(y_sd_ptr%c_sweet_data_ptr,   &
                         t,                           &
                         dtq,                         &
                         this%level-1,                &
                         rhs_sd_ptr%c_sweet_data_ptr, &
                         this%ctx,                    &
                         f_sd_ptr%c_sweet_data_ptr)

       case DEFAULT
          print *, 'Piece argument in f_comp can only be 2'! or 3'
          call exit(0)
    end select

  end subroutine sweet_f_comp_sweeper

  ! destructor

  subroutine sweet_sweeper_destroy(this, pf, level_index)
    class(sweet_sweeper_t), intent(inout) :: this
    type(pf_pfasst_t),   intent(inout),target :: pf
    integer,             intent(in)    :: level_index

    ! this is copy-pasted from LibPFASST, unsure about this
    type(pf_level_t), pointer  :: lev       !  Current level
    lev => pf%levels(level_index)           !  Assign level pointer

    ! need the following line since the "final" keyword is not supported by some (older) compilers
    ! it forces Fortran to destroy the parent class data structures
    call this%misdcQ_destroy(pf, level_index)

  end subroutine sweet_sweeper_destroy

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!             STEPPER             !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! evaluate the right-hand side

  subroutine sweet_f_eval_stepper(this, y, t, level_index, f, piece)
    class(sweet_stepper_t),    intent(inout) :: this
    class(pf_encap_t),         intent(in)    :: y
    real(pfdp),                intent(in)    :: t
    integer,                   intent(in)    :: level_index
    class(pf_encap_t),         intent(inout) :: f
    integer,                   intent(in)    :: piece

    class(sweet_data_encap_t), pointer       :: y_sd_ptr
    class(sweet_data_encap_t), pointer       :: f_sd_ptr

    y_sd_ptr  => as_sweet_data_encap(y)
    f_sd_ptr  => as_sweet_data_encap(f)

    select case (piece)

       case (1) ! misdc rhs
          call ceval_f1(y_sd_ptr%c_sweet_data_ptr, &
                        t,                         &
                        this%ctx,                  &
                        f_sd_ptr%c_sweet_data_ptr)

       case (2) ! first implicit rhs
          call ceval_f2(y_sd_ptr%c_sweet_data_ptr, &
                        t,                         &
                        this%ctx,                  &
                        f_sd_ptr%c_sweet_data_ptr)

       case (3) ! second implicit rhs
          call ceval_f3(y_sd_ptr%c_sweet_data_ptr, &
                        t,                         &
                        level_index-1,             &
                        this%ctx,                  &
                        f_sd_ptr%c_sweet_data_ptr)

       case DEFAULT
          print *, 'Piece argument in f_eval can only be 1, 2'!, or 3'
          call exit(0)
    end select

  end subroutine sweet_f_eval_stepper

  ! solve for y in the implicit system + update f

  subroutine sweet_f_comp_stepper(this, y, t, dtq, rhs, level_index, f, piece)
    class(sweet_stepper_t),   intent(inout) :: this
    class(pf_encap_t),        intent(inout) :: y, f
    real(pfdp),               intent(in)    :: t, dtq
    class(pf_encap_t),        intent(in)    :: rhs
    integer,                  intent(in)    :: level_index
    integer,                  intent(in)    :: piece

    class(sweet_data_encap_t), pointer      :: y_sd_ptr
    class(sweet_data_encap_t), pointer      :: f_sd_ptr
    class(sweet_data_encap_t), pointer      :: rhs_sd_ptr

    y_sd_ptr   => as_sweet_data_encap(y)
    f_sd_ptr   => as_sweet_data_encap(f)
    rhs_sd_ptr => as_sweet_data_encap(rhs)

    select case (piece)

         case (2) ! first implicit solve
          call ccomp_f2(y_sd_ptr%c_sweet_data_ptr,   &
                        t,                           &
                        dtq,                          &
                        rhs_sd_ptr%c_sweet_data_ptr, &
                        this%ctx,                    &
                        f_sd_ptr%c_sweet_data_ptr)

        case (3) ! second implicit solve
           call ccomp_f3(y_sd_ptr%c_sweet_data_ptr,   &
                         t,                           &
                         dtq,                          &
                         level_index-1,               &
                         rhs_sd_ptr%c_sweet_data_ptr, &
                         this%ctx,                    &
                         f_sd_ptr%c_sweet_data_ptr)

       case DEFAULT
          print *, 'Piece argument in f_comp can only be 2'! or 3'
          call exit(0)
    end select

  end subroutine sweet_f_comp_stepper

  ! apply artificial viscosity

    subroutine sweet_f_finalize_stepper(this, y, t, dtq, level_index)
    class(sweet_stepper_t), intent(inout) :: this
    class(pf_encap_t),      intent(inout) :: y
    real(pfdp),             intent(in   ) :: t
    real(pfdp),             intent(in   ) :: dtq
    integer,                intent(in   ) :: level_index

    class(sweet_data_encap_t),    pointer :: y_sd_ptr

    y_sd_ptr   => as_sweet_data_encap(y)

    call cfinalize(y_sd_ptr%c_sweet_data_ptr, &
                   t,                         &
                   dtq,                       &
                   this%ctx)

  end subroutine sweet_f_finalize_stepper

  ! destructor

  subroutine sweet_stepper_destroy(this, pf, level_index)
    class(sweet_stepper_t),         intent(inout) :: this
    type(pf_pfasst_t),      target, intent(inout) :: pf
    integer,                        intent(in   ) :: level_index

    ! this is copy-pasted from LibPFASST, unsure about this
    type(pf_level_t), pointer  :: lev       !  Current level
    lev => pf%levels(level_index)           !  Assign level pointer

    ! need the following line since the "final" keyword is not supported by some (older) compilers
    ! it forces Fortran to destroy the parent class data structures
    call this%ark_destroy(pf, level_index)

  end subroutine sweet_stepper_destroy

  subroutine sweet_sweeper_compute_dt(this, pf, level_index, t0, dt, flags)
    class(sweet_sweeper_t),         intent(inout) :: this
    type(pf_pfasst_t), target,      intent(inout) :: pf
    integer,                        intent(in   ) :: level_index
    real(pfdp),                     intent(in   ) :: t0
    real(pfdp),                     intent(inout) :: dt
    integer, optional,              intent(in   ) :: flags

    type(pf_level_t),    pointer :: lev
    lev => pf%levels(level_index)   !!  Assign level pointer
    !  Do nothing now (copy-pasted from pf_imex_sweeper)
    return
  end subroutine sweet_sweeper_compute_dt

end module feval_module

max_stars_repo_name: schreibm/sweet

max_issues_repo_name: schreibm/sweet

max_forks_repo_name: schreibm/sweet


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:23 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE READ_VALUE_INTEGER_NEXT_P__genmod
          INTERFACE
            SUBROUTINE READ_VALUE_INTEGER_NEXT_P(VALUE_TO_READ,NAME)
              INTEGER(KIND=4) :: VALUE_TO_READ
              CHARACTER(*) :: NAME
            END SUBROUTINE READ_VALUE_INTEGER_NEXT_P
          END INTERFACE
        END MODULE READ_VALUE_INTEGER_NEXT_P__genmod


text:
      module rrsw_tbl

      use shr_kind_mod, only: r8 => shr_kind_r8

!      use parkind, only : jpim, jprb

      implicit none
      save

!------------------------------------------------------------------
! rrtmg_sw lookup table arrays

! Initial version: MJIacono, AER, may2007
! Revised: MJIacono, AER, aug2007
!------------------------------------------------------------------

!  name     type     purpose
! -----  :  ----   : ----------------------------------------------
! ntbl   :  integer: Lookup table dimension
! tblint :  real   : Lookup table conversion factor
! tau_tbl:  real   : Clear-sky optical depth
! exp_tbl:  real   : Exponential lookup table for transmittance
! od_lo  :  real   : Value of tau below which expansion is used
!                  : in place of lookup table
! pade   :  real   : Pade approximation constant
! bpade  :  real   : Inverse of Pade constant
!------------------------------------------------------------------

      integer, parameter :: ntbl = 10000

      real(kind=r8), parameter :: tblint = 10000.0

      real(kind=r8), parameter :: od_lo = 0.06

      real(kind=r8) :: tau_tbl
      real(kind=r8) , dimension(0:ntbl) :: exp_tbl

      real(kind=r8), parameter :: pade = 0.278_r8
      real(kind=r8) :: bpade

      end module rrsw_tbl

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
Looking for a place to shoot some hoops?

Regular pickup games:
  Sundays @ 8:00 a.m.

Public basketball courts:
  Central Park
   1 half court  isnt it lit?
  Dairy Basketball Courts near the south end of the Dairy IM Field on UC Davis campus.
  La Playa Park
  Manor Park or Slide Hill Park  pretty close to regulation height
  Oak Grove Park
  Pioneer Park
  West Manor Park
  Willowcreek Park  half court
  http://www.dcn.davis.ca.us/go/stoneg8/ Stonegate Country Club  pseudopublic: must be a CC member (or rent from someone that is, or know someone that is and doesnt use their membership).  Hoops are generally way too low, which leads to people dunking, which exacerbates the situation.
  Mace Ranch Park

Lighted outdoor courts:
  Chestnut Park
  Hacienda Park
  Playfields Park  I think they stay lit until 12.  I may be wrong.
   At the western end of Research Park Dr.
   2 half courts.  in good shape.

Indoor courts:
  The ARC
      There are pickup games in the mornings, weekends and evenings. Many good players, not just undergrads, guys in their 40s and 50s as well. Fun to watch.
  Hickey Gym
   Neither of these are open to the general public.  In fact, evening and weekend rec. ball isnt scheduled at Hickey at all, as of April 19th.

Apartments with courts
  Allegre Apartments
   1 full court
  Arlington Farms Apartments
  The Lexington Apartments, on the southeast edge
  Silverstone Apartments
  Tanglewood Apartments
  Temescal Apartments
   1 full court.  crappy.  lots of rocks and dirt on the court.

Schools with outdoor courts
  Any Davis elementary school
   Their hoops are usually misaligned, but are really forgiving.
  Any Davis Jr. High.
   Emersons hoops are all 46 too low.  They resurfaced their courts a few years back and never raised the hoops to match.
   The new Jr. High has nets! Havent measured the hoop height yet.

20080713 13:12:26 nbsp are there any leagues or clubs? Users/atwong


text:
      subroutine cross(x,y,j)
C             ha ha eh e  po"pamm :o pe e  t   koop  hat   to k
C              epece eh    byx  "eometp  eck x   emehtob,  h opma    o
C             kotop x  a  caha b   oke
C              opma o pa eh  : call cross(x,y,j)
C             o  cah e  apametpob:
C             x,y - koop  hat  to k   epece eh
C             j - xapaktep ct ka to k   epece eh
C             j=1 -  epeceka8tc  cam    emeht
C             j=2 -  epb     emeht  epecekaet  po o  eh e btopo"o
C             j=3 - btopo    emeht  epecekaet  po o  eh e  epbo"o
C             j=4 -  epeceka8tc   po o  eh     emehtob
      common/gfgel/gf(28)
      dimension n(2)
      n(1)=0
      n(2)=0
      i=1
      if(gf(1)) 28,1,9
  1   if(gf(11)) 28,2,11
C              epece eh e  p m x
  2   d=gf(5)*gf(14)-gf(15)*gf(4)
      if(abs(d).lt.1e-8) go to 23
      gf(19)=(gf(6)*gf(15)-gf(16)*gf(5))/d
      gf(20)=(gf(16)*gf(4)-gf(6)*gf(14))/d
      gf(10)=1.
  3   continue
C              pobepka  p ha  e hoct  to k   epece eh    y"e okpy hoct
C             otpe ky  p mo
      do 7 k=1,2
      if(gf(i).ne.0.) go to 4
      if((gf(i+1)-gf(19))*(gf(19)-gf(i+6)).lt.0.) go to 6
       if((gf(i+2)-gf(20))*(gf(20)-gf(i+7))) 6,5,5
  4   if(gf(i).eq.2.) go to 5
      th=angler(gf(19)-gf(i+1),gf(20)-gf(i+2))
      if((gf(i+6)-th)*(th-gf(i+7)).lt.0.) go to 6
  5   n(k)=1
  6    if(i.gt.10) go to 24
      i=11
  7   continue
  8   return
  9   if(gf(11)) 28,10,18
  10  k=11
      l=1
      go to 12
  11  k=1
      l=11
  12  if(gf(10)) 8,30,13
  13  if(gf(10)-1.) 8,27,14
C              epece eh e okpy hoct  c  p mo
  14  s=sqrt(gf(k+3)**2+gf(k+4)**2)
      if(s.eq.0.) go to 23
      p=(gf(k+3)*gf(l+1)+gf(k+4)*gf(l+2)+gf(k+5))/s
      st=gf(k+3)/s
      ct=-gf(k+4)/s
      r=gf(l+3)
      if(abs(p)-abs(r)) 16,15,23
  15  gf(10)=1.
      gf(19)=-p*st+gf(l+1)
      gf(20)=p*ct+gf(l+2)
       go to 3
  16  gf(10)=2.
      gf(27)=21
      x=sqrt(gf(l+3)**2-p**2)
      gf(19)=x*ct-p*st+gf(l+1)
      gf(20)=x*st+p*ct+gf(l+2)
      gf(21)=-x*ct-p*st+gf(l+1)
      gf(22)=-x*st+p*ct+gf(l+2)
       go to 3
  18  if(gf(10)) 8,30,19
  19  if(gf(10)-1.) 8,27,20
C              epece eh e  byx okpy hocte       x  y"
  20  r1=abs(gf(4))
      r2=abs(gf(14))
      d1=sqrt((gf(12)-gf(2))**2+(gf(13)-gf(3))**2)
      if(d1.eq.0.) go to 23
      if(d1.gt.r1+r2) go to 23
      x1=(r1**2-r2**2+d1**2)/(2*d1)
      st=(gf(13)-gf(3))/d1
      ct=(gf(12)-gf(2))/d1
      h=sqrt(r1**2-x1**2)
      if(h) 8,21,22
  21  gf(10)=1.
      gf(19)=r1*ct+gf(2)
      gf(20)=r1*st+gf(3)
       go to 3
  22  gf(10)=2.
      gf(27)=21
      gf(19)=x1*ct-h*st+gf(2)
      gf(20)=x1*st+h*ct+gf(3)
      gf(21)=x1*ct+h*st+gf(2)
      gf(22)=x1*st-h*ct+gf(3)
       go to 3
  23  j=0
      go to 29
C              opm pobah e j - xapaktep ct k  to k   epece eh
  24  j=4-(2*n(1)+n(2))
      gf(10)=gf(10)-1
      if(gf(10)) 8,26,25
  25  j=-j
  26  x=gf(19)
      y=gf(20)
      return
  27  l=gf(27)
      gf(19)=gf(l)
      gf(20)=gf(l+1)
       go to 3
  28  call gffals(6)
      return
  29  call gffals(5)
      return
  30  call gffals(7)
      return
       end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief prepare (log) apriori and bayes values
!>    @param[in] ismpl local sample index
      SUBROUTINE prepare_jacobian(ismpl)
        use arrays
#ifdef AD
        use g_arrays
#endif
#ifdef AD_RM
        use arrays_ad
#endif
        use mod_genrl
        use mod_inverse
        use mod_linfos
        IMPLICIT NONE
        integer :: ismpl
        integer :: k
        INTEGER pe
        DOUBLE PRECISION get_optia
        INTEGER param_enabled
        EXTERNAL get_optia, param_enabled

        IF (linfos(2)>=1) WRITE(*,*) ' ... calculate (lin/log) values'
!       opti_*    : 0,1,(2) = off, lin, log
        DO k = 1, mpara
          pe = param_enabled(k,ismpl)
          IF (pe==2) THEN
            apri_input(k) = log(get_optia(k,ismpl))
            linlog_input(k) = log(main_input(k,ismpl))
          ELSE IF (pe==1) THEN
            apri_input(k) = get_optia(k,ismpl)
            linlog_input(k) = main_input(k,ismpl)
          END IF
        END DO
        RETURN
      END