text:
! This file created from test/mpi/f77/io/i_setviewcurf.f with f77tof90
! -*- Mode: Fortran; -*-
!
!  (C) 2003 by Argonne National Laboratory.
!      See COPYRIGHT in top-level directory.
!
      program main
      use mpi
      integer (kind=MPI_OFFSET_KIND) offset

      integer errs, ierr, size, rank
      integer fh, comm, status(MPI_STATUS_SIZE)
      integer buf(1024)
      integer request

      errs = 0
      call MTest_Init( ierr )

!    This test reads a header then sets the view to every "size" int,
!    using set view and current displacement.  The file is first written
!    using a combination of collective and ordered writes

      comm = MPI_COMM_WORLD
      call MPI_File_open( comm, "test.ord", MPI_MODE_WRONLY + &
      &     MPI_MODE_CREATE, MPI_INFO_NULL, fh, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Open(1)", ierr )
      endif
      call MPI_Comm_size( comm, size, ierr )
      call MPI_Comm_rank( comm, rank, ierr )
      if (size .gt. 1024) then
         if (rank .eq. 0) then
            print *, &
      &"This program must be run with no more than 1024 processes"
            call MPI_Abort( MPI_COMM_WORLD, 1, ierr )
         endif
      endif

      buf(1) = size
      call MPI_File_iwrite_all( fh, buf, 1, MPI_INTEGER, request, ierr)
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "iwrite_all", ierr )
      endif

      call MPI_Wait(request, status, ierr)

      call MPI_File_get_position( fh, offset, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Get_position", ierr )
      endif
      call MPI_File_seek_shared( fh, offset, MPI_SEEK_SET, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Seek_shared", ierr )
      endif
      buf(1) = rank
      call MPI_File_write_ordered( fh, buf, 1, MPI_INTEGER, status,ierr)
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Write_ordered", ierr )
      endif
      call MPI_File_close( fh, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Close(1)", ierr )
      endif

! Reopen the file as sequential
      call MPI_File_open( comm, "test.ord", MPI_MODE_RDONLY + &
      &     MPI_MODE_SEQUENTIAL + MPI_MODE_DELETE_ON_CLOSE, &
      &     MPI_INFO_NULL, fh, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Open(Read)", ierr )
      endif

      if (rank .eq. 0) then
         call MPI_File_read_shared( fh, buf, 1, MPI_INTEGER, status, &
      &        ierr )
         if (ierr .ne. MPI_SUCCESS) then
            errs = errs + 1
            call MTestPrintErrorMsg( "Read_all", ierr )
         endif
         if (buf(1) .ne. size) then
           errs = errs + 1
           print *, "Unexpected value for the header = ", buf(1), &
      &          ", should be ", size
        endif
      endif
      call MPI_Barrier( comm, ierr  )
!   All processes must provide the same file view for MODE_SEQUENTIAL
      call MPI_File_set_view( fh, MPI_DISPLACEMENT_CURRENT, MPI_INTEGER &
      &     ,MPI_INTEGER, "native", MPI_INFO_NULL, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Set_view", ierr )
      endif
      buf(1) = -1
      call MPI_File_read_ordered( fh, buf, 1, MPI_INTEGER, status, ierr &
      &     )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Read_all", ierr )
      endif
      if (buf(1) .ne. rank) then
         errs = errs + 1
         print *, rank, ": buf(1) = ", buf(1)
      endif

      call MPI_File_close( fh, ierr )
      if (ierr .ne. MPI_SUCCESS) then
         errs = errs + 1
         call MTestPrintErrorMsg( "Close(2)", ierr )
      endif

      call MTest_Finalize( errs )

      call MPI_Finalize( ierr )
      end

max_stars_repo_name: dominikwerder/mpich-cache

max_issues_repo_name: dominikwerder/mpich-cache

max_forks_repo_name: dominikwerder/mpich-cache


text:
PROGRAM fptest
  !--------- -------- --------- --------- --------- --------- --------- --------- -----
  !
  ! Example program 1 for using the function parser module
  !
  !--------- -------- --------- --------- --------- --------- --------- --------- -----
  USE parameters, ONLY: rn
  USE fparser,    ONLY: initf, parsef, evalf, EvalErrType, EvalErrMsg
  IMPLICIT NONE
  INTEGER,                             PARAMETER :: nfunc = 4
  CHARACTER (LEN=*), DIMENSION(nfunc), PARAMETER :: func = (/ '-1.0*x        ',  &
                                                              '-x            ',  &
                                                              'A*COS(B*x)+5  ',  &
                                                              'A*COS(B*x)+5.0' /)
  INTEGER,                             PARAMETER :: nvar = 3
  CHARACTER (LEN=*), DIMENSION(nvar),  PARAMETER :: var  = (/ 'x', &
                                                              'A', &
                                                              'B'  /)
  REAL(rn),          DIMENSION(nvar),  PARAMETER :: val  = (/  2., 3., 4. /)
  REAL(rn)                                       :: res
  INTEGER                                        :: i
  REAL(rn)                                       :: x,A,B
  !--------- -------- --------- --------- --------- --------- --------- --------- -----
  !
  CALL initf (nfunc)                      ! Initialize function parser for nfunc functions
  DO i=1,nfunc
     CALL parsef (i, func(i), var)        ! Parse and bytecompile ith function string
  END DO
  WRITE(*,*)'-----------------------------------------------------------------------------'
  WRITE(*,*)'==> Bytecode evaluation:'
  DO i=1,nfunc
     res = evalf (i, val)                 ! Interprete bytecode representation of ith function
     IF (EvalErrType > 0) WRITE(*,*)'*** Error: ',EvalErrMsg ()
     WRITE(*,*) func(i),'=',res
  END DO
  WRITE(*,*)'==> Direct evaluation:'
  x  = val(1)
  A  = val(2)
  B  = val(3)
  WRITE(*,*)'-1.0*x        =',-1.0*x
  WRITE(*,*)'-x            =',-x
  WRITE(*,*)'A*COS(B*x)+5  =',A*COS(B*x)+5
  WRITE(*,*)'A*COS(B*X)+5.0=',A*COS(B*X)+5.0
  !
END PROGRAM fptest

max_stars_repo_name: alexispaz/fparser

max_issues_repo_name: alexispaz/fparser

max_forks_repo_name: alexispaz/fparser


! Take 2 plotfiles as input and compare them zone by zone for differences.
!
! For the comparison to take place, the grids must be identical.
!

program fcompare

  use f2kcli
  use bl_space
  use bl_error_module
  use bl_constants_module
  use bl_IO_module
  use plotfile_module
  use multifab_module

  implicit none

  type(plotfile) :: pf_a, pf_b
  character (len=256) :: plotfile_a, plotfile_b
  character (len=256) :: diffvar
  integer :: unit_a, unit_b

  real(kind=dp_t), pointer :: p_a(:,:,:,:), p_b(:,:,:,:)
  real(kind=dp_t), pointer :: mp(:,:,:,:)

  integer :: lo_a(MAX_SPACEDIM), hi_a(MAX_SPACEDIM)

  integer :: flo_b(MAX_SPACEDIM), fhi_b(MAX_SPACEDIM)

  real(kind=dp_t) :: dx_a(MAX_SPACEDIM), dx_b(MAX_SPACEDIM)

  integer :: nboxes_a, nboxes_b

  integer :: n_a, n_b
  integer, allocatable :: ivar_b(:)
  integer :: save_var_a

  real(kind=dp_t), allocatable :: aerror(:), rerror(:), rerror_denom(:)

  logical, allocatable :: has_nan_a(:), has_nan_b(:)
  logical :: any_nans

  integer :: norm

  integer :: narg, farg
  character (len=256) :: fname

  integer :: i, j
  integer :: ii, jj, kk

  integer ir_a, ir_b, ng

  integer :: itest

  real(kind=dp_t) :: global_error
  real(kind=dp_t) :: pa, pb, pd, aerr, rerr

  integer :: dm
  type(box) :: bx_a, bx_b

  type(multifab), allocatable :: mf_array(:)
  type(layout) :: la
  type(boxarray) :: ba
  type(list_box) :: bl
  integer, allocatable :: ref_ratio(:)

  character(len=20), allocatable :: plot_names(:)

  logical :: do_ghost, gc_warn

  !---------------------------------------------------------------------------
  ! process the command line arguments

  narg = command_argument_count()

  ! defaults
  norm = 0
  plotfile_a = ""
  plotfile_b = ""

  diffvar = ""
  do_ghost = .false.

  allocate(plot_names(1))

  farg = 1
  do while (farg <= narg)
     call get_command_argument(farg, value = fname)

     select case (fname)

     case ('--infile1')
        farg = farg + 1
        call get_command_argument(farg, value = plotfile_a)

     case ('--infile2')
        farg = farg + 1
        call get_command_argument(farg, value = plotfile_b)

     case ('-n','--norm')
        farg = farg + 1
        call get_command_argument(farg, value = fname)
        read(fname, *) norm

     case ('-g','--ghost')
        farg = farg + 1
        do_ghost = .true.

     case ('--diffvar')
        farg = farg + 1
        call get_command_argument(farg, value = diffvar)
        plot_names(1) = trim(diffvar)

     case default
        exit

     end select
     farg = farg + 1
  enddo

  if (len_trim(plotfile_a) == 0) then
     call get_command_argument(farg, value = plotfile_a)
     farg = farg + 1
  endif

  if (len_trim(plotfile_b) == 0) then
     call get_command_argument(farg, value = plotfile_b)
     farg = farg + 1
  endif

  if (len_trim(plotfile_a) == 0 .OR. len_trim(plotfile_b) == 0) then
     print *, " "
     print *, "Compare two plotfiles, zone by zone, to machine precision"
     print *, "and report the maximum absolute and relative errors for each"
     print *, "variable."
     print *, " "
     print *, "usage:"
     print *, "   fcompare [-g|--ghost] [-n|--norm num] [--diffvar var] file1 file2"
     print *, " "
     print *, "optional arguments:"
     print *, "   -g|--ghost    : compare the ghost cells too (if stored)"
     print *, "   -n|--norm num : what norm to use (default is 0 for inf norm)"
     print *, "   --diffvar var : output a plotfile showing the differences for "
     print *, "                   variable var"
     print *, " "
     stop
  endif

  !---------------------------------------------------------------------------
  ! build the plotfiles and do initial comparisons
  !---------------------------------------------------------------------------

  unit_a = unit_new()
  call build(pf_a, plotfile_a, unit_a)

  unit_b = unit_new()
  call build(pf_b, plotfile_b, unit_b)

  dm = pf_a%dim

  ! check if they are the same dimensionality
  if (pf_a%dim /= pf_b%dim) then
     call bl_error("ERROR: plotfiles have different numbers of spatial dimensions")
  endif

  ! check if they have the same number of levels
  if (pf_a%flevel /= pf_b%flevel) then
     call bl_error("ERROR: number of levels do not match")
  endif

  ! check if the finest domains are the same size
  bx_a = plotfile_get_pd_box(pf_a, pf_a%flevel)
  bx_b = plotfile_get_pd_box(pf_b, pf_b%flevel)

  if (.not. box_equal(bx_a, bx_b)) then
     call bl_error("ERROR: grids do not match")
  endif

  ! check if they have the same number of variables
  if (pf_a%nvars /= pf_b%nvars) then
     print *, " "
     print *, "WARNING: number of variables do not match"
  endif

  allocate(aerror(pf_a%nvars))
  allocate(rerror(pf_a%nvars))
  allocate(rerror_denom(pf_a%nvars))

  allocate(has_nan_a(pf_a%nvars))
  allocate(has_nan_b(pf_a%nvars))

  any_nans = .false.

  save_var_a = -1

  ! in case the variables are not in the same order, figure out the
  ! mapping between pf_a and pf_b variables
  allocate(ivar_b(pf_a%nvars))

  do n_a = 1, pf_a%nvars

     ivar_b(n_a) = -1
     do n_b = 1, pf_b%nvars

        if (pf_a%names(n_a) == pf_b%names(n_b)) then
           ivar_b(n_a) = n_b
           exit
        endif

     enddo

     if (ivar_b(n_a) == -1) then
        print *, "WARNING: variable ", trim(pf_a%names(n_a)), &
                 " not found in plotfile 2"
     endif

     if (.not. diffvar == "") then
        if (pf_a%names(n_a) == trim(diffvar)) then
           save_var_a = n_a
        endif
     endif

  enddo

  ! also print out, as a diagnostic, those variables in plotfile 1 that
  ! are not in plotfile 2
  do n_b = 1, pf_b%nvars
     itest = -1
     do n_a = 1, pf_a%nvars

        if (pf_a%names(n_a) == pf_b%names(n_b)) then
           itest = n_a
           exit
        endif

     enddo

     if (itest == -1) then
        print *, "WARNING: variable ", trim(pf_b%names(n_b)), &
                 " not found in plotfile 1"
     endif

  enddo

  !---------------------------------------------------------------------------
  ! create a multifab to store the difference for output, if desired
  !---------------------------------------------------------------------------

  if (save_var_a > 0) then

     allocate(mf_array(plotfile_nlevels(pf_a)))
     allocate(ref_ratio(plotfile_nlevels(pf_a)-1))

     ! define ref_ratio
     do i = 1, pf_a%flevel-1
        ref_ratio(i) = pf_a%refrat(i,1)
     enddo

     ! loop over levels and grids and define the boxes needed to build the
     ! multifab
     do i = 1, pf_a%flevel

        do j = 1, nboxes(pf_a, i)
           call push_back(bl, get_box(pf_a,i,j))
        enddo

        call build(ba,bl)
        call layout_build_ba(la,ba,plotfile_get_pd_box(pf_a,1))

        ! destroy the list and boxarray so we start over next level
        call destroy(bl)
        call destroy(ba)

        ! create a new multifab with 0 ghost cells and 1 component
        call multifab_build(mf_array(i),la,1,0)

     enddo

  endif

  !---------------------------------------------------------------------------
  ! go level-by-level and patch-by-patch and compare the data
  !---------------------------------------------------------------------------

998 format(1x,a24,2x,a24,   2x,a24)
999 format(1x,70("-"))

  write (*,*) " "
  write (*,998) "variable name", "absolute error",  "relative error"
  write (*,998) "",              "(||A - B||)",     "(||A - B||/||A||)"
  write (*,999)

  gc_warn = .false.

  do i = 1, pf_a%flevel

     aerror(:) = ZERO
     rerror(:) = ZERO
     rerror_denom(:) = ZERO

     has_nan_a(:) = .false.
     has_nan_b(:) = .false.

     ! make sure the dx's agree
     dx_a = 0.0_dp_t
     dx_b = 0.0_dp_t
     dx_a(1:dm) = plotfile_get_dx(pf_a, i)
     dx_b(1:dm) = plotfile_get_dx(pf_b, i)

     if ((dx_a(1) /= dx_b(1)) .OR. &
         (pf_a%dim >= 2 .AND. dx_a(2) /= dx_b(2)) .OR. &
         (pf_a%dim == 3 .AND. dx_a(3) /= dx_b(3))) then
        call bl_error("ERROR: grid dx does not match")
     endif

     ! make sure the number of boxes agree
     nboxes_a = nboxes(pf_a, i)
     nboxes_b = nboxes(pf_b, i)

     if (nboxes_a /= nboxes_b) then
        call bl_error("ERROR: number of boxes do not match")
     endif

     do j = 1, nboxes_a

        ! make sure that the grids match
        bx_a = get_box(pf_a, i, j)

        lo_a = 1
        hi_a = 1
        lo_a(1:dm) = lwb(bx_a)
        hi_a(1:dm) = upb(bx_a)

        bx_b = get_box(pf_b, i, j)

        if (.not. box_equal(bx_a, bx_b)) then
           call bl_error("ERROR: grids do not match")
        endif

        if (.not. nghost(pf_a, i, j) == nghost(pf_b, i, j)) then
           if (.not. gc_warn) then
              call bl_warn("WARNING: grids have different numbers of ghost cells")
              gc_warn = .true.
           endif
           ng = 0
        else
           if (do_ghost) then
              ng = nghost(pf_a, i, j)
           else
              ng = 0
           endif
        endif

        ! loop over the variables.  Take plotfile_a to be the one defining
        ! the list of variables, and bind them one-by-one.  Don't assume that
        ! the variables are in the same order in plotfile_b.

        do n_a = 1, pf_a%nvars

           n_b = ivar_b(n_a)
           if (n_b == -1) cycle

           call fab_bind_comp_vec(pf_a, i, j, (/ n_a /) )
           call fab_bind_comp_vec(pf_b, i, j, (/ n_b /) )

           p_a => dataptr(pf_a, i, j)
           p_b => dataptr(pf_b, i, j)

           ! are we storing the diff?
           if (n_a == save_var_a) then
              mp => dataptr(mf_array(i), j, get_box(mf_array(i), j))

              mp(:,:,:,1) = abs(p_a(:,:,:,1) - p_b(:,:,:,1))
           endif

           ! check for NaNs -- comparisons don't work when they are present
           ! note: regardless of do_ghost, this will check the ghostcells
           ! too if they are present.
           call fab_contains_nan(p_a, volume(get_pbox(pf_a, i, j)), ir_a)
           if (ir_a == 1) has_nan_a(n_a) = .true.

           call fab_contains_nan(p_b, volume(get_pbox(pf_b, i, j)), ir_b)
           if (ir_b == 1) has_nan_b(n_a) = .true.

           if (has_nan_a(n_a) .or. has_nan_b(n_a)) cycle

           do kk = lo_a(3)-ng, hi_a(3)+ng
              do jj = lo_a(2)-ng, hi_a(2)+ng
                 do ii = lo_a(1)-ng, hi_a(1)+ng

                    pa = abs(p_a(ii,jj,kk,1))
                    pb = abs(p_b(ii,jj,kk,1))
                    pd = abs(p_a(ii,jj,kk,1) - p_b(ii,jj,kk,1))

                    if (norm == 0) then
                       aerror(n_a) = max(aerror(n_a), pd)

                       rerror(n_a) = max(rerror(n_a), pd)
                       rerror_denom(n_a) = max(rerror_denom(n_a), pa)
                    else
                       aerror(n_a) = aerror(n_a) + pd**norm

                       rerror(n_a) = rerror(n_a) + pd**norm
                       rerror_denom(n_a) = rerror_denom(n_a) + pa**norm
                    endif

                 enddo
              enddo
           enddo

           call fab_unbind(pf_a, i, j)
           call fab_unbind(pf_b, i, j)

        enddo  ! variable loop

     enddo  ! boxes loop

     ! Normalize
     if (norm > 0) then

        do n_a = 1, pf_a%nvars
           aerror(n_a) = aerror(n_a)*product(dx_a(1:pf_a%dim))
           aerror(n_a) = aerror(n_a)**(ONE/real(norm,dp_t))

           ! since we are taking the ratio of two norms, no grid normalization
           ! is needed
           rerror(n_a) = (rerror(n_a)/rerror_denom(n_a))**(ONE/real(norm,dp_t))
        enddo

     else

        do n_a = 1, pf_a%nvars
           rerror(n_a) = rerror(n_a)/rerror_denom(n_a)
        enddo

     endif

     !------------------------------------------------------------------------
     ! print out the comparison report for this level
     !------------------------------------------------------------------------

1000 format(1x,"level = ", i2)
1001 format(1x,a24,2x,g24.10,2x,g24.10)
1002 format(1x,a24,2x,a50)

     write (*,1000) i

     do n_a = 1, pf_a%nvars
        if (ivar_b(n_a) == -1) then
           write (*,1002) pf_a%names(n_a), "< variable not present in both files > "
        else if (has_nan_a(n_a) .or. has_nan_b(n_a)) then
           write (*,1002) pf_a%names(n_a), "< NaN present > "
        else
           if (aerror(n_a) > 0.0d0) then
              aerr = min(max(aerror(n_a), 1.d-99), 1.d98)
           else
              aerr = 0.0d0
           endif

           if (rerror(n_a) > 0.0d0) then
              rerr = min(max(rerror(n_a), 1.d-99), 1.d98)
           else
              rerr = 0.0d0
           endif

           write (*,1001) pf_a%names(n_a), aerr, rerr
        endif
     enddo

     if (i == 1) then
        global_error = maxval(aerror(:))
     else
        global_error = max(global_error, maxval(aerror(:)))
     endif

     any_nans = any(has_nan_a .or. has_nan_b) .or. any_nans

  enddo  ! level loop

  if (global_error == ZERO .and. .not. any_nans) then
     print *, "PLOTFILES AGREE"
  endif

  if (save_var_a > 0) then
     call fabio_ml_multifab_write_d(mf_array, ref_ratio, &
                                    "diffs", plot_names, &
                                    plotfile_get_pd_box(pf_a,1), &
                                    pf_a%plo, pf_a%phi, &
                                    pf_a%tm, &
                                    plotfile_get_dx(pf_a,1))
  endif

end program fcompare


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE BUSH ( INT_ELEM_ID, OPT, WRITE_WARN )

! Calculates, for BUSH element

!  1) KE        = element stiffness matrix if OPT(4) = 'Y'
!  2) SE1, STE1 = element stress data recovery matrices if OPT(3) = 'Y'

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  F04, F06, WRT_LOG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE CONSTANTS_1, ONLY           :  ZERO
      USE MODEL_STUF, ONLY            :  BE1, BE2, BUSH_DXA, BUSH_DXB, BUSH_DY, BUSH_DZ, EPROP, KE, OFFDIS_GA_GB, SE1, SE2
      USE SUBR_BEGEND_LEVELS, ONLY    :  BUSH_BEGEND

      USE BUSH_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BUSH'
      CHARACTER(LEN=*) , INTENT(IN)   :: WRITE_WARN        ! If 'Y" write warning messages, otherwise do not
      CHARACTER(1*BYTE), INTENT(IN)   :: OPT(6)            ! 'Y'/'N' flags for whether to calc certain elem matrices

      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BUSH_BEGEND
      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID       ! Internal element ID
      INTEGER(LONG)                   :: I,J               ! DO loop indices

      REAL(DOUBLE)                    :: DUM0(6,12)        ! Intermediate matrix in a calc
      REAL(DOUBLE)                    :: DUM1(6,12)        ! Intermediate matrix in a calc
      REAL(DOUBLE)                    :: STIFF(12,12)      ! 12x12 partition from KE (complete stiff matrix for this elem)
      REAL(DOUBLE)                    :: STRE_RCV(2)       ! Two stress recovery values
      REAL(DOUBLE)                    :: STRN_RCV(2)       ! Two strain recovery values

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Initialize:

      BUSH_DXA = OFFDIS_GA_GB(1,1)
      BUSH_DY  = OFFDIS_GA_GB(1,2)
      BUSH_DZ  = OFFDIS_GA_GB(1,3)
      BUSH_DXB = OFFDIS_GA_GB(2,1)

! Calc stiffness matrix

      IF ((OPT(3) == 'Y') .OR. (OPT(4) == 'Y')) THEN

         DO I=1,12
            DO J=1,12
               KE(I,J) = ZERO
            ENDDO
         ENDDO

         KE( 1, 1) =  EPROP(1)                             ! NB *** new 09/10/21. Big change. KE (before offsets): just EPROP vals
         KE( 2, 2) =  EPROP(2)
         KE( 3, 3) =  EPROP(3)
         KE( 4, 4) =  EPROP(4)
         KE( 5, 5) =  EPROP(5)
         KE( 6, 6) =  EPROP(6)

         KE( 7, 7) =  KE( 1, 1)
         KE( 8, 8) =  KE( 2, 2)
         KE( 9, 9) =  KE( 3, 3)
         KE(10,10) =  KE( 4, 4)
         KE(11,11) =  KE( 5, 5)
         KE(12,12) =  KE( 6, 6)

         KE( 1, 7) = -KE( 1, 1)
         KE( 2, 8) = -KE( 2, 2)
         KE( 3, 9) = -KE( 3, 3)
         KE( 4,10) = -KE( 4, 4)
         KE( 5,11) = -KE( 5, 5)
         KE( 6,12) = -KE( 6, 6)

         DO I=1,12
            DO J=I,12
               KE(J,I) = KE(I,J)
            ENDDO
         ENDDO

      ENDIF

! Calc stress and strain recovery matrices. The stress recovery coefficients on the Bulk Data entry are intended to be
! multiplied times the element ENGINEERING (NOT NODE) forces to get stresses. In order to calc stresses in the MYSTRAN fasion
! of multiplying SEi times displacements, the SEi matrices are generated with the aid of the matrix that converts element
! nodal forces to element engineering forces. The DUM0 matrix (below), if all of the STRE_RCV numbers were 1.0 is that matrix
! (see conversion code in subr OFP3_ELFE_1D. With the STRE_RCV values in DUM0, then DUM0 times KE would be the SE matrices
! that convert displacements to stresses

      IF (OPT(3) == 'Y') THEN
                                                           ! Calc SEi stress recovery matrices
         STRE_RCV(1) = EPROP(19)                           ! --- Stress/strain recovery coefficients from PBUSH entry
         STRE_RCV(2) = EPROP(20)
         STRN_RCV(1) = EPROP(21)
         STRN_RCV(2) = EPROP(22)

         DO I=1,6
            DO J=1,12
               DUM0(I,J) = ZERO
               DUM1(I,J) = ZERO
            ENDDO
         ENDDO

         DUM0(1, 7) =  STRE_RCV(1)
         DUM0(2, 8) =  STRE_RCV(1)
         DUM0(3, 9) =  STRE_RCV(1)

         DUM0(4,10) =  STRE_RCV(2)
         DUM0(5,11) =  STRE_RCV(2)
         DUM0(6,12) =  STRE_RCV(2)

         DUM0(4, 8) =  STRE_RCV(2)*BUSH_DZ
         DUM0(4, 9) = -STRE_RCV(2)*BUSH_DY
         DUM0(5, 7) = -STRE_RCV(2)*BUSH_DZ
         DUM0(5, 9) = -STRE_RCV(2)*BUSH_DXB
         DUM0(6, 7) =  STRE_RCV(2)*BUSH_DY
         DUM0(6, 8) =  STRE_RCV(2)*BUSH_DXB

         DO I=1,12
            DO J=1,12
               STIFF(I,J) = KE(I,J)
            ENDDO
         ENDDO

         CALL MATMULT_FFF ( DUM0, STIFF, 6, 12, 12, DUM1 ) ! --- Dummy 6x12 matrix that has the SEi matrix data

         DO    I=1,3                                                                                !    ---    SEi
            DO    J=1,12
               SE1(I,J,1) = DUM1(I  ,J)
               SE2(I,J,1) = DUM1(I+3,J)
            ENDDO
         ENDDO

         BE1(1, 1,1) = -STRN_RCV(1)                        ! Calc BEi strain recovery matrices
         BE1(1, 7,1) =  STRN_RCV(1)

         BE1(2, 2,1) = -STRN_RCV(1)
         BE1(2, 8,1) =  STRN_RCV(1)

         BE1(3, 3,1) = -STRN_RCV(1)
         BE1(3, 9,1) =  STRN_RCV(1)

         BE2(1, 4,1) = -STRN_RCV(2)
         BE2(1,10,1) =  STRN_RCV(2)

         BE2(2, 5,1) = -STRN_RCV(2)
         BE2(2,11,1) =  STRN_RCV(2)

         BE2(3, 6,1) = -STRN_RCV(2)
         BE2(3,12,1) =  STRN_RCV(2)

      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

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

      END SUBROUTINE BUSH


text:
MODULE mo_column_extra
 USE mo_column , ONLY: ty_column , compute_column

CONTAINS
 SUBROUTINE compute_one ( nz , q , t , nproma )
  INTEGER , INTENT(IN) :: nproma

  INTEGER , INTENT(IN) :: nz
  REAL , INTENT(INOUT) :: t ( : , : )
  REAL , INTENT(INOUT) :: q ( : , : )
  TYPE ( ty_column ) :: column

  CALL column % compute_column ( nz , q , t , nproma = nproma )
 END SUBROUTINE compute_one

END MODULE mo_column_extra

max_stars_repo_name: grnydawn/Frelpt

max_issues_repo_name: grnydawn/Frelpt

max_forks_repo_name: grnydawn/Frelpt


text:
C     ----------------------------------------------------------------------
C     This file contains the LBFGS algorithm and supporting routines
C
C     ****************
C     LBFGS SUBROUTINE
C     ****************
C
      SUBROUTINE LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
C
      INTEGER N,M,IPRINT(2),IFLAG
      DOUBLE PRECISION X(N),G(N),DIAG(N),W(N*(2*M+1)+2*M)
      DOUBLE PRECISION F,EPS,XTOL
      LOGICAL DIAGCO
C
C        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
C                          JORGE NOCEDAL
C                        *** July 1990 ***
C
C
C     This subroutine solves the unconstrained minimization problem
C
C                      min F(x),    x= (x1,x2,...,xN),
C
C      using the limited memory BFGS method. The routine is especially
C      effective on problems involving a large number of variables. In
C      a typical iteration of this method an approximation Hk to the
C      inverse of the Hessian is obtained by applying M BFGS updates to
C      a diagonal matrix Hk0, using information from the previous M steps.
C      The user specifies the number M, which determines the amount of
C      storage required by the routine. The user may also provide the
C      diagonal matrices Hk0 if not satisfied with the default choice.
C      The algorithm is described in "On the limited memory BFGS method
C      for large scale optimization", by D. Liu and J. Nocedal,
C      Mathematical Programming B 45 (1989) 503-528.
C
C      The user is required to calculate the function value F and its
C      gradient G. In order to allow the user complete control over
C      these computations, reverse  communication is used. The routine
C      must be called repeatedly under the control of the parameter
C      IFLAG.
C
C      The steplength is determined at each iteration by means of the
C      line search routine MCVSRCH, which is a slight modification of
C      the routine CSRCH written by More' and Thuente.
C
C      The calling statement is
C
C          CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
C
C      where
C
C     N       is an INTEGER variable that must be set by the user to the
C             number of variables. It is not altered by the routine.
C             Restriction: N>0.
C
C     M       is an INTEGER variable that must be set by the user to
C             the number of corrections used in the BFGS update. It
C             is not altered by the routine. Values of M less than 3 are
C             not recommended; large values of M will result in excessive
C             computing time. 3<= M <=7 is recommended. Restriction: M>0.
C
C     X       is a DOUBLE PRECISION array of length N. On initial entry
C             it must be set by the user to the values of the initial
C             estimate of the solution vector. On exit with IFLAG=0, it
C             contains the values of the variables at the best point
C             found (usually a solution).
C
C     F       is a DOUBLE PRECISION variable. Before initial entry and on
C             a re-entry with IFLAG=1, it must be set by the user to
C             contain the value of the function F at the point X.
C
C     G       is a DOUBLE PRECISION array of length N. Before initial
C             entry and on a re-entry with IFLAG=1, it must be set by
C             the user to contain the components of the gradient G at
C             the point X.
C
C     DIAGCO  is a LOGICAL variable that must be set to .TRUE. if the
C             user  wishes to provide the diagonal matrix Hk0 at each
C             iteration. Otherwise it should be set to .FALSE., in which
C             case  LBFGS will use a default value described below. If
C             DIAGCO is set to .TRUE. the routine will return at each
C             iteration of the algorithm with IFLAG=2, and the diagonal
C              matrix Hk0  must be provided in the array DIAG.
C
C
C     DIAG    is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE.,
C             then on initial entry or on re-entry with IFLAG=2, DIAG
C             it must be set by the user to contain the values of the
C             diagonal matrix Hk0.  Restriction: all elements of DIAG
C             must be positive.
C
C     IPRINT  is an INTEGER array of length two which must be set by the
C             user.
C
C             IPRINT(1) specifies the frequency of the output:
C                IPRINT(1) < 0 : no output is generated,
C                IPRINT(1) = 0 : output only at first and last iteration,
C                IPRINT(1) > 0 : output every IPRINT(1) iterations.
C
C             IPRINT(2) specifies the type of output generated:
C                IPRINT(2) = 0 : iteration count, number of function
C                                evaluations, function value, norm of the
C                                gradient, and steplength,
C                IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of
C                                variables and  gradient vector at the
C                                initial point,
C                IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of
C                                variables,
C                IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.
C
C
C     EPS     is a positive DOUBLE PRECISION variable that must be set by
C             the user, and determines the accuracy with which the solution
C             is to be found. The subroutine terminates when
C
C                         ||G|| < EPS max(1,||X||),
C
C             where ||.|| denotes the Euclidean norm.
C
C     XTOL    is a  positive DOUBLE PRECISION variable that must be set by
C             the user to an estimate of the machine precision (e.g.
C             10**(-16) on a SUN station 3/60). The line search routine will
C             terminate if the relative width of the interval of uncertainty
C             is less than XTOL.
C
C     W       is a DOUBLE PRECISION array of length N(2M+1)+2M used as
C             workspace for LBFGS. This array must not be altered by the
C             user.
C
C     IFLAG   is an INTEGER variable that must be set to 0 on initial entry
C             to the subroutine. A return with IFLAG<0 indicates an error,
C             and IFLAG=0 indicates that the routine has terminated without
C             detecting errors. On a return with IFLAG=1, the user must
C             evaluate the function F and gradient G. On a return with
C             IFLAG=2, the user must provide the diagonal matrix Hk0.
C
C             The following negative values of IFLAG, detecting an error,
C             are possible:
C
C              IFLAG=-1  The line search routine MCSRCH failed. The
C                        parameter INFO provides more detailed information
C                        (see also the documentation of MCSRCH):
C
C                       INFO = 0  IMPROPER INPUT PARAMETERS.
C
C                       INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF
C                                 UNCERTAINTY IS AT MOST XTOL.
C
C                       INFO = 3  MORE THAN 20 FUNCTION EVALUATIONS WERE
C                                 REQUIRED AT THE PRESENT ITERATION.
C
C                       INFO = 4  THE STEP IS TOO SMALL.
C
C                       INFO = 5  THE STEP IS TOO LARGE.
C
C                       INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C                                 THERE MAY NOT BE A STEP WHICH SATISFIES
C                                 THE SUFFICIENT DECREASE AND CURVATURE
C                                 CONDITIONS. TOLERANCES MAY BE TOO SMALL.
C
C
C              IFLAG=-2  The i-th diagonal element of the diagonal inverse
C                        Hessian approximation, given in DIAG, is not
C                        positive.
C
C              IFLAG=-3  Improper input parameters for LBFGS (N or M are
C                        not positive).
C
C
C
C    ON THE DRIVER:
C
C    The program that calls LBFGS must contain the declaration:
C
C                       EXTERNAL LB2
C
C    LB2 is a BLOCK DATA that defines the default values of several
C    parameters described in the COMMON section.
C
C
C
C    COMMON:
C
C     The subroutine contains one common area, which the user may wish to
C    reference:
C
         COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
C
C    MP  is an INTEGER variable with default value 6. It is used as the
C        unit number for the printing of the monitoring information
C        controlled by IPRINT.
C
C    LP  is an INTEGER variable with default value 6. It is used as the
C        unit number for the printing of error messages. This printing
C        may be suppressed by setting LP to a non-positive value.
C
C    GTOL is a DOUBLE PRECISION variable with default value 0.9, which
C        controls the accuracy of the line search routine MCSRCH. If the
C        function and gradient evaluations are inexpensive with respect
C        to the cost of the iteration (which is sometimes the case when
C        solving very large problems) it may be advantageous to set GTOL
C        to a small value. A typical small value is 0.1.  Restriction:
C        GTOL should be greater than 1.D-04.
C
C    STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
C        specify lower and uper bounds for the step in the line search.
C        Their default values are 1.D-20 and 1.D+20, respectively. These
C        values need not be modified unless the exponents are too large
C        for the machine being used, or unless the problem is extremely
C        badly scaled (in which case the exponents should be increased).
C
C
C  MACHINE DEPENDENCIES
C
C        The only variables that are machine-dependent are XTOL,
C        STPMIN and STPMAX.
C
C
C  GENERAL INFORMATION
C
C    Other routines called directly:  DAXPY, DDOT, LB1, MCSRCH
C
C    Input/Output  :  No input; diagnostic messages on unit MP and
C                     error messages on unit LP.
C
C
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
      DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN,
     .                 STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM
      INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO,
     .        BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN
      LOGICAL FINISH
C
      SAVE
      DATA ONE,ZERO/1.0D+0,0.0D+0/
C
C     INITIALIZE
C     ----------
C
      IF(IFLAG.EQ.0) GO TO 10
      GO TO (172,100) IFLAG
  10  ITER= 0
      IF(N.LE.0.OR.M.LE.0) GO TO 196
      IF(GTOL.LE.1.D-04) THEN
        IF(LP.GT.0) WRITE(LP,245)
        GTOL=9.D-01
      ENDIF
      NFUN= 1
      POINT= 0
      FINISH= .FALSE.
      IF(DIAGCO) THEN
         DO 30 I=1,N
 30      IF (DIAG(I).LE.ZERO) GO TO 195
      ELSE
         DO 40 I=1,N
 40      DIAG(I)= 1.0D0
      ENDIF
C
C     THE WORK VECTOR W IS DIVIDED AS FOLLOWS:
C     ---------------------------------------
C     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND
C         OTHER TEMPORARY INFORMATION.
C     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.
C     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED
C         IN THE FORMULA THAT COMPUTES H*G.
C     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH
C         STEPS.
C     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M
C         GRADIENT DIFFERENCES.
C
C     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A
C     CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT.
C
      ISPT= N+2*M
      IYPT= ISPT+N*M
      DO 50 I=1,N
 50   W(ISPT+I)= -G(I)*DIAG(I)
      GNORM= DSQRT(DDOT(N,G,1,G,1))
      STP1= ONE/GNORM
C
C     PARAMETERS FOR LINE SEARCH ROUTINE
C
      FTOL= 1.0D-4
      MAXFEV= 20
C
      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
     *                     GNORM,N,M,X,F,G,STP,FINISH)
C
C    --------------------
C     MAIN ITERATION LOOP
C    --------------------
C
 80   ITER= ITER+1
      INFO=0
      BOUND=ITER-1
      IF(ITER.EQ.1) GO TO 165
      IF (ITER .GT. M)BOUND=M
C
         YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1)
      IF(.NOT.DIAGCO) THEN
         YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1)
         DO 90 I=1,N
   90    DIAG(I)= YS/YY
      ELSE
         IFLAG=2
         RETURN
      ENDIF
 100  CONTINUE
      IF(DIAGCO) THEN
        DO 110 I=1,N
 110    IF (DIAG(I).LE.ZERO) GO TO 195
      ENDIF
C
C     COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
C     "Updating quasi-Newton matrices with limited storage",
C     Mathematics of Computation, Vol.24, No.151, pp. 773-782.
C     ---------------------------------------------------------
C
      CP= POINT
      IF (POINT.EQ.0) CP=M
      W(N+CP)= ONE/YS
      DO 112 I=1,N
 112  W(I)= -G(I)
      CP= POINT
      DO 125 I= 1,BOUND
         CP=CP-1
         IF (CP.EQ. -1)CP=M-1
         SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1)
         INMC=N+M+CP+1
         IYCN=IYPT+CP*N
         W(INMC)= W(N+CP+1)*SQ
         CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1)
 125  CONTINUE
