text:
!** Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
!** See https://llvm.org/LICENSE.txt for license information.
!** SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

!* Tests for runtime library MATMUL routines

program p

  parameter(NbrTests=208)

  integer*4, dimension(4,3) :: arr1
  integer*4, dimension(3) :: arr2
  integer*4, dimension(4) :: arr3
  integer*4, dimension(4,4) :: arr4
  integer*4, dimension(0:3,-1:1) :: arr5
  integer*4, dimension(-3:-1) :: arr6
  integer*4, dimension(-1:2,0:3) :: arr7
  integer*4, dimension(2:5,3) :: arr8
  integer*4, dimension(2:4) :: arr9
  integer*4, dimension(2:5) :: arr10
  integer*4, dimension(4,2:4) :: arr11
  integer*4, dimension(2:5) :: arr12

  data arr1 /0,1,2,3,4,5,6,7,8,9,10,11/
  data arr5 /0,1,2,3,4,5,6,7,8,9,10,11/
  data arr2 /0,1,2/
  data arr6 /0,1,2/
  data arr3 /0,1,2,3/
  data arr4 /0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15/
  data arr7 /0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15/
  data arr8 /0,1,2,3,4,5,6,7,8,9,10,11/
  data arr9 /0,1,2/
  data arr10 /0,1,2,3/
  data arr11 /0,1,2,3,4,5,6,7,8,9,10,11/

  integer*4 :: expect(NbrTests)
  integer*4 :: results(NbrTests)

  data expect /  &
  ! tests 1-4
      20, 23, 26, 29,  &
  ! tests 5-8
      0, 23, 26, 29,  &
  ! tests 9-12
      23, 26, 29, 0,  &
  ! tests 13-16
      0, 20, 23, 26,  &
  ! tests 17-20
      0, 5, 6, 7,  &
  ! tests 21-24
      8, 9, 10, 11,  &
  ! tests 25-28
      8, 11, 14, 17,  &
  ! tests 29-32
      0, 8, 11, 14,  &
  ! tests 33-36
      9, 10, 11, 0,  &
  ! tests 37-40
      18, 20, 22, 0,  &
  ! tests 41-44
      0, 8, 0, 14,  &
  ! tests 45-48
      9, 0, 11,  &
      0,  &
  ! tests 49-64
      0, 9, 0, 0, 0, 0, 0, 0, 0, 11, 0, 0,  &
      0, 0, 0, 0,  &
  ! tests 65-80
      0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 11, 0,  &
      0, 0, 0, 0,  &
  ! tests 81-96
      0, 9, 0, 0, 0, 0, 0, 0, 0, 11, 0, 0,  &
      0, 0, 0, 0,  &
  ! tests 97-112
      0, 0, 0, 0, 0, 0, 0, 0, 23, 0, 29, 0,  &
      0, 0, 0, 0,  &
  ! tests 113-128
      22, 20, 18, 0,  &
  ! tests 129-145
      9, 0, 11, 0,  &
  ! tests 145-160
      0, 11, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0,  &
      0, 0, 0, 0,  &
  ! tests 161-176
      0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 11, 0,  &
      0, 0, 0, 0,  &
  ! tests 177-192
      0, 0, 0, 0, 0, 0, 0, 0, 23, 0, 29, 0,  &
      0, 0, 0, 0,  &
  ! tests 193-196
      18, 20, 22, 0, &
  ! tests 173-176
     20, 23, 26, 29,  &
  ! tests  177-180
     20, 23, 26, 29,  &
  ! tests 181-184
      23, 26, 29, 0,  &
  ! tests 185-188
      0, 20, 23, 26,  &
  ! tests 189-192
      0, 5, 6, 7,  &
  ! tests 193-196
      8, 9, 10, 11,  &
  ! tests 197-200
      8, 11, 14, 17,  &
  ! tests 201-204
      0, 8, 11, 14,  &
  ! tests 205-208
      9, 10, 11, 0/

  results = -1

  ! tests 1-4
  arr3=0
  arr3 = matmul(arr1,arr2)
  call assign_result(1,4,arr3,results)
  !print *,arr3

  ! tests 5-8
  arr3=0
  arr3(2:4) = matmul(arr1(2:4,:),arr2)
  call assign_result(5,8,arr3,results)
  !print *,arr3

  ! tests 9-12
  arr3=0
  arr3(1:3) = matmul(arr1(2:4,:),arr2)
  call assign_result(9,12,arr3,results)
  !print *,arr3

  !tests 13-16
  arr3=0
  arr3(2:4) = matmul(arr1(1:3,:),arr2)
  call assign_result(13,16,arr3,results)
  !print *,arr3

  !tests 17-20
  arr3=0
  arr3(2:4) = matmul(arr1(2:4,1:2),arr2(1:2))
  call assign_result(17,20,arr3,results)
  !print *,arr3

  !tests 21-24
  arr3=0
  arr3 = matmul(arr1(:,2:3),arr2(1:2))
  call assign_result(21,24,arr3,results)
  !print *,arr3

  !tests 25-28
  arr3=0
  arr3 = matmul(arr1(:,1:2),arr2(2:3))
  call assign_result(25,28,arr3,results)
  !print *,arr3

  !tests 29-32
  arr3=0
  arr3(2:4)  = matmul(arr1(1:3,1:2),arr2(2:3))
  call assign_result(29,32,arr3,results)
  !print *,arr3

  !tests 33-36
  arr3=0
  arr3(1:3)  = matmul(arr1(2:4,2:3),arr2(1:2))
  call assign_result(33,36,arr3,results)
  !print *,arr3

  !tests 37-40
  arr3=0
  arr3(1:3) = matmul(arr1(2:4,1:3:2),arr2(1:3:2))
  call assign_result(37,40,arr3,results)
  !print *,arr3

  !tests 41-44
  arr3=0
  arr3(2:4:2)  = matmul(arr1(1:3:2,1:2),arr2(2:3))
  call assign_result(41,44,arr3,results)
  !print *,arr3

  !tests 45-48
  arr3=0
  arr3(1:3:2)  = matmul(arr1(2:4:2,2:3),arr2(1:2))
  call assign_result(45,48,arr3,results)
  !print *,arr3

  !tests 49-64
  arr4=0
  arr4(2,1:3:2)  = matmul(arr1(2:4:2,2:3),arr2(1:2))
  call assign_result(49,64,arr4,results)
  !print *,arr4

  !tests 65-80
  arr4=0
  arr4(1:3:2,3)  = matmul(arr1(2:4:2,2:3),arr2(1:2))
  call assign_result(65,80,arr4,results)
  !print *,arr4

  !tests 81-96
  arr7=0
  arr7(0,0:2:2)  = matmul(arr5(1:3:2,0:1),arr6(-3:-2))
  call assign_result(81,96,arr7,results)
  !print *,arr7

  !tests 97-112
  arr7=0
  arr7(-1:1:2,2)  = matmul(arr5(1:3:2,0:1),arr6(-2:-1))
  call assign_result(97,112,arr7,results)
  !print *,arr7

  !tests 113-116
  arr3=0
  arr3(3:1:-1) = matmul(arr1(2:4,3:1:-2),arr2(3:1:-2))
  call assign_result(113,116,arr3,results)
  !print *,arr3

  !tests 117-120
  arr3=0
  arr3(3:1:-2)  = matmul(arr1(4:2:-2,2:3),arr2(1:2))
  call assign_result(117,120,arr3,results)
  !print *,arr3

  !tests 121,136
  arr4=0
  arr4(2,3:1:-2)  = matmul(arr1(2:4:2,2:3),arr2(1:2))
  call assign_result(121,136,arr4,results)
  !print *,arr4

  !tests 137-152
  arr4=0
  arr4(3:1:-2,3)  = matmul(arr1(4:2:-2,2:3),arr2(1:2))
  call assign_result(137,152,arr4,results)
  !print *,arr4

  !tests 153-168
  arr7=0
  arr7(1:-1:-2,2)  = matmul(arr5(3:1:-2,0:1),arr6(-2:-1))
  call assign_result(153,168,arr7,results)
  !print *,arr7

  !tests 169-172
  arr3=0
  arr3(1:3) = matmul(arr1(2:4,3:1:-2),arr2(3:1:-2))
  call assign_result(169,172,arr3,results)
  !print *,arr3

  arr12 = 0

  !print *,"tests 173-176"
  arr10=0
  arr10 = arr12 + matmul(arr8,arr9)
  call assign_result(173,176,arr10,results)
  !print *,arr10

  !print *,"tests 177-180"
  arr10=0
  arr10 = arr12 +  matmul(arr11,arr9)
  call assign_result(177,180,arr10,results)
  !print *,arr10

  !print *,"tests 181-184"
  arr10=0
  arr10(2:4) = arr12(2:4) + matmul(arr8(3:5,:),arr9)
  call assign_result(181,184,arr10,results)
  !print *,arr10

  !print *,"tests 185-188"
  arr10=0
  arr10(3:5) = arr12(3:5) + matmul(arr8(2:4,:),arr9)
  call assign_result(185,188,arr10,results)
  !print *,arr10

  !print *,"tests 189-192"
  arr10=0
  arr10(3:5) = arr12(3:5) + matmul(arr8(3:5,1:2),arr9(2:3))
  call assign_result(189,192,arr10,results)
  !print *,arr10

  !print *,"tests 193-196"
  arr10=0
  arr10 = arr12 + matmul(arr8(:,2:3),arr9(2:3))
  call assign_result(193,196,arr10,results)
  !print *,arr10

  !print *,"tests 197-200"
  arr10=0
  arr10 = arr12 + matmul(arr8(:,1:2),arr9(3:4))
  call assign_result(197,200,arr10,results)
  !print *,arr10

  !print *,"tests 201-204"
  arr10=0
  arr10(3:5)  = arr12(3:5) + matmul(arr8(2:4,1:2),arr9(3:4))
  call assign_result(201,204,arr10,results)
  !print *,arr10

  !print *,"tests 205-208"
  arr10=0
  arr10(2:4)  = arr12(2:4) + matmul(arr8(3:5,2:3),arr9(2:3))
  call assign_result(205,208,arr10,results)
  !print *,arr10

  call check(results, expect, NbrTests)

end program

subroutine assign_result(s_idx, e_idx , arr, rslt)
  integer*4, dimension(1:e_idx-s_idx+1) :: arr
  integer*4, dimension(e_idx) :: rslt
  integer:: s_idx, e_idx

  rslt(s_idx:e_idx) = arr

