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
      19             : !>  [getLogSumExp](@ref pm_mathLogSumExp::getLogSumExp).
      20             : !>
      21             : !>  \author
      22             : !>  \AmirShahmoradi, Sunday 3:33 PM, September 19, 2021, Dallas, TX
      23             : 
      24             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      25             : 
      26             : #if     CK_ENABLED
      27             : #define TYPE_KIND complex(TKC)
      28             : #elif   RK_ENABLED
      29             : #define TYPE_KIND real(TKC)
      30             : #else
      31             : #error  "Unrecognized interface."
      32             : #endif
      33             :         integer(IK) :: i
      34             :         real(TKC), parameter :: LOGTINY = log(tiny(0._TKC))
      35             : #if     Sel_ENABLED
      36             :         TYPE_KIND :: logRatio
      37             : #elif   !Seq_ENABLED
      38             : #error  "Unrecognized interface."
      39             : #endif
      40             : #if     Def_ENABLED
      41             :         TYPE_KIND :: maxArray
      42     1960000 :         maxArray = maxval(real(array, TKC))
      43             : #elif   Max_ENABLED
      44      312693 :         CHECK_ASSERTION(__LINE__, real(maxArray, TKC) == maxval(real(array, TKC)), SK_"@getLogSumExp(): The condition `real(maxArray) == maxval(real(array))` must hold. real(maxArray), maxval(real(array)) = "//getStr([real(maxArray, TKC), maxval(real(array, TKC))]))
      45             : #else
      46             : #error  "Unrecognized interface."
      47             : #endif
      48        4297 :         logSumExp = 0._TKC
      49             : #if     Seq_ENABLED
      50     1588203 :         do i = 1_IK, size(array, kind = IK)
      51     1588203 :             logSumExp = logSumExp + exp(array(i) - maxArray)
      52             :         end do
      53      515429 :         logSumExp = maxArray + log(logSumExp)
      54             : #elif   Sel_ENABLED
      55           0 :         do i = 1_IK, size(array, kind = IK)
      56           0 :             logRatio = array(i) - maxArray
      57           0 :             if (real(logRatio, kind = kind(logRatio)) > LOGTINY) logSumExp = logSumExp + exp(logRatio) ! fpp
      58             :         end do
      59           0 :         logSumExp = maxArray + log(logSumExp)
      60             : #else
      61             : #error  "Unrecognized interface."
      62             : #endif
      63             : #undef  TYPE_KIND