C
      DO 130 I=1,N
 130  W(I)=DIAG(I)*W(I)
C
      DO 145 I=1,BOUND
         YR= DDOT(N,W(IYPT+CP*N+1),1,W,1)
         BETA= W(N+CP+1)*YR
         INMC=N+M+CP+1
         BETA= W(INMC)-BETA
         ISCN=ISPT+CP*N
         CALL DAXPY(N,BETA,W(ISCN+1),1,W,1)
         CP=CP+1
         IF (CP.EQ.M)CP=0
 145  CONTINUE
C
C     STORE THE NEW SEARCH DIRECTION
C     ------------------------------
C
       DO 160 I=1,N
 160   W(ISPT+POINT*N+I)= W(I)
C
C     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION
C     BY USING THE LINE SEARCH ROUTINE MCSRCH
C     ----------------------------------------------------
 165  NFEV=0
      STP=ONE
      IF (ITER.EQ.1) STP=STP1
      DO 170 I=1,N
 170  W(I)=G(I)
 172  CONTINUE
      CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,
     *            XTOL,MAXFEV,INFO,NFEV,DIAG)
      IF (INFO .EQ. -1) THEN
        IFLAG=1
        RETURN
      ENDIF
      IF (INFO .NE. 1) GO TO 190
      NFUN= NFUN + NFEV