end subroutine

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
      program t
      i=1
      do while (i.le.3)
         print *,i
         i=i+1
      end do
      end program t

ext: f

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   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.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/master/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!#define CAT(a,b) a##b
!#define XCAT(a,b) CAT(a,b)
#if defined PARADRAM

#define SAMPLER_MOD ParaDRAM_mod
#define SAMPLER_TYPE ParaDRAM_type
#define SAMPLER_PROPOSAL_ABSTRACT_MOD ParaDRAMProposalAbstract_mod

#elif defined PARADISE

#define SAMPLER_MOD ParaDISE_mod
#define SAMPLER_TYPE ParaDISE_type
#define SAMPLER_PROPOSAL_ABSTRACT_MOD ParaDISEProposalAbstract_mod

#else
#error "Unrecognized sampler in ParaDRAM_mod.inc.f90"
#endif

    use Constants_mod, only: IK, RK, CK
    use ParaMonte_mod, only: ParaMonteNumFunCall_type
    use ParaMCMC_mod, only: ParaMCMC_type, ParaMCMC_Statistics_type, ParaMCMC_Chain_type
    use SpecDRAM_mod, only: SpecDRAM_type
    use ParaMonteLogFunc_mod, only: getLogFunc_proc
    use SAMPLER_PROPOSAL_ABSTRACT_MOD, only: ProposalAbstract_type

    implicit none

#if defined PARADRAM
    character(*), parameter :: MODULE_NAME = "@ParaDRAM_mod"
#elif defined PARADISE
    character(*), parameter :: MODULE_NAME = "@ParaDISE_mod"
#endif

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    type, extends(ParaMonteNumFunCall_type) :: NumFunCall_type
        integer(IK)                         :: acceptedRejectedDelayed          ! accepted + rejected function calls, including the delayed rejections
        integer(IK)                         :: acceptedRejectedDelayedUnused    ! by all processes, used or unused
    end type NumFunCall_type

    type, extends(ParaMCMC_Statistics_type) :: Statistics_type
        type(ParaMCMC_Chain_type)           :: AdaptationBurninLoc              ! burning loc based on the minimum adaptation measure requested
        type(NumFunCall_type)               :: NumFunCall
    end type Statistics_type

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    type, extends(ParaMCMC_type)                    :: SAMPLER_TYPE
        type(SpecDRAM_type)                         :: SpecDRAM
        type(Statistics_type)                       :: Stats
       !type(RefinedChain_type)                     :: RefinedChain
       !class(ProposalAbstract_type), pointer       :: Proposal => null()
        class(ProposalAbstract_type), allocatable   :: Proposal
    contains
        procedure, pass, private                    :: getSpecFromInputFile
        procedure, pass, public                     :: runSampler
        procedure, pass, private                    :: runKernel
    end type SAMPLER_TYPE ! ParaDRAM_type or ParaDISE_type

    !interface ParaDRAM_type
    !    module procedure :: runParaDRAM
    !end interface ParaDRAM_type

    interface
    module subroutine getSpecFromInputFile( self, nd )
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: getSpecFromInputFile
#endif
        use Constants_mod, only: IK
        class(SAMPLER_TYPE), intent(inout)  :: self
        integer(IK), intent(in)             :: nd
    end subroutine getSpecFromInputFile
    end interface

    interface
    module subroutine runKernel( self, getLogFunc )
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: runKernel
#endif
        use ParaMonteLogFunc_mod, only: getLogFunc_proc
        class(SAMPLER_TYPE), intent(inout)  :: self
        procedure(getLogFunc_proc)          :: getLogFunc
    end subroutine runKernel
    end interface

    interface
    module subroutine runSampler( self                                  &
                                , ndim                                  &
                                , getLogFunc                            &
                                , inputFile                             &
                                ! ParaMonte variables
                                , sampleSize                            &
                                , randomSeed                            &
                                , description                           &
                                , outputFileName                        &
                                , outputDelimiter                       &
                                , chainFileFormat                       &
                                , variableNameList                      &
                                , restartFileFormat                     &
                                , outputColumnWidth                     &
                                , outputRealPrecision                   &
                                , silentModeRequested                   &
                                , domainLowerLimitVec                   &
                                , domainUpperLimitVec                   &
                                , parallelizationModel                  &
                                , progressReportPeriod                  &
                                , TargetAcceptanceRate                  &
                                , mpiFinalizeRequested                  &
                                , maxNumDomainCheckToWarn               &
                                , maxNumDomainCheckToStop               &
                                ! ParaMCMC variables
                                , chainSize                             &
                                , startPointVec                         &
                                , sampleRefinementCount                 &
                                , sampleRefinementMethod                &
                                , randomStartPointRequested             &
                                , randomStartPointDomainLowerLimitVec   &
                                , randomStartPointDomainUpperLimitVec   &
                                ! ParaDRAM variables
                                , scaleFactor                           &
                                , proposalModel                         &
                                , proposalStartStdVec                   &
                                , proposalStartCorMat                   &
                                , proposalStartCovMat                   &
                                , adaptiveUpdateCount                   &
                                , adaptiveUpdatePeriod                  &
                                , greedyAdaptationCount                 &
                                , delayedRejectionCount                 &
                                , burninAdaptationMeasure               &
                                , delayedRejectionScaleFactorVec        &
                                ) ! result(self)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !!DEC$ ATTRIBUTES DLLEXPORT :: SAMPLER_TYPE
        !DEC$ ATTRIBUTES DLLEXPORT :: runSampler
#endif
        use ParaMonteLogFunc_mod, only: getLogFunc_proc
        use Constants_mod, only: IK, RK

        implicit none

        ! self
        class(SAMPLER_TYPE), intent(inout)  :: self

        ! mandatory variables
        integer(IK) , intent(in)            :: ndim
        procedure(getLogFunc_proc)          :: getLogFunc

        character(*), intent(in), optional  :: inputFile

        ! ParaMonte variables
        integer(IK) , intent(in), optional  :: sampleSize
        integer(IK) , intent(in), optional  :: randomSeed
        character(*), intent(in), optional  :: description
        character(*), intent(in), optional  :: outputFileName
        character(*), intent(in), optional  :: outputDelimiter
        character(*), intent(in), optional  :: chainFileFormat
        character(*), intent(in), optional  :: variableNameList(ndim)
        character(*), intent(in), optional  :: restartFileFormat
        integer(IK) , intent(in), optional  :: outputColumnWidth
        integer(IK) , intent(in), optional  :: outputRealPrecision
        logical     , intent(in), optional  :: silentModeRequested
        real(RK)    , intent(in), optional  :: domainLowerLimitVec(ndim)
        real(RK)    , intent(in), optional  :: domainUpperLimitVec(ndim)
        character(*), intent(in), optional  :: parallelizationModel
        integer(IK) , intent(in), optional  :: progressReportPeriod
        real(RK)    , intent(in), optional  :: TargetAcceptanceRate(2)
        logical     , intent(in), optional  :: mpiFinalizeRequested
        integer(IK) , intent(in), optional  :: maxNumDomainCheckToWarn
        integer(IK) , intent(in), optional  :: maxNumDomainCheckToStop

        ! ParaMCMC variables
        integer(IK) , intent(in), optional  :: chainSize
        real(RK)    , intent(in), optional  :: startPointVec(ndim)
        integer(IK) , intent(in), optional  :: sampleRefinementCount
        character(*), intent(in), optional  :: sampleRefinementMethod
        logical     , intent(in), optional  :: randomStartPointRequested
        real(RK)    , intent(in), optional  :: randomStartPointDomainLowerLimitVec(ndim)
        real(RK)    , intent(in), optional  :: randomStartPointDomainUpperLimitVec(ndim)

        ! ParaDRAM variables
        character(*), intent(in), optional  :: scaleFactor
        character(*), intent(in), optional  :: proposalModel
        real(RK)    , intent(in), optional  :: proposalStartStdVec(ndim)
        real(RK)    , intent(in), optional  :: proposalStartCorMat(ndim,ndim)
        real(RK)    , intent(in), optional  :: proposalStartCovMat(ndim,ndim)
        integer(IK) , intent(in), optional  :: adaptiveUpdateCount
        integer(IK) , intent(in), optional  :: adaptiveUpdatePeriod
        integer(IK) , intent(in), optional  :: greedyAdaptationCount
        integer(IK) , intent(in), optional  :: delayedRejectionCount
        real(RK)    , intent(in), optional  :: burninAdaptationMeasure
        real(RK)    , intent(in), optional  :: delayedRejectionScaleFactorVec(:)

    end subroutine runSampler
    end interface

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#if defined PARADISE
#elif defined PARADRAM
#endif

#undef SAMPLER_MOD
#undef SAMPLER_TYPE
#undef SAMPLER_PROPOSAL_ABSTRACT_MOD


