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 7250 : 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 1062 : real(RK) :: diagonal
77 : integer(IK) :: i, j
78 1062 : diagonal = 1._RK
79 1062 : if (present(diag)) diagonal = diag
80 1062 : do concurrent(i=1:n, j=1:m)
81 5126 : if (i==j) then
82 1541 : eye(i,j) = diagonal
83 : else
84 979 : eye(i,j) = 0._RK
85 : end if
86 : end do
87 1062 : 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 31154 : 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 1062 : 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 31154 : real(RK) :: summ
117 : integer(IK) :: i
118 81444 : do i=1,nd
119 69480 : summ = PosDefMat(i,i) - dot_product(PosDefMat(i,1:i-1),PosDefMat(i,1:i-1))
120 50302 : if (summ <= 0._RK) then
121 12 : Diagonal(1) = -1._RK
122 12 : return
123 : end if
124 50290 : Diagonal(i) = sqrt(summ)
125 100610 : 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 31154 : 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 31154 : 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 15 : 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 165 : real(RK) :: CholeskyLower(nd,nd) ! Cholesky factor
192 69 : real(RK) :: Diagonal(nd) ! diagonal terms of the inverse matrix
193 15 : real(RK) :: summ
194 : integer(IK) :: i,j,k
195 15 : if (nd==1_IK) then
196 3 : MatInvMat(1,1) = 1._RK / MatInvMat(1,1)
197 3 : sqrtDetInvPosDefMat = MatInvMat(1,1)
198 3 : return
199 : end if
200 48 : do j=1,nd
201 120 : do i=1,j
202 108 : CholeskyLower(i,j) = MatInvMat(i,j)
203 : end do
204 : end do
205 12 : call getCholeskyFactor(nd,CholeskyLower,Diagonal)
206 12 : if (Diagonal(1)<0._RK) then
207 3 : sqrtDetInvPosDefMat = -1._RK
208 3 : return
209 : end if
210 36 : sqrtDetInvPosDefMat = 1._RK / product(Diagonal)
211 36 : do i = 1,nd
212 27 : CholeskyLower(i,i) = 1._RK / Diagonal(i)
213 63 : do j = i+1,nd
214 27 : summ = 0._RK
215 63 : do k = i,j-1
216 63 : summ = summ - CholeskyLower(j,k) * CholeskyLower(k,i)
217 : end do
218 54 : CholeskyLower(j,i) = summ / Diagonal(j)
219 : end do
220 : end do
221 36 : do i = 1,nd
222 81 : do j = i,nd
223 171 : MatInvMat(j,i) = dot_product(CholeskyLower(j:nd,j),CholeskyLower(j:nd,i))
224 : end do
225 144 : MatInvMat(i,i:nd) = MatInvMat(i:nd,i)
226 : end do
227 15 : 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 2544920 : 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 15 : 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 301832 : real(RK) :: summ
257 : integer(IK) :: i,j,k
258 301832 : if (nd==1_IK) then
259 42897 : InvMatFromCholFac(1,1) = 1._RK / Diagonal(1)**2
260 42897 : return
261 : end if
262 1812560 : InvMatFromCholFac = 0._RK
263 517873 : do j=1,nd-1
264 776814 : do i=j+1,nd
265 517879 : InvMatFromCholFac(i,j) = CholeskyLower(i,j)
266 : end do
267 : end do
268 776808 : do i = 1,nd
269 517873 : InvMatFromCholFac(i,i) = 1._RK / Diagonal(i)
270 1035750 : do j = i+1,nd
271 258941 : summ = 0._RK
272 517885 : do k = i,j-1
273 517885 : summ = summ - InvMatFromCholFac(j,k) * InvMatFromCholFac(k,i)
274 : end do
275 776814 : InvMatFromCholFac(j,i) = summ / Diagonal(j)
276 : end do
277 : end do
278 776808 : do i = 1,nd
279 1553620 : do j = i,nd
280 1812570 : InvMatFromCholFac(j,i) = dot_product(InvMatFromCholFac(j:nd,j),InvMatFromCholFac(j:nd,i))
281 1294690 : InvMatFromCholFac(i,j) = InvMatFromCholFac(j,i)
282 : end do
283 : !InvMatFromCholFac(i,i:nd) = InvMatFromCholFac(i:nd,i)
284 : end do
285 603664 : 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 90 : 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 301832 : use Constants_mod, only: RK, IK
321 : implicit none
322 : integer(IK), intent(in) :: nd
323 : real(RK) , intent(in) :: PosDefMat(nd,nd)
324 78 : real(RK) :: InvPosDefMat(nd,nd), CholeskyLower(nd,nd)
325 30 : real(RK) :: Diagonal(nd)
326 6 : 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 24 : do j=1,nd
333 60 : do i=1,j
334 54 : CholeskyLower(i,j) = PosDefMat(i,j)
335 : end do
336 : end do
337 6 : call getCholeskyFactor(nd,CholeskyLower,Diagonal)
338 6 : if (Diagonal(1)<0._RK) then
339 39 : InvPosDefMat = -1._RK ! error occurred: getCholeskyFactor() failed in getInvPosDefMat()
340 3 : return
341 : end if
342 12 : do i = 1,nd
343 9 : CholeskyLower(i,i) = 1._RK / Diagonal(i)
344 21 : do j = i+1,nd
345 9 : summ = 0._RK
346 21 : do k = i,j-1
347 21 : summ = summ - CholeskyLower(j,k) * CholeskyLower(k,i)
348 : end do
349 18 : CholeskyLower(j,i) = summ / Diagonal(j)
350 : end do
351 : end do
352 12 : do i = 1,nd
353 30 : do j = i,nd
354 48 : InvPosDefMat(j,i) = dot_product(CholeskyLower(j:nd,j),CholeskyLower(j:nd,i))
355 27 : InvPosDefMat(i,j) = InvPosDefMat(j,i)
356 : end do
357 : end do
358 12 : 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 3 : 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 6 : 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 12 : do i = 1,nd
393 36 : do j = 1,nd
394 36 : InverseMatrix(i,j) = 0._RK
395 : end do
396 12 : InverseMatrix(i,i) = 1._RK
397 : end do
398 3 : call getLU(nd,MatrixLU,Permutation,detInvMat)
399 12 : do j = 1,nd
400 9 : detInvMat = detInvMat * MatrixLU(j,j)
401 12 : call solveLinearSystem(nd,MatrixLU,Permutation,InverseMatrix(1:nd,j))
402 : end do
403 3 : detInvMat = 1._RK/detInvMat
404 3 : 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 45 : 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 3 : use Constants_mod, only: RK, IK
430 : implicit none
431 : integer(IK), intent(in) :: nd
432 : real(RK) , intent(in) :: Matrix(nd,nd)
433 42 : real(RK) :: InverseMatrix(nd,nd),LU(nd,nd)
434 3 : integer(IK) :: i,j,Permutation(nd)
435 3 : real(RK) :: parity
436 12 : do i = 1,nd
437 36 : do j = 1,nd
438 36 : InverseMatrix(i,j) = 0._RK
439 : end do
440 12 : InverseMatrix(i,i) = 1._RK
441 : end do
442 39 : LU = Matrix
443 3 : call getLU(nd,LU,Permutation,parity)
444 12 : do j = 1,nd
445 12 : call solveLinearSystem(nd,LU,Permutation,InverseMatrix(1:nd,j))
446 : end do
447 3 : 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 18 : 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 3 : 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 18 : real(RK) :: summ
481 18 : ii = 0
482 72 : do i=1,nd
483 54 : summ = InputSolution(Permutation(i))
484 54 : InputSolution(Permutation(i)) = InputSolution(i)
485 54 : if (ii /= 0) then
486 42 : summ = summ - dot_product( LU(i,ii:i-1) , InputSolution(ii:i-1) )
487 36 : else if (summ /= 0._RK) then
488 18 : ii = i
489 : end if
490 72 : InputSolution(i) = summ
491 : end do
492 72 : do i = nd,1,-1
493 126 : InputSolution(i) = ( InputSolution(i) - dot_product( LU(i,i+1:nd) , InputSolution(i+1:nd)) ) / LU(i,i)
494 : end do
495 18 : 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 9 : 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 18 : 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 36 : real(RK) :: aamax,dum,summ,vv(nd)
525 : integer(IK) :: i,imax,j,k
526 : !errorOccurred = .false.
527 9 : parity = 1._RK
528 36 : do i = 1,nd
529 27 : aamax = 0._RK
530 108 : do j=1,nd
531 108 : if ( abs(MatrixLU(i,j)) > aamax ) aamax = abs( MatrixLU(i,j) )
532 : end do
533 27 : 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 36 : vv(i) = 1._RK/aamax
542 : end do
543 36 : do j=1,nd
544 54 : do i=1,j-1
545 27 : summ = MatrixLU(i,j)
546 36 : do k=1,i-1
547 36 : summ = summ - MatrixLU(i,k) * MatrixLU(k,j)
548 : end do
549 54 : MatrixLU(i,j) = summ
550 : end do
551 27 : aamax = 0._RK
552 81 : do i = j, nd
553 54 : summ = MatrixLU(i,j)
554 90 : do k = 1, j-1
555 90 : summ = summ - MatrixLU(i,k) * MatrixLU(k,j)
556 : end do
557 54 : MatrixLU(i,j) = summ
558 54 : dum = vv(i) * abs(summ)
559 81 : if (dum >= aamax) then
560 27 : imax = i
561 27 : aamax = dum
562 : end if
563 : end do
564 27 : 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 27 : Permutation(j) = imax
574 27 : if (MatrixLU(j,j) == 0._RK) MatrixLU(j,j) = TINY
575 36 : if (j /= nd) then
576 18 : dum = 1._RK / MatrixLU(j,j)
577 45 : do i = j+1,nd
578 45 : MatrixLU(i,j) = MatrixLU(i,j) * dum
579 : end do
580 : endif
581 : end do
582 9 : 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 3 : 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 9 : 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 3 : if (colsA == rowsB) then
610 3 : rowsC = rowsA
611 3 : 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 12 : do i = 1, rowsC
621 39 : do j = 1, colsC
622 36 : C(i,j) = 0._RK
623 : end do
624 : end do
625 : ! Find product as per above formula
626 12 : do i = 1, rowsA
627 39 : do j = 1, colsB
628 117 : do k = 1, colsA
629 108 : C(i,j) = C(i,j) + A(i,k)*B(k,j)
630 : end do
631 : end do
632 : end do
633 3 : 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 3 : 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 3 : use Constants_mod, only: RK, IK
653 : implicit none
654 : integer(IK), intent(in) :: nd
655 : real(RK) , intent(in) :: Matrix(nd,nd)
656 39 : real(RK) :: DummyMat(nd,nd)
657 : real(RK) :: determinant
658 3 : integer(IK) :: Permutation(nd)
659 : integer(IK) :: j
660 39 : DummyMat = Matrix
661 3 : call getLU(nd,DummyMat,Permutation,determinant) ! This returns determinant as +-1.
662 12 : do j=1,nd
663 12 : determinant = determinant * DummyMat(j,j)
664 : end do
665 3 : 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 6 : 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 3 : use Constants_mod, only: RK, IK
687 : implicit none
688 : integer(IK), intent(in) :: nd
689 : real(RK) , intent(in) :: PosDefMat(nd,nd)
690 96 : real(RK) :: Diagonal(nd),DummyMat(nd,nd)
691 : real(RK) :: sqrtDetPosDefMat
692 : integer(IK) :: i,j
693 24 : do j=1,nd
694 60 : do i=1,j
695 54 : DummyMat(i,j) = PosDefMat(i,j)
696 : end do
697 : end do
698 6 : call getCholeskyFactor(nd,DummyMat,Diagonal)
699 6 : if (Diagonal(1)<0._RK) then
700 3 : sqrtDetPosDefMat = -1._RK
701 3 : return
702 : end if
703 12 : sqrtDetPosDefMat = product(Diagonal)
704 6 : 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 115130 : 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 6 : 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 333601 : real(RK) :: Diagonal(nd)
730 115130 : real(RK) :: summ
731 : integer(IK) :: i
732 115130 : failed = .false.
733 333592 : do i=1,nd
734 321803 : summ = PosDefMat(i,i) - dot_product(PosDefMat(i,1:i-1),PosDefMat(i,1:i-1))
735 218465 : if (summ <= 0._RK) then
736 3 : failed = .true.
737 3 : return
738 : end if
739 218462 : Diagonal(i) = sqrt(summ)
740 436927 : 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 333589 : logSqrtDetPosDefMat = sum(log(Diagonal))
743 115130 : 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 2976 : 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 115130 : use Constants_mod, only: RK, IK
767 : implicit none
768 : integer(IK), intent(in) :: nd
769 : real(RK) , intent(in) :: Matrix(nd,nd)
770 24042 : real(RK) :: Array(nd,nd),p(nd)
771 2976 : real(RK) :: dummySum
772 : integer(IK) :: i,j,k
773 2976 : isPosDef = .true.
774 7755 : do i = 1, nd
775 14361 : do j = i, nd
776 6627 : dummySum = Matrix(i,j)
777 8457 : do k = i-1,1,-1
778 8457 : dummySum = dummySum - Array(i,k) * Array(j,k)
779 : end do
780 11406 : if(i==j)then
781 4800 : if(dummySum<=0._RK) then
782 21 : isPosDef = .false.
783 21 : return
784 : end if
785 4779 : p(i) = sqrt(dummySum)
786 : else
787 1827 : Array(j,i)=dummySum/p(i)
788 : endif
789 : end do
790 : end do
791 2976 : 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 33 : 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 2976 : use Constants_mod, only: RK, IK
811 : real(RK), intent(in) :: Vector1(:),Vector2(:)
812 6 : real(RK) :: OuterProd(size(Vector1), size(Vector2))
813 27 : OuterProd = spread(Vector1,dim=2,ncopies=size(Vector2)) * spread(Vector2,dim=1,ncopies=size(Vector1))
814 6 : 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 99 : 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 3 : 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 6 : loopIndx: do i = 1, lenColIndx
862 21 : do iCol = 1, rank
863 15 : iColNew = iCol
864 15 : if (iColNew==ColIndx(i)) then
865 3 : iColNew = ColIndxMap(i)
866 12 : elseif (iColNew==ColIndxMap(i)) then
867 3 : iColNew = ColIndx(i)
868 : end if
869 63 : do iRow = 1, iCol
870 45 : iRowNew = iRow
871 45 : if (iRowNew==ColIndx(i)) then
872 6 : iRowNew = ColIndxMap(i)
873 39 : elseif (iRowNew==ColIndxMap(i)) then
874 3 : iRowNew = ColIndx(i)
875 : end if
876 60 : if (iRowNew>iColNew) then
877 3 : SortedPosDefMatUpper(iRow,iCol) = PosDefMatUpper(iColNew,iRowNew)
878 : else
879 42 : SortedPosDefMatUpper(iRow,iCol) = PosDefMatUpper(iRowNew,iColNew)
880 : end if
881 : end do
882 : end do
883 : end do loopIndx
884 3 : 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 3 : 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 3 : 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 12 : do j = 1, nd
906 21 : do i = 1,j-1
907 18 : UpperSquareMatrix(j,i) = UpperSquareMatrix(i,j)
908 : end do
909 : end do
910 3 : 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 6 : 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 15 : real(RK) :: InvS22(rankS22,rankS22), S22(rankS22,rankS22)
950 : integer(IK) :: startS22
951 3 : startS22 = rankS11 + 1
952 9 : S22 = PosDefMat(startS22:rankPDM,startS22:rankPDM)
953 3 : if (rankS22==1) then
954 3 : InvS22(1,1) = 1._RK/S22(1,1)
955 : else
956 0 : InvS22 = getInvPosDefMat(rankS22,S22)
957 : end if
958 3 : if (InvS22(1,1)<0._RK) then
959 0 : RegresCoefMat(1,1) = -1._RK
960 0 : return
961 : end if
962 3 : RegresCoefMat = matmul( PosDefMat(1:rankS11,startS22:rankPDM) , InvS22 )
963 3 : if (present(SchurComplement)) then
964 39 : SchurComplement = PosDefMat(1:rankS11,1:rankS11) - matmul( RegresCoefMat , PosDefMat(startS22:rankPDM,1:rankS11) )
965 : end if
966 6 : end subroutine getRegresCoef
967 :
968 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
969 :
970 : end module Matrix_mod ! LCOV_EXCL_LINE
|