C
C     COMPUTE THE NEW STEP AND GRADIENT CHANGE
C     -----------------------------------------
C
      NPT=POINT*N
      DO 175 I=1,N
      W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
 175  W(IYPT+NPT+I)= G(I)-W(I)
      POINT=POINT+1
      IF (POINT.EQ.M)POINT=0
C
C     TERMINATION TEST
C     ----------------
C
      GNORM= DSQRT(DDOT(N,G,1,G,1))
      XNORM= DSQRT(DDOT(N,X,1,X,1))
      XNORM= DMAX1(1.0D0,XNORM)
      IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
C
      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
     *               GNORM,N,M,X,F,G,STP,FINISH)
      IF (FINISH) THEN
         IFLAG=0
         RETURN
      ENDIF
      GO TO 80
C
C     ------------------------------------------------------------
C     END OF MAIN ITERATION LOOP. ERROR EXITS.
C     ------------------------------------------------------------
C
 190  IFLAG=-1
      IF(LP.GT.0) WRITE(LP,200) INFO
      RETURN
 195  IFLAG=-2
      IF(LP.GT.0) WRITE(LP,235) I
      RETURN
 196  IFLAG= -3
      IF(LP.GT.0) WRITE(LP,240)
C
C     FORMATS
C     -------
C
 200  FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE'
     .          ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN'
     .          ' OF LINE SEARCH: INFO= ',I2,/
     .          ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,
     .          ' OR INCORRECT TOLERANCES')
 235  FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,
     .       ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
 240  FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',
     .       ' ARE NOT POSITIVE)')
 245  FORMAT(/'  GTOL IS LESS THAN OR EQUAL TO 1.D-04',
     .       / ' IT HAS BEEN RESET TO 9.D-01')
      RETURN
      END