!##############################################################################
! PROGRAM portfoliochoice_MCM
!
! ## Portfolio choice with Monte Carlo minimization (MCM)
!
! This code is published under the GNU General Public License v3
!                         (https://www.gnu.org/licenses/gpl-3.0.en.html)
!
! Authors: Hans Fehr and Fabian Kindermann
!          contact@ce-fortran.com
!
! #VC# VERSION: 1.0  (23 January 2018)
!
!##############################################################################
program portfoliochoice_MCM

    use toolbox

    implicit none
    integer, parameter :: TT = 5        ! number of time periods
    integer, parameter :: NN = 3        ! number of stocks
    integer, parameter :: NP = 100      ! number of mc portfolios
    real*8, parameter :: gamma = 10d0   ! risk aversion
    real*8, parameter :: r_f = 0.05     ! the risk free rate

    real*8 :: omega(0:NN)
    real*8 :: mu(0:NN), sig(0:NN,0:NN)
    real*8 :: mu_p, sig_p

    real*8 :: omega_mc(0:NN,NP), mu_mc(NP), sig_mc(NP)
    integer :: ii, ik, maxi

    ! get returns and their statistical properties
    call returns

    ! simulate initial portfolios
    do ii = 1, NP
        call simulate_uniform(omega_mc(:, ii), 0d0, 1d0)
        omega_mc(:, ii) = omega_mc(:, ii)/sum(omega_mc(:, ii))
    enddo

    ! perform updating steps
    do ik = 1, 1000

        ! calculate means and returns of the portfolios
        do ii = 1, NP
            mu_mc(ii)  = dot_product(mu, omega_mc(:, ii))
            sig_mc(ii) = dot_product(omega_mc(:, ii), matmul(sig, omega_mc(:, ii)))
        enddo

        ! get portfolio that maximizes investors utility
        maxi = maxloc(mu_mc - gamma/2d0*sig_mc, 1)

        ! store maximum
        omega = omega_mc(:, maxi)

        ! simulate an updated set of portfolios around maximum
        do ii = 1, NP-1

            ! get new normalized portfolio
            call simulate_uniform(omega_mc(:, ii), 0d0, 1d0)
            omega_mc(:, ii) = omega_mc(:, ii)/sum(omega_mc(:, ii))

            ! create linear combination between old and new
            omega_mc(:, ii) = 0.9d0*omega + 0.1d0*omega_mc(:, ii)
        enddo

        ! store maximum in last entry
        omega_mc(:, NP) = omega
    enddo

    ! calculate return and variance
    mu_p   = dot_product(mu, omega)
    sig_p  = dot_product(omega, matmul(sig, omega))

    ! print results
    write(*,'(a)') 'Portfolio shares:   Riskless   Stock 1   Stock 2   Stock 3      mu_p     sig_p'
    write(*,'(a,6f10.4)')'Optimal Portfolio ', omega(0), omega(1), omega(2), omega(3), mu_p, sig_p

contains

    ! get returns and statictical properties
    subroutine returns

        implicit none
        integer :: ii, ij
        real*8 :: odat(0:TT,NN), r(TT,NN)

        ! initialize price data
        odat(0,:) = (/ 1.00d0, 2.00d0, 3.00d0 /)
        odat(1,:) = (/ 1.02d0, 2.65d0, 2.80d0 /)
        odat(2,:) = (/ 1.17d0, 2.40d0, 4.50d0 /)
        odat(3,:) = (/ 1.08d0, 2.70d0, 4.20d0 /)
        odat(4,:) = (/ 1.16d0, 2.85d0, 3.20d0 /)
        odat(5,:) = (/ 1.26d0, 2.75d0, 4.20d0 /)

        ! compute return per period
        do ii = 1, NN
            do ij = 1, TT
                r(ij, ii) = odat(ij, ii)/odat(ij-1, ii)-1d0
            enddo
        enddo

        ! compute expected value
        do ii = 1, NN
            mu(ii) = sum(r(:, ii))/dble(TT)
        enddo

        ! compute covariance matrix
        do ii = 1, NN
            do ij = 1, NN
                sig(ii, ij) = dot_product(r(:, ii)-mu(ii), r(:, ij)-mu(ij))
            enddo
        enddo
        sig = sig/dble(TT)

        ! store risk free rate and variance
        mu(0) = r_f
        sig(0, :) = 0d0
        sig(:, 0) = 0d0

    end subroutine

end program


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

      subroutine ssor(niter)

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

c---------------------------------------------------------------------
c   to perform pseudo-time stepping SSOR iterations
c   for five nonlinear pde's.
c---------------------------------------------------------------------

      implicit none
      integer niter

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m, n, lst, lend
      integer istep
      double precision  tmp, tv(5*isiz1*isiz2)
      double precision  delunm(5)

      external timer_read
      double precision timer_read

c---------------------------------------------------------------------
c   begin pseudo-time stepping iterations
c---------------------------------------------------------------------
      tmp = 1.0d+00 / ( omega * ( 2.0d+00 - omega ) )

c---------------------------------------------------------------------
c   initialize a,b,c,d to zero (guarantees that page tables have been
c   formed, if applicable on given architecture, before timestepping).
c---------------------------------------------------------------------
      do j=1,isiz2
         do i=1,isiz1
            do n=1,5
               do m=1,5
                  a(m,n,i,j) = 0.d0
                  b(m,n,i,j) = 0.d0
                  c(m,n,i,j) = 0.d0
                  d(m,n,i,j) = 0.d0
               enddo
            enddo
         enddo
      enddo
      do i = 1, t_last
         call timer_clear(i)
      end do

c---------------------------------------------------------------------
c   compute the steady-state residuals
c---------------------------------------------------------------------
      call rhs

c---------------------------------------------------------------------
c   compute the L2 norms of newton iteration residuals
c---------------------------------------------------------------------
      call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
     >             ist, iend, jst, jend,
     >             rsd, rsdnm )

c      if ( ipr .eq. 1 ) then
c         write (*,*) '          Initial residual norms'
c         write (*,*)
c         write (*,1007) ( rsdnm(m), m = 1, 5 )
c         write (*,'(/a)') 'Iteration RMS-residual of 5th PDE'
c      end if

      do i = 1, t_last
         call timer_clear(i)
      end do
      call timer_start(1)

c---------------------------------------------------------------------
c   the timestep loop
c---------------------------------------------------------------------
      do istep = 1, niter

c         if ( ( mod ( istep, inorm ) .eq. 0 ) .and.
c     >          ipr .eq. 1 ) then
c             write ( *, 1001 ) istep
c         end if
         if (mod ( istep, 20) .eq. 0 .or.
     >         istep .eq. itmax .or.
     >         istep .eq. 1) then
            if (niter .gt. 1) write( *, 200) istep
 200        format(' Time step ', i4)
         endif

c---------------------------------------------------------------------
c   perform SSOR iteration
c---------------------------------------------------------------------
         if (timeron) call timer_start(t_rhs)
         do k = 2, nz - 1
            do j = jst, jend
               do i = ist, iend
                  do m = 1, 5
                     rsd(m,i,j,k) = dt * rsd(m,i,j,k)
                  end do
               end do
            end do
         end do
         if (timeron) call timer_stop(t_rhs)

         lst = ist + jst
         lend = iend + jend

         do k = 2, nz -1
c---------------------------------------------------------------------
c   form the lower triangular part of the jacobian matrix
c---------------------------------------------------------------------
            if (timeron) call timer_start(t_jacld)
            call jacld(k)
            if (timeron) call timer_stop(t_jacld)

c---------------------------------------------------------------------
c   perform the lower triangular solution
c---------------------------------------------------------------------
            if (timeron) call timer_start(t_blts)
            call blts( isiz1, isiz2, isiz3,
     >                 nx, ny, nz, k,
     >                 omega,
     >                 rsd,
     >                 a, b, c, d,
     >                 ist, iend, jst, jend,
     >                 lst, lend )
            if (timeron) call timer_stop(t_blts)
         end do

         do k = nz - 1, 2, -1
c---------------------------------------------------------------------
c   form the strictly upper triangular part of the jacobian matrix
c---------------------------------------------------------------------
            if (timeron) call timer_start(t_jacu)
            call jacu(k)
            if (timeron) call timer_stop(t_jacu)

c---------------------------------------------------------------------
c   perform the upper triangular solution
c---------------------------------------------------------------------
            if (timeron) call timer_start(t_buts)
            call buts( isiz1, isiz2, isiz3,
     >                 nx, ny, nz, k,
     >                 omega,
     >                 rsd, tv,
     >                 d, a, b, c,
     >                 ist, iend, jst, jend,
     >                 lst, lend )
            if (timeron) call timer_stop(t_buts)
         end do

c---------------------------------------------------------------------
c   update the variables
c---------------------------------------------------------------------

         if (timeron) call timer_start(t_add)
         do k = 2, nz-1
            do j = jst, jend
               do i = ist, iend
                  do m = 1, 5
                     u( m, i, j, k ) = u( m, i, j, k )
     >                    + tmp * rsd( m, i, j, k )
                  end do
               end do
            end do
         end do
         if (timeron) call timer_stop(t_add)

c---------------------------------------------------------------------
c   compute the max-norms of newton iteration corrections
c---------------------------------------------------------------------
         if ( mod ( istep, inorm ) .eq. 0 ) then
            if (timeron) call timer_start(t_l2norm)
            call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
     >                   ist, iend, jst, jend,
     >                   rsd, delunm )
            if (timeron) call timer_stop(t_l2norm)
c            if ( ipr .eq. 1 ) then
c                write (*,1006) ( delunm(m), m = 1, 5 )
c            else if ( ipr .eq. 2 ) then
c                write (*,'(i5,f15.6)') istep,delunm(5)
c            end if
         end if

c---------------------------------------------------------------------
c   compute the steady-state residuals
c---------------------------------------------------------------------
         if (timeron) call timer_start(t_rhs)
         call rhs
         if (timeron) call timer_stop(t_rhs)

c---------------------------------------------------------------------
c   compute the max-norms of newton iteration residuals
c---------------------------------------------------------------------
         if ( ( mod ( istep, inorm ) .eq. 0 ) .or.
     >        ( istep .eq. itmax ) ) then
            if (timeron) call timer_start(t_l2norm)
            call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
     >                   ist, iend, jst, jend,
     >                   rsd, rsdnm )
            if (timeron) call timer_stop(t_l2norm)
c            if ( ipr .eq. 1 ) then
c                write (*,1007) ( rsdnm(m), m = 1, 5 )
c            end if
         end if

c---------------------------------------------------------------------
c   check the newton-iteration residuals against the tolerance levels
c---------------------------------------------------------------------
         if ( ( rsdnm(1) .lt. tolrsd(1) ) .and.
     >        ( rsdnm(2) .lt. tolrsd(2) ) .and.
     >        ( rsdnm(3) .lt. tolrsd(3) ) .and.
     >        ( rsdnm(4) .lt. tolrsd(4) ) .and.
     >        ( rsdnm(5) .lt. tolrsd(5) ) ) then
c            if (ipr .eq. 1 ) then
c               write (*,1004) istep
c            end if
            return
         end if

      end do

      call timer_stop(1)
      maxtime= timer_read(1)

      return

 1001 format (1x/5x,'pseudo-time SSOR iteration no.=',i4/)
 1004 format (1x/1x,'convergence was achieved after ',i4,
     >   ' pseudo-time steps' )
 1006 format (1x/1x,'RMS-norm of SSOR-iteration correction ',
     > 'for first pde  = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for second pde = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for third pde  = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for fourth pde = ',1pe12.5/,
     > 1x,'RMS-norm of SSOR-iteration correction ',
     > 'for fifth pde  = ',1pe12.5)
 1007 format (1x/1x,'RMS-norm of steady-state residual for ',
     > 'first pde  = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'second pde = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'third pde  = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'fourth pde = ',1pe12.5/,
     > 1x,'RMS-norm of steady-state residual for ',
     > 'fifth pde  = ',1pe12.5)

      end

max_stars_repo_name: rahlk/Parallel-Programming

