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 tests of the module [Matrix_mod](@ref matrix_mod).
44 : !> \author Amir Shahmoradi
45 :
46 : module Test_Matrix_mod
47 :
48 : use Matrix_mod ! LCOV_EXCL_LINE
49 : use Err_mod, only: Err_type
50 : use Test_mod, only: Test_type
51 : implicit none
52 :
53 : private
54 : public :: test_Matrix
55 :
56 : type(Test_type) :: Test
57 :
58 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59 :
60 : contains
61 :
62 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 :
64 1 : subroutine test_Matrix()
65 :
66 : implicit none
67 :
68 1 : Test = Test_type(moduleName=MODULE_NAME)
69 1 : call Test%run(test_eye_1, "test_eye_1")
70 1 : call Test%run(test_eye_2, "test_eye_2")
71 1 : call Test%run(test_isPosDef_1, "test_isPosDef_1")
72 1 : call Test%run(test_isPosDef_2, "test_isPosDef_2")
73 1 : call Test%run(test_getInvMat_1, "test_getInvMat_1")
74 1 : call Test%run(test_getOuterProd_1, "test_getOuterProd_1")
75 1 : call Test%run(test_getInvMatDet_1, "test_getInvMatDet_1")
76 1 : call Test%run(test_sortPosDefMat_1, "test_sortPosDefMat_1")
77 1 : call Test%run(test_getRegresCoef_1, "test_getRegresCoef_1")
78 1 : call Test%run(test_multiplyMatrix_1, "test_multiplyMatrix_1")
79 1 : call Test%run(test_getDeterminant_1, "test_getDeterminant_1")
80 1 : call Test%run(test_getInvPosDefMat_1, "test_getInvPosDefMat_1")
81 1 : call Test%run(test_getInvPosDefMat_2, "test_getInvPosDefMat_2")
82 1 : call Test%run(test_getCholeskyFactor_1, "test_getCholeskyFactor_1")
83 1 : call Test%run(test_getCholeskyFactor_2, "test_getCholeskyFactor_2")
84 1 : call Test%run(test_getInvMatFromCholFac, "test_getInvMatFromCholFac")
85 1 : call Test%run(test_getSqrtDetPosDefMat_1, "test_getSqrtDetPosDefMat_1")
86 1 : call Test%run(test_getSqrtDetPosDefMat_2, "test_getSqrtDetPosDefMat_2")
87 1 : call Test%run(test_getInvPosDefMatSqrtDet_1, "test_getInvPosDefMatSqrtDet_1")
88 1 : call Test%run(test_getInvPosDefMatSqrtDet_2, "test_getInvPosDefMatSqrtDet_2")
89 1 : call Test%run(test_getInvPosDefMatSqrtDet_3, "test_getInvPosDefMatSqrtDet_3")
90 1 : call Test%run(test_getLogSqrtDetPosDefMat_1, "test_getLogSqrtDetPosDefMat_1")
91 1 : call Test%run(test_getLogSqrtDetPosDefMat_2, "test_getLogSqrtDetPosDefMat_2")
92 1 : call Test%run(test_symmetrizeUpperSquareMatrix_1, "test_symmetrizeUpperSquareMatrix_1")
93 1 : call Test%finalize()
94 :
95 1 : end subroutine test_Matrix
96 :
97 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
98 :
99 : !> \brief
100 : !> Test whether [getEye()](@ref matrix_mod::geteye) return a correct output eye matrix.
101 1 : function test_eye_1() result(assertion)
102 :
103 1 : use Constants_mod, only: IK, RK
104 : implicit none
105 : logical :: assertion
106 : integer(IK) , parameter :: n = 3_IK, m = 4_IK
107 : real(RK) , parameter :: Eye_ref(n, m) = reshape( [ 1._RK, 0._RK, 0._RK &
108 : , 0._RK, 1._RK, 0._RK &
109 : , 0._RK, 0._RK, 1._RK &
110 : , 0._RK, 0._RK, 0._RK ], shape = shape(eye_ref) )
111 : real(RK), allocatable :: Eye(:,:)
112 : real(RK), allocatable :: Difference(:,:)
113 1 : Eye = getEye(n,m)
114 :
115 17 : Difference = abs(Eye - Eye_ref)
116 17 : assertion = all(Difference==0._RK)
117 :
118 1 : if (Test%isDebugMode .and. .not. assertion) then
119 : ! LCOV_EXCL_START
120 : write(Test%outputUnit,"(*(g0,:,', '))")
121 : write(Test%outputUnit,"(*(g0,:,', '))") "Eye_ref = ", Eye_ref
122 : write(Test%outputUnit,"(*(g0,:,', '))") "Eye = ", Eye
123 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference = ", Difference
124 : write(Test%outputUnit,"(*(g0,:,', '))")
125 : end if
126 : ! LCOV_EXCL_STOP
127 :
128 1 : end function test_eye_1
129 :
130 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 :
132 : !> \brief
133 : !> Test whether [getEye()](@ref matrix_mod::geteye) return a correct output eye matrix with the input optional argument.
134 1 : function test_eye_2() result(assertion)
135 :
136 1 : use Constants_mod, only: IK, RK
137 : implicit none
138 : logical :: assertion
139 : integer(IK) , parameter :: n = 3_IK, m = 3_IK
140 : real(RK) , parameter :: Eye_ref(n, m) = reshape( [ 2._RK, 0._RK, 0._RK &
141 : , 0._RK, 2._RK, 0._RK &
142 : , 0._RK, 0._RK, 2._RK ], shape = shape(eye_ref) )
143 : real(RK), allocatable :: Eye(:,:)
144 : real(RK), allocatable :: Difference(:,:)
145 1 : Eye = getEye(n,m, diag = 2._RK)
146 :
147 13 : Difference = abs(Eye - Eye_ref)
148 13 : assertion = all(Difference==0._RK)
149 :
150 1 : if (Test%isDebugMode .and. .not. assertion) then
151 : ! LCOV_EXCL_START
152 : write(Test%outputUnit,"(*(g0,:,', '))")
153 : write(Test%outputUnit,"(*(g0,:,', '))") "Eye_ref = ", Eye_ref
154 : write(Test%outputUnit,"(*(g0,:,', '))") "Eye = ", Eye
155 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference = ", Difference
156 : write(Test%outputUnit,"(*(g0,:,', '))")
157 : end if
158 : ! LCOV_EXCL_STOP
159 :
160 1 : end function test_eye_2
161 :
162 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163 :
164 1 : function test_getCholeskyFactor_1() result(assertion)
165 :
166 1 : use Constants_mod, only: IK, RK
167 : implicit none
168 : logical :: assertion
169 : integer(IK) , parameter :: nd = 3_IK
170 : real(RK) , parameter :: tolerance = 1.e-12_RK
171 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
172 : , 0._RK, 2._RK, 0._RK &
173 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
174 : real(RK) , parameter :: CholeskyLower_ref(nd,nd) = reshape( [ 1.000000000000000_RK, 0.000000000000000_RK, 1.000000000000000_RK &
175 : , 0.000000000000000_RK, 2.000000000000000_RK, 0.000000000000000_RK &
176 : , 1.000000000000000_RK, 0.000000000000000_RK, 3.000000000000000_RK ] &
177 : , shape = shape(CholeskyLower_ref) )
178 : real(RK) , parameter :: CholeskyDiagonal_ref(nd) = [ 1.000000000000000_RK, 1.414213562373095_RK, 1.414213562373095_RK ]
179 : real(RK) :: CholeskyLower(nd,nd), CholeskyDiagonal(nd)
180 : real(RK), allocatable :: CholeskyLower_diff(:,:), CholeskyDiagonal_diff(:)
181 1 : CholeskyLower = PosDefMat
182 1 : call getCholeskyFactor(nd = nd, PosDefMat = CholeskyLower, Diagonal = CholeskyDiagonal)
183 :
184 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
185 1 : if (allocated(CholeskyLower_diff)) deallocate(CholeskyLower_diff); allocate(CholeskyLower_diff, mold = PosDefMat)
186 14 : CholeskyLower_diff = abs(PosDefMat - CholeskyLower_ref)
187 :
188 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
189 1 : if (allocated(CholeskyDiagonal_diff)) deallocate(CholeskyDiagonal_diff); allocate(CholeskyDiagonal_diff, mold = CholeskyDiagonal)
190 5 : CholeskyDiagonal_diff = abs(CholeskyDiagonal - CholeskyDiagonal_ref)
191 16 : assertion = all(CholeskyLower_diff < tolerance) .and. all(CholeskyDiagonal_diff < tolerance)
192 :
193 1 : if (Test%isDebugMode .and. .not. assertion) then
194 : ! LCOV_EXCL_START
195 : write(Test%outputUnit,"(*(g0,:,', '))")
196 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower_ref = ", CholeskyLower_ref
197 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower = ", CholeskyLower
198 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower_diff = ", CholeskyLower_diff
199 : write(Test%outputUnit,"(*(g0,:,', '))")
200 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyDiagonal_ref = ", CholeskyDiagonal_ref
201 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyDiagonal = ", CholeskyDiagonal
202 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyDiagonal_diff = ", CholeskyDiagonal_diff
203 : write(Test%outputUnit,"(*(g0,:,', '))")
204 : end if
205 : ! LCOV_EXCL_STOP
206 :
207 1 : end function test_getCholeskyFactor_1
208 :
209 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210 :
211 1 : function test_getCholeskyFactor_2() result(assertion)
212 1 : use Constants_mod, only: IK, RK
213 : implicit none
214 : logical :: assertion
215 : integer(IK) , parameter :: nd = 3_IK
216 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, -1._RK &
217 : , 0._RK, 2._RK, -0._RK &
218 : , 1._RK, 0._RK, -3._RK ], shape = shape(PosDefMat) )
219 : real(RK) :: CholeskyLower(nd,nd), CholeskyDiagonal(nd)
220 1 : CholeskyLower = PosDefMat
221 1 : call getCholeskyFactor(nd = nd, PosDefMat = CholeskyLower, Diagonal = CholeskyDiagonal)
222 1 : assertion = CholeskyDiagonal(1) < 0._RK
223 1 : if (Test%isDebugMode .and. .not. assertion) then
224 : ! LCOV_EXCL_START
225 : write(Test%outputUnit,"(*(g0,:,', '))")
226 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower = ", CholeskyLower
227 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyDiagonal = ", CholeskyDiagonal
228 : write(Test%outputUnit,"(*(g0,:,', '))")
229 : end if
230 : ! LCOV_EXCL_STOP
231 1 : end function test_getCholeskyFactor_2
232 :
233 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
234 :
235 1 : function test_getInvPosDefMatSqrtDet_1() result(assertion)
236 :
237 1 : use Constants_mod, only: IK, RK
238 : implicit none
239 :
240 : logical :: assertion
241 : integer(IK) , parameter :: nd = 3_IK
242 : real(RK) , parameter :: tolerance = 1.e-12_RK
243 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
244 : , 0._RK, 2._RK, 0._RK &
245 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
246 : real(RK) , parameter :: MatInvMat_ref(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
247 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
248 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
249 : , shape = shape(MatInvMat_ref) )
250 : real(RK) , parameter :: CholeskyDiagonal_ref(nd) = [ 1.000000000000000_RK, 1.414213562373095_RK, 1.414213562373095_RK ]
251 : real(RK) , parameter :: sqrtDetInvPosDefMat_ref = 0.5_RK
252 1 : real(RK) :: MatInvMat(nd,nd), sqrtDetInvPosDefMat
253 1 : real(RK), allocatable :: MatInvMat_diff(:,:), sqrtDetInvPosDefMat_diff
254 :
255 1 : MatInvMat = PosDefMat
256 :
257 1 : call getInvPosDefMatSqrtDet(nd = nd, MatInvMat = MatInvMat, sqrtDetInvPosDefMat = sqrtDetInvPosDefMat)
258 :
259 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
260 1 : if (allocated(MatInvMat_diff)) deallocate(MatInvMat_diff); allocate(MatInvMat_diff, mold = MatInvMat)
261 :
262 14 : MatInvMat_diff = abs(MatInvMat - MatInvMat_ref)
263 1 : sqrtDetInvPosDefMat_diff = abs(sqrtDetInvPosDefMat - sqrtDetInvPosDefMat_ref)
264 :
265 13 : assertion = all(MatInvMat_diff < tolerance) .and. sqrtDetInvPosDefMat_diff < tolerance
266 :
267 1 : if (Test%isDebugMode .and. .not. assertion) then
268 : ! LCOV_EXCL_START
269 : write(Test%outputUnit,"(*(g0,:,', '))")
270 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_ref = ", MatInvMat_ref
271 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat = ", MatInvMat
272 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_diff = ", MatInvMat_diff
273 : write(Test%outputUnit,"(*(g0,:,', '))")
274 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat_ref = ", sqrtDetInvPosDefMat_ref
275 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat_diff = ", sqrtDetInvPosDefMat
276 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat = ", sqrtDetInvPosDefMat_diff
277 : write(Test%outputUnit,"(*(g0,:,', '))")
278 : end if
279 : ! LCOV_EXCL_STOP
280 :
281 1 : end function test_getInvPosDefMatSqrtDet_1
282 :
283 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
284 :
285 1 : function test_getInvPosDefMatSqrtDet_2() result(assertion)
286 1 : use Constants_mod, only: IK, RK
287 : implicit none
288 : logical :: assertion
289 : integer(IK) , parameter :: nd = 3_IK
290 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, -1._RK &
291 : , 0._RK, 2._RK, -0._RK &
292 : , 1._RK, 0._RK, -3._RK ], shape = shape(PosDefMat) )
293 1 : real(RK) :: MatInvMat(nd,nd), sqrtDetInvPosDefMat
294 :
295 1 : MatInvMat = PosDefMat
296 :
297 1 : call getInvPosDefMatSqrtDet(nd = nd, MatInvMat = MatInvMat, sqrtDetInvPosDefMat = sqrtDetInvPosDefMat)
298 :
299 1 : assertion = sqrtDetInvPosDefMat < 0._RK
300 :
301 1 : if (Test%isDebugMode .and. .not. assertion) then
302 : ! LCOV_EXCL_START
303 : write(Test%outputUnit,"(*(g0,:,', '))")
304 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat = ", MatInvMat
305 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat = ", sqrtDetInvPosDefMat
306 : write(Test%outputUnit,"(*(g0,:,', '))")
307 : end if
308 : ! LCOV_EXCL_STOP
309 :
310 1 : end function test_getInvPosDefMatSqrtDet_2
311 :
312 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
313 :
314 : !> \brief
315 : !> Test with an 1-dimensional input matrix.
316 1 : function test_getInvPosDefMatSqrtDet_3() result(assertion)
317 1 : use Constants_mod, only: IK, RK
318 : implicit none
319 : logical :: assertion
320 : integer(IK) , parameter :: nd = 1_IK
321 : real(RK) , parameter :: tolerance = 1.e-12_RK
322 : real(RK) , parameter :: MatInvMat_ref(nd,nd) = reshape( [ 0.5_RK ], shape = shape(MatInvMat_ref) )
323 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 2._RK ], shape = shape(PosDefMat) )
324 : real(RK) , parameter :: sqrtDetInvPosDefMat_ref = 0.5_RK
325 3 : real(RK) :: MatInvMat(nd,nd), sqrtDetInvPosDefMat
326 1 : real(RK), allocatable :: MatInvMat_diff(:,:), sqrtDetInvPosDefMat_diff
327 :
328 1 : MatInvMat = PosDefMat
329 :
330 1 : call getInvPosDefMatSqrtDet(nd = nd, MatInvMat = MatInvMat, sqrtDetInvPosDefMat = sqrtDetInvPosDefMat)
331 1 : if (allocated(MatInvMat_diff)) deallocate(MatInvMat_diff); allocate(MatInvMat_diff, mold = MatInvMat)
332 :
333 4 : MatInvMat_diff = abs(MatInvMat - MatInvMat_ref)
334 1 : sqrtDetInvPosDefMat_diff = abs(sqrtDetInvPosDefMat - sqrtDetInvPosDefMat_ref)
335 :
336 3 : assertion = all(MatInvMat_diff < tolerance) .and. sqrtDetInvPosDefMat_diff < tolerance
337 :
338 1 : if (Test%isDebugMode .and. .not. assertion) then
339 : ! LCOV_EXCL_START
340 : write(Test%outputUnit,"(*(g0,:,', '))")
341 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_ref = ", MatInvMat_ref
342 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat = ", MatInvMat
343 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_diff = ", MatInvMat_diff
344 : write(Test%outputUnit,"(*(g0,:,', '))")
345 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat_ref = ", sqrtDetInvPosDefMat_ref
346 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat_diff = ", sqrtDetInvPosDefMat
347 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvPosDefMat = ", sqrtDetInvPosDefMat_diff
348 : write(Test%outputUnit,"(*(g0,:,', '))")
349 : end if
350 : ! LCOV_EXCL_STOP
351 :
352 :
353 1 : end function test_getInvPosDefMatSqrtDet_3
354 :
355 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356 :
357 1 : function test_getInvMatFromCholFac() result(assertion)
358 :
359 1 : use Constants_mod, only: IK, RK
360 : implicit none
361 :
362 : logical :: assertion
363 : integer(IK) , parameter :: nd = 3_IK
364 : real(RK) , parameter :: tolerance = 1.e-12_RK
365 : real(RK) , parameter :: CholeskyLower(nd,nd) = reshape( [ 1.000000000000000_RK, 0.000000000000000_RK, 1.000000000000000_RK &
366 : , 0.000000000000000_RK, 2.000000000000000_RK, 0.000000000000000_RK &
367 : , 1.000000000000000_RK, 0.000000000000000_RK, 3.000000000000000_RK ] &
368 : , shape = shape(CholeskyLower) )
369 : real(RK) , parameter :: CholeskyDiagonal(nd) = [ 1.000000000000000_RK, 1.414213562373095_RK, 1.414213562373095_RK ]
370 : real(RK) , parameter :: CholeskyDiagonal_ref(nd) = [ 1.000000000000000_RK, 1.414213562373095_RK, 1.414213562373095_RK ]
371 : real(RK) , parameter :: InvMatFromCholFac_ref(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
372 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
373 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
374 : , shape = shape(InvMatFromCholFac_ref) )
375 : real(RK) , parameter :: sqrtDetInvPosDefMat_ref = 0.5_RK
376 : real(RK) :: InvMatFromCholFac(nd,nd)
377 : real(RK), allocatable :: InvMatFromCholFac_diff(:,:)
378 :
379 1 : InvMatFromCholFac = getInvMatFromCholFac(nd = nd, CholeskyLower = CholeskyLower, Diagonal = CholeskyDiagonal)
380 :
381 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
382 1 : if (allocated(InvMatFromCholFac_diff)) deallocate(InvMatFromCholFac_diff); allocate(InvMatFromCholFac_diff, mold = InvMatFromCholFac)
383 :
384 14 : InvMatFromCholFac_diff = abs(InvMatFromCholFac - InvMatFromCholFac_ref)
385 :
386 13 : assertion = all(InvMatFromCholFac_diff < tolerance)
387 :
388 1 : if (Test%isDebugMode .and. .not. assertion) then
389 : ! LCOV_EXCL_START
390 : write(Test%outputUnit,"(*(g0,:,', '))")
391 : write(Test%outputUnit,"(*(g0,:,', '))") "InvMatFromCholFac_ref = ", InvMatFromCholFac_ref
392 : write(Test%outputUnit,"(*(g0,:,', '))") "InvMatFromCholFac = ", InvMatFromCholFac
393 : write(Test%outputUnit,"(*(g0,:,', '))") "InvMatFromCholFac_diff = ", InvMatFromCholFac_diff
394 : write(Test%outputUnit,"(*(g0,:,', '))")
395 : end if
396 : ! LCOV_EXCL_STOP
397 :
398 1 : end function test_getInvMatFromCholFac
399 :
400 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
401 :
402 1 : function test_getInvPosDefMat_1() result(assertion)
403 :
404 1 : use Constants_mod, only: IK, RK
405 : implicit none
406 :
407 : logical :: assertion
408 : integer(IK) , parameter :: nd = 3_IK
409 : real(RK) , parameter :: tolerance = 1.e-12_RK
410 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
411 : , 0._RK, 2._RK, 0._RK &
412 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
413 : real(RK) , parameter :: MatInvMat_ref(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
414 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
415 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
416 : , shape = shape(MatInvMat_ref) )
417 : real(RK) :: MatInvMat(nd,nd)
418 : real(RK), allocatable :: MatInvMat_diff(:,:)
419 :
420 1 : MatInvMat = getInvPosDefMat(nd = nd, PosDefMat = PosDefMat)
421 :
422 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
423 1 : if (allocated(MatInvMat_diff)) deallocate(MatInvMat_diff); allocate(MatInvMat_diff, mold = MatInvMat)
424 14 : MatInvMat_diff = abs(MatInvMat - MatInvMat_ref)
425 :
426 13 : assertion = all(MatInvMat_diff < tolerance)
427 :
428 1 : if (Test%isDebugMode .and. .not. assertion) then
429 : ! LCOV_EXCL_START
430 : write(Test%outputUnit,"(*(g0,:,', '))")
431 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_ref = ", MatInvMat_ref
432 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat = ", MatInvMat
433 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_diff = ", MatInvMat_diff
434 : write(Test%outputUnit,"(*(g0,:,', '))")
435 : end if
436 : ! LCOV_EXCL_STOP
437 :
438 1 : end function test_getInvPosDefMat_1
439 :
440 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
441 :
442 1 : function test_getInvPosDefMat_2() result(assertion)
443 :
444 1 : use Constants_mod, only: IK, RK
445 : implicit none
446 :
447 : logical :: assertion
448 : integer(IK) , parameter :: nd = 3_IK
449 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, -1._RK &
450 : , 0._RK, 2._RK, -0._RK &
451 : , 1._RK, 0._RK, -3._RK ], shape = shape(PosDefMat) )
452 : real(RK) :: MatInvMat(nd,nd)
453 :
454 1 : MatInvMat = getInvPosDefMat(nd = nd, PosDefMat = PosDefMat)
455 :
456 1 : assertion = MatInvMat(1,1) < 0._RK
457 :
458 1 : if (Test%isDebugMode .and. .not. assertion) then
459 : ! LCOV_EXCL_START
460 : write(Test%outputUnit,"(*(g0,:,', '))")
461 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat = ", MatInvMat
462 : write(Test%outputUnit,"(*(g0,:,', '))")
463 : end if
464 : ! LCOV_EXCL_STOP
465 :
466 1 : end function test_getInvPosDefMat_2
467 :
468 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
469 :
470 1 : function test_getInvMatDet_1() result(assertion)
471 :
472 1 : use Constants_mod, only: IK, RK
473 : implicit none
474 :
475 : logical :: assertion
476 : integer(IK) , parameter :: nd = 3_IK
477 : real(RK) , parameter :: tolerance = 1.e-12_RK
478 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
479 : , 0._RK, 2._RK, 0._RK &
480 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
481 : real(RK) , parameter :: MatInvMat_ref(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
482 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
483 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
484 : , shape = shape(MatInvMat_ref) )
485 : real(RK) , parameter :: MatrixLU_ref(nd,nd) = reshape( [ 1.000000000000000_RK, 0.000000000000000_RK, 1.000000000000000_RK &
486 : , 0.000000000000000_RK, 2.000000000000000_RK, 0.000000000000000_RK &
487 : , 1.000000000000000_RK, 0.000000000000000_RK, 2.000000000000000_RK ] &
488 : , shape = shape(MatInvMat_ref) )
489 : real(RK) , parameter :: detInvMat_ref = 0.25_RK
490 1 : real(RK) :: MatInvMat(nd,nd), detInvMat, detInvMat_diff
491 : real(RK), allocatable :: MatrixLU(:,:), MatInvMat_diff(:,:), MatrixLU_diff(:,:)
492 :
493 13 : MatrixLU = PosDefMat
494 :
495 1 : call getInvMatDet(nd = nd, MatrixLU = MatrixLU, InverseMatrix = MatInvMat, detInvMat = detInvMat)
496 :
497 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
498 1 : if (allocated(MatrixLU_diff)) deallocate(MatrixLU_diff); allocate(MatrixLU_diff, mold = MatrixLU)
499 14 : MatrixLU_diff = abs(MatrixLU - MatrixLU_ref)
500 :
501 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
502 1 : if (allocated(MatInvMat_diff)) deallocate(MatInvMat_diff); allocate(MatInvMat_diff, mold = MatInvMat)
503 14 : MatInvMat_diff = abs(MatInvMat - MatInvMat_ref)
504 :
505 1 : detInvMat_diff = abs(detInvMat - detInvMat_ref)
506 :
507 13 : assertion = all(MatInvMat_diff < tolerance) .and. detInvMat_diff < tolerance
508 :
509 1 : if (Test%isDebugMode .and. .not. assertion) then
510 : ! LCOV_EXCL_START
511 : write(Test%outputUnit,"(*(g0,:,', '))")
512 : write(Test%outputUnit,"(*(g0,:,', '))") "MatrixLU_ref = ", MatrixLU_ref
513 : write(Test%outputUnit,"(*(g0,:,', '))") "MatrixLU = ", MatrixLU
514 : write(Test%outputUnit,"(*(g0,:,', '))") "MatrixLU_diff = ", MatrixLU_diff
515 : write(Test%outputUnit,"(*(g0,:,', '))")
516 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_ref = ", MatInvMat_ref
517 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat = ", MatInvMat
518 : write(Test%outputUnit,"(*(g0,:,', '))") "MatInvMat_diff = ", MatInvMat_diff
519 : write(Test%outputUnit,"(*(g0,:,', '))")
520 : write(Test%outputUnit,"(*(g0,:,', '))") "detInvMat_ref = ", detInvMat_ref
521 : write(Test%outputUnit,"(*(g0,:,', '))") "detInvMat = ", detInvMat
522 : write(Test%outputUnit,"(*(g0,:,', '))") "detInvMat_diff = ", detInvMat_diff
523 : write(Test%outputUnit,"(*(g0,:,', '))")
524 : end if
525 : ! LCOV_EXCL_STOP
526 :
527 1 : end function test_getInvMatDet_1
528 :
529 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
530 :
531 1 : function test_getInvMat_1() result(assertion)
532 :
533 1 : use Constants_mod, only: IK, RK
534 : implicit none
535 :
536 : logical :: assertion
537 : integer(IK) , parameter :: nd = 3_IK
538 : real(RK) , parameter :: tolerance = 1.e-12_RK
539 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
540 : , 0._RK, 2._RK, 0._RK &
541 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
542 : real(RK) , parameter :: InverseMatrix_ref(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
543 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
544 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
545 : , shape = shape(InverseMatrix_ref) )
546 : real(RK) :: InverseMatrix(nd,nd)
547 : real(RK), allocatable :: InverseMatrix_diff(:,:)
548 :
549 1 : InverseMatrix = getInvMat(nd = nd, Matrix = PosDefMat)
550 :
551 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
552 1 : if (allocated(InverseMatrix_diff)) deallocate(InverseMatrix_diff); allocate(InverseMatrix_diff, mold = InverseMatrix)
553 :
554 14 : InverseMatrix_diff = abs(InverseMatrix - InverseMatrix_ref)
555 :
556 13 : assertion = all(InverseMatrix_diff < tolerance)
557 :
558 1 : if (Test%isDebugMode .and. .not. assertion) then
559 : ! LCOV_EXCL_START
560 : write(Test%outputUnit,"(*(g0,:,', '))")
561 : write(Test%outputUnit,"(*(g0,:,', '))") "InverseMatrix_ref = ", InverseMatrix_ref
562 : write(Test%outputUnit,"(*(g0,:,', '))") "InverseMatrix = ", InverseMatrix
563 : write(Test%outputUnit,"(*(g0,:,', '))") "InverseMatrix_diff = ", InverseMatrix_diff
564 : write(Test%outputUnit,"(*(g0,:,', '))")
565 : end if
566 : ! LCOV_EXCL_STOP
567 :
568 1 : end function test_getInvMat_1
569 :
570 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
571 :
572 1 : function test_multiplyMatrix_1() result(assertion)
573 :
574 1 : use Constants_mod, only: IK, RK
575 : implicit none
576 :
577 : logical :: assertion
578 : integer(IK) , parameter :: nd = 3_IK
579 : real(RK) , parameter :: tolerance = 1.e-12_RK
580 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
581 : , 0._RK, 2._RK, 0._RK &
582 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
583 : real(RK) :: MatrixProduct(nd,nd)
584 : real(RK), allocatable :: MatrixProduct_diff(:,:), MatrixProduct_ref(:,:)
585 :
586 13 : MatrixProduct_ref = matmul(PosDefMat,PosDefMat)
587 :
588 1 : call multiplyMatrix(A = PosDefMat, rowsA = nd, colsA = nd, B = PosDefMat, rowsB = nd, colsB = nd, C = MatrixProduct)
589 :
590 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
591 1 : if (allocated(MatrixProduct_diff)) deallocate(MatrixProduct_diff); allocate(MatrixProduct_diff, mold = MatrixProduct)
592 :
593 14 : MatrixProduct_diff = abs(MatrixProduct - MatrixProduct_ref)
594 :
595 13 : assertion = all(MatrixProduct_diff < tolerance)
596 :
597 1 : if (Test%isDebugMode .and. .not. assertion) then
598 : ! LCOV_EXCL_START
599 : write(Test%outputUnit,"(*(g0,:,', '))")
600 : write(Test%outputUnit,"(*(g0,:,', '))") "MatrixProduct_ref = ", MatrixProduct_ref
601 : write(Test%outputUnit,"(*(g0,:,', '))") "MatrixProduct = ", MatrixProduct
602 : write(Test%outputUnit,"(*(g0,:,', '))") "MatrixProduct_diff = ", MatrixProduct_diff
603 : write(Test%outputUnit,"(*(g0,:,', '))")
604 : end if
605 : ! LCOV_EXCL_STOP
606 :
607 1 : end function test_multiplyMatrix_1
608 :
609 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
610 :
611 1 : function test_getDeterminant_1() result(assertion)
612 :
613 1 : use Constants_mod, only: IK, RK
614 : implicit none
615 :
616 : logical :: assertion
617 : integer(IK) , parameter :: nd = 3_IK
618 : real(RK) , parameter :: tolerance = 1.e-12_RK
619 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ -1._RK, -0._RK, -1._RK &
620 : , -0._RK, -2._RK, -0._RK &
621 : , -1._RK, -0._RK, -3._RK ], shape = shape(PosDefMat) )
622 : real(RK) , parameter :: determinant_ref = -4._RK
623 1 : real(RK) :: determinant, determinant_diff
624 :
625 2 : determinant = getDeterminant(nd = nd, Matrix = PosDefMat)
626 :
627 1 : determinant_diff = abs(determinant - determinant_ref)
628 :
629 1 : assertion = determinant_diff < tolerance
630 :
631 1 : if (Test%isDebugMode .and. .not. assertion) then
632 : ! LCOV_EXCL_START
633 : write(Test%outputUnit,"(*(g0,:,', '))")
634 : write(Test%outputUnit,"(*(g0,:,', '))") "determinant_ref = ", determinant_ref
635 : write(Test%outputUnit,"(*(g0,:,', '))") "determinant = ", determinant
636 : write(Test%outputUnit,"(*(g0,:,', '))") "determinant_diff = ", determinant_diff
637 : write(Test%outputUnit,"(*(g0,:,', '))")
638 : end if
639 : ! LCOV_EXCL_STOP
640 :
641 1 : end function test_getDeterminant_1
642 :
643 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
644 :
645 1 : function test_getSqrtDetPosDefMat_1() result(assertion)
646 :
647 1 : use Constants_mod, only: IK, RK
648 : implicit none
649 :
650 : logical :: assertion
651 : integer(IK) , parameter :: nd = 3_IK
652 : real(RK) , parameter :: tolerance = 1.e-12_RK
653 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
654 : , 0._RK, 2._RK, 0._RK &
655 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
656 : real(RK) , parameter :: sqrtDetPosDefMat_ref = 2._RK
657 1 : real(RK) :: sqrtDetPosDefMat, sqrtDetPosDefMat_diff
658 :
659 1 : sqrtDetPosDefMat = getSqrtDetPosDefMat(nd = nd, PosDefMat = PosDefMat)
660 :
661 1 : sqrtDetPosDefMat_diff = abs(sqrtDetPosDefMat - sqrtDetPosDefMat_ref)
662 :
663 1 : assertion = sqrtDetPosDefMat_diff < tolerance
664 :
665 1 : if (Test%isDebugMode .and. .not. assertion) then
666 : ! LCOV_EXCL_START
667 : write(Test%outputUnit,"(*(g0,:,', '))")
668 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetPosDefMat_ref = ", sqrtDetPosDefMat_ref
669 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetPosDefMat = ", sqrtDetPosDefMat
670 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetPosDefMat_diff = ", sqrtDetPosDefMat_diff
671 : write(Test%outputUnit,"(*(g0,:,', '))")
672 : end if
673 : ! LCOV_EXCL_STOP
674 :
675 1 : end function test_getSqrtDetPosDefMat_1
676 :
677 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
678 :
679 1 : function test_getSqrtDetPosDefMat_2() result(assertion)
680 :
681 1 : use Constants_mod, only: IK, RK
682 : implicit none
683 :
684 : logical :: assertion
685 : integer(IK) , parameter :: nd = 3_IK
686 : real(RK) , parameter :: tolerance = 1.e-12_RK
687 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ -1._RK, -0._RK, -1._RK &
688 : , -0._RK, -2._RK, -0._RK &
689 : , -1._RK, -0._RK, -3._RK ], shape = shape(PosDefMat) )
690 1 : real(RK) :: sqrtDetPosDefMat
691 :
692 1 : sqrtDetPosDefMat = getSqrtDetPosDefMat(nd = nd, PosDefMat = PosDefMat)
693 :
694 1 : assertion = sqrtDetPosDefMat < 0._RK
695 :
696 1 : if (Test%isDebugMode .and. .not. assertion) then
697 : ! LCOV_EXCL_START
698 : write(Test%outputUnit,"(*(g0,:,', '))")
699 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetPosDefMat = ", sqrtDetPosDefMat
700 : write(Test%outputUnit,"(*(g0,:,', '))")
701 : end if
702 : ! LCOV_EXCL_STOP
703 :
704 1 : end function test_getSqrtDetPosDefMat_2
705 :
706 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
707 :
708 1 : function test_getLogSqrtDetPosDefMat_1() result(assertion)
709 :
710 1 : use Constants_mod, only: IK, RK
711 : implicit none
712 :
713 : logical :: assertion
714 : integer(IK) , parameter :: nd = 3_IK
715 : real(RK) , parameter :: tolerance = 1.e-12_RK
716 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
717 : , 0._RK, 2._RK, 0._RK &
718 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
719 : real(RK) , parameter :: logSqrtDetPosDefMat_ref = log(2._RK)
720 1 : real(RK) :: logSqrtDetPosDefMat, logSqrtDetPosDefMat_diff
721 : real(RK), allocatable :: Matrix(:,:)
722 : logical :: failed
723 :
724 13 : Matrix = PosDefMat
725 :
726 1 : call getLogSqrtDetPosDefMat(nd = nd, PosDefMat = Matrix, logSqrtDetPosDefMat = logSqrtDetPosDefMat, failed = failed)
727 :
728 1 : assertion = .not. failed
729 1 : if (.not. assertion) return
730 :
731 1 : logSqrtDetPosDefMat_diff = abs(logSqrtDetPosDefMat - logSqrtDetPosDefMat_ref)
732 :
733 1 : assertion = logSqrtDetPosDefMat_diff < tolerance
734 :
735 1 : if (Test%isDebugMode .and. .not. assertion) then
736 : ! LCOV_EXCL_START
737 : write(Test%outputUnit,"(*(g0,:,', '))")
738 : write(Test%outputUnit,"(*(g0,:,', '))") "logSqrtDetPosDefMat_ref = ", logSqrtDetPosDefMat_ref
739 : write(Test%outputUnit,"(*(g0,:,', '))") "logSqrtDetPosDefMat = ", logSqrtDetPosDefMat
740 : write(Test%outputUnit,"(*(g0,:,', '))") "logSqrtDetPosDefMat_diff = ", logSqrtDetPosDefMat_diff
741 : write(Test%outputUnit,"(*(g0,:,', '))")
742 : end if
743 : ! LCOV_EXCL_STOP
744 :
745 1 : end function test_getLogSqrtDetPosDefMat_1
746 :
747 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
748 :
749 1 : function test_getLogSqrtDetPosDefMat_2() result(assertion)
750 :
751 1 : use Constants_mod, only: IK, RK
752 : implicit none
753 :
754 : logical :: assertion
755 : integer(IK) , parameter :: nd = 3_IK
756 : real(RK) , parameter :: tolerance = 1.e-12_RK
757 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
758 : , 0._RK, 2._RK, 0._RK &
759 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
760 : real(RK) , parameter :: logSqrtDetPosDefMat_ref = log(2._RK)
761 1 : real(RK) :: logSqrtDetPosDefMat
762 : real(RK), allocatable :: Matrix(:,:)
763 : logical :: failed
764 :
765 13 : Matrix = -PosDefMat
766 :
767 1 : call getLogSqrtDetPosDefMat(nd = nd, PosDefMat = Matrix, logSqrtDetPosDefMat = logSqrtDetPosDefMat, failed = failed)
768 :
769 1 : assertion = failed
770 1 : if (.not. assertion) return
771 :
772 1 : end function test_getLogSqrtDetPosDefMat_2
773 :
774 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
775 :
776 1 : function test_isPosDef_1() result(assertion)
777 1 : use Constants_mod, only: IK, RK
778 : implicit none
779 : logical :: assertion
780 : integer(IK) , parameter :: nd = 3_IK
781 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
782 : , 0._RK, 2._RK, 0._RK &
783 : , 1._RK, 0._RK, 3._RK ], shape = shape(PosDefMat) )
784 :
785 1 : assertion = isPosDef(nd = nd, Matrix = PosDefMat)
786 2 : end function test_isPosDef_1
787 :
788 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
789 :
790 1 : function test_isPosDef_2() result(assertion)
791 1 : use Constants_mod, only: IK, RK
792 : implicit none
793 : logical :: assertion
794 : integer(IK) , parameter :: nd = 3_IK
795 : real(RK) , parameter :: PosDefMat(nd,nd) = reshape( [ -1._RK, -0._RK, -1._RK &
796 : , -0._RK, -2._RK, -0._RK &
797 : , -1._RK, -0._RK, -3._RK ], shape = shape(PosDefMat) )
798 :
799 1 : assertion = .not. isPosDef(nd = nd, Matrix = PosDefMat)
800 2 : end function test_isPosDef_2
801 :
802 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
803 :
804 1 : function test_symmetrizeUpperSquareMatrix_1() result(assertion)
805 1 : use Constants_mod, only: IK, RK
806 : implicit none
807 : logical :: assertion
808 : integer(IK) , parameter :: nd = 3_IK
809 : real(RK) , parameter :: UpperSquareMatrix_ref(nd,nd) = reshape( [ -1._RK, -0._RK, -1._RK &
810 : , -0._RK, -2._RK, -0._RK &
811 : , -1._RK, -0._RK, -3._RK ], shape = shape(UpperSquareMatrix_ref) )
812 : real(RK) :: UpperSquareMatrix(nd,nd)
813 : integer(IK) :: i, j
814 :
815 1 : UpperSquareMatrix = 0._RK
816 4 : do j = 1, nd
817 10 : do i = 1, j
818 9 : UpperSquareMatrix(i,j) = UpperSquareMatrix_ref(i,j)
819 : end do
820 : end do
821 :
822 1 : call symmetrizeUpperSquareMatrix(nd,UpperSquareMatrix)
823 13 : assertion = all( abs(UpperSquareMatrix-UpperSquareMatrix_ref) < 1.e-14_RK )
824 :
825 1 : if (Test%isDebugMode .and. .not. assertion) then
826 : ! LCOV_EXCL_START
827 : write(Test%outputUnit,"(*(g0,:,', '))")
828 : write(Test%outputUnit,"(*(g0,:,', '))") "UpperSquareMatrix_ref = ", UpperSquareMatrix_ref
829 : write(Test%outputUnit,"(*(g0,:,', '))") "UpperSquareMatrix = ", UpperSquareMatrix
830 : write(Test%outputUnit,"(*(g0,:,', '))")
831 : end if
832 : ! LCOV_EXCL_STOP
833 :
834 1 : end function test_symmetrizeUpperSquareMatrix_1
835 :
836 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
837 :
838 1 : function test_getOuterProd_1() result(assertion)
839 :
840 1 : use Constants_mod, only: IK, RK
841 : implicit none
842 :
843 : logical :: assertion
844 : real(RK) , parameter :: tolerance = 1.e-12_RK
845 : real(RK) , parameter :: Vector2(*) = [4._RK, 5._RK]
846 : real(RK) , parameter :: Vector1(*) = [1._RK, 2._RK, 3._RK]
847 : real(RK) , parameter :: OuterProduct_ref(3,2) = reshape( [ 4._RK, 8._RK, 12._RK &
848 : , 5._RK, 10._RK, 15._RK ], shape = shape(OuterProduct_ref) )
849 : real(RK), allocatable :: OuterProduct(:,:), OuterProduct_diff(:,:)
850 :
851 1 : OuterProduct = getOuterProd(Vector1 = Vector1, Vector2 = Vector2)
852 :
853 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
854 1 : if (allocated(OuterProduct_diff)) deallocate(OuterProduct_diff); allocate(OuterProduct_diff, mold = OuterProduct)
855 :
856 10 : OuterProduct_diff = abs(OuterProduct - OuterProduct_ref)
857 :
858 9 : assertion = all(OuterProduct_diff < tolerance)
859 :
860 1 : if (Test%isDebugMode .and. .not. assertion) then
861 : ! LCOV_EXCL_START
862 : write(Test%outputUnit,"(*(g0,:,', '))")
863 : write(Test%outputUnit,"(*(g0,:,', '))") "OuterProduct_ref = ", OuterProduct_ref
864 : write(Test%outputUnit,"(*(g0,:,', '))") "OuterProduct = ", OuterProduct
865 : write(Test%outputUnit,"(*(g0,:,', '))") "OuterProduct_diff = ", OuterProduct_diff
866 : write(Test%outputUnit,"(*(g0,:,', '))")
867 : end if
868 : ! LCOV_EXCL_STOP
869 :
870 1 : end function test_getOuterProd_1
871 :
872 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
873 :
874 1 : function test_sortPosDefMat_1() result(assertion)
875 :
876 1 : use Constants_mod, only: RK, IK
877 :
878 : implicit none
879 : logical :: assertion
880 : integer(IK), parameter :: ColIndx(*) = [4]
881 : integer(IK), parameter :: ColIndxMap(*) = [5]
882 : integer(IK) :: rank
883 : real(RK), allocatable :: PosDefMat(:,:), OutPosDefMat(:,:), OutPosDefMat_ref(:,:)
884 : integer(IK) :: i, j
885 :
886 1 : assertion = .true.
887 :
888 1 : rank = 5
889 1 : allocate(PosDefMat(rank,rank))
890 6 : do j = 1,rank
891 21 : do i = 1,j
892 20 : PosDefMat(i,j) = i*10 + j
893 : end do
894 : end do
895 :
896 : ! switch the variable 4 with 5, such that the output remains a positive-definite matrix.
897 :
898 1 : OutPosDefMat = sortPosDefMat(rank, PosDefMat, 1_IK, ColIndx, ColIndxMap)
899 :
900 : OutPosDefMat_ref = reshape( [ 11._RK, 12._RK, 13._RK, 15._RK, 14._RK &
901 : , 0._RK, 22._RK, 23._RK, 25._RK, 24._RK &
902 : , 0._RK, 0._RK, 33._RK, 35._RK, 34._RK &
903 : , 0._RK, 0._RK, 0._RK, 55._RK, 45._RK &
904 : , 0._RK, 0._RK, 0._RK, 0._RK, 44._RK ] &
905 3 : , shape = [rank,rank] )
906 62 : OutPosDefMat_ref = transpose(OutPosDefMat_ref)
907 6 : do j = 1, rank
908 21 : do i = 1, j
909 20 : assertion = assertion .and. OutPosDefMat_ref(i,j) == OutPosDefMat(i,j)
910 : end do
911 : end do
912 :
913 1 : if (Test%isDebugMode .and. .not. assertion) then
914 : ! LCOV_EXCL_START
915 :
916 : write(Test%outputUnit,"(*(g0))")
917 : write(Test%outputUnit,"(*(g0))") "OutPosDefMat_ref:"
918 : do i = 1,rank
919 : write(Test%outputUnit,"(*(F7.1))") (OutPosDefMat_ref(i,j),j=1,rank)
920 : end do
921 : write(Test%outputUnit,"(*(g0))")
922 :
923 : write(Test%outputUnit,"(*(g0))")
924 : write(Test%outputUnit,"(*(g0))") "Output Positive-Definite Matrix with variables 4 and 5 swapped:"
925 : do i = 1,rank
926 : write(Test%outputUnit,"(*(F7.1))") (OutPosDefMat(i,j),j=1,rank)
927 : end do
928 : write(Test%outputUnit,"(*(g0))")
929 :
930 : write(Test%outputUnit,"(*(g0))")
931 : write(Test%outputUnit,"(*(g0))") "Original Matrix:"
932 : do i = 1,rank
933 : write(Test%outputUnit,"(*(F7.1))") (PosDefMat(i,j),j=1,rank)
934 : end do
935 :
936 : end if
937 : ! LCOV_EXCL_STOP
938 :
939 1 : end function test_sortPosDefMat_1
940 :
941 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
942 :
943 1 : function test_getRegresCoef_1() result(assertion)
944 :
945 1 : use Constants_mod, only: RK, IK
946 :
947 : implicit none
948 : logical :: assertion
949 1 : real(RK) :: normalizedDifference
950 : integer , parameter :: rankPDM = 4, rankS11 = 3, rankS22 = 1
951 : real(RK), parameter :: PosDefMat(rankPDM,rankPDM) = reshape( &
952 : [ 4.414182620515998_RK, 1.173760167060120_RK, 0.757607629189287_RK, 5.075277296976230_RK &
953 : , 1.173760167060120_RK, 0.866956750091570_RK, 0.310654936099342_RK, 1.621274787164182_RK &
954 : , 0.757607629189287_RK, 0.310654936099342_RK, 0.955157221699132_RK, 1.254186231887444_RK &
955 : , 5.075277296976230_RK, 1.621274787164182_RK, 1.254186231887444_RK, 6.407791808961157_RK ] &
956 : , shape=shape(PosDefMat) )
957 : real(RK), parameter :: RegresCoefMatRef(rankS11,rankS22) = reshape( &
958 : [ 0.792047783119072 &
959 : , 0.253016145889269 &
960 : , 0.195728305363088 ] &
961 : , shape=shape(RegresCoefMatRef) )
962 : real(RK), parameter :: SchurComplementRef(rankS11,rankS11) = reshape( &
963 : [ 0.3943204887314210_RK, -0.110366933940115_RK, -0.235767795395625_RK &
964 : , -0.110366933940115_RK, 0.456748052015843_RK, -0.006674430520204_RK &
965 : , -0.235767795395625_RK, -0.006674430520204_RK, 0.709677475922085_RK ] &
966 : , shape=shape(SchurComplementRef) )
967 : real(RK) :: SchurComplement(rankS11,rankS11)
968 : real(RK) :: RegresCoefMat(rankS11,rankS22)
969 : integer :: i,j
970 :
971 1 : assertion = .true.
972 :
973 : call getRegresCoef ( rankPDM = rankPDM &
974 : , rankS11 = rankS11 &
975 : , rankS22 = rankS22 &
976 : , PosDefMat = PosDefMat &
977 : , RegresCoefMat = RegresCoefMat &
978 : , SchurComplement = SchurComplement &
979 1 : )
980 :
981 4 : do i = 1,rankS11
982 : !write(Test%outputUnit,"(*(F22.15))") (RegresCoefMat(i,j),j=1,rankS22), (RegresCoefMatRef(i,j),j=1,rankS22)
983 7 : do j = 1,rankS22
984 3 : normalizedDifference = abs(RegresCoefMatRef(i,j) - RegresCoefMat(i,j)) / (RegresCoefMatRef(i,j) + RegresCoefMat(i,j))
985 6 : assertion = assertion .and. normalizedDifference < 1.e-5_RK
986 : end do
987 : end do
988 :
989 1 : if (Test%isDebugMode .and. .not. assertion) then
990 : ! LCOV_EXCL_START
991 : write(Test%outputUnit,"(*(g0))")
992 : write(Test%outputUnit,"(*(g0))") "Original Covariance Matrix:"
993 : do i = 1,rankPDM
994 : write(Test%outputUnit,"(*(F22.15))") (PosDefMat(i,j),j=1,rankPDM)
995 : end do
996 : write(Test%outputUnit,"(*(g0))")
997 : write(Test%outputUnit,"(*(g0))") "RegresCoefMat, RegresCoefMatRef, difference, normalizedDifference:"
998 : write(Test%outputUnit,"(*(g0))")
999 : write(Test%outputUnit,"(*(g0))") "SchurComplement, SchurComplementRef, difference, normalizedDifference:"
1000 : write(Test%outputUnit,"(*(g0))")
1001 : do i = 1,rankS11
1002 : do j = 1,rankS22
1003 : write(Test%outputUnit,"(*(F22.15))" ) RegresCoefMat(i,j) &
1004 : , RegresCoefMatRef(i,j) &
1005 : , abs(RegresCoefMat(i,j)-RegresCoefMatRef(i,j)) &
1006 : , normalizedDifference
1007 : end do
1008 : end do
1009 : write(Test%outputUnit,"(*(g0))")
1010 : end if
1011 : ! LCOV_EXCL_STOP
1012 :
1013 4 : do i = 1,rankS11
1014 13 : do j = 1,rankS11
1015 9 : normalizedDifference = abs(SchurComplementRef(i,j) - SchurComplement(i,j)) / (SchurComplementRef(i,j) + SchurComplement(i,j))
1016 12 : assertion = assertion .and. normalizedDifference < 1.e-6_RK
1017 : end do
1018 : end do
1019 :
1020 1 : if (Test%isDebugMode .and. .not. assertion) then
1021 : ! LCOV_EXCL_START
1022 : do i = 1,rankS11
1023 : do j = 1,rankS11
1024 : write(Test%outputUnit,"(*(F22.15))" ) SchurComplement(i,j) &
1025 : , SchurComplementRef(i,j) &
1026 : , abs(SchurComplement(i,j)-SchurComplementRef(i,j)) &
1027 : , normalizedDifference
1028 : end do
1029 : end do
1030 : end if
1031 : ! LCOV_EXCL_STOP
1032 :
1033 1 : end function test_getRegresCoef_1
1034 :
1035 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1036 :
1037 : end module Test_Matrix_mod ! LCOV_EXCL_LINE
|
