The ParaMonte Documentation Website
Current view: top level - kernel - ParaDXXX_mod@Kernel_smod@nextMove.inc.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 54 58 93.1 %
Date: 2021-01-08 13:07:16 Functions: 0 0 -
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : #if defined SINGLCHAIN_PARALLELISM
      44             : #define LOOP_NEXT_MOVE loopNextMoveSinglChain
      45             : #else
      46             : #define LOOP_NEXT_MOVE loopNextMove
      47             : #endif
      48             : 
      49             : #if defined SINGLCHAIN_PARALLELISM
      50      281394 :                 proposalFoundSinglChainMode = 0_IK
      51             : #if defined CAF_ENABLED
      52             :                 ! This syncing is necessary since the co_LogFuncState has to be fetched from the 
      53             :                 ! first image by all other images before it is updated again below by the first image.
      54             :                 if (self%Image%isLeader) then
      55             :                     call self%Timer%toc()
      56             :                     sync images(*)
      57             :                     call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
      58             :                 else
      59             :                     sync images(1)
      60             :                 end if
      61             : #endif
      62             : #endif
      63             : 
      64      749243 :                 LOOP_NEXT_MOVE : do counterDRS = 0, self%SpecDRAM%DelayedRejectionCount%val
      65             : 
      66             : #if defined SINGLCHAIN_PARALLELISM
      67             :                     ! self%Stats%NumFunCall%acceptedRejectedDelayedUnused is relevant only on the first image, despite being updated by all images
      68      318840 :                     self%Stats%NumFunCall%acceptedRejectedDelayedUnused = self%Stats%NumFunCall%acceptedRejectedDelayedUnused + self%Image%count
      69             : #endif
      70             : 
      71             :                     co_LogFuncState(1:nd,counterDRS) = self%Proposal%getNew ( nd            = nd & ! LCOV_EXCL_LINE
      72             :                                                                             , counterDRS    = counterDRS & ! LCOV_EXCL_LINE
      73             :                                                                             , StateOld      = co_LogFuncState(1:nd,counterDRS-1) & ! LCOV_EXCL_LINE
      74     1150833 :                                                                             )
      75             : #if (defined MATLAB_ENABLED || defined PYTHON_ENABLED || defined R_ENABLED) && !defined CAF_ENABLED && !defined MPI_ENABLED
      76             :                     if(ProposalErr%occurred) then; self%Err%occurred = .true.; self%Err%msg = ProposalErr%msg; return; end if
      77             : #elif defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED
      78             :                     ! This block is exclusively used to test the deterministic restart functionality of ParaDXXX samplers. 
      79             :                     ! This block must not be activated under any other circumstances.
      80             :                     ! This block must be executed by all images.
      81      416105 :                     self%testSamplingCounter = self%testSamplingCounter + 1_IK
      82      416105 :                     self%Err%occurred = self%Err%occurred .or. self%testSamplingCounter > self%testSamplingCountTarget
      83      416105 :                     if (self%Err%occurred) then
      84           0 :                         if (self%testSamplingCounter > self%testSamplingCountTarget) self%Err%msg = "The simulation was interrupted at the requested sampling count for restart testing purposes."
      85             : #if defined MPI_ENABLED || defined CAF_ENABLED
      86           0 :                         block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
      87             : #endif
      88           0 :                         return
      89             :                     end if
      90             : #endif
      91             : 
      92             :                     ! The following random function call is only needed fresh runs to evaluate the acceptance of a proposal.
      93             :                     ! However, it is taken out of the subsequent loops to achieve 100% deterministic reproducibility when
      94             :                     ! a simulation is restarted.
      95             : 
      96      416105 :                     call random_number(uniformRnd)
      97             : 
      98      487001 :                     if (self%isFreshRun .or. numFunCallAcceptedPlusOne==self%Chain%Count%compact) then
      99             : 
     100      327672 :                         call self%Timer%toc()
     101      327672 :                         co_LogFuncState(0,counterDRS) = getLogFunc(nd,co_LogFuncState(1:nd,counterDRS))
     102      327672 :                         call self%Timer%toc(); self%Stats%avgTimePerFunCalInSec = self%Stats%avgTimePerFunCalInSec + self%Timer%Time%delta
     103             : 
     104             :                         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     105             :                         ! accept or reject the proposed state
     106             :                         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     107             : 
     108      327672 :                         if ( co_LogFuncState(0,counterDRS) >= co_LogFuncState(0,-1) ) then ! accept the proposed state
     109             : 
     110       63104 :                             co_AccRate(counterDRS) = 1._RK
     111       63104 :                             co_AccRate(-1) = real(counterDRS,kind=RK)
     112      224200 :                             co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,counterDRS)
     113             : #if !defined SINGLCHAIN_PARALLELISM
     114        7419 :                             exit LOOP_NEXT_MOVE
     115             : #endif
     116             : 
     117      264568 :                             elseif ( co_LogFuncState(0,counterDRS) < maxLogFuncRejectedProposal ) then  ! reject the proposed state. This step should be reachable only when delayedRejectionCount > 0
     118             : 
     119       15021 :                                 co_AccRate(counterDRS) = 0._RK  ! proposal rejected XXX is this correct? could co_AccRatebe absolute zero?
     120             : 
     121             :                             else    ! accept with probability co_AccRate
     122             : 
     123      249547 :                             if ( counterDRS == 0_IK ) then ! This should be equivalent to maxLogFuncRejectedProposal == NEGINF_RK
     124      245418 :                                 logFuncDiff = co_LogFuncState(0,counterDRS) - co_LogFuncState(0,-1)
     125      245418 :                                 if (logFuncDiff < NEGLOGINF_RK) then ! xxx should the condition for LOGHUGE_RK be also added?
     126       18383 :                                     co_AccRate(counterDRS) = 0._RK
     127             :                                 else
     128      227035 :                                     co_AccRate(counterDRS) = exp( logFuncDiff )
     129             :                                 end if
     130             :                             else    ! ensure no arithmetic overflow/underflow. ATT: co_LogFuncState(0,-1) > co_LogFuncState(0,counterDRS) > maxLogFuncRejectedProposal
     131             :                                 co_AccRate(counterDRS)  = getLogSubExp( co_LogFuncState(0,counterDRS)   , maxLogFuncRejectedProposal ) & ! LCOV_EXCL_LINE
     132        4129 :                                                         - getLogSubExp( co_LogFuncState(0,-1)           , maxLogFuncRejectedProposal )
     133        4129 :                                 if (co_AccRate(counterDRS) < NEGLOGINF_RK) then
     134           0 :                                     co_AccRate(counterDRS) = 0._RK
     135             :                                 else
     136        4129 :                                     co_AccRate(counterDRS) = exp(co_AccRate(counterDRS))
     137             :                                 end if
     138             :                             end if
     139             : 
     140      249547 :                             if (uniformRnd<co_AccRate(counterDRS)) then ! accept the proposed state
     141       61855 :                                 co_AccRate(-1) = real(counterDRS,kind=RK)
     142      219463 :                                 co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,counterDRS)
     143             : #if !defined SINGLCHAIN_PARALLELISM
     144        7586 :                                 exit LOOP_NEXT_MOVE
     145             : #endif
     146             :                             end if
     147             : 
     148             :                         end if
     149             : 
     150      312667 :                         maxLogFuncRejectedProposal = max( maxLogFuncRejectedProposal, co_LogFuncState(0,counterDRS) )
     151             : 
     152             :                     else ! if dryrun
     153             : 
     154             : #if defined SINGLCHAIN_PARALLELISM
     155       86844 :                         if (self%Image%id==self%Chain%ProcessID(numFunCallAcceptedPlusOne) .and. &
     156       86844 :                             currentStateWeight+self%Image%id-1_IK==self%Chain%Weight(self%Stats%NumFunCall%accepted) .and. &
     157       43422 :                             counterDRS==self%Chain%DelRejStage(numFunCallAcceptedPlusOne)) then
     158        6972 :                             co_AccRate(-1) = real(counterDRS,kind=RK)
     159        6972 :                             co_LogFuncState(   0,counterDRS) = self%Chain%LogFunc   (numFunCallAcceptedPlusOne)
     160       20916 :                             co_LogFuncState(1:nd,counterDRS) = self%Chain%State(1:nd,numFunCallAcceptedPlusOne)
     161       27888 :                             co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,counterDRS)
     162             :                         end if
     163             : #else
     164       45011 :                         if (currentStateWeight==self%Chain%Weight(self%Stats%NumFunCall%accepted) .and. counterDRS==self%Chain%DelRejStage(numFunCallAcceptedPlusOne)) then
     165       11364 :                             co_AccRate(-1) = real(counterDRS,kind=RK)
     166       11364 :                             co_LogFuncState(   0,-1) = self%Chain%LogFunc   (numFunCallAcceptedPlusOne)
     167       34092 :                             co_LogFuncState(1:nd,-1) = self%Chain%State(1:nd,numFunCallAcceptedPlusOne)
     168       11364 :                             exit LOOP_NEXT_MOVE
     169             :                         end if
     170             : #endif
     171             :                     end if
     172             : 
     173             : #if defined SINGLCHAIN_PARALLELISM
     174             :                     ! broadcast the sampling status from the all images to all images.
     175             :                     ! This is needed in singlChain parallelism to avoid unnecessary delayed rejection if a proposal is already found. 
     176             :                     ! Note that this is not necessary when the delayed rejection is deactivated.
     177      581082 :                     if (delayedRejectionRequested) then
     178       64272 :                         if (co_AccRate(-1) > -0.5_RK) proposalFoundSinglChainMode = 1_IK ! -0.5 instead of -1. avoids possible real roundoff errors.
     179       64272 :                         call self%Timer%toc()
     180             : #if defined CAF_ENABLED
     181             :                         call co_sum(proposalFoundSinglChainMode)
     182             : #elif defined MPI_ENABLED
     183             :                         call mpi_allreduce  ( proposalFoundSinglChainMode           & ! send buffer
     184             :                                             , proposalFoundSinglChainModeReduced    & ! recv buffer
     185             :                                             , 1                                     & ! buffer size
     186             :                                             , mpi_integer                           & ! datatype
     187             :                                             , mpi_sum                               & ! mpi reduction operation
     188             :                                             , mpi_comm_world                        & ! comm
     189             :                                             , ierrMPI                               & ! ierr
     190       64272 :                                             )
     191       64272 :                         proposalFoundSinglChainMode = proposalFoundSinglChainModeReduced
     192             : #endif
     193       64272 :                         call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     194       64272 :                         if (proposalFoundSinglChainMode>0_IK) exit LOOP_NEXT_MOVE
     195             :                     end if
     196             : #endif
     197             : 
     198             :                 end do LOOP_NEXT_MOVE
     199             : 
     200             : #if defined SINGLCHAIN_PARALLELISM && defined MPI_ENABLED
     201             :                 ! This SINGLCHAIN_PARALLELISM preprocessing avoids unnecessary communication in serial or multichain-parallel modes.
     202             :                 ! gather all accr on the leader image. 
     203      281394 :                 call self%Timer%toc()
     204             :                 call mpi_gather ( co_AccRate(:)                 & ! send buffer
     205             :                                 , delayedRejectionCountPlusTwo  & ! send count
     206             :                                 , mpi_double_precision          & ! send datatype
     207             :                                 , AccRateMatrix(:,:)            & ! recv buffer
     208             :                                 , delayedRejectionCountPlusTwo  & ! recv count
     209             :                                 , mpi_double_precision          & ! recv datatype
     210             :                                 , 0                             & ! root rank
     211             :                                 , mpi_comm_world                & ! comm
     212             :                                 , ierrMPI                       & ! mpi error flag
     213      281394 :                                 )
     214      281394 :                 call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     215             : #elif defined SINGLCHAIN_PARALLELISM && defined CAF_ENABLED
     216             :                 if (self%Image%isLeader) then
     217             :                     call self%Timer%toc()
     218             :                     sync images(*)
     219             :                     call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     220             :                 else
     221             :                     sync images(1)
     222             :                 end if
     223             : #endif
     224             : 
     225             : #undef LOOP_NEXT_MOVE

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021