max_issues_repo_name: rahlk/Parallel-Programming

max_forks_repo_name: rahlk/Parallel-Programming


text:
C
C
*+ plot_cursor_get

       subroutine plot_cursor_get(uv,s)
C      --------------------------------
C
C Read cursor position and return it
C
C Returned:
C  cursor position -- map coordinates
       real*4    uv(2)
C
C Updated:
C  error status
       integer   s
C
C A position is read from the current graphics map display and returned
C in uv.  The standard parameters %map-u, %map-v %map-pos, %u, %v,
C %uv-pos and %gr-char are also set.
C
*-
       include '../include/plt_basic_defn.inc'
       include '../include/plt_image_defn.inc'
       include '/mrao/include/constants.inc'
       include '/mrao/include/chrlib_functions.inc'

C graphic position
       real*4       gx, gy
       character    ch_typed*1, rnumber*20, value
       integer      i1, i2, len_v

C cursor position in map coordinates
       real*4       pos_u, pos_v

C functions and flags
       integer      pgcurse
       logical      cursor_available, cmd_dblev

       if (s.ne.0) return

       if (.not.plot_open) then
         call cmd_wrerr('CURSOR-GET',
     *                  'Plot not open - unable to access cursor')
         return
       end if
       if (.not.frame_init) then
         call cmd_wrerr('CURSOR-GET',
     *                  'Frame not initialised for image')
         return
       end if
       call pgqcur(cursor_available)
       if (.not.cursor_available) then
         call cmd_wrerr('CURSOR-GET',
     *                   'Device does not support input')
         return
       end if

C find initial position for cursor
       if (cmd_dblev(1)) then
         call io_wrout('.. move cursor: type any character to end')
       endif

       pos_u = (frame_win(1)+frame_win(2))/2.0
       pos_v = (frame_win(3)+frame_win(4))/2.0

C read cursor position
       s = pgcurse(pos_u,pos_v,ch_typed) - 1
       if (s.ne.0) then
         call cmd_wrerr('CURSOR-GET','Error reading CURSOR')
         return
       end if

