Line data    Source code

       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!                                                                                                                            !!!!
       4             : !!!!    ParaMonte: Parallel Monte Carlo and Machine Learning Library.                                                           !!!!
       5             : !!!!                                                                                                                            !!!!
       6             : !!!!    Copyright (C) 2012-present, The Computational Data Science Lab                                                          !!!!
       7             : !!!!                                                                                                                            !!!!
       8             : !!!!    This file is part of the ParaMonte library.                                                                             !!!!
       9             : !!!!                                                                                                                            !!!!
      10             : !!!!    LICENSE                                                                                                                 !!!!
      11             : !!!!                                                                                                                            !!!!
      12             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md                                                          !!!!
      13             : !!!!                                                                                                                            !!!!
      14             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      15             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      16             : 
      17             : !>  \brief
      18             : !>  This include file contains procedure implementation of [pm_ziggurat](@ref pm_ziggurat).
      19             : !>
      20             : !>  \finmain
      21             : !>
      22             : !>  \author
      23             : !>  \AmirShahmoradi, April 25, 2015, 2:21 PM, National Institute for Fusion Studies, The University of Texas at Austin
      24             : 
      25             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      26             : 
      27             :         !%%%%%%%%%%%%%
      28             : #if     getZig_ENABLED
      29             :         !%%%%%%%%%%%%%
      30             : 
      31             :         real(RKC)   :: ub, area
      32             :         real(RKC)   , parameter :: SQRT_HUGE = sqrt(huge(0._RKC))
      33             :         real(RKC)   , parameter :: lb = 10 * epsilon(0._RKC)
      34           3 :         CHECK_ASSERTION(__LINE__, 1_IK < nlay, SK_"@getZig(): The condition `1 <= nlay` must hold. shape(nlay) = "//getStr(nlay))
      35           3 :         zig(1, nlay) = 0._RKC ! The highest rectangle point (origin).
      36           3 :         zig(2, nlay) = getFunc(zig(1, nlay)) ! maxFunc (corresponding to the highest rectangle point)
      37             :         ! Define the upper bound of the search limit.
      38           3 :         ub = lb
      39           9 :         CHECK_ASSERTION(__LINE__, 0 <= getZ(lb), SK_"@getZig(): Internal library error occurred. The condition `0 <= getZ(lb)` must hold. lb, getZ(lb) = "//getStr([lb, getZ(lb)]))
      40             :         do
      41         272 :             ub = 2._RKC * ub
      42         272 :             if (0._RKC <= getZ(ub)) cycle ! \todo should we compare the sign change from zlb to zub here instead?
      43           9 :             CHECK_ASSERTION(__LINE__, ub < SQRT_HUGE, SK_"@getZig(): The specified `getFunc()` does not appear to be monotonically decreasing. The condition `ub < SQRT_HUGE` must hold. ub, SQRT_HUGE = "//getStr([ub, SQRT_HUGE]))
      44         269 :             exit
      45             :         end do
      46           3 :         zig(1, 1) = getRoot(getZ, lb, ub, abstol = abstol) ! r: first/lowest/largest rectangle-corner point.
      47           3 :         zig(1, 0) = area / zig(2, 1)
      48           3 :         zig(2, 0) = 0._RKC
      49             : 
      50             :     contains
      51             : 
      52             :         ! This function returns the difference between the first computed area (from the lowest layer) and the last computed area from the highest layer.
      53         343 :         impure function getZ(r) result(z)
      54             :             real(RKC)   , intent(in)    :: r
      55             :             real(RKC)                   :: z
      56             :             real(RKC)                   :: func
      57             :             integer(IK)                 :: ilay
      58         343 :             zig(1, 1) = r
      59         343 :             zig(2, 1) = getFunc(r)
      60         343 :             area = getZigArea(r) ! v: partition area. Store it in the first (always zero) element to save space.
      61        7516 :             do ilay = 2, nlay - 1
      62        7408 :                 func = zig(2, ilay - 1) + area / zig(1, ilay - 1)
      63        7408 :                 if (zig(2, nlay) < func) then
      64         235 :                     z = SQRT_HUGE
      65         235 :                     abserr = z
      66         235 :                     return
      67             :                 end if
      68        7173 :                 zig(1, ilay) = getFuncInv(func)
      69        7281 :                 zig(2, ilay) = getFunc(zig(1, ilay))
      70             :             end do
      71         108 :             ilay = nlay - 1_IK
      72         108 :             z = area - zig(1, ilay) * (zig(2, nlay) - zig(2, ilay))
      73         108 :             abserr = abs(z)
      74         108 :         end function
      75             : 
      76             : !        !%%%%%%%%%%%%%%%%%%
      77             : !#elif   getZig_ENABLED && 0
      78             : !        !%%%%%%%%%%%%%%%%%%
      79             : !
      80             : !        real(RKC)   :: ub, area
      81             : !        real(RKC)   , parameter :: SQRT_HUGE = sqrt(huge(0._RKC))
      82             : !        real(RKC)   , parameter :: lb = 10 * epsilon(0._RKC)
      83             : !        integer(IK) :: nlayTwice
      84             : !        nlayTwice = 2 * nlay
      85             : !        CHECK_ASSERTION(__LINE__, 1_IK < nlay, SK_"@getZig(): The condition `1 <= nlay` must hold. shape(nlay) = "//getStr(nlay))
      86             : !        zig(nlay) = 0._RKC ! The highest rectangle point (origin).
      87             : !        zig(nlayTwice) = getFunc(zig(nlay)) ! maxFunc (corresponding to the highest rectangle point)
      88             : !        ! Define the upper bound of the search limit.
      89             : !        ub = lb
      90             : !        if (getZ(lb) < 0._RKC) error stop "@getZig(): Internal error occurred. Please report to the ParaMonte library developers."
      91             : !        do
      92             : !            ub = 2._RKC * ub
      93             : !            if (getZ(ub) < 0._RKC) exit
      94             : !        end do
      95             : !        zig(1) = getRoot(getZ, lb, ub) ! r: first/lowest/largest rectangle-corner point.
      96             : !        zig(0) = area / zig(1 + nlay)
      97             : !        !zig(nlayTwice + 1) = getRoot(getZ, lb, ub) ! r: first/lowest/largest rectangle-corner point.
      98             : !        !zig(2 : nlay) = zig(2 : nlay) / zig(1 : nlay - 1)
      99             : !        !zig(1) = zig(1) * zig(1 + nlay) / area ! x1 / x0 = x1 / (area / f(r)) = x1 * f(r) / area.
     100             : !        !!zig(1 : nlay) = zig(1 : nlay) * 2._RKC**31
     101             : !
     102             : !    contains
     103             : !
     104             : !        impure function getZ(r) result(z)
     105             : !            real(RKC)   , intent(in)    :: r
     106             : !            real(RKC)                   :: z
     107             : !            real(RKC)                   :: func
     108             : !            integer(IK)                 :: ilay
     109             : !            zig(1) = r
     110             : !            zig(1 + nlay) = getFunc(r)
     111             : !            area = getZigArea(r) ! v: partition area. Store it in the first (always zero) element to save space.
     112             : !            do ilay = 2_IK, nlay - 1_IK
     113             : !                func = area / zig(ilay - 1) + zig(ilay - 1 + nlay)
     114             : !                if (zig(nlayTwice) < func) then
     115             : !                    z = SQRT_HUGE
     116             : !                    abserr = z
     117             : !                    return
     118             : !                end if
     119             : !                zig(ilay) = getFuncInv(func)
     120             : !                zig(ilay + nlay) = getFunc(zig(ilay))
     121             : !            end do
     122             : !            ilay = nlay - 1_IK
     123             : !            z = area + zig(ilay) * (zig(ilay + nlay) - 1._RKC)
     124             : !            abserr = z
     125             : !        end function
     126             : !
     127             : !        !%%%%%%%%%%%%%%%%
     128             : !#elif   getZigCRD_ENABLED
     129             : !        !%%%%%%%%%%%%%%%%
     130             : !
     131             : !        integer(IK) :: lenZig ! ilay, nlay,
     132             : !        lenZig = size(zig, 1, IK)
     133             : !        CHECK_ASSERTION(__LINE__, mod(lenZig, 2_IK) == 1_IK, SK_"@getZigCRD(): The condition `mod(size(zig), 2) == 1` must hold. shape(size(zig)) = "//getStr(lenZig))
     134             : !        crd(1,:) = zig(1 : lenZig / 2)
     135             : !        crd(2,:) = zig(1 + lenZig / 2 : lenZig - 1)
     136             : !        !nlay = lenZig / 2_IK
     137             : !        !crd(1, 1) = zig(lenZig)
     138             : !        !crd(2, 1) = zig(nlay + 1)
     139             : !        !do ilay = 2_IK, nlay
     140             : !        !    crd(1, ilay) = crd(1, ilay - 1) * zig(ilay) ! * 0.5_RKC**31
     141             : !        !    crd(2, ilay) = zig(nlay + ilay)
     142             : !        !end do
     143             : 
     144             : #else
     145             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     146             : #error  "Unlayognized interface."
     147             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     148             : #endif