C
C     LAST LINE OF SUBROUTINE LBFGS
C
C
      SUBROUTINE LB1(IPRINT,ITER,NFUN,
     *                     GNORM,N,M,X,F,G,STP,FINISH)
C
C     -------------------------------------------------------------
C     THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND
C     AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT.
C     -------------------------------------------------------------
C
      INTEGER IPRINT(2),ITER,NFUN,LP,MP,N,M
      DOUBLE PRECISION X(N),G(N),F,GNORM,STP,GTOL,STPMIN,STPMAX
      LOGICAL FINISH
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
C
      IF (ITER.EQ.0)THEN
           WRITE(MP,10)
           WRITE(MP,20) N,M
           WRITE(MP,30)F,GNORM
                 IF (IPRINT(2).GE.1)THEN
                     WRITE(MP,40)
                     WRITE(MP,50) (X(I),I=1,N)
                     WRITE(MP,60)
                     WRITE(MP,50) (G(I),I=1,N)
                  ENDIF
           WRITE(MP,10)
           WRITE(MP,70)
      ELSE
          IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN
              IF (IPRINT(1).NE.0)THEN
                   IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN
                         IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70)
                         WRITE(MP,80)ITER,NFUN,F,GNORM,STP
                   ELSE
                         RETURN
                   ENDIF
              ELSE
                   IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70)
                   WRITE(MP,80)ITER,NFUN,F,GNORM,STP
              ENDIF
              IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN
                    IF (FINISH)THEN
                        WRITE(MP,90)
                    ELSE
                        WRITE(MP,40)
                    ENDIF
                      WRITE(MP,50)(X(I),I=1,N)
                  IF (IPRINT(2).EQ.3)THEN
                      WRITE(MP,60)
                      WRITE(MP,50)(G(I),I=1,N)
                  ENDIF
              ENDIF
            IF (FINISH) WRITE(MP,100)
      ENDIF