C set value for graphic position and character
       gx = (pos_u-frame_win(1))*(frame_vp(2)-frame_vp(1))/
     *      (frame_win(2)-frame_win(1)) + frame_vp(1)
       gy = (pos_v-frame_win(3))*(frame_vp(4)-frame_vp(3))/
     *      (frame_win(4)-frame_win(3)) + frame_vp(3)
       call cmd_setlocal( '%gr-char', ch_typed, s )
       call chr_chrtoc(gx,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal( '%gr-x', rnumber(i1:i2), s )
       call chr_chrtoc(gy,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal( '%gr-y', rnumber(i1:i2), s )

       call chr_chitoc(nint(pos_u),rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%u',rnumber(i1:i2),s)
       value = rnumber(i1:i2)
       call chr_chitoc(nint(pos_v),rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%v',rnumber(i1:i2),s)
       len_v = chr_lenb(value)
       value = value(1:len_v)//','//rnumber(i1:i2)
       len_v = chr_lenb(value)
       call cmd_setlocal( '%uv-pos', value(1:len_v), s )

       call chr_chrtoc(pos_u,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%map-u',rnumber(i1:i2),s)
       value = rnumber(i1:i2)
       call chr_chrtoc(pos_v,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%map-v',rnumber(i1:i2),s)
       len_v = chr_lenb(value)
       value = value(1:len_v)//','//rnumber(i1:i2)
       len_v = chr_lenb(value)
       call cmd_setlocal( '%map-pos', value(1:len_v), s )

       if (s.eq.0) then
         uv(1) = pos_u
         uv(2) = pos_v
       endif
       call cmd_err( s,'PLOT_CURSOR_GET',' ')

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
constant a = 1
constant b = 0
constant c = a/b

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
module viscosity_module

  use amrex_fort_module, only : rt => amrex_real
  use network
  use eos_type_module
  use eos_module

  implicit none

contains

  subroutine viscous_coeff(eos_state, visc_coeff)

    use extern_probin_module, only: const_viscosity

    type (eos_t), intent(in) :: eos_state
    real (rt), intent(inout) :: visc_coeff

    visc_coeff = const_viscosity

  end subroutine viscous_coeff

end module viscosity_module

max_stars_repo_name: doreenfan/Microphysics

max_issues_repo_name: Youhichka/Microphysics

max_forks_repo_name: Youhichka/Microphysics


text:
C$Procedure      PODEGC ( Pod, end group, character )

      SUBROUTINE PODEGC ( POD )
      IMPLICIT NONE

C$ Abstract
C
C     End the active group of a pod, restoring the previous group
C     unchanged.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     PODS
C
C$ Keywords
C
C     ARRAY
C     CELLS
C     PODS
C
C$ Declarations

      INTEGER               LBCELL
      PARAMETER           ( LBCELL = -5 )

      CHARACTER*(*)         POD ( LBCELL : * )

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     POD       I,O  Pod.
C
C$ Detailed_Input
C
C     POD       on input, is an arbitrary pod.
C
C$ Detailed_Output
C
C     POD       on output, is the same pod after the active group
C               has been removed and the previous group has been
C               restored unchanged.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If the active group is the only group in the pod, the
C        cardinality of the POD is set to zero.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     The active group of a pod may be removed by any of the
C     following routines: PODEG (end group), PODCG (close group),
C     or PODRG (replace group).
C
C     PODEG effectively returns the pod to its state before the
C     active group was created. The contents of the active group
C     are simply lost.
C
C     PODCG appends the contents of the active group to the previous
C     group to obtain the new active group, reducing the number of
C     groups in the pod by one.
C
C     PODRG also reduces the number of groups, but by replacing the
C     previous group with the active group, as though the previous
C     group had never existed.
C
C$ Examples
C
C     Let NAMES be a character POD containing the following groups:
C
C        Group 1:  NEWTON
C                  GALILEO
C                  KEPLER
C
C        Group 2:  EINSTEIN
C                  BOHR
C                  HEISENBERG
C
C        Group 3:  FEYNMAN
C                  BARDEEN
C
C     Following the call
C
C        CALL PODEGC ( NAMES )
C
C     the active group (Group 2) contains EINSTEIN, BOHR, and
C     HEISENBERG. Following the call
C
C        CALL PODCGC ( NAMES )
C
C     the active group (again, Group 2) contains EINSTEIN, BOHR,
C     HEISENBERG, FEYNMAN, and BARDEEN. Following the call
C
C        CALL PODRGC ( NAMES )
C
C     the active group (also Group 2) contains FEYNMAN and BARDEEN.
C
C$ Restrictions
C
C     1) In any pod, only the active group should be accessed,
C        and its location should always be determined by PODBE
C        or PODON. Never assume that the active group begins
C        at POD(1).
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber     (JPL)
C     I.M. Underwood (JPL)
C
C$ Version
C
C-    Beta Version 1.0.0, 15-JUL-1989 (WLT) (IMU)
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN

C
C     Local parameters
C
      INTEGER               GRPOFF
      PARAMETER           ( GRPOFF = -2 )

C
C     Local variables
C
      INTEGER               NUMBER
      INTEGER               OFFSET

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'PODEGC' )
      END IF

C
C     At any given time, the offset of the active group is stored
C     in location GRPOFF of the control area, so POD(GRPOFF) tells
C     us the location of the element preceding the active group.
C
C     This element is a backward pointer, containing the offset of
C     the previous group; and so on, with turtles all the way down.
C     For example, consider a pod with three groups
C
C         G.  <10>
C         1.  Bob
C         2.  Carol
C         3.  Ted
C         4.  Alice
C         5.  <0>
C         6.  Fred
C         7.  Wilma
C         8.  Barney
C         9.  Bettey
C        10.  <5>
C        11.  Ricky
C        12.  Lucy
C        13.  Fred
C        14.  Ethel
C
C     When the second group was created, the offset of the first
C     group (zero) was appended to the pod; the location of this
C     offset became the offset for the second group. When the
C     third group was created, the offset of the second group (5)
C     was appended; the location of this offset became the offset for
C     the third group. The offset for the third group is located
C     in element GRPOFF.
C
C     To remove a group then, all that is necessary is to look at
C     element GRPOFF to get the offset of the current group; go to
C     that location to get the offset of the previous group; and
C     move that offset into element GRPOFF. The new cardinality,
C     of course, should be one less than the original offset (9).
C
C     If the pod contains only one group, it can't be removed, but
C     it can be emptied by setting the cardinality of the pod to
C     zero.
C
C
      CALL PODONC ( POD, OFFSET, NUMBER )

      IF ( OFFSET .EQ. 0 ) THEN
         CALL SCARDC ( 0, POD )

      ELSE
         POD(GRPOFF) = POD(OFFSET)
         CALL SCARDC ( OFFSET - 1, POD )
      END IF

      CALL CHKOUT ( 'PODEGC' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Copyright (c) 1997 Neil N. Carlson
!!
!! This file is part of MFE1 which is released under the MIT license.  See the
!! file LICENSE or visit http://opensource.org/licenses/MIT for details.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module mfe_constants

  use kind_parameters, only: r8
  implicit none
  private

  integer, parameter, public :: wp = r8                 ! Use 8-byte reals

  integer, parameter, public :: NEQNS = 2               ! Number equations in the PDE system

  !!! DO NOT MODIFY THESE PARAMETERS !!!
  integer, parameter, public :: DIMEN = 1               ! Spatial dimension
  integer, parameter, public :: NVARS = NEQNS + DIMEN   ! Total number of variables per node
  integer, parameter, public :: NVERT = DIMEN + 1       ! Number of vertices per cell
  integer, parameter, public :: NFACE = NVERT           ! Number of faces per cell

end module mfe_constants


text:
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                          .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... file shsgc.f
c
c     this file contains code and documentation for subroutines
c     shsgc and shsgci
c
c ... files which must be loaded with shsgc.f
c
c     sphcom.f, hrfft.f, gaqd.f
c
c     subroutine shsgc(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab,
c    +                 wshsgc,lshsgc,work,lwork,ierror)
c
c     subroutine shsgc performs the spherical harmonic synthesis
c     on the arrays a and b and stores the result in the array g.
c     the synthesis is performed on an equally spaced longitude grid
c     and a gaussian colatitude grid.  the associated legendre functions
c     are recomputed rather than stored as they are in subroutine
c     shsgs.  the synthesis is described below at output parameter
c     g.
c
c     input parameters
c
c     nlat   the number of points in the gaussian colatitude grid on the
c            full sphere. these lie in the interval (0,pi) and are compu
c            in radians in theta(1),...,theta(nlat) by subroutine gaqd.
c            if nlat is odd the equator will be included as the grid poi
c            theta((nlat+1)/2).  if nlat is even the equator will be
c            excluded as a grid point and will lie half way between
c            theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3.
c            note: on the half sphere, the number of grid points in the
c            colatitudinal direction is nlat/2 if nlat is even or
c            (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than or equal to 4. the efficiency of the computation is
c            improved when nlon is a product of small prime numbers.
c
c     isym   = 0  no symmetries exist about the equator. the synthesis
c                 is performed on the entire sphere.  i.e. on the
c                 array g(i,j) for i=1,...,nlat and j=1,...,nlon.
c                 (see description of g below)
c
c            = 1  g is antisymmetric about the equator. the synthesis
c                 is performed on the northern hemisphere only.  i.e.
c                 if nlat is odd the synthesis is performed on the
c                 array g(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon.
c                 if nlat is even the synthesis is performed on the
c                 array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon.
c
c
c            = 2  g is symmetric about the equator. the synthesis is
c                 performed on the northern hemisphere only.  i.e.
c                 if nlat is odd the synthesis is performed on the
c                 array g(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon.
c                 if nlat is even the synthesis is performed on the
c                 array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon.
c
c     nt     the number of syntheses.  in the program that calls shsgc,
c            the arrays g,a and b can be three dimensional in which
c            case multiple synthesis will be performed.  the third
c            index is the synthesis index which assumes the values
c            k=1,...,nt.  for a single synthesis set nt=1. the
c            discription of the remaining parameters is simplified
c            by assuming that nt=1 or that the arrays g,a and b
c            have only two dimensions.
c
c     idg    the first dimension of the array g as it appears in the
c            program that calls shsgc. if isym equals zero then idg
c            must be at least nlat.  if isym is nonzero then idg must
c            be at least nlat/2 if nlat is even or at least (nlat+1)/2
c            if nlat is odd.
c
c     jdg    the second dimension of the array g as it appears in the
c            program that calls shsgc. jdg must be at least nlon.
c
c     mdab   the first dimension of the arrays a and b as it appears
c            in the program that calls shsgc. mdab must be at least
c            min0((nlon+2)/2,nlat) if nlon is even or at least
c            min0((nlon+1)/2,nlat) if nlon is odd
c
c     ndab   the second dimension of the arrays a and b as it appears
c            in the program that calls shsgc. ndab must be at least nlat
c
c     a,b    two or three dimensional arrays (see the input parameter
c            nt) that contain the coefficients in the spherical harmonic
c            expansion of g(i,j) given below at the definition of the
c            output parameter g.  a(m,n) and b(m,n) are defined for
c            indices m=1,...,mmax and n=m,...,nlat where mmax is the
c            maximum (plus one) longitudinal wave number given by
c            mmax = min0(nlat,(nlon+2)/2) if nlon is even or
c            mmax = min0(nlat,(nlon+1)/2) if nlon is odd.
c
c     wshsgc an array which must be initialized by subroutine shsgci.
c            once initialized, wshsgc can be used repeatedly by shsgc
c            as long as nlat and nlon remain unchanged.  wshsgc must
c            not be altered between calls of shsgc.
c
c     lshsgc the dimension of the array wshsgc as it appears in the
c            program that calls shsgc. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshsgc must be at least
c
c               nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/2+nlon+15
c
c     work   a work array that does not have to be saved.
c
c     lwork  the dimension of the array work as it appears in the
c            program that calls shsgc. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            if isym is zero then lwork must be at least
c
c                      nlat*(nlon*nt+max0(3*l2,nlon))
c
c            if isym is not zero then lwork must be at least
c
c                      l2*(nlon*nt+max0(3*nlat,nlon))
c
c     **************************************************************
c
c     output parameters
c
c     g      a two or three dimensional array (see input parameter nt)
c            that contains the discrete function which is synthesized.
c            g(i,j) contains the value of the function at the gaussian
c            colatitude point theta(i) and longitude point
c            phi(j) = (j-1)*2*pi/nlon. the index ranges are defined
c            above at the input parameter isym.  for isym=0, g(i,j)
c            is given by the the equations listed below.  symmetric
c            versions are used when isym is greater than zero.
c
c     the normalized associated legendre functions are given by
c
c     pbar(m,n,theta) = sqrt((2*n+1)*factorial(n-m)/(2*factorial(n+m)))
c                       *sin(theta)**m/(2**n*factorial(n)) times the
c                       (n+m)th derivative of (x**2-1)**n with respect
c                       to x=cos(theta)
c
c
c     define the maximum (plus one) longitudinal wave number
c     as   mmax = min0(nlat,(nlon+2)/2) if nlon is even or
c          mmax = min0(nlat,(nlon+1)/2) if nlon is odd.
c
c     then g(i,j) = the sum from n=0 to n=nlat-1 of
c
c                   .5*pbar(0,n,theta(i))*a(1,n+1)
c
c              plus the sum from m=1 to m=mmax-1 of
c
c                   the sum from n=m to n=nlat-1 of
c
c              pbar(m,n,theta(i))*(a(m+1,n+1)*cos(m*phi(j))
c                                    -b(m+1,n+1)*sin(m*phi(j)))
c
c     ierror = 0  no errors
c            = 1  error in the specification of nlat
c            = 2  error in the specification of nlon
c            = 3  error in the specification of isym
c            = 4  error in the specification of nt
c            = 5  error in the specification of idg
c            = 6  error in the specification of jdg
c            = 7  error in the specification of mdab
c            = 8  error in the specification of ndab
c            = 9  error in the specification of lwshig
c            = 10 error in the specification of lwork
c
c
c ****************************************************************
c
c     subroutine shsgci(nlat,nlon,wshsgc,lshsgc,dwork,ldwork,ierror)
c
c     subroutine shsgci initializes the array wshsgc which can then
c     be used repeatedly by subroutines shsgc. it precomputes
c     and stores in wshsgc quantities such as gaussian weights,
c     legendre polynomial coefficients, and fft trigonometric tables.
c
c     input parameters
c
c     nlat   the number of points in the gaussian colatitude grid on the
c            full sphere. these lie in the interval (0,pi) and are compu
c            in radians in theta(1),...,theta(nlat) by subroutine gaqd.
c            if nlat is odd the equator will be included as the grid poi
c            theta((nlat+1)/2).  if nlat is even the equator will be
c            excluded as a grid point and will lie half way between
c            theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3.
c            note: on the half sphere, the number of grid points in the
c            colatitudinal direction is nlat/2 if nlat is even or
c            (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than or equal to 4. the efficiency of the computation is
c            improved when nlon is a product of small prime numbers.
c
c     wshsgc an array which must be initialized by subroutine shsgci.
c            once initialized, wshsgc can be used repeatedly by shsgc
c            as long as nlat and nlon remain unchanged.  wshsgc must
c            not be altered between calls of shsgc.
c
c     lshsgc the dimension of the array wshsgc as it appears in the
c            program that calls shsgc. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshsgc must be at least
c
c                  nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/2+nlon+15
c
c     dwork  a double precision work array that does not have to be saved.
c
c     ldwork the dimension of the array dwork as it appears in the
c            program that calls shsgci. ldwork must be at least
c
c                 nlat*(nlat+4)
c
c     output parameter
c
c     wshsgc an array which must be initialized before calling shsgc.
c            once initialized, wshsgc can be used repeatedly by shsgc
c            as long as nlat and nlon remain unchanged.  wshsgc must not
c            altered between calls of shsgc.
c
c     ierror = 0  no errors
c            = 1  error in the specification of nlat
c            = 2  error in the specification of nlon
c            = 3  error in the specification of lshsgc
c            = 4  error in the specification of ldwork
c            = 5  failure in gaqd to compute gaussian points
c                 (due to failure in eigenvalue routine)
c
c
c ****************************************************************
      subroutine shsgc(nlat,nlon,mode,nt,g,idg,jdg,a,b,mdab,ndab,
     1                    wshsgc,lshsgc,work,lwork,ierror)
c     subroutine shsgc performs the spherical harmonic synthesis on
c     a gaussian grid using the coefficients in array(s) a,b and returns
c     the results in array(s) g.  the legendre polynomials are computed
c     as needed in this version.
c
      dimension g(idg,jdg,1),a(mdab,ndab,1),b(mdab,ndab,1),
     1          wshsgc(lshsgc),work(lwork)
c     check input parameters
      ierror = 1
      if (nlat.lt.3) return
      ierror = 2
      if (nlon.lt.4) return
      ierror = 3
      if (mode.lt.0 .or.mode.gt.2) return
      ierror = 4
      if (nt.lt.1) return
c     set limit for m iin a(m,n),b(m,n) computation
      l = min0((nlon+2)/2,nlat)
c     set gaussian point nearest equator pointer
      late = (nlat+mod(nlat,2))/2
c     set number of grid points for analysis/synthesis
      lat = nlat
      if (mode.ne.0) lat = late
      ierror = 5
      if (idg.lt.lat) return
      ierror = 6
      if (jdg.lt.nlon) return
      ierror = 7
      if(mdab .lt. l) return
      ierror = 8
      if(ndab .lt. nlat) return
      l1 = l
      l2 = late
      ierror = 9
c     check permanent work space length
      if (lshsgc .lt. nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/2+nlon+15)return
      ierror = 10
c     check temporary work space length
      if (mode.eq.0) then
      if(lwork.lt.nlat*(nlon*nt+max0(3*l2,nlon)))return
      else
c     mode.ne.0
      if(lwork.lt.l2*(nlon*nt+max0(3*nlat,nlon))) return
      end if
      ierror = 0
c     starting address  fft values
      ifft = nlat+2*nlat*late+3*(l*(l-1)/2+(nlat-l)*(l-1))+1
c     set pointers for internal storage of g and legendre polys
      ipmn = lat*nlon*nt+1
      call shsgc1(nlat,nlon,l,lat,mode,g,idg,jdg,nt,a,b,mdab,ndab,
     1wshsgc,wshsgc(ifft),late,work(ipmn),work)
      return
      end
      subroutine shsgc1(nlat,nlon,l,lat,mode,gs,idg,jdg,nt,a,b,mdab,
     1                  ndab,w,wfft,late,pmn,g)
      dimension gs(idg,jdg,nt),a(mdab,ndab,nt),b(mdab,ndab,nt)
      dimension w(1),pmn(nlat,late,3),g(lat,nlon,nt),wfft(1)
c     reconstruct fourier coefficients in g on gaussian grid
c     using coefficients in a,b
c     set m+1 limit for b coefficient calculation
      lm1 = l
      if (nlon .eq. l+l-2) lm1 = l-1
c     initialize to zero
      do 100 k=1,nt
      do 100 j=1,nlon
      do 100 i=1,lat
      g(i,j,k) = 0.0
  100 continue
      if (mode.eq.0) then
c     set first column in g
      m = 0
c     compute pmn for all i and n=m,...,l-1
      call legin(mode,l,nlat,m,w,pmn,km)
      do 101 k=1,nt
c     n even
      do 102 np1=1,nlat,2
      do 102 i=1,late
      g(i,1,k) = g(i,1,k)+a(1,np1,k)*pmn(np1,i,km)
  102 continue
c     n odd
      nl2 = nlat/2
      do 103 np1=2,nlat,2
      do 103 i=1,nl2
      is = nlat-i+1
      g(is,1,k) = g(is,1,k)+a(1,np1,k)*pmn(np1,i,km)
  103 continue
c     restore m=0 coefficents (reverse implicit even/odd reduction)
      do 112 i=1,nl2
      is = nlat-i+1
      t1 = g(i,1,k)
      t3 = g(is,1,k)
      g(i,1,k) = t1+t3
      g(is,1,k) = t1-t3
  112 continue
  101 continue
c     sweep  columns of g for which b is available
      do 104 mp1=2,lm1
      m = mp1-1
      mp2 = m+2
c     compute pmn for all i and n=m,...,l-1
      call legin(mode,l,nlat,m,w,pmn,km)
      do 105 k=1,nt
c     for n-m even store (g(i,p,k)+g(nlat-i+1,p,k))/2 in g(i,p,k) p=2*m,
c     for i=1,...,late
      do 106 np1=mp1,nlat,2
      do 107 i=1,late
      g(i,2*m,k) = g(i,2*m,k)+a(mp1,np1,k)*pmn(np1,i,km)
      g(i,2*m+1,k) = g(i,2*m+1,k)+b(mp1,np1,k)*pmn(np1,i,km)
  107 continue
  106 continue
c     for n-m odd store g(i,p,k)-g(nlat-i+1,p,k) in g(nlat-i+1,p,k)
c     for i=1,...,nlat/2 (p=2*m,p=2*m+1)
      do 108 np1=mp2,nlat,2
      do 109 i=1,nl2
      is = nlat-i+1
      g(is,2*m,k) = g(is,2*m,k)+a(mp1,np1,k)*pmn(np1,i,km)
      g(is,2*m+1,k) = g(is,2*m+1,k)+b(mp1,np1,k)*pmn(np1,i,km)
  109 continue
  108 continue
c     now set fourier coefficients using even-odd reduction above
      do 110 i=1,nl2
      is = nlat-i+1
      t1 = g(i,2*m,k)
      t2 = g(i,2*m+1,k)
      t3 = g(is,2*m,k)
      t4 = g(is,2*m+1,k)
      g(i,2*m,k) = t1+t3
      g(i,2*m+1,k) = t2+t4
      g(is,2*m,k) = t1-t3
      g(is,2*m+1,k) = t2-t4
  110 continue
  105 continue
  104 continue
c     set last column (using a only)
      if (nlon.eq. l+l-2) then
      m = l-1
      call legin(mode,l,nlat,m,w,pmn,km)
      do 111 k=1,nt
c     n-m even
      do 131 np1=l,nlat,2
      do 131 i=1,late
      g(i,nlon,k) = g(i,nlon,k)+2.0*a(l,np1,k)*pmn(np1,i,km)
  131 continue
      lp1 = l+1
c     n-m odd
      do 132 np1=lp1,nlat,2
      do 132 i=1,nl2
      is = nlat-i+1
      g(is,nlon,k) = g(is,nlon,k)+2.0*a(l,np1,k)*pmn(np1,i,km)
  132 continue
      do 133 i=1,nl2
      is = nlat-i+1
      t1 = g(i,nlon,k)
      t3 = g(is,nlon,k)
      g(i,nlon,k)= t1+t3
      g(is,nlon,k)= t1-t3
  133 continue
  111 continue
      end if
      else
c     half sphere (mode.ne.0)
c     set first column in g
      m = 0
      meo = 1
      if (mode.eq.1) meo = 2
      ms = m+meo
c     compute pmn for all i and n=m,...,l-1
      call legin(mode,l,nlat,m,w,pmn,km)
      do 113 k=1,nt
      do 113 np1=ms,nlat,2
      do 113 i=1,late
      g(i,1,k) = g(i,1,k)+a(1,np1,k)*pmn(np1,i,km)
  113 continue
c     sweep interior columns of g
      do 114 mp1=2,lm1
      m = mp1-1
      ms = m+meo
c     compute pmn for all i and n=m,...,l-1
      call legin(mode,l,nlat,m,w,pmn,km)
      do 115 k=1,nt
      do 115 np1=ms,nlat,2
      do 115 i=1,late
      g(i,2*m,k) = g(i,2*m,k)+a(mp1,np1,k)*pmn(np1,i,km)
      g(i,2*m+1,k) = g(i,2*m+1,k)+b(mp1,np1,k)*pmn(np1,i,km)
  115 continue
  114 continue
      if (nlon.eq.l+l-2) then
c     set last column
      m = l-1
      call legin(mode,l,nlat,m,w,pmn,km)
      ns = l
      if (mode.eq.1) ns = l+1
      do 116 k=1,nt
      do 116 i=1,late
      do 116 np1=ns,nlat,2
      g(i,nlon,k) = g(i,nlon,k)+2.0*a(l,np1,k)*pmn(np1,i,km)
  116 continue
      end if
      end if
c     do inverse fourier transform
      do 120 k=1,nt
      call hrfftb(lat,nlon,g(1,1,k),lat,wfft,pmn)
  120 continue
c     scale output in gs
      do 122 k=1,nt
      do 122 j=1,nlon
      do 122 i=1,lat
      gs(i,j,k) = 0.5*g(i,j,k)
  122 continue
      return
      end
      subroutine shsgci(nlat,nlon,wshsgc,lshsgc,dwork,ldwork,ierror)
c     this subroutine must be called before calling shsgc with
c     fixed nlat,nlon. it precomputes quantites such as the gaussian
c     points and weights, m=0,m=1 legendre polynomials, recursion
c     recursion coefficients.
      dimension wshsgc(lshsgc)
      double precision dwork(ldwork)
      ierror = 1
      if (nlat.lt.3) return
      ierror = 2
      if (nlon.lt.4) return
c     set triangular truncation limit for spherical harmonic basis
      l = min0((nlon+2)/2,nlat)
c     set equator or nearest point (if excluded) pointer
      late = (nlat+mod(nlat,2))/2
      l1 = l
      l2 = late
      ierror = 3
c     check permanent work space length
      if (lshsgc .lt. nlat*(2*l2+3*l1-2)+3*l1*(1-l1)/2+nlon+15)return
      ierror = 4
      if (ldwork .lt. nlat*(nlat+4)) return
      ierror = 0
c     set pointers
      i1 = 1
      i2 = i1+nlat
      i3 = i2+nlat*late
      i4 = i3+nlat*late
      i5 = i4+l*(l-1)/2 +(nlat-l)*(l-1)
      i6 = i5+l*(l-1)/2 +(nlat-l)*(l-1)
      i7 = i6+l*(l-1)/2 +(nlat-l)*(l-1)
c     set indices in temp work for double precision gaussian wts and pts
      idth = 1
      idwts = idth+nlat
      iw = idwts+nlat
      call shsgci1(nlat,nlon,l,late,wshsgc(i1),wshsgc(i2),wshsgc(i3),
     1wshsgc(i4),wshsgc(i5),wshsgc(i6),wshsgc(i7),dwork(idth),
     2dwork(idwts),dwork(iw),ierror)
      if (ierror.ne.0) ierror = 5
      return
      end
      subroutine shsgci1(nlat,nlon,l,late,wts,p0n,p1n,abel,bbel,cbel,
     1                  wfft,dtheta,dwts,work,ier)
      dimension wts(nlat),p0n(nlat,late),p1n(nlat,late),abel(1),bbel(1),
     1 cbel(1),wfft(1),dtheta(nlat),dwts(nlat)
      double precision pb,dtheta,dwts,work(*)
c     compute the nlat  gaussian points and weights, the
c     m=0,1 legendre polys for gaussian points and all n,
c     and the legendre recursion coefficients
c     define index function used in storing
c     arrays for recursion coefficients (functions of (m,n))
c     the index function indx(m,n) is defined so that
c     the pairs (m,n) map to [1,2,...,indx(l-1,l-1)] with no
c     "holes" as m varies from 2 to n and n varies from 2 to l-1.
c     (m=0,1 are set from p0n,p1n for all n)
c     define for 2.le.n.le.l-1
      indx(m,n) = (n-1)*(n-2)/2+m-1
c     define index function for l.le.n.le.nlat
      imndx(m,n) = l*(l-1)/2+(n-l-1)*(l-1)+m-1
c     preset quantites for fourier transform
      call hrffti(nlon,wfft)
c     compute double precision gaussian points and weights
c     lw = 4*nlat*(nlat+1)+2
      lw = nlat*(nlat+2)
      call gaqd(nlat,dtheta,dwts,work,lw,ier)
      if (ier.ne.0) return
c     store gaussian weights single precision to save computation
c     in inner loops in analysis
      do 100 i=1,nlat
      wts(i) = dwts(i)
  100 continue
c     initialize p0n,p1n using double precision dnlfk,dnlft
      do 101 np1=1,nlat
      do 101 i=1,late
      p0n(np1,i) = 0.0
      p1n(np1,i) = 0.0
  101 continue
c     compute m=n=0 legendre polynomials for all theta(i)
      np1 = 1
      n = 0
      m = 0
      call dnlfk(m,n,work)
      do 103 i=1,late
      call dnlft(m,n,dtheta(i),work,pb)
      p0n(1,i) = pb
  103 continue
c     compute p0n,p1n for all theta(i) when n.gt.0
      do 104 np1=2,nlat
      n = np1-1
      m = 0
      call dnlfk(m,n,work)
      do 105 i=1,late
      call dnlft(m,n,dtheta(i),work,pb)
      p0n(np1,i) = pb
  105 continue
c     compute m=1 legendre polynomials for all n and theta(i)
      m = 1
      call dnlfk(m,n,work)
      do 106 i=1,late
      call dnlft(m,n,dtheta(i),work,pb)
      p1n(np1,i) = pb
  106 continue
  104 continue
c     compute and store swarztrauber recursion coefficients
c     for 2.le.m.le.n and 2.le.n.le.nlat in abel,bbel,cbel
      do 107 n=2,nlat
      mlim = min0(n,l)
      do 107 m=2,mlim
      imn = indx(m,n)
      if (n.ge.l) imn = imndx(m,n)
      abel(imn)=sqrt(float((2*n+1)*(m+n-2)*(m+n-3))/
     1               float(((2*n-3)*(m+n-1)*(m+n))))
      bbel(imn)=sqrt(float((2*n+1)*(n-m-1)*(n-m))/
     1               float(((2*n-3)*(m+n-1)*(m+n))))
      cbel(imn)=sqrt(float((n-m+1)*(n-m+2))/
     1               float(((n+m-1)*(n+m))))
  107 continue
      return
      end


text:
! The MIT License (MIT)
!
! Copyright (c) 2016 odlomax
!
! 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.

! node for image quadtree
module m_image_tree_node

   use m_kind_parameters
   use m_maths
   use m_particle

   implicit none

   private
   public :: image_tree_node
   public :: proj_particle

   ! module parameters
   integer(kind=int_kind),parameter :: n_dim_im=2                  ! dimensions of image

   ! projected particle type
   type :: proj_particle

      real(kind=rel_kind) :: r(n_dim_im)                        ! projected particle position
      real(kind=rel_kind) :: h                                  ! particle smoothing length

      contains

      procedure,non_overridable :: initialise=>initialise_proj_particle

   end type

   ! define quadtree node type
   type :: image_tree_node

      real(kind=rel_kind) :: aabb(n_dim_im,2)                  ! axis aligned bounding box (grid volume)
      integer(kind=int_kind) :: n_leaf                         ! number of leaves on tree
      integer(kind=int_kind) :: level                          ! level of tree
      integer(kind=int_kind) :: max_level                      ! maximum level of tree
      real(kind=rel_kind) :: sigma                             ! column density
      real(kind=rel_kind),allocatable :: t_sigma(:)            ! array of temperature differential column densities
      real(kind=rel_kind),allocatable :: t_moments(:)          ! temperature moments
      real(kind=rel_kind),allocatable :: i_lambda(:)           ! intensity array
      real(kind=rel_kind),allocatable :: j_lambda(:)           ! intensity array
      type(proj_particle),pointer,contiguous :: particle_array(:) ! particle array
      type(image_tree_node),pointer,contiguous :: children(:)  ! child node array

      contains

      procedure,non_overridable :: initialise=>initialise_image_tree_node
      procedure,non_overridable :: destroy
      procedure,non_overridable :: get_leaf_aabb
      procedure,non_overridable :: set_leaf_i_lambda
      procedure,non_overridable :: get_node_pointer
      procedure,non_overridable :: write_leaf
      procedure,non_overridable :: rebuild_tree

   end type

   contains

   pure subroutine initialise_proj_particle(self,sph_particle,unit_vector)

      ! argument declarations
      class(proj_particle),intent(inout) :: self               ! projected particle class
      type(particle),intent(in) :: sph_particle                ! sph particle
      real(kind=rel_kind),intent(in) :: unit_vector(n_dim,n_dim_im)  ! unit vectors of image plane

      ! variable declarations
      integer(kind=int_kind) :: i                              ! counter

      ! set position
      do i=1,n_dim_im
         self%r(i)=dot_product(sph_particle%r,unit_vector(:,i))
      end do

      ! set smoothing length
      self%h=sph_particle%h

      return

   end subroutine

   pure recursive subroutine initialise_image_tree_node(self,particle_array,aabb,level)

      ! argument declarations
      class(image_tree_node),intent(inout) :: self             ! quadtree node
      type(proj_particle),intent(inout),target :: particle_array(:)! array of particles
      real(kind=rel_kind),intent(in) :: aabb(n_dim_im,2)       ! axis aligned bounding box
      integer(kind=int_kind),intent(in),optional :: level      ! level on tree

      ! set level
      if (present(level)) then
         self%level=level
      else
         self%level=0
      end if

      ! associate particle array
      self%particle_array=>particle_array

      ! set aabb
      self%aabb=aabb

      ! build more sells if minimum smoothing length is smaller than cell
      if (size(self%particle_array)>0) then

         if (minval(self%particle_array%h)<maxval(self%aabb(:,2)-self%aabb(:,1))) then

            ! recursively build children
            allocate(self%children(2**n_dim_im))
            call build_children(self%particle_array,self%children,self%aabb,self%level)
            self%n_leaf=sum(self%children%n_leaf)
            self%max_level=maxval(self%children%max_level)
            return

         end if

      end if

      self%children=>null()
      self%n_leaf=1
      self%max_level=self%level

      return

   end subroutine

   pure recursive subroutine destroy(self)

      ! argument declarations
      class(image_tree_node),intent(inout) :: self              ! quadtree node

      ! variable declarations
      integer(kind=int_kind) :: i                               ! counter

      if (allocated(self%i_lambda)) deallocate(self%i_lambda)
      if (allocated(self%j_lambda)) deallocate(self%j_lambda)
      if (allocated(self%t_sigma)) deallocate(self%t_sigma)
      if (allocated(self%t_moments)) deallocate(self%t_moments)

      if (associated(self%children)) then

         do i=1,size(self%children)
            call self%children(i)%destroy()
         end do
         deallocate(self%children)

      end if

      return

   end subroutine

   ! get array of aabb of all leaves in tree
   pure recursive function get_leaf_aabb(self) result(aabb)

      ! argument declarations
      class(image_tree_node),intent(in) :: self                ! quadtree node

      ! result declaration
      real(kind=rel_kind) :: aabb(n_dim_im,2,self%n_leaf)      ! array of axis-aligned bounding boxes

      ! variable declarations
      integer(kind=int_kind) :: i,j                            ! counter

      if (associated(self%children)) then

         j=1
         do i=1,size(self%children)
            aabb(:,:,j:j+self%children(i)%n_leaf-1)=self%children(i)%get_leaf_aabb()
            j=j+self%children(i)%n_leaf
         end do

      else

         aabb(:,:,1)=self%aabb

      end if

      return

   end function

   ! set assign i_lambda array to all leaves in tree
   pure recursive subroutine set_leaf_i_lambda(self,i_lambda,j_lambda,sigma,t_sigma,t_moments)

      ! argument declarations
      class(image_tree_node),intent(inout) :: self          ! quadtree node
      real(kind=rel_kind),intent(in),contiguous :: i_lambda(:,:)  ! i_lambda array. second dimension must equal self%n_leaf
      real(kind=rel_kind),intent(in),contiguous :: j_lambda(:,:)  ! j_lambda_array
      real(kind=rel_kind),intent(in),contiguous :: sigma(:)       ! column density array
      real(kind=rel_kind),intent(in),contiguous,optional :: t_sigma(:,:)   ! temperature differential column density
      real(kind=rel_kind),intent(in),contiguous,optional :: t_moments(:,:) ! temperature moments

      ! variable declarations
      integer(kind=int_kind) :: i,j                         ! counter

      allocate(self%i_lambda(size(i_lambda,1)))
      allocate(self%j_lambda(size(j_lambda,1)))
      if (present(t_sigma)) allocate(self%t_sigma(size(t_sigma,1)))
      if (present(t_moments)) allocate(self%t_moments(size(t_moments,1)))

      if (associated(self%children)) then

         ! recurse down tree and average i_lambda over children cells
         j=1
         self%i_lambda=0._rel_kind
         self%j_lambda=0._rel_kind
         self%sigma=0._rel_kind
         if (present(t_sigma)) self%t_sigma=0._rel_kind
         if (present(t_moments)) self%t_moments=0._rel_kind
         do i=1,size(self%children)
            call self%children(i)%set_leaf_i_lambda(i_lambda(:,j:j+self%children(i)%n_leaf-1),&
               &j_lambda(:,j:j+self%children(i)%n_leaf-1),sigma(j:j+self%children(i)%n_leaf-1),&
               &t_sigma(:,j:j+self%children(i)%n_leaf-1),t_moments(:,j:j+self%children(i)%n_leaf-1))
            j=j+self%children(i)%n_leaf
            self%i_lambda=self%i_lambda+self%children(i)%i_lambda
            self%j_lambda=self%j_lambda+self%children(i)%j_lambda
            self%sigma=self%sigma+self%children(i)%sigma
            if (present(t_sigma)) self%t_sigma=self%t_sigma+self%children(i)%t_sigma
            if (present(t_moments)) self%t_moments=self%children(i)%t_moments*self%children(i)%sigma
         end do
         self%i_lambda=self%i_lambda/real(size(self%children),rel_kind)
         self%j_lambda=self%j_lambda/real(size(self%children),rel_kind)
         self%sigma=self%sigma/real(size(self%children),rel_kind)
         if (present(t_sigma)) self%t_sigma=self%t_sigma/real(size(self%children),rel_kind)
         if (present(t_moments)) self%t_moments=self%t_moments/sum(self%children%sigma)
      else

         self%i_lambda=i_lambda(:,1)
         self%j_lambda=j_lambda(:,1)
         self%sigma=sigma(1)
         if (present(t_sigma)) self%t_sigma=t_sigma(:,1)
         if (present(t_moments)) self%t_moments=t_moments(:,1)

      end if

      return

   end subroutine

   ! returns pointer to the node containing r
   recursive function get_node_pointer(self,r,max_level_in) result(node_pointer)

      ! argument declarations
      class(image_tree_node),intent(in),target :: self   ! quadtree node
      real(kind=rel_kind),intent(in) :: r(n_dim_im)      ! position
      integer(kind=int_kind),intent(in),optional :: max_level_in  ! maximum recursion level

      ! result declaration
      class(image_tree_node),pointer :: node_pointer     ! i_lambda

      ! variable declarations
      integer(kind=int_kind) :: j                        ! counter
      integer(kind=int_kind) :: max_level                ! maximum recursion level

      ! set max level to max_level_in or really big number
      if (present(max_level_in)) then
         max_level=max_level_in
      else
         max_level=huge(1)
      end if

      if (associated(self%children).and.self%level<max_level) then

         ! find child cell which contains r
         do j=1,size(self%children)

            if (all(r>=self%children(j)%aabb(:,1).and.r<self%children(j)%aabb(:,2))) then
               node_pointer=>self%children(j)%get_node_pointer(r,max_level_in)
               return
            end if

         end do

         ! if r is not in child cells, pointer set to null
         node_pointer=>null()

      else

         ! found correct cell
         node_pointer=>self

      end if

      return

   end function

   ! rebuild branch nodes from modified leaves
   pure recursive subroutine rebuild_tree(self)

      ! argument declarations
      class(image_tree_node),intent(inout) :: self                ! quadtree node

      ! variable declarations
      integer(kind=int_kind) :: i                                 ! counter

      if (associated(self%children)) then

         self%i_lambda=0._rel_kind
         self%j_lambda=0._rel_kind
         do i=1,size(self%children)

            call self%children(i)%rebuild_tree()
            self%i_lambda=self%i_lambda+self%children(i)%i_lambda
            self%j_lambda=self%j_lambda+self%children(i)%j_lambda

         end do
         self%i_lambda=self%i_lambda/real(size(self%children),rel_kind)
         self%j_lambda=self%j_lambda/real(size(self%children),rel_kind)

      end if

      return

   end subroutine

   recursive subroutine write_leaf(self,file_id)

      ! argument declarations
      class(image_tree_node),intent(in) :: self                   ! quadtree node
      integer(kind=int_kind),intent(in) :: file_id                ! file identifier

      ! variable declarations
      integer(kind=int_kind) :: i                                 ! counter
      character(kind=chr_kind,len=string_length) :: format_string ! character format string

      if (associated(self%children)) then

         do i=1,size(self%children)
            call self%children(i)%write_leaf(file_id)
         end do

      else

         if (allocated(self%t_sigma).and.allocated(self%t_moments)) then

            write(format_string,"(A,I0,A)") "(",size(self%aabb)+1+size(self%t_sigma)+&
               &size(self%t_moments)+2*size(self%i_lambda),"(E25.17))"
            write(file_id,trim(format_string)) self%aabb,self%sigma,self%t_sigma,&
               &nth_rt(self%t_moments,(/(i,i=1,size(self%t_moments))/)),self%i_lambda,self%j_lambda

         else

            write(format_string,"(A,I0,A)") "(",size(self%aabb)+1+2*size(self%i_lambda),"(E25.17))"
            write(file_id,trim(format_string)) self%aabb,self%sigma,self%i_lambda,self%j_lambda

         end if

      end if

      return

   end subroutine

   ! partition particle array about boundary along dimension k_dim
   pure subroutine partition_particles(particle_array,k_part,k_dim,boundary)

      ! argument declarations
      type(proj_particle),intent(inout),contiguous :: particle_array(:) ! array of particles
      integer(kind=int_kind),intent(out) :: k_part              ! partition partilce
      integer(kind=int_kind),intent(in) :: k_dim                ! splitting dimension
      real(kind=rel_kind),intent(in) :: boundary                ! split array about this boundary

      ! variable declarations
      type(proj_particle) :: temp_particle                      ! temporary particle
      integer(kind=int_kind) :: i                               ! counter

      k_part=0

      ! k_part = number of particles to the left of boundary
      do i=1,size(particle_array)

         if (particle_array(i)%r(k_dim)<boundary) then

            k_part=k_part+1
            temp_particle=particle_array(k_part)
            particle_array(k_part)=particle_array(i)
            particle_array(i)=temp_particle

         end if

      end do

      return

   end subroutine

   ! build children nodes
   pure recursive subroutine build_children(particle_array,children,aabb,tree_level,centre_in,level_in,lhs_in)

      ! argument declarations
      type(proj_particle),intent(inout),contiguous,target :: particle_array(:)   ! array of particles
      type(image_tree_node),intent(inout),contiguous :: children(:)         ! array of child nodes
      real(kind=rel_kind),intent(in) :: aabb(n_dim_im,2)                    ! parent axis aligned bounding box
      integer(kind=int_kind),intent(in) :: tree_level                       ! level on tree
      real(kind=rel_kind),intent(in),optional :: centre_in(n_dim_im)        ! centre of aabb
      integer(kind=int_kind),intent(in),optional :: level_in                ! level of recursion (1 to n_dim_im)
      logical(kind=log_kind),intent(inout),optional :: lhs_in(n_dim_im)     ! (left hand side) trace which cell we're in

      ! variable declarations
      integer(kind=int_kind) :: k_part                                      ! partition element
      integer(kind=int_kind) :: level                                       ! level of recursion
      real(kind=rel_kind) :: new_aabb(n_dim_im,2)                           ! new axis aligned bounding box
      real(kind=rel_kind) :: centre(n_dim_im)                               ! centre of aabb
      logical(kind=log_kind) :: lhs(n_dim_im)

      ! optional arguments aren't set on first call
      ! set aabb centre
      if (present(centre_in)) then
         centre=centre_in
      else
         centre=0.5_rel_kind*(aabb(:,1)+aabb(:,2))
      end if
      ! set level
      if (present(level_in)) then
         level=level_in
      else
         level=1
      end if
      ! set lhs
      if (present(lhs_in)) lhs=lhs_in

      ! if level is less than n_dim_im, split up the volume
      ! when level exceeds n_dim_im, assign particles to child nodes
      if (level<=n_dim_im) then

         ! partition particle array
         call partition_particles(particle_array,k_part,level,centre(level))

         ! split up particle array and children until children array has only one element
         lhs(level)=.true.
         call build_children(particle_array(:k_part),&
            &children(:size(children)/2),aabb,tree_level,centre,level+1,lhs)

         lhs(level)=.false.
         call build_children(particle_array(k_part+1:),&
            &children(size(children)/2+1:),aabb,tree_level,centre,level+1,lhs)

      else

         ! children array now only has a single element
         new_aabb(:,1)=merge(aabb(:,1),centre,lhs)
         new_aabb(:,2)=merge(centre,aabb(:,2),lhs)
         call children(1)%initialise(particle_array,new_aabb,tree_level+1)

      end if

      return

   end subroutine

end module

max_stars_repo_name: odlomax/spamcart-dev

max_issues_repo_name: odlomax/spamcart-dev

max_forks_repo_name: odlomax/spamcart-dev


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!!  MFE1 --
!!
!!  This is a Fortran 90 implementation of the piecewise-linear
!!  Gradient-Weighted Moving Finite Element Method (GWMFE) for
!!  time-dependent systems of partial differential equations in
!!  one space dimension.
!!
!!  Version 0.3, 12 August 1997
!!
!!  Neil N. Carlson, Dept of Math, Purdue University
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Copyright (c) 1997, Neil N. Carlson
!!
!! 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.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program mfe1

  use mfe_constants
  use mfe_types
  use common_io
  use system_io
  use output
  use initialize
  use mfe_ode_solver
  use mfe_procs, only: eval_udot
  implicit none

  real(kind=wp) :: t
  real(kind=wp), dimension(6) :: rvar
  integer :: mode, rtype, j, errc, debug_unit, nstep
  integer, dimension(3) :: ivar
  type(NodeVar), dimension(:), pointer :: u, udot

  character(len=16) :: string

  real :: cpusec, cpusec0

  call cpu_time (cpusec0)

  ! Open the input and output files.

  call new_unit (input_unit)
  open (unit=input_unit, file="mfein", position="rewind", action="read", status="old")

  call new_unit (log_unit)
  open (unit=log_unit, file="mfelog", position="rewind", action="write", status="replace")

  call new_unit (out_unit)
  open (unit=out_unit, file="mfegrf", position="rewind", action="write", status="replace")

  call read_soln (u, udot)
  call read_data ()

  t = tout(1)
  mode = START_SOLN
  call eval_udot (u, t, udot, errc)
  call write_soln (u, t)
  j = 2

  if (errc /= 0) then
    call abort ( (/log_unit,STDERR/), "Bad initial solution.")
  end if

  if (mstep <= 0) then
    stop
  end if

  if (debug /= 0) then
    call new_unit (debug_unit)
    open (unit=debug_unit, file="bdfout", action="write", status="replace")
    call set_solver_messages (debug, debug_unit)
  end if

  rvar(1) = h
  rvar(2) = hlb
  rvar(3) = hub
  rvar(4) = ntol
  rvar(5) = margin
  rvar(6) = vtol

  ivar(1) = mtry
  ivar(2) = mitr
  ivar(3) = mvec

  write (unit=string, fmt="(es11.3)" ) t
  call info ( (/ log_unit, STDOUT /), "BEGINNING SOLUTION AT T = " // string )

  if (ofreq <= 0) then
    ofreq = mstep
  end if

  do

    call bdf2_solver (mode, rvar, ivar, tout(j), ofreq, rtype, u, udot, t)

    mode = RESUME_SOLN

    call cpu_time (cpusec)
    cpusec = cpusec - cpusec0
    call write_soln (u, t)

    select case (rtype)

      case (SOLN_AT_TOUT)       ! Integrated to TOUT.
        call write_status( (/ log_unit, STDOUT /), cpusec, t )
        j = j + 1

      case (SOLN_AT_STEP)       ! Integrated OFREQ more steps.
        call write_status( (/ log_unit, STDOUT /), cpusec )

      case (FAIL_ON_STEP)
        call write_status( (/ log_unit, STDOUT /), cpusec )
        call abort ( (/ log_unit, STDERR /), "Repeated failure at a step.")

      case (SMALL_H_FAIL)
        call write_status( (/ log_unit, STDOUT /), cpusec )
        call abort ( (/ log_unit, STDERR /), "Next time step is too small.")

      case (FAIL_ON_START)
        call write_status( (/ log_unit, STDOUT /), cpusec )
        call abort ( (/ log_unit, STDERR /), "Starting procedure failed.")

      case default
        call write_status( (/ log_unit, STDOUT /), cpusec )
        call abort ( (/ log_unit, STDERR /), "Unknown return type!")

    end select

    if (j > size(tout)) then
      call info( (/ log_unit, STDOUT /), "Integrated to final TOUT.  Done." )
      exit
    end if

    call bdf2_inquire( nstep=nstep )
    if (nstep >= mstep) then
      call info( (/ log_unit, STDOUT /), "Maximum number of steps taken.  Done." )
      exit
    end if

  end do

end program mfe1