The ParaMonte Documentation Website
Current view: top level - kernel - Matrix_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: Serial Kernel - Code Coverage Report Lines: 276 291 94.8 %
Date: 2021-01-08 13:03:42 Functions: 18 19 94.7 %
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             : !>  \brief This module contains mathematical procedures.
      44             : !>  \author Amir Shahmoradi
      45             : 
      46             : module Matrix_mod
      47             : 
      48             :     implicit none
      49             : 
      50             :     character(*), parameter :: MODULE_NAME = "@Matrix_mod"
      51             : 
      52             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      53             : 
      54             : contains
      55             : 
      56             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      57             : 
      58             :     !> \brief
      59             :     !> Return the an eye matrix of rank (nd, md).
      60             :     !>
      61             :     !> \param[in]   n       :   The number of rows of the eye matrix.
      62             :     !> \param[in]   m       :   The number of columns of the eye matrix.
      63             :     !> \param[in]   diag    :   The value to appear on the diagonal elements of the output eye matrix (**optional**, default = `1.`).
      64             :     !>
      65             :     !> \return
      66             :     !> `eye` : A matrix of rank n * m whose diagonal elements are set to the input `diag` value or, if not provided, to `1.`.
      67        2432 :     pure function getEye(n,m,diag) result(eye)
      68             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      69             :         !DEC$ ATTRIBUTES DLLEXPORT :: getEye
      70             : #endif
      71             :         use Constants_mod, only: RK, IK
      72             :         implicit none
      73             :         integer(IK), intent(in)             :: n, m
      74             :         real(RK)   , intent(in), optional   :: diag
      75             :         real(RK)                            :: eye(n,m)
      76         356 :         real(RK)                            :: diagonal
      77             :         integer(IK)                         :: i, j
      78         356 :         diagonal = 1._RK
      79         356 :         if (present(diag)) diagonal = diag
      80         356 :         do concurrent(i=1:n, j=1:m)
      81        1720 :             if (i==j) then
      82         517 :                 eye(i,j) = diagonal
      83             :             else
      84         329 :                 eye(i,j) = 0._RK
      85             :             end if
      86             :         end do
      87         356 :     end function getEye
      88             : 
      89             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      90             : 
      91             :     !> \brief
      92             :     !> Return the Cholesky factorization of the input positive-definite matrix.
      93             :     !>
      94             :     !> \param[in]       nd          :   The size of the input square matrix - `nd` by `nd`.
      95             :     !> \param[in,out]   PosDefMat   :   The input square matrix.
      96             :     !> \param[out]      Diagonal    :   The Diagonal elements of the Cholesky factorization.
      97             :     !>
      98             :     !> \remark
      99             :     !> If `nd = 1`, `PosDefMat` will not be touched, only `sqrt(PosDefMat)` will be output to `Diagonal`.
     100             :     !>
     101             :     !> \warning
     102             :     !> If Cholesky factorization fails, `Diagonal(1)` will be set to `-1` to indicate error on return.
     103             :     !>
     104             :     !> \details
     105             :     !> Returns in the lower triangle of `PosDefMat`, the Cholesky factorization L of \f$\texttt{PosDefMat} = L.L^T\f$.
     106             :     !> On input, the upper upper triangle of `PosDefMat` should be given, which remains intact on output.
     107       24107 :     pure subroutine getCholeskyFactor(nd,PosDefMat,Diagonal)
     108             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     109             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCholeskyFactor
     110             : #endif
     111         356 :         use Constants_mod, only: RK, IK
     112             :         implicit none
     113             :         integer(IK), intent(in)    :: nd
     114             :         real(RK)   , intent(inout) :: PosDefMat(nd,nd) ! Upper triangle + diagonal is input matrix, lower is output.
     115             :         real(RK)   , intent(out)   :: Diagonal(nd)
     116       24107 :         real(RK)                   :: summ
     117             :         integer(IK)                :: i
     118       61065 :         do i=1,nd
     119       49827 :             summ = PosDefMat(i,i) - dot_product(PosDefMat(i,1:i-1),PosDefMat(i,1:i-1))
     120       36962 :             if (summ <= 0._RK) then
     121           4 :                 Diagonal(1) = -1._RK
     122           4 :                 return
     123             :             end if
     124       36958 :             Diagonal(i) = sqrt(summ)
     125       73926 :             PosDefMat(i+1:nd,i) = ( PosDefMat(i,i+1:nd) - matmul(PosDefMat(i+1:nd,1:i-1),PosDefMat(i,1:i-1)) ) / Diagonal(i)
     126             :         end do
     127       24107 :     end subroutine getCholeskyFactor
     128             : 
     129             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     130             : 
     131             :     !> \brief
     132             :     !> Solve the linear equation system of the form: \f$ \texttt{PosDefMat} \times \texttt{InputSolution} = \texttt{Intercept} \f$
     133             :     !>
     134             :     !> \param[in]       nd              :   The size of the input square matrix - `nd` by `nd`.
     135             :     !> \param[in]       PosDefMat       :   The input square matrix.
     136             :     !> \param[in]       Diagonal        :   The Diagonal elements of the Cholesky factorization.
     137             :     !> \param[in]       Intercept       :   The intercept.
     138             :     !> \param[in,out]   InputSolution   :   The input right-hand-side which becomes the solution on return.
     139             :     !>
     140             :     !> \remark
     141             :     !> `PosDefMat` and `Diagonal` are the output of subroutine [getCholeskyFactor](@ref getcholeskyfactor)
     142             :     !> (i.g., only the lower triangle of `PosDefMat` is used).
     143           0 :     subroutine solveLinearPosDefSystem(nd,PosDefMat,Diagonal,Intercept,InputSolution)
     144             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     145             :         !DEC$ ATTRIBUTES DLLEXPORT :: solveLinearPosDefSystem
     146             : #endif
     147       24107 :         use Constants_mod, only: RK, IK
     148             :         implicit none
     149             :         integer(IK), intent(in)    :: nd
     150             :         real(RK)   , intent(in)    :: PosDefMat(nd,nd),Diagonal(nd),Intercept(nd)
     151             :         real(RK)   , intent(inout) :: InputSolution(nd)
     152             :         integer(IK)                :: i
     153           0 :         do i=1,nd
     154           0 :             InputSolution(i) = ( Intercept(i) - dot_product(PosDefMat(i,1:i-1),InputSolution(1:i-1)) ) / Diagonal(i)
     155             :         end do
     156           0 :         do i = nd,1,-1
     157           0 :             InputSolution(i) = ( InputSolution(i) - dot_product(PosDefMat(i+1:nd,i),InputSolution(i+1:nd)) ) / Diagonal(i)
     158             :         end do
     159           0 :     end subroutine solveLinearPosDefSystem
     160             : 
     161             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     162             : 
     163             :     !> \brief
     164             :     !> Return the inverse matrix of a symmetric-positive-definite input matrix, which is given in the upper triangle of `MatInvMat`.
     165             :     !> On output `MatInvMat` is completely overwritten by in the inverse of the matrix. Also, return the square root of determinant
     166             :     !> of the inverse matrix.
     167             :     !>
     168             :     !> \param[in]       nd                  :   The size of the input square matrix - `nd` by `nd`.
     169             :     !> \param[in,out]   MatInvMat           :   The input square matrix. On input, the upper half must be covariance matrix.
     170             :     !>                                          On output, it is completely overwritten by the inverse matrix.
     171             :     !> \param[out]      sqrtDetInvPosDefMat :   The square root of the determinant of the inverse matrix.
     172             :     !>
     173             :     !> \warning
     174             :     !> Do not call this routine when `nd = 1`. There is no need to call this routine when `nd = 1`.
     175             :     !>
     176             :     !> \warning
     177             :     !> If the input matrix is not positive-definite, the output `sqrtDetInvPosDefMat`
     178             :     !> will be set to `-1` on return to indicate the occurrence of an error.
     179             :     !>
     180             :     !> \author
     181             :     !> Amir Shahmoradi, Apr 21, 2017, 1:54 AM, ICES, UT Austin
     182           5 :     pure subroutine getInvPosDefMatSqrtDet(nd,MatInvMat,sqrtDetInvPosDefMat)
     183             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     184             :         !DEC$ ATTRIBUTES DLLEXPORT :: getInvPosDefMatSqrtDet
     185             : #endif
     186             :         use Constants_mod, only: RK, IK ! LCOV_EXCL_LINE
     187             :         implicit none
     188             :         integer(IK), intent(in)    :: nd
     189             :         real(RK)   , intent(inout) :: MatInvMat(nd,nd)           ! input: upper half is covariance matrix, output: inverse matrix
     190             :         real(RK)   , intent(out)   :: sqrtDetInvPosDefMat        !
     191          55 :         real(RK)                   :: CholeskyLower(nd,nd)       ! Cholesky factor
     192          23 :         real(RK)                   :: Diagonal(nd)               ! diagonal terms of the inverse matrix
     193           5 :         real(RK)                   :: summ
     194             :         integer(IK)                :: i,j,k
     195           5 :         if (nd==1_IK) then
     196           1 :             MatInvMat(1,1) = 1._RK / MatInvMat(1,1)
     197           1 :             sqrtDetInvPosDefMat = MatInvMat(1,1)
     198           1 :             return
     199             :         end if
     200          16 :         do j=1,nd
     201          40 :             do i=1,j
     202          36 :                 CholeskyLower(i,j) = MatInvMat(i,j)
     203             :             end do
     204             :         end do
     205           4 :         call getCholeskyFactor(nd,CholeskyLower,Diagonal)
     206           4 :         if (Diagonal(1)<0._RK) then
     207           1 :             sqrtDetInvPosDefMat = -1._RK
     208           1 :             return
     209             :         end if
     210          12 :         sqrtDetInvPosDefMat = 1._RK / product(Diagonal)
     211          12 :         do i = 1,nd
     212           9 :             CholeskyLower(i,i) = 1._RK / Diagonal(i)
     213          21 :             do j = i+1,nd
     214           9 :                 summ = 0._RK
     215          21 :                 do k = i,j-1
     216          21 :                     summ = summ - CholeskyLower(j,k) * CholeskyLower(k,i)
     217             :                 end do
     218          18 :                 CholeskyLower(j,i) = summ / Diagonal(j)
     219             :             end do
     220             :         end do
     221          12 :         do i = 1,nd
     222          27 :             do j = i,nd
     223          57 :                 MatInvMat(j,i) = dot_product(CholeskyLower(j:nd,j),CholeskyLower(j:nd,i))
     224             :             end do
     225          48 :             MatInvMat(i,i:nd) = MatInvMat(i:nd,i)
     226             :         end do
     227           5 :   end subroutine getInvPosDefMatSqrtDet
     228             : 
     229             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     230             : 
     231             :     !> \brief
     232             :     !> Return the inverse matrix of a symmetric-positive-definite matrix, whose Cholesky Lower triangle is given in the lower part
     233             :     !> of `CholeskyLower` and and its diagonal elements in `Diagonal`.
     234             :     !>
     235             :     !> \param[in]       nd              :   The size of the input square matrix - `nd` by `nd`.
     236             :     !> \param[in]       CholeskyLower   :   The Cholesky factorization of the matrix.
     237             :     !> \param[in]       Diagonal        :   The diagonal elements of the Cholesky factorization of the matrix.
     238             :     !>
     239             :     !> \return
     240             :     !> `InvMatFromCholFac` : The full inverse matrix.
     241             :     !>
     242             :     !> \warning
     243             :     !> Do not call this routine when `nd = 1`. For `nd = 1`: `InvMatFromCholFac = 1._RK / Diagonal(1)^2`.
     244             :     !>
     245             :     !> \author
     246             :     !> Amir Shahmoradi, Apr 21, 2017, 1:54 AM, ICES, UT Austin
     247     1059270 :     pure function getInvMatFromCholFac(nd,CholeskyLower,Diagonal) result(InvMatFromCholFac)
     248             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     249             :         !DEC$ ATTRIBUTES DLLEXPORT :: getInvMatFromCholFac
     250             : #endif
     251           5 :         use Constants_mod, only: RK, IK
     252             :         implicit none
     253             :         integer(IK), intent(in) :: nd
     254             :         real(RK)   , intent(in) :: CholeskyLower(nd,nd),Diagonal(nd)
     255             :         real(RK)                :: InvMatFromCholFac(nd,nd)
     256      123498 :         real(RK)                :: summ
     257             :         integer(IK)             :: i,j,k
     258      123498 :         if (nd==1_IK) then
     259       13055 :             InvMatFromCholFac(1,1) = 1._RK / Diagonal(1)**2
     260       13055 :             return
     261             :         end if
     262      773107 :         InvMatFromCholFac = 0._RK
     263      220887 :         do j=1,nd-1
     264      331332 :             do i=j+1,nd
     265      220889 :                 InvMatFromCholFac(i,j) = CholeskyLower(i,j)
     266             :             end do
     267             :         end do
     268      331330 :         do i = 1,nd
     269      220887 :             InvMatFromCholFac(i,i) = 1._RK / Diagonal(i)
     270      441775 :             do j = i+1,nd
     271      110445 :                 summ = 0._RK
     272      220891 :                 do k = i,j-1
     273      220891 :                     summ = summ - InvMatFromCholFac(j,k) * InvMatFromCholFac(k,i)
     274             :                 end do
     275      331332 :                 InvMatFromCholFac(j,i) = summ / Diagonal(j)
     276             :             end do
     277             :         end do
     278      331330 :         do i = 1,nd
     279      662662 :             do j = i,nd
     280      773110 :                 InvMatFromCholFac(j,i) = dot_product(InvMatFromCholFac(j:nd,j),InvMatFromCholFac(j:nd,i))
     281      552219 :                 InvMatFromCholFac(i,j) = InvMatFromCholFac(j,i)
     282             :             end do
     283             :             !InvMatFromCholFac(i,i:nd) = InvMatFromCholFac(i:nd,i)
     284             :         end do
     285      246996 :     end function getInvMatFromCholFac
     286             : 
     287             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     288             : 
     289             :     !> \brief
     290             :     !> Return the inverse matrix of an input symmetric-positive-definite matrix `PosDefMat`.
     291             :     !>
     292             :     !> \param[in]   nd          :   The size of the input square matrix - `nd` by `nd`.
     293             :     !> \param[in]   PosDefMat   :   The input symmetric-positive-definite matrix.
     294             :     !>
     295             :     !> \return
     296             :     !> `InvPosDefMat` : The full inverse matrix.
     297             :     !>
     298             :     !> \warning
     299             :     !> If the input matrix is not positive definite, the function will return with `InvPosDefMat(1,1) = -1._RK` to indicate
     300             :     !> the occurrence of an error. This is done to preserve the purity of the function.
     301             :     !>
     302             :     !> \warning
     303             :     !> Do not call this routine when `nd = 1`. For `nd = 1`: `InvPosDefMat = 1._RK / InvPosDefMat`.
     304             :     !>
     305             :     !> \remark
     306             :     !> On input, only the upper half of `PosDefMat` needs to be given.
     307             :     !>
     308             :     !> \remark
     309             :     !> This routine has the same functionality as the subroutine [getInvPosDefMatSqrtDet()](@ref getinvposdefmatsqrtdet),
     310             :     !> but returns the result as a function output.
     311             :     !> According to the timing tests, `compare_InvMatRoutines_1()`, the function version appears to be 15-30% faster than
     312             :     !> the subroutine version [getInvPosDefMatSqrtDet()](@ref getinvposdefmatsqrtdet).
     313             :     !>
     314             :     !> \author
     315             :     !> Amir Shahmoradi, Apr 21, 2017, 1:54 AM, ICES, UT Austin
     316          30 :     pure function getInvPosDefMat(nd,PosDefMat) result(InvPosDefMat)
     317             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     318             :         !DEC$ ATTRIBUTES DLLEXPORT :: getInvPosDefMat
     319             : #endif
     320      123498 :         use Constants_mod, only: RK, IK
     321             :         implicit none
     322             :         integer(IK), intent(in) :: nd
     323             :         real(RK)   , intent(in) :: PosDefMat(nd,nd)
     324          26 :         real(RK)                :: InvPosDefMat(nd,nd), CholeskyLower(nd,nd)
     325          10 :         real(RK)                :: Diagonal(nd)
     326           2 :         real(RK)                :: summ
     327             :         integer(IK)             :: i,j,k
     328             :        !if (nd==1) then
     329             :        !    InvPosDefMat(1,1) = 1._RK / sqrt(PosDefMat(1,1))
     330             :        !    return
     331             :        !end if
     332           8 :         do j=1,nd
     333          20 :             do i=1,j
     334          18 :                 CholeskyLower(i,j) = PosDefMat(i,j)
     335             :             end do
     336             :         end do
     337           2 :         call getCholeskyFactor(nd,CholeskyLower,Diagonal)
     338           2 :         if (Diagonal(1)<0._RK) then
     339          13 :             InvPosDefMat = -1._RK   ! error occurred: getCholeskyFactor() failed in getInvPosDefMat()
     340           1 :             return
     341             :         end if
     342           4 :         do i = 1,nd
     343           3 :             CholeskyLower(i,i) = 1._RK / Diagonal(i)
     344           7 :             do j = i+1,nd
     345           3 :                 summ = 0._RK
     346           7 :                 do k = i,j-1
     347           7 :                     summ = summ - CholeskyLower(j,k) * CholeskyLower(k,i)
     348             :                 end do
     349           6 :                 CholeskyLower(j,i) = summ / Diagonal(j)
     350             :             end do
     351             :         end do
     352           4 :         do i = 1,nd
     353          10 :             do j = i,nd
     354          16 :                 InvPosDefMat(j,i) = dot_product(CholeskyLower(j:nd,j),CholeskyLower(j:nd,i))
     355           9 :                 InvPosDefMat(i,j) = InvPosDefMat(j,i)
     356             :             end do
     357             :         end do
     358           4 :   end function getInvPosDefMat
     359             : 
     360             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     361             : 
     362             :     !> \brief
     363             :     !> Return the inverse matrix `InverseMatrix` of a `nd*nd` input matrix `MatrixLU`, and its determinant, using `LU` decomposition.
     364             :     !>
     365             :     !> \param[in]       nd              : The size of the square matrix and its `LU` decomposition - `nd` by `nd`.
     366             :     !> \param[inout]    MatrixLU        : The target symmetric-positive-definite matrix. On input it is the matrix, on output it is the LU decomposition.
     367             :     !> \param[out]      InverseMatrix   : The input symmetric-positive-definite matrix.
     368             :     !> \param[out]      detInvMat       : The determinant of the inverse of the symmetric-positive-definite matrix.
     369             :     !>
     370             :     !> \warning
     371             :     !> Do not call this routine when `nd = 1`. For `nd = 1`: `InverseMatrix = detInvMat = 1._RK / InverseMatrix`.
     372             :     !>
     373             :     !> \remark
     374             :     !> This routine has the same functionality as the function [getInvMat()](@ref getinvmat),
     375             :     !> but returns the result as a function output.
     376             :     !> According to the timing tests, `compare_InvMatRoutines_1()`, the function version of this code [getInvMat()](@ref getinvmat)
     377             :     !> appears to be 15-30% faster than this subroutine version.
     378             :     !>
     379             :     !> \author
     380             :     !> Amir Shahmoradi, Oct 18, 2009, 1:54 AM, MTU
     381           1 :     subroutine getInvMatDet(nd,MatrixLU,InverseMatrix,detInvMat)
     382             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     383             :         !DEC$ ATTRIBUTES DLLEXPORT :: getInvMatDet
     384             : #endif
     385           2 :         use Constants_mod, only: RK, IK
     386             :         implicit none
     387             :         integer(IK), intent(in)    :: nd
     388             :         real(RK)   , intent(inout) :: MatrixLU(nd,nd) ! On input it is the matrix, on output it is the LU decomposition.
     389             :         real(RK)   , intent(out)   :: InverseMatrix(nd,nd)
     390             :         real(RK)   , intent(out)   :: detInvMat
     391             :         integer(IK)                :: i,j,Permutation(nd) ! LCOV_EXCL_LINE
     392           4 :         do i = 1,nd
     393          12 :             do j = 1,nd
     394          12 :                 InverseMatrix(i,j) = 0._RK
     395             :             end do
     396           4 :             InverseMatrix(i,i) = 1._RK
     397             :         end do
     398           1 :         call getLU(nd,MatrixLU,Permutation,detInvMat)
     399           4 :         do j = 1,nd
     400           3 :             detInvMat = detInvMat * MatrixLU(j,j)
     401           4 :             call solveLinearSystem(nd,MatrixLU,Permutation,InverseMatrix(1:nd,j))
     402             :         end do
     403           1 :         detInvMat = 1._RK/detInvMat
     404           1 :   end subroutine getInvMatDet
     405             : 
     406             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     407             : 
     408             :     !> \brief
     409             :     !> Return the inverse matrix `InverseMatrix` of a `nd*nd` input matrix `Matrix`, and its determinant, using `LU` decomposition.
     410             :     !>
     411             :     !> \param[in]       nd      :   The size of the square matrix and its `LU` decomposition - `nd` by `nd`.
     412             :     !> \param[inout]    Matrix  :   The target symmetric-positive-definite matrix.
     413             :     !>
     414             :     !> \return
     415             :     !> `InverseMatrix` : The full inverse matrix.
     416             :     !>
     417             :     !> \remark
     418             :     !> This routine has the same functionality as the subroutine [getInvMatDet()](@ref getinvmatdet),
     419             :     !> but returns the result as a function output.
     420             :     !> According to the timing tests, `compare_InvMatRoutines_1()`, this function version appears
     421             :     !> to be 15-30% faster than the subroutine equivalent [getInvMatDet()](@ref getinvmatdet).
     422             :     !>
     423             :     !> \author
     424             :     !> Amir Shahmoradi, Apr 8, 2017, 1:54 PM, MTU
     425          15 :     function getInvMat(nd,Matrix) result(InverseMatrix)
     426             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     427             :         !DEC$ ATTRIBUTES DLLEXPORT :: getInvMat
     428             : #endif
     429           1 :         use Constants_mod, only: RK, IK
     430             :         implicit none
     431             :         integer(IK), intent(in) :: nd
     432             :         real(RK)   , intent(in) :: Matrix(nd,nd)
     433          14 :         real(RK)                :: InverseMatrix(nd,nd),LU(nd,nd)
     434           1 :         integer(IK)             :: i,j,Permutation(nd)
     435           1 :         real(RK)                :: parity
     436           4 :         do i = 1,nd
     437          12 :             do j = 1,nd
     438          12 :                 InverseMatrix(i,j) = 0._RK
     439             :             end do
     440           4 :             InverseMatrix(i,i) = 1._RK
     441             :         end do
     442          13 :         LU = Matrix
     443           1 :         call getLU(nd,LU,Permutation,parity)
     444           4 :         do j = 1,nd
     445           4 :             call solveLinearSystem(nd,LU,Permutation,InverseMatrix(1:nd,j))
     446             :         end do
     447           1 :     end function getInvMat
     448             : 
     449             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     450             : 
     451             :     !> \brief
     452             :     !> Solve the set of `nd` linear equations `Matrix * X = InputSolution`.
     453             :     !>
     454             :     !> \param[in]       nd              :   The size of the square matrix and its `LU` decomposition - `nd` by `nd`.
     455             :     !> \param[in]       LU              :   The LU factorization of the matrix, determined by the routine [getLU()](@ref getlu).
     456             :     !> \param[in]       Permutation     :   The permutation vector returned by [getLU()](@ref getlu).
     457             :     !> \param[inout]    InputSolution   :   The right-hand-side vector of length `nd`.
     458             :     !>                                      On output, it is completely overwritten by the solution to the system.
     459             :     !>
     460             :     !> \remark
     461             :     !> The input `nd`, `LU`, and `Permutation` parameters are not modified by this routine and,
     462             :     !> can be left in place for successive calls with different right-hand sides `InputSolution`.
     463             :     !>
     464             :     !> \remark
     465             :     !> This routine takes into account the possibility that `InputSolution` will begin with many zero elements,
     466             :     !> and is therefore efficient for matrix inversion.
     467             :     !>
     468             :     !> \author
     469             :     !> Amir Shahmoradi, Apr 8, 2017, 1:54 PM, MTU
     470           6 :     subroutine solveLinearSystem(nd, LU, Permutation, InputSolution)
     471             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     472             :         !DEC$ ATTRIBUTES DLLEXPORT :: solveLinearSystem
     473             : #endif
     474           1 :         use Constants_mod, only: RK, IK
     475             :         integer(IK), intent(in)    :: nd
     476             :         integer(IK), intent(in)    :: Permutation(nd)
     477             :         real(RK)   , intent(in)    :: LU(nd,nd)
     478             :         real(RK)   , intent(inout) :: InputSolution(nd)
     479             :         integer(IK)                :: i,ii
     480           6 :         real(RK)                   :: summ
     481           6 :         ii = 0
     482          24 :         do i=1,nd
     483          18 :             summ = InputSolution(Permutation(i))
     484          18 :             InputSolution(Permutation(i)) = InputSolution(i)
     485          18 :             if (ii /= 0) then
     486          14 :                 summ = summ - dot_product( LU(i,ii:i-1) , InputSolution(ii:i-1) )
     487          12 :             else if (summ /= 0._RK) then
     488           6 :                 ii = i
     489             :             end if
     490          24 :             InputSolution(i) = summ
     491             :         end do
     492          24 :         do i = nd,1,-1
     493          42 :             InputSolution(i) = ( InputSolution(i) - dot_product( LU(i,i+1:nd) , InputSolution(i+1:nd)) ) / LU(i,i)
     494             :         end do
     495           6 :     end subroutine solveLinearSystem
     496             : 
     497             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     498             : 
     499             :     !> \brief
     500             :     !> Return the LU decomposition of the input matrix `MatrixLU(nd,nd)`.
     501             :     !>
     502             :     !> \param[in]       nd          :   The size of the square matrix and its `LU` decomposition - `nd` by `nd`.
     503             :     !> \param[inout]    MatrixLU    :   The input matrix. On output, it is completely overwritten by its LU decomposition.
     504             :     !> \param[out]      Permutation :   An output vector of length `nd` that records the row permutation effected by the partial pivoting.
     505             :     !> \param[out]      Parity      :   An output real as `+-1` depending on whether the number of row interchanges was even or odd, respectively.
     506             :     !>
     507             :     !> \remark
     508             :     !> This routine is used in combination with [solveLinearSystem](@ref solvelinearsystem) to solve linear equations or invert a matrix.
     509             :     !>
     510             :     !> \author
     511             :     !> Amir Shahmoradi, Apr 21, 2017, 1:43 PM, ICES, UT Austin
     512           3 :     subroutine getLU(nd,MatrixLU,Permutation,parity) ! ,errorOccurred)
     513             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     514             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLU
     515             : #endif
     516           6 :         use Constants_mod, only: RK, IK
     517             :         use, intrinsic :: iso_fortran_env, only: output_unit
     518             :         integer(IK) , intent(in)    :: nd
     519             :         integer(IK) , intent(out)   :: Permutation(nd)
     520             :         real(RK)    , intent(inout) :: MatrixLU(nd,nd)
     521             :         real(RK)    , intent(out)   :: parity
     522             :         !logical     , intent(out)   :: errorOccurred
     523             :         real(RK)    , parameter     :: TINY = 1.e-20_RK
     524          12 :         real(RK)                    :: aamax,dum,summ,vv(nd)
     525             :         integer(IK)                 :: i,imax,j,k
     526             :         !errorOccurred = .false.
     527           3 :         parity = 1._RK
     528          12 :         do i = 1,nd
     529           9 :             aamax = 0._RK
     530          36 :             do j=1,nd
     531          36 :                 if ( abs(MatrixLU(i,j)) > aamax ) aamax = abs( MatrixLU(i,j) )
     532             :             end do
     533           9 :             if (aamax == 0._RK) then
     534             :             ! LCOV_EXCL_START
     535             :                 write(output_unit,"(A)") "Statistics@getLU() failed. Singular matrix detected."
     536             :                 error stop
     537             :                 !errorOccurred = .true.
     538             :                 !return
     539             :             ! LCOV_EXCL_STOP
     540             :             end if
     541          12 :             vv(i) = 1._RK/aamax
     542             :         end do
     543          12 :         do j=1,nd
     544          18 :             do i=1,j-1
     545           9 :                 summ = MatrixLU(i,j)
     546          12 :                 do k=1,i-1
     547          12 :                     summ = summ - MatrixLU(i,k) * MatrixLU(k,j)
     548             :                 end do
     549          18 :                 MatrixLU(i,j) = summ
     550             :             end do
     551           9 :             aamax = 0._RK
     552          27 :             do i = j, nd
     553          18 :                 summ = MatrixLU(i,j)
     554          30 :                 do k = 1, j-1
     555          30 :                     summ = summ - MatrixLU(i,k) * MatrixLU(k,j)
     556             :                 end do
     557          18 :                 MatrixLU(i,j) = summ
     558          18 :                 dum = vv(i) * abs(summ)
     559          27 :                 if (dum >= aamax) then
     560           9 :                     imax = i
     561           9 :                     aamax = dum
     562             :                 end if
     563             :             end do
     564           9 :             if (j /= imax)then
     565           0 :                 do k=1,nd
     566           0 :                     dum = MatrixLU(imax,k)
     567           0 :                     MatrixLU(imax,k) = MatrixLU(j,k)
     568           0 :                     MatrixLU(j,k) = dum
     569             :                 end do
     570           0 :                 parity = - parity
     571           0 :                 vv(imax) = vv(j)
     572             :             end if
     573           9 :             Permutation(j) = imax
     574           9 :             if (MatrixLU(j,j) == 0._RK) MatrixLU(j,j) = TINY
     575          12 :             if (j /= nd) then
     576           6 :                 dum = 1._RK / MatrixLU(j,j)
     577          15 :                 do i = j+1,nd
     578          15 :                     MatrixLU(i,j) = MatrixLU(i,j) * dum
     579             :                 end do
     580             :             endif
     581             :         end do
     582           3 :   end subroutine getLU
     583             : 
     584             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     585             : 
     586             :     !> \brief
     587             :     !> Return the product of two matrices.
     588             :     !> The product of two matrices is defined by `c(i,j) = a(i,1)*b(1,j) + a(i,2)*b(2,j) + ... + a(i,n)*b(n,j)`.
     589             :     !>
     590             :     !> \remark
     591             :     !> There is not reason to use this routine as the Fortran intrinsic already provides an optimized implementation of matrix
     592             :     !> multiplication via `matmul()`. It is present here only for archival and legacy reasons.
     593             :     !>
     594             :     !> \author
     595             :     !> Amir Shahmoradi, Oct 20, 2009, 10:56 PM, MTU
     596           1 :     subroutine multiplyMatrix(A,rowsA,colsA,B,rowsB,colsB,C)
     597             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     598             :         !DEC$ ATTRIBUTES DLLEXPORT :: multiplyMatrix
     599             : #endif
     600           3 :         use Constants_mod, only: RK, IK
     601             :         use, intrinsic :: iso_fortran_env, only: output_unit
     602             :         !use Misc  , only: abortProgram
     603             :         implicit none
     604             :         integer(IK), intent(in)  :: rowsA, colsA, rowsB, colsB !Matrix Dimensions
     605             :         real(RK)   , intent(in)  :: A(rowsA,colsA)
     606             :         real(RK)   , intent(in)  :: B(rowsB,colsB)
     607             :         real(RK)   , intent(out) :: C(rowsA,colsB)
     608             :         integer(IK)              :: i,j,k,rowsC,colsC !Counters
     609           1 :         if (colsA == rowsB) then
     610           1 :             rowsC = rowsA
     611           1 :             colsC = colsB
     612             :         ! LCOV_EXCL_START
     613             :         else
     614             :             write(output_unit,"(A)") "Matrix@multiplyMatrix() failed. dimensions of matrices do not match."
     615             :             stop
     616             :             !call abortProgram( output_unit , 1 , 1 , 'Statistics@multiplyMatrix() failed. dimensions of matrices do not match.' )
     617             :         ! LCOV_EXCL_STOP
     618             :         end if
     619             :         ! Initialize product matrix to 0
     620           4 :         do i = 1, rowsC
     621          13 :             do j = 1, colsC
     622          12 :                 C(i,j) = 0._RK
     623             :             end do
     624             :         end do
     625             :         ! Find product as per above formula
     626           4 :         do i = 1, rowsA
     627          13 :             do j = 1, colsB
     628          39 :                 do k = 1, colsA
     629          36 :                     C(i,j) = C(i,j) + A(i,k)*B(k,j)
     630             :                 end do
     631             :             end do
     632             :         end do
     633           1 :     end subroutine multiplyMatrix
     634             : 
     635             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     636             : 
     637             :     !> \brief
     638             :     !> Return the determinant of a given `nd * nd` matrix via LU factorization.
     639             :     !>
     640             :     !> \param[in]       nd      :   The size of the square matrix - `nd` by `nd`.
     641             :     !> \param[in]       Matrix  :   The input matrix.
     642             :     !>
     643             :     !> \return
     644             :     !> `determinant` : The determinant of the matrix.
     645             :     !>
     646             :     !> \author
     647             :     !> Amir Shahmoradi, Oct 18, 2009, 4:10 PM, MTU
     648           1 :     function getDeterminant(nd,Matrix) result(determinant)
     649             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     650             :         !DEC$ ATTRIBUTES DLLEXPORT :: getDeterminant
     651             : #endif
     652           1 :         use Constants_mod, only: RK, IK
     653             :         implicit none
     654             :         integer(IK), intent(in) :: nd
     655             :         real(RK)   , intent(in) :: Matrix(nd,nd)
     656          13 :         real(RK)                :: DummyMat(nd,nd)
     657             :         real(RK)                :: determinant
     658           1 :         integer(IK)             :: Permutation(nd)
     659             :         integer(IK)             :: j
     660          13 :         DummyMat = Matrix
     661           1 :         call getLU(nd,DummyMat,Permutation,determinant) ! This returns determinant as +-1.
     662           4 :         do j=1,nd
     663           4 :             determinant = determinant * DummyMat(j,j)
     664             :         end do
     665           1 :     end function getDeterminant
     666             : 
     667             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     668             : 
     669             :     !> Return the square root of the determinant of a given positive-definite `nd * nd` matrix `PosDefMat` using Cholesky factorization.
     670             :     !>
     671             :     !> \param[in]       nd          :   The size of the square matrix - `nd` by `nd`.
     672             :     !> \param[in]       PosDefMat   :   The input matrix.
     673             :     !>
     674             :     !> \return
     675             :     !> `sqrtDetPosDefMat` : The square root of the determinant of the matrix.
     676             :     !>
     677             :     !> \warning
     678             :     !> If the input matrix is not positive definite, then `sqrtDetPosDefMat = -1._RK` upon return.
     679             :     !>
     680             :     !> \author
     681             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     682           2 :     function getSqrtDetPosDefMat(nd,PosDefMat) result(sqrtDetPosDefMat)
     683             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     684             :         !DEC$ ATTRIBUTES DLLEXPORT :: getSqrtDetPosDefMat
     685             : #endif
     686           1 :         use Constants_mod, only: RK, IK
     687             :         implicit none
     688             :         integer(IK), intent(in) :: nd
     689             :         real(RK)   , intent(in) :: PosDefMat(nd,nd)
     690          32 :         real(RK)                :: Diagonal(nd),DummyMat(nd,nd)
     691             :         real(RK)                :: sqrtDetPosDefMat
     692             :         integer(IK)             :: i,j
     693           8 :         do j=1,nd
     694          20 :             do i=1,j
     695          18 :                 DummyMat(i,j) = PosDefMat(i,j)
     696             :             end do
     697             :         end do
     698           2 :         call getCholeskyFactor(nd,DummyMat,Diagonal)
     699           2 :         if (Diagonal(1)<0._RK) then
     700           1 :             sqrtDetPosDefMat = -1._RK
     701           1 :             return
     702             :         end if
     703           4 :         sqrtDetPosDefMat = product(Diagonal)
     704           2 :     end function getSqrtDetPosDefMat
     705             : 
     706             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     707             : 
     708             :     !> \brief
     709             :     !> Return the natural logarithm of the square root of the determinant of a given positive-definite `nd * nd` matrix `PosDefMat`
     710             :     !> using Cholesky factorization.
     711             :     !>
     712             :     !> \param[in]       nd                  :   The size of the square matrix - `nd` by `nd`.
     713             :     !> \param[inout]    PosDefMat           :   The input matrix. On input, the upper triangle should be given, which remains intact on output.
     714             :     !> \param[out]      logSqrtDetPosDefMat :   The natural logarithm of the square root of the determinant of `PosDefMat`.
     715             :     !> \param[out]      failed              :   A logical value. If `.true.`, the determinant computation has failed.
     716             :     !>
     717             :     !> \author
     718             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     719       60512 :     pure subroutine getLogSqrtDetPosDefMat(nd,PosDefMat,logSqrtDetPosDefMat,failed)
     720             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     721             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogSqrtDetPosDefMat
     722             : #endif
     723           2 :         use Constants_mod, only: RK, IK
     724             :         implicit none
     725             :         integer(IK), intent(in)     :: nd
     726             :         real(RK)   , intent(inout)  :: PosDefMat(nd,nd) ! Upper triangle + diagonal is input matrix, lower is output.
     727             :         real(RK)   , intent(out)    :: logSqrtDetPosDefMat
     728             :         logical    , intent(out)    :: failed
     729      170332 :         real(RK)                    :: Diagonal(nd)
     730       60512 :         real(RK)                    :: summ
     731             :         integer(IK)                 :: i
     732       60512 :         failed = .false.
     733      170329 :         do i=1,nd
     734      159125 :             summ = PosDefMat(i,i) - dot_product(PosDefMat(i,1:i-1),PosDefMat(i,1:i-1))
     735      109818 :             if (summ <= 0._RK) then
     736           1 :                 failed = .true.
     737           1 :                 return
     738             :             end if
     739      109817 :             Diagonal(i) = sqrt(summ)
     740      219635 :             PosDefMat(i+1:nd,i) = ( PosDefMat(i,i+1:nd) - matmul(PosDefMat(i+1:nd,1:i-1),PosDefMat(i,1:i-1)) ) / Diagonal(i)
     741             :         end do
     742      170328 :         logSqrtDetPosDefMat = sum(log(Diagonal))
     743       60512 :     end subroutine getLogSqrtDetPosDefMat
     744             : 
     745             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     746             : 
     747             :     !> \brief
     748             :     !> Return `.false.` value for `isPosDef`, if the Cholesky decomposition of the input matrix fails
     749             :     !> (i.e. matrix is not positive definite), otherwise return `.true.`.
     750             :     !>
     751             :     !> \param[in]       nd      :   The size of the square matrix - `nd` by `nd`.
     752             :     !> \param[inout]    Matrix  :   The input matrix. Note that only the upper triangle of the matrix is used.
     753             :     !>
     754             :     !> \return
     755             :     !> `isPosDef` : A logical value indicating whether the input matrix is positive-definite.
     756             :     !>
     757             :     !> \warning
     758             :     !> Do not call this routine when `nd = 1`. In such a case, if `Matrix(1,1) <= 0`, then the matrix is not positive-definite.
     759             :     !>
     760             :     !> \author
     761             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     762         992 :     pure logical function isPosDef(nd,Matrix)
     763             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     764             :         !DEC$ ATTRIBUTES DLLEXPORT :: isPosDef
     765             : #endif
     766       60512 :         use Constants_mod, only: RK, IK
     767             :         implicit none
     768             :         integer(IK), intent(in) :: nd
     769             :         real(RK)   , intent(in) :: Matrix(nd,nd)
     770        8014 :         real(RK)                :: Array(nd,nd),p(nd)
     771         992 :         real(RK)                :: dummySum
     772             :         integer(IK)             :: i,j,k
     773         992 :         isPosDef = .true.
     774        2585 :         do i = 1, nd
     775        4787 :             do j = i, nd
     776        2209 :                 dummySum = Matrix(i,j)
     777        2819 :                 do k = i-1,1,-1
     778        2819 :                     dummySum = dummySum - Array(i,k) * Array(j,k)
     779             :                 end do
     780        3802 :                 if(i==j)then
     781        1600 :                     if(dummySum<=0._RK) then
     782           7 :                         isPosDef = .false.
     783           7 :                         return
     784             :                     end if
     785        1593 :                     p(i) = sqrt(dummySum)
     786             :                 else
     787         609 :                     Array(j,i)=dummySum/p(i)
     788             :                 endif
     789             :             end do
     790             :         end do
     791         992 :     end function isPosDef
     792             : 
     793             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     794             : 
     795             :     !> \brief
     796             :     !> Return the outer product of the two input matrices.
     797             :     !>
     798             :     !> \param[in]       Vector1 : The first input vector.
     799             :     !> \param[in]       Vector2 : The second input vector.
     800             :     !>
     801             :     !> \return
     802             :     !> `OuterProd` : A matrix of size `( size(Vector1), size(Vector2) )`.
     803             :     !>
     804             :     !> \author
     805             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     806          11 :     pure function getOuterProd(Vector1,Vector2) result(OuterProd)
     807             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     808             :         !DEC$ ATTRIBUTES DLLEXPORT :: getOuterProd
     809             : #endif
     810         992 :         use Constants_mod, only: RK, IK
     811             :         real(RK), intent(in) :: Vector1(:),Vector2(:)
     812           2 :         real(RK)             :: OuterProd(size(Vector1), size(Vector2))
     813           9 :         OuterProd = spread(Vector1,dim=2,ncopies=size(Vector2)) * spread(Vector2,dim=1,ncopies=size(Vector1))
     814           2 :     end function getOuterProd
     815             : 
     816             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     817             : 
     818             :     !> \brief
     819             :     !> Return an ordered matrix of the input matrix by rearranging the columns corresponding to `ColIndx` to into the corresponding
     820             :     !> columns in `ColIndxMap`, while keeping the rest of the matrix structure intact,
     821             :     !> such that if it is positive-definite, it remains positive-definite.
     822             :     !> For example,
     823             :     !> ```
     824             :     !> 11.0    12.0    13.0    14.0    15.0
     825             :     !>  0.0    22.0    23.0    24.0    25.0
     826             :     !>  0.0     0.0    33.0    34.0    35.0
     827             :     !>  0.0     0.0     0.0    44.0    45.0
     828             :     !>  0.0 *******     0.0     0.0    55.0
     829             :     !> ```
     830             :     !> goes to the following by a `ColIndx=[4]`, `ColIndxMap=[5]`,
     831             :     !> ```
     832             :     !> 11.0    12.0    13.0    15.0    14.0
     833             :     !>  0.0    22.0    23.0    25.0    24.0
     834             :     !>  0.0     0.0    33.0    35.0    34.0
     835             :     !>  0.0     0.0     0.0    55.0    45.0
     836             :     !>  0.0     0.0     0.0     0.0    44.0
     837             :     !> ```
     838             :     !>
     839             :     !> \param[in]       rank            :   The number of columns (or rows) of the input square matrix `PosDefMatUpper`.
     840             :     !> \param[in]       PosDefMatUpper  :   The input matrix.
     841             :     !> \param[in]       lenColIndx      :   The length of the input `ColIndx` vector.
     842             :     !> \param[in]       ColIndx         :   An input array of indices indicating the order by which
     843             :     !>                                      the columns in the input matrix must be rearranged.
     844             :     !> \param[in]       ColIndxMap      :   An input array of length `lenColIndxs` that indicated the
     845             :     !>                                      new column index of the corresponding column indices in `ColIndx`.
     846             :     !>
     847             :     !> \return
     848             :     !> `SortedPosDefMatUpper` : An ordered matrix of size `( rank, rank )`.
     849             :     !>
     850             :     !> \author
     851             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     852          33 :     pure function sortPosDefMat(rank,PosDefMatUpper,lenColIndx,ColIndx,ColIndxMap) result(SortedPosDefMatUpper)
     853             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     854             :         !DEC$ ATTRIBUTES DLLEXPORT :: sortPosDefMat
     855             : #endif
     856           1 :         use Constants_mod, only: RK, IK
     857             :         integer(IK), intent(in) :: rank, lenColIndx, ColIndx(lenColIndx), ColIndxMap(lenColIndx)
     858             :         real(RK), intent(in)    :: PosDefMatUpper(rank,rank)
     859             :         real(RK)                :: SortedPosDefMatUpper(rank,rank)
     860             :         integer(IK)             :: iRow, iCol, iRowNew, iColNew, i
     861           2 :         loopIndx: do i = 1, lenColIndx
     862           7 :             do iCol = 1, rank
     863           5 :                 iColNew = iCol
     864           5 :                 if (iColNew==ColIndx(i)) then
     865           1 :                     iColNew = ColIndxMap(i)
     866           4 :                 elseif (iColNew==ColIndxMap(i)) then
     867           1 :                     iColNew = ColIndx(i)
     868             :                 end if
     869          21 :                 do iRow = 1, iCol
     870          15 :                     iRowNew = iRow
     871          15 :                     if (iRowNew==ColIndx(i)) then
     872           2 :                         iRowNew = ColIndxMap(i)
     873          13 :                     elseif (iRowNew==ColIndxMap(i)) then
     874           1 :                         iRowNew = ColIndx(i)
     875             :                     end if
     876          20 :                     if (iRowNew>iColNew) then
     877           1 :                         SortedPosDefMatUpper(iRow,iCol) = PosDefMatUpper(iColNew,iRowNew)
     878             :                     else
     879          14 :                         SortedPosDefMatUpper(iRow,iCol) = PosDefMatUpper(iRowNew,iColNew)
     880             :                     end if
     881             :                 end do
     882             :             end do
     883             :         end do loopIndx
     884           1 :     end function sortPosDefMat
     885             : 
     886             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     887             : 
     888             :     !> \brief
     889             :     !> Symmetrize an input upper-triangular matrix by copying the upper to the lower triangle.
     890             :     !>
     891             :     !> \param[in]       nd                  :   The size of the input matrix `PosDefMatUpper`.
     892             :     !> \param[inout]    UpperSquareMatrix   :   The input upper, and output full matrix. On output, the matrix is symmetrized.
     893             :     !>
     894             :     !> \author
     895             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     896           1 :     pure subroutine symmetrizeUpperSquareMatrix(nd,UpperSquareMatrix)
     897             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     898             :         !DEC$ ATTRIBUTES DLLEXPORT :: symmetrizeUpperSquareMatrix
     899             : #endif
     900           1 :         use Constants_mod, only: RK, IK
     901             :         implicit none
     902             :         integer(IK) , intent(in)    :: nd
     903             :         real(RK)    , intent(inout) :: UpperSquareMatrix(nd,nd)
     904             :         integer(IK)                 :: i, j
     905           4 :         do j = 1, nd
     906           7 :             do i = 1,j-1
     907           6 :                 UpperSquareMatrix(j,i) = UpperSquareMatrix(i,j)
     908             :             end do
     909             :         end do
     910           1 :     end subroutine symmetrizeUpperSquareMatrix
     911             : 
     912             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     913             : 
     914             :     !> \brief
     915             :     !> Return the the Regression Coefficient Matrix, whose dimension is `rankS11` by `rankS22`, as well as optionally the
     916             :     !> Schur complement of the `S22` block of rank `rankS22` of the input matrix `PosDefMat` of rank `rankPDM`.
     917             :     !> For example, if,
     918             :     !> \f{equation}{
     919             :     !> \texttt{PosDefMat} = \begin{pmatrix}
     920             :     !>                          S11 & S12 \\
     921             :     !>                          S21 & S22
     922             :     !>                      \end{pmatrix}
     923             :     !> \f}
     924             :     !> then, the Regression Coefficient Matrix given `S22` is: \f$ S12 \times S22^{-1} \f$, whose rank is `rankS11` by `rankS22`.
     925             :     !> The Schur complement of `S22` is:
     926             :     !> \f{equation}
     927             :     !>     \texttt{SchurComplement} = S11 - S12 \times S22^{-1} \times S21 ~,
     928             :     !> \f}
     929             :     !> whose rank is `rankS11` by `rankS11`.
     930             :     !>
     931             :     !> \warning
     932             :     !> On input, the full matrix PosDefMat must be given.
     933             :     !>
     934             :     !> \remark
     935             :     !> For clarity, note that `rankS11 + rankS22 = rankPDM`.
     936             :     !>
     937             :     !> \author
     938             :     !> Amir Shahmoradi, Apr 21, 2017, 4:10 PM, ICES, UT
     939           2 :     pure subroutine getRegresCoef(rankPDM,rankS11,rankS22,PosDefMat,RegresCoefMat,SchurComplement)
     940             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     941             :         !DEC$ ATTRIBUTES DLLEXPORT :: getRegresCoef
     942             : #endif
     943             :         use Constants_mod, only: RK, IK ! LCOV_EXCL_LINE
     944             :         implicit none
     945             :         integer(IK), intent(in)             :: rankPDM, rankS11, rankS22
     946             :         real(RK), intent(in)                :: PosDefMat(rankPDM,rankPDM)
     947             :         real(RK), intent(out)               :: RegresCoefMat(rankS11,rankS22)
     948             :         real(RK), intent(out), optional     :: SchurComplement(rankS11,rankS11)
     949           5 :         real(RK)                            :: InvS22(rankS22,rankS22), S22(rankS22,rankS22)
     950             :         integer(IK)                         :: startS22
     951           1 :         startS22 = rankS11 + 1
     952           3 :         S22 = PosDefMat(startS22:rankPDM,startS22:rankPDM)
     953           1 :         if (rankS22==1) then
     954           1 :             InvS22(1,1) = 1._RK/S22(1,1)
     955             :         else
     956           0 :             InvS22 = getInvPosDefMat(rankS22,S22)
     957             :         end if
     958           1 :         if (InvS22(1,1)<0._RK) then
     959           0 :             RegresCoefMat(1,1) = -1._RK
     960           0 :             return
     961             :         end if
     962           1 :         RegresCoefMat = matmul( PosDefMat(1:rankS11,startS22:rankPDM) , InvS22 )
     963           1 :         if (present(SchurComplement)) then
     964          13 :             SchurComplement = PosDefMat(1:rankS11,1:rankS11) - matmul( RegresCoefMat , PosDefMat(startS22:rankPDM,1:rankS11) )
     965             :         end if
     966           2 :     end subroutine getRegresCoef
     967             : 
     968             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     969             : 
     970             : end module Matrix_mod ! LCOV_EXCL_LINE

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