C
 10   FORMAT('*************************************************')
 20   FORMAT('  N=',I5,'   NUMBER OF CORRECTIONS=',I2,
     .       /,  '       INITIAL VALUES')
 30   FORMAT(' F= ',1PD10.3,'   GNORM= ',1PD10.3)
 40   FORMAT(' VECTOR X= ')
 50   FORMAT(6(2X,1PD10.3))
 60   FORMAT(' GRADIENT VECTOR G= ')
 70   FORMAT(/'   I   NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
 80   FORMAT(2(I4,1X),3X,3(1PD10.3,2X))
 90   FORMAT(' FINAL POINT X= ')
 100  FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.',
     .       /' IFLAG = 0')
C
      RETURN
      END
C     ******
C
C
C   ----------------------------------------------------------
C     DATA
C   ----------------------------------------------------------
C
      BLOCK DATA LB2
      INTEGER LP,MP
      DOUBLE PRECISION GTOL,STPMIN,STPMAX
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
      DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,9.0D-01,1.0D-20,1.0D+20/
      END
C    ------------------------------------------------------------------
C
C     **************************
C     LINE SEARCH ROUTINE MCSRCH
C     **************************
C
      SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA)
      INTEGER N,MAXFEV,INFO,NFEV
      DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
      DOUBLE PRECISION X(N),G(N),S(N),WA(N)
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
      SAVE
C
C                     SUBROUTINE MCSRCH
C
C     A slight modification of the subroutine CSRCH of More' and Thuente.
C     The changes are to allow reverse communication, and do not affect
C     the performance of the routine.
C
C     THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES
C     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
C
C     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
C     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
C     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
C     MINIMIZER OF THE MODIFIED FUNCTION
C
C          F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
C
C     IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
C     HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
C     THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
C     CONTAINS A MINIMIZER OF F(X+STP*S).
C
C     THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
C     THE SUFFICIENT DECREASE CONDITION
C
C           F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
C
C     AND THE CURVATURE CONDITION
C
C           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
C
C     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
C     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
C     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
C     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
C     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
C     SATISFIES THE SUFFICIENT DECREASE CONDITION.
C
C     THE SUBROUTINE STATEMENT IS
C
C        SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF VARIABLES.
C
C       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
C         X + STP*S.
C
C       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
C         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
C
C       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
C         OF F AT X + STP*S.
C
C       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
C         SEARCH DIRECTION.
C
C       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
C         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
C         STP CONTAINS THE FINAL ESTIMATE.
C
C       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse
C         communication implementation GTOL is defined in a COMMON
C         statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE
C         CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
C         SATISFIED.
C
C       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
C         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C         IS AT MOST XTOL.
C
C       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
C         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse
C         communication implementatin they are defined in a COMMON
C         statement).
C
C       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
C         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
C         MAXFEV BY THE END OF AN ITERATION.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C
C         INFO = 0  IMPROPER INPUT PARAMETERS.
C
C         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
C
C         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
C                   DIRECTIONAL DERIVATIVE CONDITION HOLD.
C
C         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C                   IS AT MOST XTOL.
C
C         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
C
C         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
C
C         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
C
C         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C                   THERE MAY NOT BE A STEP WHICH SATISFIES THE
C                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
C                   TOLERANCES MAY BE TOO SMALL.
C
C       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
C         CALLS TO FCN.
C
C       WA IS A WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       MCSTEP
C
C       FORTRAN-SUPPLIED...ABS,MAX,MIN
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C     JORGE J. MORE', DAVID J. THUENTE
C
C     **********
      INTEGER INFOC,J
      LOGICAL BRACKT,STAGE1
      DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO
      DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
      IF(INFO.EQ.-1) GO TO 45
      INFOC = 1
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.
     *    GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO
     *    .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN
C
C     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
C     AND CHECK THAT S IS A DESCENT DIRECTION.
C
      DGINIT = ZERO
      DO 10 J = 1, N
         DGINIT = DGINIT + G(J)*S(J)
   10    CONTINUE
      IF (DGINIT .GE. ZERO) then
         write(LP,15)
   15    FORMAT(/'  THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION')
         RETURN
         ENDIF
C
C     INITIALIZE LOCAL VARIABLES.
C
      BRACKT = .FALSE.
      STAGE1 = .TRUE.
      NFEV = 0
      FINIT = F
      DGTEST = FTOL*DGINIT
      WIDTH = STPMAX - STPMIN
      WIDTH1 = WIDTH/P5
      DO 20 J = 1, N
         WA(J) = X(J)
   20    CONTINUE
C
C     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
C     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
C     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
C     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
C     THE INTERVAL OF UNCERTAINTY.
C     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
C     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
C
      STX = ZERO
      FX = FINIT
      DGX = DGINIT
      STY = ZERO
      FY = FINIT
      DGY = DGINIT
C
C     START OF ITERATION.
C
   30 CONTINUE
C
C        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
C        TO THE PRESENT INTERVAL OF UNCERTAINTY.
C
         IF (BRACKT) THEN
            STMIN = MIN(STX,STY)
            STMAX = MAX(STX,STY)
         ELSE
            STMIN = STX
            STMAX = STP + XTRAPF*(STP - STX)
            END IF
C
C        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
C
         STP = MAX(STP,STPMIN)
         STP = MIN(STP,STPMAX)
C
C        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
C        STP BE THE LOWEST POINT OBTAINED SO FAR.
C
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
     *      .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0
     *      .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
C
C        EVALUATE THE FUNCTION AND GRADIENT AT STP
C        AND COMPUTE THE DIRECTIONAL DERIVATIVE.
C        We return to main program to obtain F and G.
C
         DO 40 J = 1, N
            X(J) = WA(J) + STP*S(J)
   40       CONTINUE
         INFO=-1
         RETURN
C
   45    INFO=0
         NFEV = NFEV + 1
         DG = ZERO
         DO 50 J = 1, N
            DG = DG + G(J)*S(J)
   50       CONTINUE
         FTEST1 = FINIT + STP*DGTEST
