Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains the implementations of procedures of [test_pm_complexDiv](@ref test_pm_complexDiv).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Wednesday 12:20 AM, October 13, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : use pm_val2str, only: getStr
28 : use pm_kind, only: SK, IK, LK
29 : use pm_distUnif, only: setUnifRand
30 : use pm_complexCompareAll, only: operator(<=)
31 : use pm_complexAbs, only: abs
32 :
33 : #if getDiv_ENABLED
34 :
35 : integer(IK) , parameter :: NP = 5_IK
36 : complex(CKC) , parameter :: SQRT_HUGE = sqrt(huge(0._CKC)), EPS = 100 * cmplx(epsilon(0._CKC), epsilon(0._CKC), CKC)
37 : complex(CKC) , parameter :: LOWER = -SQRT_HUGE, UPPER = +SQRT_HUGE, ZERO = (0._CKC, 0._CKC)
38 : complex(CKC) :: dividend(NP,NP), divisor(NP,NP), Quotient(NP,NP), Quotient_ref(NP,NP)
39 : integer(IK) :: i, j, itest
40 :
41 4 : assertion = .true._LK
42 :
43 44 : do itest = 1, 10
44 44 : call runTests()
45 : end do
46 :
47 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 :
49 : contains
50 :
51 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 :
53 40 : subroutine runTests()
54 :
55 1240 : call setUnifRand(dividend, LOWER, UPPER)
56 : do
57 1240 : call setUnifRand(divisor, LOWER, UPPER)
58 1240 : if (all(divisor /= ZERO)) exit
59 : end do
60 1240 : Quotient_ref = getDiv_ref(dividend, divisor)
61 :
62 240 : do i = 1, NP
63 1240 : do j = 1, NP
64 1200 : Quotient(i,j) = getDiv(dividend(i,j), divisor(i,j))
65 : end do
66 : end do
67 40 : call report(int(__LINE__, IK))
68 :
69 240 : do j = 1, NP
70 2240 : Quotient(:,j) = getDiv(dividend(:,j), divisor(:,j))
71 : end do
72 40 : call report(int(__LINE__, IK))
73 :
74 2440 : Quotient(:,:) = getDiv(dividend(:,:), divisor(:,:))
75 40 : call report(int(__LINE__, IK))
76 :
77 : ! Finite-range tests.
78 :
79 1240 : call setUnifRand(dividend)
80 : do
81 1240 : call setUnifRand(divisor)
82 1240 : if (all(divisor /= ZERO)) exit
83 : end do
84 1240 : Quotient_ref = dividend / divisor
85 2440 : Quotient(:,:) = getDiv(dividend(:,:), divisor(:,:))
86 40 : call report(int(__LINE__, IK))
87 :
88 40 : end subroutine
89 :
90 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 :
92 1000 : pure elemental function getDiv_ref(dividend, divisor) result(quotient)
93 : complex(CKC), intent(in) :: dividend, divisor
94 : complex(CKC) :: quotient
95 : real(CKC) :: r, d
96 1000 : if (abs(divisor%re) < abs(divisor%im)) then
97 0 : r = divisor%re / divisor%im
98 0 : d = divisor%im + r * divisor%re
99 0 : quotient%re = (dividend%re * r + dividend%im) / d
100 0 : quotient%im = (dividend%im * r - dividend%re) / d
101 : else
102 1000 : r = divisor%im / divisor%re
103 1000 : d = divisor%re + r * divisor%im
104 1000 : quotient%re = (dividend%re + dividend%im * r) / d
105 1000 : quotient%im = (dividend%im - dividend%re * r) / d
106 : end if
107 1000 : end function getDiv_ref
108 :
109 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110 :
111 160 : subroutine report(line)
112 : integer(IK), intent(in) :: line
113 : integer(IK) :: i, j
114 : complex(CKC) :: diff
115 960 : do i = 1, NP
116 4960 : do j = 1, NP
117 4000 : diff = abs(Quotient(i,j) - Quotient_ref(i,j))
118 4000 : assertion = assertion .and. diff <= EPS
119 4000 : if (test%traceable .and. .not. assertion) then
120 : ! LCOV_EXCL_START
121 : write(test%disp%unit,"(*(g0,:,', '))")
122 : write(test%disp%unit,"(*(g0,:,', '))") "i, j", i, j
123 : write(test%disp%unit,"(*(g0,:,', '))") "dividend(i,j)", dividend(i,j)
124 : write(test%disp%unit,"(*(g0,:,', '))") "divisor(i,j)", divisor(i,j)
125 : write(test%disp%unit,"(*(g0,:,', '))") "Quotient(i,j)", Quotient(i,j)
126 : write(test%disp%unit,"(*(g0,:,', '))") "Quotient_ref(i,j)", Quotient_ref(i,j)
127 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
128 : write(test%disp%unit,"(*(g0,:,', '))")
129 : ! LCOV_EXCL_STOP
130 : end if
131 4800 : call test%assert(assertion, SK_"The ratio of the complex values must be correctly computed without overflow.", line)
132 : end do
133 : end do
134 160 : end subroutine
135 :
136 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137 :
138 : #else
139 :
140 : #error "Unrecognized interface."
141 :
142 : #endif
143 :
144 : #undef COMPARE
|