C
C        TEST FOR CONVERGENCE.
C
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
     *      .OR. INFOC .EQ. 0) INFO = 6
         IF (STP .EQ. STPMAX .AND.
     *       F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
         IF (STP .EQ. STPMIN .AND.
     *       (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
         IF (NFEV .GE. MAXFEV) INFO = 3
         IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
         IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1
C
C        CHECK FOR TERMINATION.
C
         IF (INFO .NE. 0) RETURN
C
C        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
C        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
C
         IF (STAGE1 .AND. F .LE. FTEST1 .AND.
     *       DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
C
C        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
C        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
C        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
C        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
C        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
C
         IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
C
C           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
C
            FM = F - STP*DGTEST
            FXM = FX - STX*DGTEST
            FYM = FY - STY*DGTEST
            DGM = DG - DGTEST
            DGXM = DGX - DGTEST
            DGYM = DGY - DGTEST
C
C           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
C           AND TO COMPUTE THE NEW STEP.
C
            CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,
     *                 BRACKT,STMIN,STMAX,INFOC)
C
C           RESET THE FUNCTION AND GRADIENT VALUES FOR F.
C
            FX = FXM + STX*DGTEST
            FY = FYM + STY*DGTEST
            DGX = DGXM + DGTEST
            DGY = DGYM + DGTEST
         ELSE
C
C           CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
C           AND TO COMPUTE THE NEW STEP.
C
            CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,
     *                 BRACKT,STMIN,STMAX,INFOC)
            END IF
C
C        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
C        INTERVAL OF UNCERTAINTY.
C
         IF (BRACKT) THEN
            IF (ABS(STY-STX) .GE. P66*WIDTH1)
     *         STP = STX + P5*(STY - STX)
            WIDTH1 = WIDTH
            WIDTH = ABS(STY-STX)
            END IF
C
C        END OF ITERATION.
C
         GO TO 30
C
C     LAST LINE OF SUBROUTINE MCSRCH.
C
      END
      SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
     *                 STPMIN,STPMAX,INFO)
      INTEGER INFO
      DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX
      LOGICAL BRACKT,BOUND
C
C     SUBROUTINE MCSTEP
C
C     THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR
C     A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
C     A MINIMIZER OF THE FUNCTION.
C
C     THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
C     VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
C     ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
C     DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
C     MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
C     WITH ENDPOINTS STX AND STY.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
C                        STPMIN,STPMAX,INFO)
C
C     WHERE
C
C       STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
C         SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
C         OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
C         SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
C
C       STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
C         THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
C         UPDATED APPROPRIATELY.
C
C       STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
C         IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
C         BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
C
C       BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER
C         HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
C         THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
C         IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
C
C       STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
C         AND UPPER BOUNDS FOR THE STEP.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C         IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
C         ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
C         INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C     JORGE J. MORE', DAVID J. THUENTE
C
      DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
      INFO = 0
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR.
     *    STP .GE. MAX(STX,STY))) .OR.
     *    DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN
C
C     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
C
      SGND = DP*(DX/ABS(DX))
C
C     FIRST CASE. A HIGHER FUNCTION VALUE.
C     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
C     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
C     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
C
      IF (FP .GT. FX) THEN
         INFO = 1
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .LT. STX) GAMMA = -GAMMA
         P = (GAMMA - DX) + THETA
         Q = ((GAMMA - DX) + GAMMA) + DP
         R = P/Q
         STPC = STX + R*(STP - STX)
         STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX)
         IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN
            STPF = STPC
         ELSE
           STPF = STPC + (STPQ - STPC)/2
           END IF
         BRACKT = .TRUE.
C
C     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
C     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
C     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
C     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
C
      ELSE IF (SGND .LT. 0.0) THEN
         INFO = 2
         BOUND = .FALSE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = ((GAMMA - DP) + GAMMA) + DX
         R = P/Q
         STPC = STP + R*(STX - STP)
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN
            STPF = STPC
         ELSE
            STPF = STPQ
            END IF
         BRACKT = .TRUE.
C
C     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
C     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
C     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
C     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
C     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
C     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
C     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
C     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
C
      ELSE IF (ABS(DP) .LT. ABS(DX)) THEN
         INFO = 3
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
C
C        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
C        TO INFINITY IN THE DIRECTION OF THE STEP.
C
         GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = (GAMMA + (DX - DP)) + GAMMA
         R = P/Q
         IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN
            STPC = STP + R*(STX - STP)
         ELSE IF (STP .GT. STX) THEN
            STPC = STPMAX
         ELSE
            STPC = STPMIN
            END IF
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (BRACKT) THEN
            IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
         ELSE
            IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
            END IF
C
C     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
C     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
C     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
C     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
C
      ELSE
         INFO = 4
         BOUND = .FALSE.
         IF (BRACKT) THEN
            THETA = 3*(FP - FY)/(STY - STP) + DY + DP
            S = MAX(ABS(THETA),ABS(DY),ABS(DP))
            GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S))
            IF (STP .GT. STY) GAMMA = -GAMMA
            P = (GAMMA - DP) + THETA
            Q = ((GAMMA - DP) + GAMMA) + DY
            R = P/Q
            STPC = STP + R*(STY - STP)
            STPF = STPC
         ELSE IF (STP .GT. STX) THEN
            STPF = STPMAX
         ELSE
            STPF = STPMIN
            END IF
         END IF
C
C     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
C     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
C
      IF (FP .GT. FX) THEN
         STY = STP
         FY = FP
         DY = DP
      ELSE
         IF (SGND .LT. 0.0) THEN
            STY = STX
            FY = FX
            DY = DX
            END IF
         STX = STP
         FX = FP
         DX = DP
         END IF
C
C     COMPUTE THE NEW STEP AND SAFEGUARD IT.
C
      STPF = MIN(STPMAX,STPF)
      STPF = MAX(STPMIN,STPF)
      STP = STPF
      IF (BRACKT .AND. BOUND) THEN
         IF (STY .GT. STX) THEN
            STP = MIN(STX+0.66*(STY-STX),STP)
         ELSE
            STP = MAX(STX+0.66*(STY-STX),STP)
            END IF
         END IF
      RETURN
C
C     LAST LINE OF SUBROUTINE MCSTEP.
C
      END

ext: f


text:
*----------------------------------------------------------------------*
      logical function iocc_bound(cbound,iocc,dagi,jocc,dagj)
*----------------------------------------------------------------------*
*     compare occupations element-wise with relation cbound
*     ('>','>=','<','<='). .TRUE. only if .TRUE. for all elements
*----------------------------------------------------------------------*

      implicit none
      include 'opdim.h'

      character, intent(in) ::
     &     cbound*(*)
      logical, intent(in) ::
     &     dagi, dagj
      integer, intent(in) ::
     &     iocc(ngastp,2), jocc(ngastp,2)

      integer ::
     &     ica, ica_i, ica_j, ihpv

      iocc_bound = .true.
      select case(trim(cbound))
      case('<','LT','lt')
        do ica = 1,2
          ica_i = ica
          ica_j = ica
          if (dagi) ica_i = 3-ica
          if (dagj) ica_j = 3-ica
          do ihpv = 1,ngastp
            iocc_bound = iocc_bound.and.
     &           iocc(ihpv,ica_i).lt.jocc(ihpv,ica_j)
          end do
        end do
      case('<=','LE','le')
        do ica = 1,2
          ica_i = ica
          ica_j = ica
          if (dagi) ica_i = 3-ica
          if (dagj) ica_j = 3-ica
          do ihpv = 1,ngastp
            iocc_bound = iocc_bound.and.
     &           iocc(ihpv,ica_i).le.jocc(ihpv,ica_j)
          end do
        end do
      case('>','GT','gt')
        do ica = 1,2
          ica_i = ica
          ica_j = ica
          if (dagi) ica_i = 3-ica
          if (dagj) ica_j = 3-ica
          do ihpv = 1,ngastp
            iocc_bound = iocc_bound.and.
     &           iocc(ihpv,ica_i).gt.jocc(ihpv,ica_j)
          end do
        end do
      case('>=','GE','ge')
        do ica = 1,2
          ica_i = ica
          ica_j = ica
          if (dagi) ica_i = 3-ica
          if (dagj) ica_j = 3-ica
          do ihpv = 1,ngastp
            iocc_bound = iocc_bound.and.
     &           iocc(ihpv,ica_i).ge.jocc(ihpv,ica_j)
          end do
        end do
      case default
        call quit(1,'iocc_bound','Unknow relation: "'//cbound//'"')
      end select

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
! Copyright (C) 2015-2017, Simone Marsili
! All rights reserved.
! License: BSD 3 clause

program fmpl
  use kinds
  use constants,     only: long_string
  use command_line,  only: read_args
  use data,          only: initialize_data,read_data,data_reweight
  use model,         only: initialize_model, model_set_myv,fix_gauge
  use scrs,          only: print_mat, compute_scores
  use dvmlm_wrapper, only: dvmlm_minimize
  implicit none
  ! input variables
  character(long_string) :: data_file, scores_file
  real(kflt)             :: w_id, lambda
  integer                :: ignore_pivot, accuracy
  character(1)           :: scores_format
  ! data dimensions
  integer :: nd ! n. of samples
  integer :: nv ! n. of variables
  integer :: ns ! n. of classes
  real(kflt) :: neff
  ! main arrays for the run
  integer,    allocatable :: data_samples(:,:)  ! nv x nd
  real(kflt), allocatable :: w(:)
  real(kflt), allocatable :: prm(:,:)           ! (ns + ns x ns x nv) x nv
  real(kflt), allocatable :: grd(:)             ! (ns + ns x ns x nv) x nv
  real(kflt), allocatable :: fields(:,:)        ! ns x nv
  real(kflt), allocatable :: couplings(:,:,:,:) ! ns x ns x nv x nv
  real(kflt), allocatable :: scores(:,:)        ! nv x nv
  !
  integer                :: udata,uscrs
  integer                :: err,iv
  integer                :: niter,neval
  real(kflt_single)      :: finish,start,start_min,tpv,elapsed_time,expected_time
  character(4) :: time_unit

  ! get command line
  call read_args(data_file,w_id,lambda,ignore_pivot,accuracy,scores_format,err)
  if(err /= 0) then
     write(0,100)
     stop
  end if
  ! print an header for the run
  write(0,101) trim(data_file), lambda, accuracy
  if (w_id > 0.0_kflt) write(0,102) w_id

  ! open data file
  open(newunit=udata,file=data_file,status='old',action='read',iostat=err)
  if(err /= 0) then
     write(0,'("error ! cannot access ",a,": file not found")') trim(data_file)
     stop
  end if

  ! open scores file
  scores_file = trim(data_file)//'.scores'
  open(newunit=uscrs,file=scores_file,status='replace',action='write',iostat=err)

  call initialize_data(udata,nd,nv,neff)
  allocate(data_samples(nv,nd),stat=err)
  allocate(w(nd),stat=err)

  write(0,'(/a)') "Reading data.."
  call read_data(udata,w_id,ns,neff,data_samples)
  write(0,103) nd, nv, ns
  flush(0)
  if(w_id > 1.E-10_kflt) then
     write(0,'(a)') "Computing weights.."
     call data_reweight(data_samples,w_id,neff,w)
     write(0,'(a,f8.1)') "Neff: ",neff
  else
     neff = nd
     w = 1.0_kflt / neff
  end if

  ! allocate parameters and gradient
  allocate(prm(ns + ns*ns*nv,nv),stat=err)
  allocate(grd(ns + ns*ns*nv),stat=err)
  allocate(scores(nv,nv),stat=err)
  call initialize_model(nv,ns,lambda)

  write(0,'(/,a)') 'Running..'
  call cpu_time(start_min)
  tpv = 0.0_kflt
  ! loop over features
  do iv = 1,nv
     call model_set_myv(nd,nv,iv,data_samples,w,prm(:,iv),grd,err)
     niter = 0
     call cpu_time(start)
     call dvmlm_minimize(nv,ns,nd,size(prm(:,iv)),data_samples,w,prm(:,iv),grd,accuracy,niter,neval)
     call cpu_time(finish)
     elapsed_time = finish - start_min
     tpv = elapsed_time / real(iv)
     expected_time = tpv * (nv - iv)
     time_unit = "secs"
     if (elapsed_time + expected_time > 2*60.) then
        elapsed_time = elapsed_time / 60.0_kflt
        expected_time = expected_time / 60.0_kflt
        time_unit = "mins"
     end if
     if (mod(iv,int(3.0/tpv)+1)==0) &
          write(0,'(i4,"/",i4," completed in ",f5.1,1x,a,"; ",f5.1," to end")')&
          iv, nv, elapsed_time, time_unit, expected_time
     if (iv == nv) &
          write(0,'("*** ",i4,"/",i4," completed in ",f9.3,1x,a," ***")')&
          iv, nv, elapsed_time, "secs"
     !write(0,'(a,i5,a,2i5,a,f8.3,a)') ' variable ', iv, &
     !'  converged (niter,neval) ', niter, neval, ' in ', finish-start, ' secs'
     !call fix_gauge(nv,ns,prm(:ns,iv),prm(ns+1:,iv))
  end do
  flush(0)

  ! reorder prm array into fields and couplings
  allocate(fields(ns,nv),couplings(ns,ns,nv,nv),stat=err)
  call reshape_prm(nv,ns,prm,fields,couplings)
  deallocate(prm,grd)

  ! compute scores and print
  call compute_scores(nv,ns,couplings,ignore_pivot,scores)
  call print_mat(scores,uscrs,scores_format,err)
  if (err /= 0) stop

101 format(&
         '# fmpl                                                             '/&
         '#                                                                  '/&
         '# input data file: ',a,'                                           '/&
         '# regularization: ',f6.4,'                                         '/&
         '# accuracy level: ',i1,'                                           ')

102 format(&
         '# reweight samples; % id threshold: ',f5.1,'                       ')

103 format(&
         'Sample size         : ',i6,'                                       '/&
         'Dimensionality      : ',i6,'                                       '/&
         'Classes per variable: ',i6,'                                       ')

100 format(/&
         'fmpl                                                           '/&
         '                                                               '/&
         'Usage:                                                         '/&
         '    fmpl [options] -i <file>                                   '/&
         '                                                               '/&
         'Options:                                                       '/&
         '-h, --help                                                     '/&
         '    Display this help and exit.                                '/&
         '                                                               '/&
         '-l, --lambda <regularization_parameter> float, optional        '/&
         '    Controls L_2 regularization strength.                      '/&
         '    [default: 0.01]                                            '/&
         '                                                               '/&
         '-w, --reweight <percentage_identity> float, optional           '/&
         '    Reweight data using a percentage identity threshold.       '/&
         '    [default: no reweight.]                                    '/&
         '                                                               '/&
         '-a, --accuracy <accuracy_level> integer, optional              '/&
         '    Larger values correspond to stricter convergence criteria  '/&
         '    and longer covergence time. Possible values are {0, 1, 2}. '/&
         '    [default: 1.]                                              '/&
         '                                                               '/&
         '--ignore_pivot <pivot_state> integer, optional                 '/&
         '    Ignore the contribution of <pivot_state> to final scores.  '/&
         '    [default: include all states.]                             '/&
         '                                                               '/&
         '--scores_format <scores_matrix_format> - string, optional      '/&
         '    Possible values are {"rcv", "array", "coordinate"}.        '/&
         '    "rcv": (row_index, column_index, value) format             '/&
         '    "coordinate", "array": Matrix Market (MM) formats. See:    '/&
         '    https://people.sc.fsu.edu/~jburkardt/data/mm/mm.html       '/&
         '    [default: "rcv"].                                          '/)

contains

  subroutine reshape_prm(nv,ns,prm,fields,couplings)
    ! reorder prm array into fields and couplings
    implicit none
    integer, intent(in) :: nv,ns
    real(kflt), intent(in) :: prm(:,:)
    real(kflt), intent(out) :: fields(ns,nv)
    real(kflt), intent(out) :: couplings(ns,ns,nv,nv)
    integer :: err
    integer :: iv,jv,is,js,k1,k2,k

    k = 0
    do iv = 1,nv
       fields(:,iv) = prm(:ns,iv)
       k = ns
       do jv = 1,nv
          do js = 1,ns
             do is = 1,ns
                k = k + 1
                couplings(js,is,jv,iv) = prm(k,iv)
             end do
          end do
       end do
    end do
  end subroutine reshape_prm

end program fmpl

max_stars_repo_name: simomarsili/fmpl

max_issues_repo_name: simomarsili/fmpl

max_forks_repo_name: simomarsili/fmpl


text:
The Wastewater Treatment Plant is located in the rural outskirts of Davis on County Road 28H. Human and miscellaneous waste products are conveyed to the plant through the mighty powers of the local sewers sewer system. Upon reaching its destination, the raw sewage is treated, and the resulting sludge is pumped out to the neighboring Davis Wetlands, which is just east of the plant.

Treatment Plant Exposure is a huge problem to the local wildlife.  Many animals have died and will continue to die due to it.  In the first two weeks of June, 2005, a good Samaritan brought in  Ducks ducks he found not moving near the area to the Vet Med Teaching Hospital VMTH, along with a few deceased corpses to be disposed of.  Some died quickly, while the others were treated at a relatively great expense. Problems include dehydration, avian cholera, botulism, heavy metal toxicity, and all kinds of microbial infections, especially gastric and respiratory ones. Poor things.

The old waste water treatment plant was located at 24662 County Road 102 in an unincorporated area of Yolo County. A solar farm known as PVUSA was built on the location in 1986.

See also: UC Davis Wastewater Treatment Plant

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
        subroutine idd_pairsamps(n,l,ind,l2,ind2,marker)
c
c       calculates the indices of the l2 pairs of integers
c       to which the l individual integers from ind belong.
c       The integers in ind may range from 1 to n.
c
c       input:
c       n -- upper bound on the integers in ind
c            (the number 1 must be a lower bound);
c            n must be even
c       l -- length of ind
c       ind -- integers selected from 1 to n
c
c       output:
c       l2 -- length of ind2
c       ind2 -- indices in the range from 1 to n/2 of the pairs
c               of integers to which the entries of ind belong
c
c       work:
c       marker -- must be at least n/2 integer elements long
c
c       _N.B._: n must be even.
c
        implicit none
        integer l,n,ind(l),ind2(l),marker(n/2),l2,k
c
c
c       Unmark all pairs.
c
        do k = 1,n/2
          marker(k) = 0
        enddo ! k
c
c
c       Mark the required pairs.
c
        do k = 1,l
          marker((ind(k)+1)/2) = marker((ind(k)+1)/2)+1
        enddo ! k
c
c
c       Record the required pairs in indpair.
c
        l2 = 0
c
        do k = 1,n/2
c
          if(marker(k) .ne. 0) then
            l2 = l2+1
            ind2(l2) = k
          endif
c
        enddo ! k
c
c
        return
        end
c
c
c
c


text:
      module linkedListMod
!
!     This module supports a linked list class.
!
!
!     Define class.
!
      implicit none
      TYPE::real_value
        real::value
        type(real_value),pointer::prev,next
      end TYPE real_value
!
      TYPE::linkedListReal
        type(real_value),pointer::head=>null()
        type(real_value),pointer::tail=>null()
        type(real_value),pointer::current=>null()
        logical::init=.false.
      end TYPE linkedListReal
!
!
!     Procedure Subroutines and Functions
!
      CONTAINS
!
      subroutine linkedListReal_incrementCurrent(myLinkedList,incrementBy)
!
!     This subroutine is used to increment the current pointer of a
!     linkedListReal object. If current has not yet been set to point at the
!     linked list, then put the current pointer to the head of the list before
!     applying the increment. If current is at the end of the list, do NOT
!     increment further.
!
      implicit none
      type(linkedListReal)::myLinkedList
      integer,optional,intent(in)::incrementBy
      integer::myIncrementBy,i
!
!     Set up myIncrementBy.
!
      myIncrementBy = 1
      if(PRESENT(incrementBy)) myIncrementBy = incrementBy
!
      if(.not.ASSOCIATED(myLinkedList%current)) myLinkedList%current=>myLinkedList%head
      if(myIncrementBy.gt.0) then
        do i = 1,myIncrementBy
          if(ASSOCIATED(myLinkedList%current%next)) myLinkedList%current=>myLinkedList%current%next
        endDo
      elseIf(myIncrementBy.lt.0) then
        do i = 1,ABS(myIncrementBy),-1
          if(ASSOCIATED(myLinkedList%current%next)) myLinkedList%current=>myLinkedList%current%next
        endDo
      endIf
!
      return
      end subroutine linkedListReal_incrementCurrent

      function linkedListReal_getCurrent(myLinkedList)
!
!     This function returns an intrinsic real value equal to the current node
!     value. The current position in the linked list always begins at the head
!     of the list and increments each time this function is used. When it
!     reaches the tail of the list, it just stays at the tail.
!
      implicit none
      type(linkedListReal)::myLinkedList
      real::linkedListReal_getCurrent
!
      if(.not.ASSOCIATED(myLinkedList%current)) myLinkedList%current=>myLinkedList%head
      linkedListReal_getCurrent = myLinkedList%current%value
      call linkedListReal_incrementCurrent(myLinkedList)
!
      return
      end function linkedListReal_getCurrent

      function linkedListReal_hasNext(myLinkedList)
!
!     This function returns a logical (TRUE or FALSE) indicating whether or not
!     the current member of the linked list object has a "next" node to go to.
!
      implicit none
      type(linkedListReal)::myLinkedList
      logical::linkedListReal_hasNext
!
      if(.not.ASSOCIATED(myLinkedList%current)) myLinkedList%current=>myLinkedList%head
      linkedListReal_hasNext = ASSOCIATED(myLinkedList%current%next)
!
      return
      end function linkedListReal_hasNext

      subroutine linkedListReal_push(myLinkedList,newValue)
!
!     This subroutine executes push onto the linked list called myLinkedList.
!
!
      implicit none
      type(linkedListReal)::myLinkedList
      real,intent(in)::newValue
      type(real_value),pointer::ptr=>null()
      integer::istat
!
!     Put newValue at the end of the linked list...
!
      if(.not.ASSOCIATED(myLinkedList%head)) then
        ALLOCATE(myLinkedList%head,STAT=istat)
        myLinkedList%tail => myLinkedList%head
        NULLIFY(myLinkedList%tail%prev,myLinkedList%tail%next)
        myLinkedList%tail%value = newValue
      else
        ALLOCATE(myLinkedList%tail%next,STAT=istat)
        ptr => myLinkedList%tail
        myLinkedList%tail => myLinkedList%tail%next
        myLinkedList%tail%prev => ptr
        NULLIFY(myLinkedList%tail%next)
        myLinkedList%tail%value = newValue
      endIf
!
      return
      end subroutine linkedListReal_push

      end module linkedListMod

max_stars_repo_name: MQCPack/summerCodingWorkshop

max_issues_repo_name: MQCPack/summerCodingWorkshop

max_forks_repo_name: MQCPack/summerCodingWorkshop


text:
!-----------------------------------------------------------------------
!    MODULE:       el_eos_modlue
!    Author:       S. W. Bruenn
!    Date:         10/26/02
!-----------------------------------------------------------------------

MODULE el_eos_module

USE wlKindModule, ONLY: double => dp

!-----------------------------------------------------------------------
!  NSUBE, QSUBE - Electron density and concentration
!
!  NEPLUS - Positron density (Not set in our simple ultrarelativistic EOS)
!
!  ACOEF, BCOEF, CCOEF, DBDT - Temporary coefficients
!
!  MUSUBE - Electron chemical potential
!
!  EPRESS, EU, FSUBE, ES - Electron pressure, internal energy, Helmholtz
!   free energy, and entropy
!
!  PPRESS, PS, PU, PF - Photon pressure, internal energy, Helmholtz free
!   energy, and entropy
!-----------------------------------------------------------------------

REAL(KIND=double)                               :: NSUBE, QSUBE
REAL(KIND=double)                               :: NEPLUS
REAL(KIND=double)                               :: ACOEF, BCOEF, CCOEF, DBDT
REAL(KIND=double)                               :: MUSUBE
REAL(KIND=double)                               :: EPRESS, EU, FSUBE, ES
REAL(KIND=double)                               :: PPRESS, PS, PU, PF
REAL(KIND=double)                               :: DEMUDT, DEMUDN, DEMUDY
REAL(KIND=double)                               :: DEPDT, DEPDN, DEPDY
REAL(KIND=double)                               :: DESDT, DESDN, DESDY
REAL(KIND=double)                               :: DEUDT, DEUDN, DEUDY
REAL(KIND=double)                               :: DPPDN, DPPDT, DPPDY
REAL(KIND=double)                               :: DPSDN, DPSDT, DPSDY
REAL(KIND=double)                               :: DPUDN, DPUDT, DPUDY

END MODULE el_eos_module

max_stars_repo_name: srichers/weaklib


text:
*DECK DCFODE
      SUBROUTINE DCFODE (METH, ELCO, TESCO)
C***BEGIN PROLOGUE  DCFODE
C***SUBSIDIARY
C***PURPOSE  Set ODE integrator coefficients.
C***TYPE      DOUBLE PRECISION (SCFODE-S, DCFODE-D)
C***AUTHOR  Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C  DCFODE is called by the integrator routine to set coefficients
C  needed there.  The coefficients for the current method, as
C  given by the value of METH, are set for all orders and saved.
C  The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
C  (A smaller value of the maximum order is also allowed.)
C  DCFODE is called once at the beginning of the problem,
C  and is not called again unless and until METH is changed.
C
C  The ELCO array contains the basic method coefficients.
C  The coefficients el(i), 1 .le. i .le. nq+1, for the method of
C  order nq are stored in ELCO(i,nq).  They are given by a genetrating
C  polynomial, i.e.,
C      l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
C  For the implicit Adams methods, l(x) is given by
C      dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1),    l(-1) = 0.
C  For the BDF methods, l(x) is given by
C      l(x) = (x+1)*(x+2)* ... *(x+nq)/K,
C  where         K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
C
C  The TESCO array contains test constants used for the
C  local error test and the selection of step size and/or order.
C  At order nq, TESCO(k,nq) is used for the selection of step
C  size at order nq - 1 if k = 1, at order nq if k = 2, and at order
C  nq + 1 if k = 3.
C
C***SEE ALSO  DLSODE
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   791129  DATE WRITTEN
C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
C   890503  Minor cosmetic changes.  (FNF)
C   930809  Renamed to allow single/double precision versions. (ACH)
C***END PROLOGUE  DCFODE
C**End
      INTEGER METH
      INTEGER I, IB, NQ, NQM1, NQP1
      DOUBLE PRECISION ELCO, TESCO
      DOUBLE PRECISION AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ,
     1   RQFAC, RQ1FAC, TSIGN, XPIN
      DIMENSION ELCO(13,12), TESCO(3,12)
      DIMENSION PC(12)
C
C***FIRST EXECUTABLE STATEMENT  DCFODE
      GO TO (100, 200), METH
C
 100  ELCO(1,1) = 1.0D0
      ELCO(2,1) = 1.0D0
      TESCO(1,1) = 0.0D0
      TESCO(2,1) = 2.0D0
      TESCO(1,2) = 1.0D0
      TESCO(3,12) = 0.0D0
      PC(1) = 1.0D0
      RQFAC = 1.0D0
      DO 140 NQ = 2,12
C-----------------------------------------------------------------------
C The PC array will contain the coefficients of the polynomial
C     p(x) = (x+1)*(x+2)*...*(x+nq-1).
C Initially, p(x) = 1.
C-----------------------------------------------------------------------
        RQ1FAC = RQFAC
        RQFAC = RQFAC/NQ
        NQM1 = NQ - 1
        FNQM1 = NQM1
        NQP1 = NQ + 1
C Form coefficients of p(x)*(x+nq-1). ----------------------------------
        PC(NQ) = 0.0D0
        DO 110 IB = 1,NQM1
          I = NQP1 - IB
 110      PC(I) = PC(I-1) + FNQM1*PC(I)
        PC(1) = FNQM1*PC(1)
C Compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
        PINT = PC(1)
        XPIN = PC(1)/2.0D0
        TSIGN = 1.0D0
        DO 120 I = 2,NQ
          TSIGN = -TSIGN
          PINT = PINT + TSIGN*PC(I)/I
 120      XPIN = XPIN + TSIGN*PC(I)/(I+1)
C Store coefficients in ELCO and TESCO. --------------------------------
        ELCO(1,NQ) = PINT*RQ1FAC
        ELCO(2,NQ) = 1.0D0
        DO 130 I = 2,NQ
 130      ELCO(I+1,NQ) = RQ1FAC*PC(I)/I
        AGAMQ = RQFAC*XPIN
        RAGQ = 1.0D0/AGAMQ
        TESCO(2,NQ) = RAGQ
        IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1
        TESCO(3,NQM1) = RAGQ
 140    CONTINUE
      RETURN
C
 200  PC(1) = 1.0D0
      RQ1FAC = 1.0D0
      DO 230 NQ = 1,5
C-----------------------------------------------------------------------
C The PC array will contain the coefficients of the polynomial
C     p(x) = (x+1)*(x+2)*...*(x+nq).
C Initially, p(x) = 1.
C-----------------------------------------------------------------------
        FNQ = NQ
        NQP1 = NQ + 1
C Form coefficients of p(x)*(x+nq). ------------------------------------
        PC(NQP1) = 0.0D0
        DO 210 IB = 1,NQ
          I = NQ + 2 - IB
 210      PC(I) = PC(I-1) + FNQ*PC(I)
        PC(1) = FNQ*PC(1)
C Store coefficients in ELCO and TESCO. --------------------------------
        DO 220 I = 1,NQP1
 220      ELCO(I,NQ) = PC(I)/PC(2)
        ELCO(2,NQ) = 1.0D0
        TESCO(1,NQ) = RQ1FAC
        TESCO(2,NQ) = NQP1/ELCO(1,NQ)
        TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ)
        RQ1FAC = RQ1FAC/FNQ
 230    CONTINUE
      RETURN
C----------------------- END OF SUBROUTINE DCFODE ----------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg