Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains procedure implementations of [pm_polynomial](@ref pm_polynomial).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Sunday 3:33 PM, September 19, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Define zero.
28 : #if CK_ENABLED || CK_CK_ENABLED
29 : #define TYPE_KIND complex(TKC)
30 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC)
31 : #elif RK_ENABLED || RK_RK_ENABLED || RK_CK_ENABLED
32 : #define TYPE_KIND real(TKC)
33 : real(TKC), parameter :: ZERO = 0._TKC, ONE = 1._TKC
34 : #else
35 : #error "Unrecognized interface."
36 : #endif
37 : #if setPolyRoot_ENABLED && Jen_ENABLED
38 : ! Bypass the non-existence of `SIND` and `COSD` for compilers (other than Intel/GNU): !(INTEL_ENABLED || GNU_ENABLED)
39 : #if 1
40 : #define COSD(x) cos(x * acos(real(-1, kind(x))) / real(180, kind(x)))
41 : #define SIND(x) sin(x * acos(real(-1, kind(x))) / real(180, kind(x)))
42 : #endif
43 : ! Enable robust calculations.
44 : ! \todo `hypot` argument still needs work.
45 : #if 0
46 : use pm_complexDiv, only: getDiv
47 : #define GET_RATIO(x, y) getDiv(x, y)
48 : #define GET_ABS(x) hypot(x%re, x%im)
49 : #else
50 : #define GET_RATIO(x, y) x / y
51 : #define GET_ABS(x) abs(x)
52 : #endif
53 : #endif
54 :
55 : !%%%%%%%%%%%%%%%%%
56 : #if getPolyVal_ENABLED
57 : !%%%%%%%%%%%%%%%%%
58 :
59 : integer(IK) :: i, degree
60 77 : degree = size(coef, 1, IK) - 1_IK
61 77 : CHECK_ASSERTION(__LINE__, 0_IK <= degree, SK_"@setPolyRoot(): The condition `0_IK < size(coef)` must hold. size(coef) = "//getStr(degree + 1_IK))
62 646 : poly = coef(degree)
63 664 : do i = degree - 1_IK, 0_IK, -1_IK
64 5875 : poly = poly * x + coef(i)
65 : end do
66 :
67 : !%%%%%%%%%%%%%%%%%
68 : #elif getPolyAdd_ENABLED
69 : !%%%%%%%%%%%%%%%%%
70 :
71 : integer(IK) :: lenLhs, lenRhs, lenMax
72 : lenLhs = size(lhs, 1, IK)
73 : lenRhs = size(rhs, 1, IK)
74 2 : do
75 24 : if (lenLhs == 0_IK) exit
76 18 : if (lhs(lenLhs) /= ZERO) exit
77 18 : lenLhs = lenLhs - 1_IK
78 : end do
79 0 : do
80 22 : if (lenRhs == 0_IK) exit
81 18 : if (rhs(lenRhs) /= ZERO) exit
82 18 : lenRhs = lenRhs - 1_IK
83 : end do
84 22 : if (0_IK < lenLhs .and. 0_IK < lenRhs) then
85 14 : lenMax = max(lenLhs, lenRhs)
86 14 : call setPolyAdd(add(1 : lenMax), lhs(1 : lenLhs), rhs(1 : lenRhs))
87 14 : add(lenMax + 1 :) = ZERO
88 8 : elseif (0_IK < lenLhs) then
89 6 : add(1 : lenLhs) = lhs(1 : lenLhs)
90 2 : add(lenLhs + 1:) = ZERO
91 : else
92 16 : add(1 : lenRhs) = rhs(1 : lenRhs)
93 6 : add(lenRhs + 1:) = ZERO
94 : end if
95 :
96 : !%%%%%%%%%%%%%%%%%
97 : #elif setPolyAdd_ENABLED
98 : !%%%%%%%%%%%%%%%%%
99 :
100 : integer(IK) :: lenRhs, lenLhs
101 22 : lenLhs = size(lhs, 1, IK)
102 22 : lenRhs = size(rhs, 1, IK)
103 22 : CHECK_ASSERTION(__LINE__, 0_IK < lenLhs, SK_"@setPolyAdd(): The condition `0 < size(lhs)` must hold. size(lhs) = "//getStr(lenLhs))
104 22 : CHECK_ASSERTION(__LINE__, 0_IK < lenRhs, SK_"@setPolyAdd(): The condition `0 < size(rhs)` must hold. size(rhs) = "//getStr(lenRhs))
105 88 : CHECK_ASSERTION(__LINE__, size(add, 1, IK) == max(lenLhs, lenRhs), \
106 : SK_"@setPolyAdd(): The condition `size(add) == max(lenLhs, lenRhs)` must hold. size(add), size(lhs), size(rhs) = "//\
107 : getStr([size(add, 1, IK), lenLhs, lenRhs]))
108 22 : CHECK_ASSERTION(__LINE__, lhs(lenLhs) /= ZERO, SK_"@setPolyAdd(): The condition `lhs(size(lhs)) /= 0.` must hold. lhs = "//getStr(lhs))
109 22 : CHECK_ASSERTION(__LINE__, rhs(lenRhs) /= ZERO, SK_"@setPolyAdd(): The condition `rhs(size(rhs)) /= 0.` must hold. rhs = "//getStr(rhs))
110 22 : if (lenLhs < lenRhs) then
111 14 : add(1 : lenLhs) = lhs(1 : lenLhs) + rhs(1 : lenLhs)
112 22 : add(lenLhs + 1 : lenRhs) = rhs(lenLhs + 1 : lenRhs)
113 : else
114 52 : add(1 : lenRhs) = lhs(1 : lenRhs) + rhs(1 : lenRhs)
115 36 : add(lenRhs + 1 : lenLhs) = lhs(lenRhs + 1 : lenLhs)
116 : end if
117 :
118 : !%%%%%%%%%%%%%%%%%
119 : #elif getPolySub_ENABLED
120 : !%%%%%%%%%%%%%%%%%
121 :
122 : integer(IK) :: lenLhs, lenRhs, lenMax
123 : lenLhs = size(lhs, 1, IK)
124 : lenRhs = size(rhs, 1, IK)
125 0 : do
126 14 : if (lenLhs == 0_IK) exit
127 10 : if (lhs(lenLhs) /= ZERO) exit
128 10 : lenLhs = lenLhs - 1_IK
129 : end do
130 0 : do
131 14 : if (lenRhs == 0_IK) exit
132 10 : if (rhs(lenRhs) /= ZERO) exit
133 10 : lenRhs = lenRhs - 1_IK
134 : end do
135 14 : if (0_IK < lenLhs .and. 0_IK < lenRhs) then
136 8 : lenMax = max(lenLhs, lenRhs)
137 8 : call setPolySub(sub(1 : lenMax), lhs(1 : lenLhs), rhs(1 : lenRhs))
138 8 : sub(lenMax + 1 :) = ZERO
139 6 : elseif (0_IK < lenLhs) then
140 6 : sub(1 : lenLhs) = lhs(1 : lenLhs)
141 2 : sub(lenLhs + 1:) = ZERO
142 : else
143 8 : sub(1 : lenRhs) = rhs(1 : lenRhs)
144 4 : sub(lenRhs + 1:) = ZERO
145 : end if
146 :
147 : !%%%%%%%%%%%%%%%%%
148 : #elif setPolySub_ENABLED
149 : !%%%%%%%%%%%%%%%%%
150 :
151 : integer(IK) :: lenRhs, lenLhs
152 16 : lenLhs = size(lhs, 1, IK)
153 16 : lenRhs = size(rhs, 1, IK)
154 16 : CHECK_ASSERTION(__LINE__, 0_IK < lenLhs, SK_"@setPolySub(): The condition `0 < size(lhs)` must hold. size(lhs) = "//getStr(lenLhs))
155 16 : CHECK_ASSERTION(__LINE__, 0_IK < lenRhs, SK_"@setPolySub(): The condition `0 < size(rhs)` must hold. size(rhs) = "//getStr(lenRhs))
156 64 : CHECK_ASSERTION(__LINE__, size(sub, 1, IK) == max(lenLhs, lenRhs), \
157 : SK_"@setPolySub(): The condition `size(sub) == max(lenLhs, lenRhs)` must hold. size(sub), size(lhs), size(rhs) = "//\
158 : getStr([size(sub, 1, IK), lenLhs, lenRhs]))
159 16 : CHECK_ASSERTION(__LINE__, lhs(lenLhs) /= ZERO, SK_"@setPolySub(): The condition `lhs(size(lhs)) /= 0.` must hold. lhs = "//getStr(lhs))
160 16 : CHECK_ASSERTION(__LINE__, rhs(lenRhs) /= ZERO, SK_"@setPolySub(): The condition `rhs(size(rhs)) /= 0.` must hold. rhs = "//getStr(rhs))
161 16 : if (lenLhs < lenRhs) then
162 0 : sub(1 : lenLhs) = lhs(1 : lenLhs) - rhs(1 : lenLhs)
163 0 : sub(lenLhs + 1 : lenRhs) = -rhs(lenLhs + 1 : lenRhs)
164 : else
165 52 : sub(1 : lenRhs) = lhs(1 : lenRhs) + rhs(1 : lenRhs)
166 36 : sub(lenRhs + 1 : lenLhs) = lhs(lenRhs + 1 : lenLhs)
167 : end if
168 :
169 : !%%%%%%%%%%%%%%%%%%
170 : #elif getPolyMul_ENABLED
171 : !%%%%%%%%%%%%%%%%%%
172 :
173 : integer(IK) :: lenLhs, lenRhs
174 : lenLhs = size(lhs, 1, IK)
175 : lenRhs = size(rhs, 1, IK)
176 0 : do
177 30 : if (lenLhs == 0_IK) exit
178 26 : if (lhs(lenLhs) /= ZERO) exit
179 26 : lenLhs = lenLhs - 1_IK
180 : end do
181 0 : do
182 30 : if (lenRhs == 0_IK) exit
183 26 : if (rhs(lenRhs) /= ZERO) exit
184 26 : lenRhs = lenRhs - 1_IK
185 : end do
186 30 : if (0_IK < lenLhs .and. 0_IK < lenRhs) then
187 24 : call setPolyMul(mul(1 : lenLhs + lenRhs - 1), lhs(1 : lenLhs), rhs(1 : lenRhs))
188 : else
189 6 : mul = ZERO
190 : end if
191 :
192 : !%%%%%%%%%%%%%%%%%
193 : #elif setPolyMul_ENABLED
194 : !%%%%%%%%%%%%%%%%%
195 :
196 : integer(IK) :: i, j
197 : integer(IK) :: lhsdeg, rhsdeg
198 32 : lhsdeg = size(lhs, 1, IK) - 1_IK
199 32 : rhsdeg = size(rhs, 1, IK) - 1_IK
200 32 : CHECK_ASSERTION(__LINE__, 0_IK <= lhsdeg, SK_"@setPolyMul(): The condition `0 < size(lhs)` must hold. size(lhs) = "//getStr(lhsdeg + 1_IK))
201 32 : CHECK_ASSERTION(__LINE__, 0_IK <= rhsdeg, SK_"@setPolyMul(): The condition `0 < size(rhs)` must hold. size(rhs) = "//getStr(rhsdeg + 1_IK))
202 128 : CHECK_ASSERTION(__LINE__, size(mul, 1, IK) == lhsdeg + rhsdeg + 1_IK, \
203 : SK_"@setPolyMul(): The condition `size(mul) == size(lhs) + size(rhs) - 1` must hold. size(mul), size(lhs), size(rhs) = "//\
204 : getStr([size(mul, 1, IK), lhsdeg + 1_IK, rhsdeg + 1_IK]))
205 32 : CHECK_ASSERTION(__LINE__, lhs(lhsdeg) /= ZERO, SK_"@setPolyMul(): The condition `lhs(size(lhs)) /= 0.` must hold. lhs = "//getStr(lhs))
206 32 : CHECK_ASSERTION(__LINE__, rhs(rhsdeg) /= ZERO, SK_"@setPolyMul(): The condition `rhs(size(rhs)) /= 0.` must hold. rhs = "//getStr(rhs))
207 168 : mul = ZERO
208 124 : do i = 0_IK, lhsdeg
209 340 : do j = 0_IK, rhsdeg
210 308 : mul(i + j) = mul(i + j) + lhs(i) * rhs(j)
211 : end do
212 : end do
213 :
214 : !%%%%%%%%%%%%%%%%%
215 : #elif setPolyDiv_ENABLED
216 : !%%%%%%%%%%%%%%%%%
217 :
218 : TYPE_KIND :: remainder, normfac
219 : TYPE_KIND, allocatable :: DividendTerm(:), DivisorTerm(:), QuotientTerm(:)
220 : integer(IK) :: i, degDividend, degDivisor, lenDivisor, lenDividend
221 8 : lenDividend = size(dividend, 1, IK)
222 8 : lenDivisor = size(divisor, 1, IK)
223 8 : CHECK_ASSERTION(__LINE__, 0_IK < lenDivisor, SK_"@setPolyDiv(): The condition `0 < size(divisor)` must hold. size(divisor) = "//getStr(lenDivisor))
224 8 : CHECK_ASSERTION(__LINE__, 0_IK < lenDividend, SK_"@setPolyDiv(): The condition `0 < size(dividend)` must hold. size(dividend) = "//getStr(lenDividend))
225 24 : CHECK_ASSERTION(__LINE__, size(quorem, 1, IK) == lenDividend, SK_"@setPolyDiv(): The condition `size(quorem) == size(dividend)` must hold. size(quorem), size(dividend) = "//getStr([size(quorem, 1, IK), lenDividend]))
226 8 : CHECK_ASSERTION(__LINE__, dividend(lenDividend) /= ZERO, SK_"@setPolyDiv(): The condition `dividend(size(dividend)) /= 0.` must hold. dividend = "//getStr(dividend))
227 8 : CHECK_ASSERTION(__LINE__, divisor(lenDivisor) /= ZERO, SK_"@setPolyDiv(): The condition `divisor(size(divisor)) /= 0.` must hold. divisor = "//getStr(divisor))
228 8 : if (lenDivisor <= lenDividend) then
229 8 : if (lenDivisor > 2_IK) then
230 18 : allocate(DivisorTerm(lenDividend), QuotientTerm(lenDividend), source = ZERO)
231 8 : DivisorTerm(1:lenDivisor) = divisor
232 2 : degDividend = lenDividend - 1_IK
233 2 : degDivisor = lenDivisor - 1_IK
234 12 : DividendTerm = dividend
235 2 : lenQuo = 0_IK
236 : do
237 40 : DivisorTerm = eoshift(DivisorTerm, degDivisor - degDividend)
238 4 : QuotientTerm(degDividend - degDivisor + 1) = DividendTerm(degDividend + 1) / DivisorTerm(degDividend + 1)
239 20 : DividendTerm = DividendTerm - DivisorTerm * QuotientTerm(degDividend - degDivisor + 1)
240 4 : lenQuo = max(lenQuo, degDividend - degDivisor)
241 : do
242 4 : degDividend = degDividend - 1_IK
243 4 : if (DividendTerm(degDividend + 1) /= ZERO) exit
244 : end do
245 16 : DivisorTerm(1 : lenDivisor) = divisor
246 8 : DivisorTerm(lenDivisor + 1 : lenDividend) = ZERO
247 4 : if (degDividend < degDivisor) exit
248 : end do
249 2 : lenQuo = lenQuo + 1_IK
250 6 : quorem(1 : lenQuo) = QuotientTerm(1 : lenQuo)
251 6 : quorem(lenQuo + 1 : lenDividend) = DividendTerm(1 : degDividend + 1)
252 : ! \bug gfortran as of version 12 fails to automatically deallocate heap allocations.
253 2 : deallocate(DividendTerm, DivisorTerm, QuotientTerm)
254 6 : elseif (lenDivisor == 2_IK) then
255 6 : if (divisor(2) /= ONE) then
256 3 : normfac = ONE / divisor(2)
257 3 : remainder = dividend(lenDividend) * normfac
258 11 : do i = lenDividend - 1_IK, 1_IK, -1_IK
259 8 : quorem(i) = remainder
260 11 : remainder = (dividend(i) - divisor(1) * remainder) * normfac
261 : end do
262 : else
263 3 : remainder = dividend(lenDividend)
264 11 : do i = lenDividend - 1_IK, 1_IK, -1_IK
265 8 : quorem(i) = remainder
266 11 : remainder = dividend(i) - divisor(1) * remainder
267 : end do
268 : end if
269 6 : quorem(lenDividend) = remainder
270 6 : lenQuo = lenDividend - 1_IK
271 : else ! lenDivisor == 1_IK
272 0 : normfac = ONE / divisor(1)
273 0 : quorem = dividend * normfac
274 : end if
275 : else
276 0 : lenQuo = 0_IK
277 0 : quorem = dividend
278 : end if
279 :
280 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 : #elif (getPolyDiff_ENABLED || setPolyDiff_ENABLED) && Def_ENABLED
282 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283 :
284 : integer(IK) :: i
285 : integer(IK), parameter :: order = 1
286 : #if setPolyDiff_ENABLED
287 0 : CHECK_ASSERTION(__LINE__, size(diff, 1, IK) == max(0_IK, size(coef, 1, IK) - order),\
288 : SK_"@setPolyDiff(): The condition `size(diff) == max(0, size(coef) - order)` must hold. size(diff), size(coef), order = "//\
289 : getStr([size(diff, 1, IK), size(coef, 1, IK), order]))
290 : #elif !getPolyDiff_ENABLED
291 : #error "Unrecognized interface."
292 : #endif
293 : do concurrent(i = order : size(coef, 1, IK) - 1_IK)
294 0 : diff(i) = coef(i) * i
295 : end do
296 :
297 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298 : #elif (getPolyDiff_ENABLED || setPolyDiff_ENABLED) && Ord_ENABLED
299 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
300 :
301 : integer(IK) :: i, lenCoef
302 40 : lenCoef = size(coef, 1, IK)
303 : #if setPolyDiff_ENABLED
304 160 : CHECK_ASSERTION(__LINE__, size(diff, 1, IK) == max(0_IK, lenCoef - order), \
305 : SK_"@setPolyDiff(): The condition `size(diff) == max(0, size(coef) - order)` must hold. size(diff), size(coef), order = "//\
306 : getStr([size(diff, 1, IK), lenCoef, order]))
307 : #elif !getPolyDiff_ENABLED
308 : #error "Unrecognized interface."
309 : #endif
310 80 : if (0_IK < order) then
311 35 : do concurrent(i = order : lenCoef - 1_IK)
312 201 : diff(i) = coef(i) * real(i, TKC)
313 : end do
314 : else
315 12 : CHECK_ASSERTION(__LINE__, 0_IK <= order, SK_"@setPolyDiff(): The condition `0 <= order` must hold. order = "//getStr(order))
316 37 : diff = coef
317 : end if
318 :
319 : !%%%%%%%%%%%%%%%%%
320 : #elif getPolyStr_ENABLED
321 : !%%%%%%%%%%%%%%%%%
322 :
323 : !> The maximum possible length of the string output from the functions under the generic interface `getStr`.
324 : !> Note that `NUMERIC_MAXLEN = 127` is very generous given that,
325 : !> + The maximum length of a complex of kind `real128` including signs, exponentiation, comma, and parentheses is `93` characters.
326 : !> `(-1.189731495357231765085759326628007016E+4932,-1.189731495357231765085759326628007016E+4932)`
327 : !> + The maximum length of a real of kind `real128` including signs and exponentiation is `44` characters.
328 : !> `-1.18973149535723176508575932662800702E+4932`
329 : !> + The maximum length of an integer of kind `int64` including signs `20` characters.
330 : !> `-9223372036854775807`
331 : character(127, SK) :: iomsg
332 : integer(IK) :: i, iostat, lenCoef, strlen
333 : integer(IK) , parameter :: NUMERIC_MAXLEN = 127_IK
334 : character(*,SKC), parameter :: PRODSYM = SKC_""
335 : character(*,SKC), parameter :: VARNAME = SKC_"x"
336 : character(*,SKC), parameter :: POWSYM = SKC_"^"
337 : character(*,SKC), parameter :: FIXED = PRODSYM//VARNAME//POWSYM
338 : #if CK_ENABLED
339 : #define GET_POLY_TERM(i) SKC_"(", coef(i)%re, SKC_",", coef(i)%im, SKC_")", FIXED, i
340 : #define POLY_TERM_CONST GET_POLY_TERM(0)
341 : #define GET_SIGNSYM(i) SKC_" + "
342 : #elif RK_ENABLED
343 : character(*,SKC), parameter :: SIGNSYM(-1:1) = [SKC_" - ", SKC_" ", SKC_" + "]
344 : #define POLY_TERM_CONST coef(0), FIXED, 0
345 : #define GET_POLY_TERM(i) abs(coef(i)), FIXED, i
346 : #define GET_SIGNSYM(i) SIGNSYM(int(sign(1.1_TKC, coef(i))))
347 : #else
348 : #error "Unrecognized interface."
349 : #endif
350 414 : lenCoef = size(coef, 1, IK)
351 414 : if (0 < lenCoef) then
352 333 : strlen = NUMERIC_MAXLEN * lenCoef
353 333 : allocate(character(strlen,SKC) :: str)
354 0 : loopExpandStr: do
355 1175 : write(str, "(ss,*(g0))", iostat = iostat, iomsg = iomsg) POLY_TERM_CONST, (GET_SIGNSYM(i), GET_POLY_TERM(i), i = 1_IK, lenCoef - 1_IK)
356 : !write(*,*); write(*,*) is_iostat_eor(iostat), iostat, iomsg; write(*,*)
357 333 : if (iostat == 0_IK) then
358 333 : str = getTrimmedTZ(str)
359 : return
360 0 : elseif (is_iostat_eor(iostat)) then
361 : call setResized(str) ! LCOV_EXCL_LINE
362 : cycle loopExpandStr ! LCOV_EXCL_LINE
363 : else
364 : error stop MODULE_NAME//& ! LCOV_EXCL_LINE
365 : SK_"@getPolyStr(): A runtime error occurred while converting the polynomial coefficients to a polynomial expression: "//trim(iomsg) ! LCOV_EXCL_LINE
366 : end if
367 : end do loopExpandStr
368 : else
369 81 : str = SKC_""
370 : end if
371 : #undef POLY_TERM_CONST
372 : #undef GET_POLY_TERM
373 : #undef GET_SIGNSYM
374 :
375 : !%%%%%%%%%%%%%%%%%%
376 : #elif getPolyRoot_ENABLED
377 : !%%%%%%%%%%%%%%%%%%
378 :
379 : integer(IK) :: count
380 : #if Def_ENABLED
381 : type(eigen_type), parameter :: method = eigen_type()
382 : #elif !(Eig_ENABLED || Jen_ENABLED || Lag_ENABLED)
383 : #error "Unrecognized interface."
384 : #endif
385 18 : allocate(root(size(coef, 1, IK) - 1_IK))
386 18 : call setPolyRoot(root, count, coef, method)
387 192 : root = root(1 : count)
388 :
389 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
390 : #elif setPolyRoot_ENABLED && Eig_ENABLED
391 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
392 :
393 : real(TKC), parameter :: EPS = epsilon(0._TKC), BASE = radix(1._TKC), BASESQ = BASE * BASE
394 : #if RK_CK_ENABLED
395 : #define GET_ABS(x) abs(x)
396 : real(TKC) :: p, q, t, w, x, y, zz
397 : integer(IK) :: k, l, m, ll, low, mm, mp2, na, niter, enm2
398 : logical(LK) :: notlas
399 : #elif CK_CK_ENABLED
400 : #define GET_ABS(x) abs(x%re) + abs(x%im)
401 : #else
402 : #error "Unrecognized interface."
403 : #endif
404 : integer(IK) :: i, j, degree
405 : logical(LK) :: normReductionEnabled
406 0 : TYPE_KIND :: workspace(size(root, 1, IK), size(root, 1, IK))
407 : real(TKC) :: r, s, c, scaleFac, g
408 29 : degree = size(coef, 1, IK) - 1_IK
409 29 : CHECK_ASSERTION(__LINE__, 0_IK < degree, SK_"@setPolyRoot(): The condition `1 < size(coef)` must hold. size(coef) = "//getStr(degree + 1_IK))
410 87 : CHECK_ASSERTION(__LINE__, size(root, 1, IK) == degree, SK_"@setPolyRoot(): The condition `size(root) == size(coef) - 1` must hold. size(root), size(coef) - 1 = "//getStr([size(root, 1, IK), degree]))
411 232 : CHECK_ASSERTION(__LINE__, all(shape(workspace, kind = IK) == degree), SK_"@setPolyRoot(): The condition `all(shape(workspace) == size(coef) - 1)` must hold. shape(workspace), size(coef) = "//getStr([shape(workspace, kind = IK), degree + 1_IK]))
412 29 : CHECK_ASSERTION(__LINE__, coef(degree) /= ZERO, SK_"@setPolyRoot(): The condition `coef(size(coef)) /= 0.` must hold. degree = "//getStr(coef(degree)))
413 :
414 : ! Gracefully capture the error and return if the highest polynomial coefficient is zero.
415 :
416 29 : if (coef(degree) == ZERO) then
417 0 : count = 0_IK ! degree + 1_IK
418 0 : return
419 : end if
420 :
421 : ! Build the count row of the upper Hessenberg Companion matrix.
422 :
423 : do concurrent(i = 1_IK : degree)
424 238 : workspace(1, i) = -coef(degree - i) / coef(degree)
425 : end do
426 :
427 : ! Extract any exact ZERO roots and set degree = degree of remaining polynomial.
428 :
429 29 : do j = degree, 1_IK, -1_IK
430 29 : if (workspace(1, j) /= ZERO) exit
431 0 : root(j) = cmplx(ZERO, kind = TKC) ! xxx
432 29 : degree = degree - 1_IK
433 : end do
434 :
435 : ! The special case of `degree == 0` or `degree == 1`.
436 :
437 29 : if (degree == 0_IK) then ! Gracefully capture the error and return if the degree of polynomial is zero.
438 0 : count = 0_IK
439 0 : return
440 29 : elseif (degree == 1_IK) then
441 0 : root(1) = workspace(1,1)
442 0 : count = 1_IK
443 0 : return
444 : end if
445 :
446 : ! Build rows 2 through `degree` of the Companion matrix.
447 :
448 209 : do i = 2_IK, degree
449 1782 : do j = 1_IK, degree
450 1782 : workspace(i, j) = ZERO
451 : end do
452 209 : workspace(i, i - 1) = ONE
453 : end do
454 :
455 : ! Balance the matrix.
456 : ! This is an adaption of the EISPACK subroutine `balanc` to the special case of a Companion matrix.
457 : ! The EISPACK balance is a translation of the Algol procedure balance, num. math. 13, 293-304(1969) by parlett and reinsch. handbook for auto. comp., vol.ii-linear algebra, 315-326(1971).
458 :
459 : ! Iterative loop for norm reduction
460 :
461 : loopNormReduction: do
462 :
463 : normReductionEnabled = .false.
464 :
465 836 : do i = 1_IK, degree
466 :
467 : ! r = sum of magnitudes in row i skipping diagonal.
468 : ! c = sum of magnitudes in col i skipping diagonal.
469 :
470 744 : if (i == 1_IK) then
471 92 : r = GET_ABS(workspace(1, 2))
472 652 : do j = 3, degree
473 652 : r = r + GET_ABS(workspace(1, j))
474 : end do
475 92 : c = GET_ABS(workspace(2, 1))
476 : else
477 652 : r = GET_ABS(workspace(i, i - 1))
478 652 : c = GET_ABS(workspace(1, i))
479 652 : if (i /= degree) c = c + GET_ABS(workspace(i + 1, i))
480 : end if
481 :
482 : ! Determine column scale factor, `scaleFac`.
483 :
484 744 : g = r / BASE
485 264 : scaleFac = real(ONE, TKC)
486 744 : s = c + r
487 468 : do
488 1212 : if (c < g) then
489 468 : scaleFac = scaleFac * BASE
490 468 : c = c * BASESQ
491 : cycle
492 : end if
493 : exit
494 : end do
495 744 : g = r * BASE
496 628 : do
497 1372 : if (c >= g) then
498 628 : scaleFac = scaleFac / BASE
499 628 : c = c / BASESQ
500 : cycle
501 : end if
502 : exit
503 : end do
504 :
505 : ! Will the factor `scaleFac` have a significant effect?
506 :
507 836 : if ((c + r) / scaleFac < .95_TKC * s) then ! yes, so do the scaling.
508 :
509 312 : g = real(ONE, TKC) / scaleFac
510 : normReductionEnabled = .true.
511 :
512 : ! scale row `i`
513 :
514 312 : if (i == 1_IK) then
515 324 : do j = 1_IK, degree
516 324 : workspace(1, j) = workspace(1, j) * g
517 : end do
518 : else
519 274 : workspace(i, i - 1) = workspace(i, i - 1) * g
520 : end if
521 :
522 : ! Scale column `i`.
523 :
524 312 : workspace(1, i) = workspace(1, i) * scaleFac
525 312 : if (i /= degree) workspace(i + 1, i) = workspace(i + 1, i) * scaleFac
526 :
527 : end if
528 :
529 : end do
530 :
531 92 : if (normReductionEnabled) cycle loopNormReduction
532 63 : exit loopNormReduction
533 :
534 : end do loopNormReduction
535 :
536 : #if RK_CK_ENABLED
537 : ! ***************** QR eigenvalue algorithm ***********************
538 : ! This is the EISPACK subroutine hqr that uses the QR
539 : ! algorithm to compute all eigenvalues of an upper
540 : ! hessenberg matrix. original algol code was due to Martin,
541 : ! Peters, and Wilkinson, numer. math., 14, 219-231(1970).
542 10 : count = degree
543 : low = 1_IK
544 3 : t = ZERO
545 : ! ********** search for next eigenvalues **********
546 : loopNextEigen: do
547 67 : if (count < low) exit loopNextEigen
548 : niter = 0_IK
549 54 : na = count - 1_IK
550 54 : enm2 = na - 1_IK
551 : ! ********** look for single small sub-diagonal element for `l = count step -1` until `low` do -- **********
552 : do
553 1191 : do ll = low, count
554 1191 : l = count + low - ll
555 1191 : if (l == low) exit
556 1191 : if (abs(workspace(l, l - 1)) <= EPS * (abs(workspace(l - 1, l - 1)) + abs(workspace(l, l)))) exit
557 : end do
558 : ! ********** form shift **********
559 261 : x = workspace(count, count)
560 261 : if (l == count) exit
561 227 : y = workspace(na,na)
562 227 : w = workspace(count,na) * workspace(na,count)
563 227 : if (l == na) then ! two roots found
564 20 : p = (y - x) * .5_TKC
565 20 : q = p * p + w
566 20 : zz = sqrt(abs(q))
567 20 : x = x + t
568 20 : if (q < ZERO) then ! complex pair
569 15 : root(na) = cmplx(x + p, zz, TKC)
570 15 : root(count) = cmplx(x + p, -zz, TKC)
571 : else ! pair of reals
572 5 : zz = p + sign(zz, p)
573 5 : root(na) = cmplx(x + zz, ZERO, TKC)
574 5 : root(count) = root(na)
575 5 : if (zz /= ZERO) root(count) = cmplx(x - w / zz, ZERO, TKC)
576 : end if
577 20 : count = enm2
578 20 : cycle loopNextEigen
579 : end if
580 : ! increased iterations for quad-precision.
581 : !if (niter == 30_IK) exit loopNextEigen ! failed. No convergence to an eigenvalue after 30 iterations.
582 207 : if (niter == 90_IK) exit loopNextEigen ! failed. No convergence to an eigenvalue after 30 iterations.
583 207 : if (niter == 10_IK .or. niter == 20_IK) then ! form exceptional shift.
584 7 : t = t + x
585 42 : do i = low, count
586 42 : workspace(i,i) = workspace(i,i) - x
587 : end do
588 7 : s = abs(workspace(count, na)) + abs(workspace(na, enm2))
589 7 : x = .75_TKC * s
590 2 : y = x
591 7 : w = -.4375_TKC * s * s
592 : end if
593 207 : niter = niter + 1_IK
594 : ! Look for two consecutive small sub-diagonal elements. for `m = count - 2 step -1` until l do.
595 647 : do mm = l, enm2
596 647 : m = enm2 + l - mm
597 647 : zz = workspace(m, m)
598 647 : r = x - zz
599 647 : s = y - zz
600 647 : p = (r * s - w) / workspace(m + 1,m) + workspace(m, m + 1)
601 647 : q = workspace(m + 1, m + 1) - zz - r - s
602 647 : r = workspace(m + 2, m + 1)
603 647 : s = abs(p) + abs(q) + abs(r)
604 647 : p = p / s
605 647 : q = q / s
606 647 : r = r / s
607 647 : if (m == l) exit
608 647 : if (abs(workspace(m, m - 1)) * (abs(q) + abs(r)) <= EPS * abs(p) * (abs(workspace(m - 1, m - 1)) + abs(zz) + abs(workspace(m + 1, m + 1)))) exit
609 : end do
610 207 : mp2 = m + 2_IK
611 854 : do i = mp2, count
612 647 : workspace(i, i - 2) = ZERO
613 647 : if (i == mp2) cycle
614 854 : workspace(i, i - 3) = ZERO
615 : end do
616 : ! Double QR step involving rows l to count and columns m to count.
617 1095 : do k = m, na
618 : notlas = k /= na
619 854 : if (k /= m) then
620 647 : p = workspace(k, k - 1)
621 647 : q = workspace(k + 1, k - 1)
622 260 : r = ZERO
623 647 : if (notlas) r = workspace(k + 2, k - 1)
624 647 : x = abs(p) + abs(q) + abs(r)
625 647 : if (x == ZERO) cycle
626 647 : p = p / x
627 647 : q = q / x
628 647 : r = r / x
629 : end if
630 854 : s = sign(sqrt(p * p + q * q + r * r), p)
631 854 : if (k /= m) then
632 647 : workspace(k, k - 1) = -s * x
633 207 : elseif (l /= m) then
634 53 : workspace(k, k - 1) = -workspace(k, k - 1)
635 : end if
636 854 : p = p + s
637 854 : x = p / s
638 854 : y = q / s
639 854 : zz = r / s
640 854 : q = q / p
641 508 : r = r / p
642 : ! row modification.
643 4287 : do j = k, count
644 3433 : p = workspace(k, j) + q * workspace(k + 1, j)
645 3433 : if (notlas) then
646 3019 : p = p + r * workspace(k + 2, j)
647 3019 : workspace(k + 2, j) = workspace(k + 2, j) - p * zz
648 : end if
649 3433 : workspace(k + 1, j) = workspace(k + 1, j) - p * y
650 4287 : workspace(k, j) = workspace(k, j) - p * x
651 : end do
652 854 : j = min(count, k + 3)
653 : ! column modification.
654 5742 : do i = l, j
655 4681 : p = x * workspace(i,k) + y * workspace(i, k + 1)
656 4681 : if (notlas) then
657 3564 : p = p + zz * workspace(i, k + 2)
658 3564 : workspace(i, k + 2) = workspace(i, k + 2) - p * r
659 : end if
660 4681 : workspace(i, k + 1) = workspace(i, k + 1) - p * q
661 5535 : workspace(i,k) = workspace(i,k) - p
662 : end do
663 : end do
664 : end do
665 : ! ********** ONE root found **********
666 34 : root(count) = cmplx(x + t, ZERO, TKC)
667 34 : count = na
668 : end do loopNextEigen
669 13 : count = count + 1_IK
670 : #elif CK_CK_ENABLED
671 16 : call dcomqr(1_IK, degree, workspace(:, 1 : degree), root(1 : degree), count)
672 : #else
673 : #error "Unrecognized interface."
674 : #endif
675 : ! This block reverses the ordering of the identified roots, such that `root(1:count)` contains it all.
676 : ! \todo
677 : ! This ordering reversal could be merged with the rest of the algorithm to avoid the extra final reversal copy below.
678 : ! However, given that most polynomials are of low degree, the final reversal copy cost should generally be quite negligible.
679 : ! Because of the complexity of the rest of the algorithm, a final copy was the preferred approach over merging the reversal with the algorithm, for now.
680 : block
681 : complex(TKC) :: temp
682 : integer(IK) :: lenRoot
683 : lenRoot = size(root, 1, IK)
684 238 : do i = lenRoot, count, -1_IK
685 209 : temp = root(i)
686 209 : root(i) = root(lenRoot - i + 1_IK)
687 238 : root(lenRoot - i + 1_IK) = temp
688 : end do
689 29 : count = lenRoot - count + 1_IK
690 : end block
691 : #if CK_CK_ENABLED
692 : contains
693 : ! This subroutine is a translation of a unitary analogue of the algol
694 : ! procedure comlr, num. math. 12, 369-376(1968) by Martin and Wilkinson.
695 : ! Handbook for auto. comp., vol.ii-linear algebra, 396-403(1971).
696 : ! the unitary analogue substitutes the QR algorithm of Francis
697 : ! (comp. jour. 4, 332-345(1962)) for the LR algorithm.
698 : ! This subroutine finds the eigenvalues of a complex
699 : ! upper hessenberg matrix by the QR method.
700 : ! `low` and `igh` are integers determined by the balancing subroutine `cbal`.
701 : ! if `cbal` has not been used, set `low = 1, igh = order`.
702 : ! `hessen` contains the complex upper hessenberg matrix.
703 : ! The lower triangle of `hessen` below the subdiagonal contains
704 : ! information about the unitary transformations used in the reduction by `corth`, if performed.
705 : ! ON OUTPUT:
706 : ! The upper hessenberg portions of `hessen` are destroyed.
707 : ! Therefore, they must be saved before calling `comqr`
708 : ! if subsequent calculation of eigenvectors is to be performed.
709 : !
710 : ! `root` contains the eigenvalues.
711 : ! If an error exit is made, the eigenvalues should be correct for indices `first+1,...,order`,
712 : ! `first` is set to,
713 : ! zero for normal return,
714 : ! j if the j-th eigenvalue has not been
715 : ! determined after 30 iterations.
716 : !
717 : ! arithmetic is real except for the replacement of the algol
718 : ! procedure cdiv by complex division and use of the subroutines
719 : ! sqrt and cmplx in computing complex square roots.
720 16 : pure subroutine dcomqr(low, igh, hessen, eigen, first)
721 : complex(TKC) , intent(out) :: eigen(:)
722 : complex(TKC) , intent(inout) :: hessen(:,:) ! hessen%im(nm, order),hessen%re(nm,order)
723 : integer(IK) , intent(in) :: low, igh
724 : integer(IK) , intent(out) :: first
725 : complex(TKC) :: s, t, x, y, z3, zz
726 : integer(IK) :: enm1, i, niter, j, l, ll, lp1, order
727 : real(TKC) :: norm
728 16 : order = size(eigen, 1, IK) ! The order of the Hessenberg matrix.
729 16 : if (low /= igh) then ! 180
730 : ! create real subdiagonal elements
731 16 : l = low + 1_IK
732 135 : do i = l, igh
733 135 : if (hessen(i, i - 1)%im /= 0._TKC) then
734 0 : norm = abs(hessen(i, i - 1))
735 0 : y = hessen(i, i - 1) / norm
736 0 : hessen(i, i - 1) = cmplx(norm, 0._TKC, TKC)
737 0 : do j = i, igh
738 0 : s%im = y%re * hessen(i, j)%im - y%im * hessen(i, j)%re
739 0 : hessen(i, j)%re = y%re * hessen(i, j)%re + y%im * hessen(i, j)%im
740 0 : hessen(i, j)%im = s%im
741 : end do
742 0 : do j = low, min(i + 1_IK, igh)
743 0 : s%im = y%re * hessen(j, i)%im + y%im * hessen(j, i)%re
744 0 : hessen(j, i)%re = y%re * hessen(j, i)%re - y%im * hessen(j, i)%im
745 0 : hessen(j, i)%im = s%im
746 : end do
747 : end if
748 : end do
749 : end if
750 : ! Store roots isolated by cbal.
751 151 : do i = 1_IK, order
752 135 : if (low <= i .and. i <= igh) cycle
753 151 : eigen(i) = hessen(i, i)
754 : end do
755 16 : first = igh
756 11 : t = (0._TKC, 0._TKC)
757 : ! Search for next eigenvalue.
758 135 : loopNextEigen: do ! 220
759 151 : if (first < low) then ! convergence achieved.
760 16 : first = low
761 16 : return
762 : end if
763 : niter = 0_IK
764 135 : enm1 = first - 1_IK
765 : ! Look for single small sub-diagonal element for `l = first` step `-1` until `low`.
766 : loop240: do
767 3514 : do ll = low, first
768 3514 : l = first + low - ll
769 3514 : if (l == low) exit
770 : if (abs(hessen(l, l - 1)%re) <= EPS * & ! LCOV_EXCL_LINE
771 : ( abs(hessen(l - 1, l - 1)%re) & ! LCOV_EXCL_LINE
772 : + abs(hessen(l - 1, l - 1)%im) & ! LCOV_EXCL_LINE
773 : + abs(hessen(l, l)%re) & ! LCOV_EXCL_LINE
774 : + abs(hessen(l, l)%im) & ! LCOV_EXCL_LINE
775 661 : ) ) exit
776 : end do
777 : ! Form shift. 300
778 661 : if (l == first) then ! 660 : a root found.
779 : ! A root found ! 660
780 135 : eigen(first) = hessen(first, first) + t
781 135 : first = enm1
782 : cycle loopNextEigen
783 : end if
784 : ! amir: increased maximum iteration to allow for quad-precision convergence.
785 526 : if (niter == 90_IK) then ! No convergence to an eigenvalue after 30 iterations.
786 0 : first = first + 1_IK
787 0 : return
788 : end if
789 526 : if (niter == 10_IK .or. niter == 20_IK) then ! 320 Form exceptional shift.
790 11 : s = cmplx(abs(hessen(first, enm1)%re) + abs(hessen(enm1, first - 2)%re), 0._TKC, TKC)
791 : else
792 515 : s = hessen(first, first)
793 515 : x%re = hessen(enm1, first)%re * hessen(first, enm1)%re
794 515 : x%im = hessen(enm1, first)%im * hessen(first, enm1)%re
795 515 : if (x%re /= 0._TKC .or. x%im /= 0._TKC) then ! 340 ! xx z3 and zz can be replaced with a single variable `z`.
796 499 : y = 0.5_TKC * (hessen(enm1, enm1) - s)
797 499 : z3 = sqrt(cmplx(y%re**2 - y%im**2 + x%re, 2._TKC * y%re * y%im + x%im, TKC))
798 393 : zz = z3
799 499 : if (y%re * zz%re + y%im * zz%im < 0._TKC) zz = -zz ! 310
800 499 : z3 = x / (y + zz)
801 499 : s = s - z3
802 : end if
803 : end if
804 4402 : do i = low, first ! 340
805 4402 : hessen(i,i) = hessen(i,i) - s
806 : end do
807 526 : t = t + s
808 526 : niter = niter + 1_IK
809 : ! Reduce to triangle (rows).
810 526 : lp1 = l + 1_IK
811 3379 : do i = lp1, first ! 500
812 2853 : s%re = hessen(i, i - 1)%re
813 2853 : hessen(i, i - 1)%re = 0._TKC
814 2853 : norm = sqrt(s%re**2 + hessen(i - 1, i - 1)%re**2 + hessen(i - 1, i - 1)%im**2)
815 2853 : x = hessen(i - 1, i - 1) / norm
816 2853 : eigen(i - 1) = x
817 2853 : hessen(i - 1, i - 1) = cmplx(norm, 0._TKC, TKC)
818 2853 : hessen(i, i - 1)%im = s%re / norm
819 15106 : do j = i, first ! 490
820 11727 : y = hessen(i - 1, j)
821 11727 : zz = hessen(i, j)
822 11727 : hessen(i - 1, j)%re = x%re * y%re + x%im * y%im + hessen(i, i - 1)%im * zz%re
823 11727 : hessen(i - 1, j)%im = x%re * y%im - x%im * y%re + hessen(i, i - 1)%im * zz%im
824 11727 : hessen(i, j)%re = x%re * zz%re - x%im * zz%im - hessen(i, i - 1)%im * y%re
825 14580 : hessen(i, j)%im = x%re * zz%im + x%im * zz%re - hessen(i, i - 1)%im * y%im
826 : end do
827 : end do
828 526 : s%im = hessen(first, first)%im
829 526 : if (s%im /= 0._TKC) then
830 503 : norm = abs(cmplx(hessen(first, first)%re, s%im, TKC))
831 503 : s%re = hessen(first, first)%re / norm
832 503 : s%im = s%im / norm
833 503 : hessen(first, first) = cmplx(norm, 0._TKC, TKC)
834 : end if
835 : ! Inverse operation (columns).
836 3379 : do j = lp1, first ! 540
837 2853 : x = eigen(j - 1)
838 17959 : do i = l, j
839 14580 : y = cmplx(hessen(i, j - 1)%re, 0._TKC, TKC)
840 14580 : zz = hessen(i, j)
841 14580 : if (i /= j) then ! 560
842 11727 : y%im = hessen(i, j - 1)%im
843 11727 : hessen(i, j - 1)%im = x%re * y%im + x%im * y%re + hessen(j, j - 1)%im * zz%im
844 : end if
845 14580 : hessen(i, j - 1)%re = x%re * y%re - x%im * y%im + hessen(j, j - 1)%im * zz%re
846 14580 : hessen(i, j)%re = x%re * zz%re + x%im * zz%im - hessen(j, j - 1)%im * y%re
847 20197 : hessen(i, j)%im = x%re * zz%im - x%im * zz%re - hessen(j, j - 1)%im * y%im
848 : end do
849 : end do
850 526 : if (s%im == 0._TKC) cycle loop240
851 3745 : do i = l, first
852 3242 : y = hessen(i,first)
853 3242 : hessen(i,first)%re = s%re * y%re - s%im * y%im
854 3768 : hessen(i,first)%im = s%re * y%im + s%im * y%re
855 : end do
856 : end do loop240
857 : end do loopNextEigen
858 : end subroutine
859 : #endif
860 : #undef GET_ABS
861 :
862 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
863 : #elif setPolyRoot_ENABLED && Jen_ENABLED && CK_CK_ENABLED
864 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
865 :
866 : ! the program has been written to reduce the chance of overflow occurring.
867 : ! if it does occur, there is still a possibility that the zerofinder will work
868 : ! provided the overflowed quantity is replaced by a large number.
869 :
870 : ! Global variables.
871 : real(TKC) , parameter :: BASE = radix(0._TKC), EPS = epsilon(0._TKC), TIN = tiny(0._TKC), INF = huge(0._TKC)
872 : real(TKC) , parameter :: SINR = SIND(94._TKC), COSR = COSD(94._TKC) ! rotation by 94 degrees.
873 : real(TKC) , parameter :: MRE = 2._TKC * sqrt(2._TKC) * EPS ! error bound on complex multiplication.
874 : real(TKC) , parameter :: ARE = EPS ! error bound on complex addition.
875 : real(TKC) , parameter :: HALF_SQRT2 = 0.5_TKC * sqrt(2._TKC) ! cos(45d)
876 : complex(TKC), allocatable :: h(:), qp(:), qh(:), sh(:), workspace(:)
877 : complex(TKC) :: s, t, pv
878 : integer(IK) :: degree
879 : integer(IK) :: nn
880 :
881 : ! Local variables.
882 : complex(TKC) :: temp
883 : real(TKC) , parameter :: HIGH = sqrt(INF), LOW = TIN / EPS, INV_LOG_BASE = 1._TKC / log(BASE)
884 : real(TKC) , parameter :: NEG_HALF_INV_LOG_BASE = -.5_TKC * INV_LOG_BASE
885 : real(TKC) :: xx, yy, xxx, bnd, max, min
886 : real(TKC) :: x, xm, f, dx, df, xni ! cauchy, noshft routines
887 : integer(IK) :: cnt1, cnt2, i, n, j, jj, nm1
888 : logical(LK) :: converged
889 :
890 14 : nn = size(coef, 1, IK)
891 14 : degree = nn - 1_IK
892 5 : xx = +HALF_SQRT2
893 5 : yy = -HALF_SQRT2
894 :
895 14 : CHECK_ASSERTION(__LINE__, 1_IK < nn, SK_"@setPolyRoot(): The condition `1 < size(coef)` must hold. size(coef) = "//getStr(nn))
896 14 : CHECK_ASSERTION(__LINE__, coef(degree) /= ZERO, SK_"@setPolyRoot(): The condition `coef(size(coef)) /= 0.` must hold. coef = "//getStr(coef))
897 42 : CHECK_ASSERTION(__LINE__, size(root, 1, IK) == degree, SK_"@setPolyRoot(): The condition `size(root) == size(coef) - 1` must hold. size(root), size(coef) - 1 = "//getStr([size(root, 1, IK), degree]))
898 :
899 : ! Algorithm fails if the leading coefficient is zero.
900 14 : if (coef(degree) == ZERO) then
901 0 : count = 0_IK
902 0 : return
903 : end if
904 :
905 : ! Allocate arrays.
906 : #if __GFORTRAN__
907 14 : if (allocated(workspace)) deallocate(workspace)
908 14 : if (allocated(qp)) deallocate(qp)
909 14 : if (allocated(qh)) deallocate(qh)
910 14 : if (allocated(sh)) deallocate(sh)
911 14 : if (allocated(h)) deallocate(h)
912 : #endif
913 14 : allocate(h(nn), qp(nn), qh(nn), sh(nn), workspace(nn)) ! xxx this should be moved down once cleaning is complete.
914 :
915 : ! Extract any exact ZERO roots and set count = degree of remaining polynomial.
916 14 : do count = 1_IK, degree ! 10
917 14 : if (coef(count - 1) /= ZERO) exit
918 0 : root(count) = ZERO
919 14 : nn = nn - 1_IK
920 : end do
921 :
922 : ! Make a copy of the coefficients.
923 155 : do i = 1_IK, nn
924 141 : workspace(i) = coef(nn - i)
925 155 : sh(i)%re = GET_ABS(workspace(i)) ! modulus of coefficients of `p`.
926 : end do
927 :
928 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
929 : ! Compute a scale factor to multiply the coefficients of the polynomial.
930 : ! The scaling is done to avoid overflow and to avoid undetected underflow interfering with the convergence criterion.
931 : ! The factor is a power of the BASE.
932 :
933 : ! Find largest and smallest moduli of coefficients.
934 5 : min = INF
935 5 : max = 0._TKC
936 155 : do i = 1_IK, nn
937 141 : if (sh(i)%re > max) max = sh(i)%re
938 155 : if (sh(i)%re /= 0._TKC .and. sh(i)%re < min) min = sh(i)%re
939 : end do
940 :
941 : ! Scale only if there are very large or very small components.
942 14 : if (LOW <= min .and. max <= HIGH) then
943 5 : bnd = 1._TKC
944 : else
945 0 : bnd = LOW / min
946 0 : if (bnd <= 1._TKC) then
947 0 : bnd = BASE ** int((log(max) + log(min)) * NEG_HALF_INV_LOG_BASE + .5_TKC)
948 0 : elseif (max < INF / bnd) then
949 0 : bnd = 1._TKC
950 : else
951 0 : bnd = BASE ** int(log(bnd) * INV_LOG_BASE + .5_TKC)
952 : end if
953 0 : if (bnd /= 1._TKC) then
954 : ! Scale the polynomial.
955 0 : do i = 1_IK, nn
956 0 : workspace(i) = bnd * workspace(i)
957 : end do
958 : end if
959 : end if
960 :
961 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
962 :
963 : ! Start the algorithm for one zero.
964 : loopIter: do ! 40
965 :
966 127 : if (nn <= 2) then
967 : ! Calculate the final zero and return.
968 14 : count = degree
969 14 : root(degree) = GET_RATIO(-workspace(2), workspace(1))
970 14 : return
971 : end if
972 :
973 : ! Calculate bnd, a lower bound on the modulus of the zeros.
974 913 : do i = 1, nn
975 913 : sh(i)%re = GET_ABS(workspace(i)) ! hypot
976 : end do
977 :
978 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
979 : ! call cauchy(nn, sh, bnd)
980 : ! Cauchy computes a lower bound on the moduli of the zeros of a polynomial -
981 : ! the real component of `sh` is the modulus of the coefficients.
982 : ! Compute upper estimate of bound.
983 113 : n = nn - 1_IK
984 113 : sh(nn)%re = -sh(nn)%re
985 113 : x = exp((log(-sh(nn)%re) - log(sh(1)%re)) / n)
986 113 : if (sh(n)%re /= 0._TKC) then
987 : ! if newton step at the origin is better, use it.
988 113 : xm = -sh(nn)%re / sh(n)%re
989 113 : if (xm < x) x = xm
990 : end if
991 : ! Chop the interval (0,x) until `f <= 0`.
992 0 : do ! 10
993 73 : xm = x * .1_TKC
994 40 : f = sh(1)%re
995 800 : do i = 2_IK, nn
996 800 : f = f * xm + sh(i)%re
997 : end do
998 113 : if (f > 0._TKC) then
999 0 : x = xm
1000 : cycle ! 10
1001 : end if
1002 : exit
1003 : end do
1004 40 : dx = x
1005 : ! Do newton iteration until x converges to two decimal places.
1006 329 : do ! 30
1007 442 : if (abs(dx / x) > .005_TKC) then
1008 329 : sh(1)%im = sh(1)%re
1009 2366 : do i = 2_IK, nn
1010 2366 : sh(i)%im = sh(i - 1)%im * x + sh(i)%re
1011 : end do
1012 329 : f = sh(nn)%im
1013 329 : df = sh(1)%im
1014 2037 : do i = 2_IK, n
1015 2037 : df = df * x + sh(i)%im
1016 : end do
1017 329 : dx = f / df
1018 329 : x = x - dx
1019 : cycle
1020 : end if
1021 : exit
1022 : end do
1023 40 : bnd = x
1024 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1025 :
1026 : ! The outer loop to control two major passes with different sequences of shifts.
1027 113 : do cnt1 = 1_IK, 2_IK
1028 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1029 : ! First stage calculation, no shift.
1030 : ! Compute the derivative polynomial as the initial `h`
1031 : ! polynomial and compute `l1` no-shift `h` polynomials.
1032 : ! call noshft(5_IK)
1033 33 : n = nn - 1_IK
1034 113 : nm1 = n - 1_IK
1035 800 : do i = 1_IK, n
1036 687 : xni = nn - i
1037 800 : h(i) = xni * workspace(i) / real(n, TKC)
1038 : end do
1039 678 : do jj = 1_IK, 5_IK
1040 678 : if (GET_ABS(h(n)) > 10._TKC * EPS * GET_ABS(workspace(n))) then ! hypot
1041 565 : t = GET_RATIO(-workspace(nn), h(n))
1042 3435 : do i = 1, nm1
1043 2870 : j = nn - i
1044 2870 : temp = h(j - 1)
1045 2870 : h(j)%re = t%re * temp%re - t%im * temp%im + workspace(j)%re
1046 3435 : h(j)%im = t%re * temp%im + t%im * temp%re + workspace(j)%im
1047 : end do
1048 565 : h(1) = workspace(1)
1049 : else
1050 : ! If the constant term is essentially zero, shift `h` coefficients.
1051 0 : do i = 1, nm1
1052 0 : j = nn - i
1053 0 : h(j) = h(j - 1)
1054 : end do
1055 0 : h(1) = ZERO
1056 : end if
1057 : end do
1058 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1059 : ! Inner loop to select a shift.
1060 115 : do cnt2 = 1_IK, 9_IK
1061 : ! Shift is chosen with modulus `bnd` and amplitude rotated by 94 degrees from the previous shift.
1062 115 : xxx = COSR * xx - SINR * yy
1063 115 : yy = SINR * xx + COSR * yy
1064 41 : xx = xxx
1065 115 : s%re = bnd * xx
1066 115 : s%im = bnd * yy
1067 : ! Second stage calculation, fixed shift.
1068 115 : count = degree - nn + 2_IK
1069 115 : call fxshft(10_IK * cnt2, root(count), converged)
1070 115 : if (converged) then
1071 : ! The second stage jumps directly to the third stage iteration.
1072 : ! If successful the zero is stored and the polynomial deflated.
1073 113 : nn = nn - 1_IK
1074 800 : do i = 1_IK, nn
1075 800 : workspace(i) = qp(i)
1076 : end do
1077 : cycle loopIter ! 40
1078 : end if
1079 : ! if the iteration is unsuccessful another shift is chosen.
1080 : end do
1081 : ! if 9 shifts fail, the outer loop is repeated with another sequence of shifts.
1082 : end do
1083 : ! The root finder has failed on two major passes. Return empty handed.
1084 0 : count = count - 1_IK
1085 0 : return
1086 : end do loopIter
1087 :
1088 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1089 :
1090 : contains
1091 :
1092 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1093 :
1094 230 : subroutine fxshft(l2, z, converged)
1095 : ! Computes l2 fixed-shift h polynomials and tests for convergence.
1096 : ! Initiates a variable-shift iteration and returns with the approximate zero if successful.
1097 : ! l2 - limit of fixed shift steps.
1098 : ! z - the output approximate zero if `converged == .true._LK`.
1099 : ! converged - logical indicating convergence of stage 3 iteration.
1100 : integer(IK) , intent(in) :: l2
1101 : complex(TKC), intent(out) :: z
1102 : logical(LK) , intent(out) :: converged
1103 : complex(TKC) :: ot, svs
1104 : logical(LK) :: test, pasd, bool
1105 : integer(LK) :: i, j, n
1106 115 : n = nn - 1_IK
1107 : ! Evaluate `p` at `s`.
1108 82 : call polyev(nn, s, workspace, qp, pv)
1109 : test = .true._LK
1110 : pasd = .false._LK
1111 : ! calculate first t = -p(s)/h(s).
1112 115 : call calct(bool)
1113 : ! main loop for one second stage step.
1114 342 : do j = 1_IK, l2
1115 333 : ot = t
1116 : ! Compute next `h` polynomial and new `t`.
1117 333 : call nexth(bool)
1118 333 : call calct(bool)
1119 333 : z = s + t
1120 : ! Test for convergence unless stage 3 has failed once or this is the last h polynomial.
1121 342 : if (.not. (bool .or. .not. test .or. j == l2)) then
1122 261 : if (GET_ABS(t - ot) < .5_TKC * GET_ABS(z)) then
1123 232 : if (pasd) then
1124 : ! The weak convergence test has been passed twice, start the third stage iteration,
1125 : ! after saving the current h polynomial and shift.
1126 806 : do i = 1, n
1127 806 : sh(i) = h(i)
1128 : end do
1129 115 : svs = s
1130 115 : call vrshft(10_IK, z, converged)
1131 115 : if (converged) return
1132 : ! The iteration failed to converge. turn off testing and restore `h, s, pv, t`.
1133 : test = .false._LK
1134 70 : do i = 1_IK, n
1135 70 : h(i) = sh(i)
1136 : end do
1137 9 : s = svs
1138 7 : call polyev(nn, s, workspace, qp, pv)
1139 9 : call calct(bool)
1140 9 : cycle
1141 : end if
1142 : pasd = .true._LK
1143 : else
1144 : pasd = .false._LK
1145 : end if
1146 : end if
1147 : end do
1148 : ! Attempt an iteration with final `h` polynomial from the second stage.
1149 9 : call vrshft(10_IK, z, converged)
1150 : end subroutine fxshft
1151 :
1152 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1153 :
1154 124 : subroutine vrshft(limss3, z, converged)
1155 : ! carries out the third stage iteration.
1156 : ! limss3: limit of steps in stage 3.
1157 : ! z: on entry contains the initial iterate, if the iteration converges it contains the final iterate on exit.
1158 : ! converged: `.true.` if iteration converges.
1159 : integer(IK) , intent(in) :: limss3
1160 : complex(TKC), intent(inout) :: z
1161 : logical(LK) , intent(out) :: converged
1162 : real(TKC) :: mp, ms, omp, relstp, r1, r2, tp
1163 : logical(LK) :: b, bool
1164 : integer(IK) :: i, j
1165 124 : converged = .false._LK
1166 : b = .false._LK
1167 124 : s = z
1168 : ! main loop for stage three.
1169 483 : do i = 1_IK, limss3
1170 : ! Evaluate p at s and test for convergence.
1171 397 : call polyev(nn, s, workspace, qp, pv)
1172 482 : mp = GET_ABS(pv)
1173 482 : ms = GET_ABS(s)
1174 482 : if (mp <= 20._TKC * getErrHorner(nn, qp, ms, mp)) then
1175 : ! Polynomial value is smaller in value than a bound on the error in evaluating p, terminate the iteration.
1176 113 : converged = .true._LK
1177 113 : z = s
1178 123 : return
1179 : end if
1180 369 : if (i /= 1_IK) then
1181 253 : if (.not. (b .or. mp < omp .or. relstp >= .05_TKC)) then
1182 : ! Iteration has stalled. probably a cluster of zeros.
1183 : ! Do 5 fixed shift steps into the cluster to force one zero to dominate.
1184 3 : tp = relstp
1185 : b = .true._LK
1186 6 : if (relstp < EPS) tp = EPS
1187 6 : r1 = sqrt(tp)
1188 6 : r2 = s%re * (1._TKC + r1) - s%im * r1
1189 6 : s%im = s%re * r1 + s%im * (1._TKC + r1)
1190 6 : s%re = r2
1191 5 : call polyev(nn, s, workspace, qp, pv)
1192 36 : do j = 1_IK, 5_IK
1193 30 : call calct(bool)
1194 36 : call nexth(bool)
1195 : end do
1196 3 : omp = INF
1197 3 : goto 20
1198 : end if
1199 : ! Exit if polynomial value increases significantly.
1200 247 : if (omp < .1_TKC * mp) return
1201 : end if
1202 179 : omp = mp
1203 : ! Calculate the next iterate.
1204 359 : 20 call calct(bool)
1205 359 : call nexth(bool)
1206 359 : call calct(bool)
1207 1201 : if (.not. bool) then
1208 359 : relstp = GET_ABS(t) / GET_ABS(s)
1209 359 : s = s + t
1210 : end if
1211 : end do
1212 : end subroutine vrshft
1213 :
1214 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1215 :
1216 1205 : subroutine calct(bool)
1217 : ! Computes t = -p(s) / h(s).
1218 : ! bool: logical(LK) , set true if h(s) is essentially zero.
1219 : logical(LK) , intent(out) :: bool
1220 : complex(TKC) :: hv
1221 : integer(IK) :: n
1222 1205 : n = nn - 1_IK
1223 : ! Evaluate h(s).
1224 1205 : call polyev(n, s, h, qh, hv)
1225 1205 : bool = GET_ABS(hv) <= 10._TKC * ARE * GET_ABS(h(n))
1226 1205 : if (.not. bool) then
1227 1205 : t = GET_RATIO(-pv, hv)
1228 1205 : return
1229 : end if
1230 0 : t = ZERO
1231 : end subroutine calct
1232 :
1233 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1234 :
1235 722 : subroutine nexth(bool)
1236 : ! Compute the next shifted `h` polynomial.
1237 : ! bool: logical(LK) , if .true._LK h(s) is essentially zero
1238 : logical(LK) , intent(in) :: bool
1239 : complex(TKC) :: temp
1240 : integer(IK) :: j,n
1241 722 : n = nn - 1_IK
1242 722 : if (.not. bool) then
1243 4210 : do j = 2_IK, n
1244 3488 : temp = qh(j - 1)
1245 3488 : h(j)%re = t%re * temp%re - t%im * temp%im + qp(j)%re
1246 4210 : h(j)%im = t%re * temp%im + t%im * temp%re + qp(j)%im
1247 : end do
1248 722 : h(1) = qp(1)
1249 722 : return
1250 : end if
1251 : ! If `h(s) == zero`, replace `h` with `qh`.
1252 0 : do j = 2_IK, n
1253 0 : h(j) = qh(j - 1)
1254 : end do
1255 0 : h(1) = ZERO
1256 : end subroutine nexth
1257 :
1258 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1259 :
1260 1817 : pure subroutine polyev(nn, s, workspace, q, pv)
1261 : ! Evaluate a polynomial p at s by the horner recurrence placing the partial sums in q and the computed value in pv.
1262 : integer(IK) , intent(in) :: nn
1263 : complex(TKC), intent(in) :: s
1264 : complex(TKC), intent(in) , contiguous :: workspace(:)
1265 : complex(TKC), intent(out) , contiguous :: q(:)
1266 : complex(TKC), intent(out) :: pv
1267 : real(TKC) :: tempo
1268 : integer(IK) :: i
1269 1817 : q(1) = workspace(1)
1270 1817 : pv = q(1)
1271 11296 : do i = 2_IK, nn
1272 9479 : tempo = pv%re * s%re - pv%im * s%im + workspace(i)%re
1273 9479 : pv%im = pv%re * s%im + pv%im * s%re + workspace(i)%im
1274 9479 : pv%re = tempo
1275 11296 : q(i) = pv
1276 : end do
1277 1817 : end subroutine polyev
1278 :
1279 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1280 :
1281 482 : pure function getErrHorner(nn, q, ms, mp) result(errHorner)
1282 : ! Bound the error in evaluating the polynomial by the horner recurrence.
1283 : ! q: the partial sums.
1284 : ! ms: modulus of the point.
1285 : ! mp: modulus of polynomial value.
1286 : integer(IK) , intent(in) :: nn
1287 : complex(TKC), intent(in) :: q(:)
1288 : real(TKC) , intent(in) :: ms
1289 : real(TKC) , intent(in) :: mp
1290 : real(TKC) :: errHorner
1291 : real(TKC) :: e
1292 : integer :: i
1293 482 : e = GET_ABS(q(1)) * MRE / (ARE + MRE)
1294 3805 : do i = 1_IK, nn
1295 3805 : e = e * ms + GET_ABS(q(i))
1296 : end do
1297 482 : errHorner = e * (ARE + MRE) - mp * MRE
1298 482 : end function
1299 :
1300 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1301 : #elif setPolyRoot_ENABLED && Jen_ENABLED && RK_CK_ENABLED
1302 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1303 :
1304 : ! Global variables.
1305 :
1306 : integer(IK) , parameter :: MAX_ITER = 30_IK
1307 : integer , parameter :: RKR = TKC ! SELECTED_REAL_KIND(int(precision(0._TKC) / 3.))
1308 : real(RKR) , parameter :: BASE = radix(0._RKR)
1309 : real(RKR) , parameter :: EPS = epsilon(0._TKC) ! important to ask for the `TKC` precision.
1310 : real(RKR) , parameter :: TIN = tiny(0._RKR)
1311 : real(RKR) , parameter :: INF = huge(0._RKR)
1312 : real(RKR) , parameter :: MRE = EPS ! error bound on complex multiplication.
1313 : real(RKR) , parameter :: ARE = EPS ! The error bound on `+` operation.
1314 : real(TKC) :: u, v, a, b, c, d, a1, a3, a7, e, f, g, h ! , a2, a6
1315 : real(TKC) , allocatable :: qp(:), k(:), qk(:), svk(:), workspace(:)
1316 : complex(TKC) :: s, sz, lz
1317 : integer(IK) :: n, nn
1318 :
1319 : ! Local variables.
1320 :
1321 : real(RKR) :: max, min, xx, yy, xxx, x, sc, bnd, xm, ff, df, dx
1322 : real(RKR) , parameter :: SINR = SIND(94._RKR), COSR = COSD(94._RKR) ! rotation by 94 degrees.
1323 : real(RKR) , parameter :: LOW = TIN / EPS
1324 : real(RKR) , allocatable :: pt(:)
1325 : real(TKC) , allocatable :: temp(:)
1326 : real(TKC) :: t, aa, bb, cc, factor
1327 : integer(IK) :: cnt, nz, i, j, jj, nm1
1328 : logical(LK) :: zerok, scalingEnabled
1329 :
1330 9 : nn = size(coef, 1, IK)
1331 9 : count = nn - 1_IK ! degree of the polynomial.
1332 9 : n = count
1333 :
1334 9 : CHECK_ASSERTION(__LINE__, 1_IK < nn, SK_"@setPolyRoot(): The condition `1 < size(coef)` must hold. size(coef) = "//getStr(nn))
1335 9 : CHECK_ASSERTION(__LINE__, coef(count) /= 0._TKC, SK_"@setPolyRoot(): The condition `coef(size(coef)) /= 0.` must hold. coef = "//getStr(coef))
1336 27 : CHECK_ASSERTION(__LINE__, size(root, 1, IK) == count, SK_"@setPolyRoot(): The condition `size(root) == size(coef) - 1` must hold. size(root), size(coef) - 1 = "//getStr([size(root, 1, IK), count]))
1337 :
1338 : ! Initialization of constants for shift rotation.
1339 2 : xx = sqrt(.5_RKR)
1340 2 : yy = -xx
1341 :
1342 : ! Algorithm fails if the leading coefficient is zero.
1343 9 : if (coef(count) == 0._TKC) then
1344 0 : count = 0_IK
1345 0 : return
1346 : end if
1347 :
1348 : ! Extract any exact ZERO roots and set count = degree of remaining polynomial.
1349 9 : do j = 1_IK, count
1350 9 : if (coef(j - 1) /= 0._TKC) exit
1351 0 : root(j) = (0._TKC, 0._TKC)
1352 0 : nn = nn - 1_IK
1353 9 : n = n - 1_IK
1354 : !count = count - 1_IK
1355 : end do
1356 :
1357 : ! bypass gfortran bug for heap allocations.
1358 : #if __GFORTRAN__
1359 9 : if (allocated(k)) deallocate(k)
1360 9 : if (allocated(pt)) deallocate(pt)
1361 9 : if (allocated(qp)) deallocate(qp)
1362 9 : if (allocated(qk)) deallocate(qk)
1363 9 : if (allocated(svk)) deallocate(svk)
1364 9 : if (allocated(temp)) deallocate(temp)
1365 9 : if (allocated(workspace)) deallocate(workspace)
1366 : #endif
1367 9 : allocate(k(nn), pt(nn), qp(nn), qk(nn), svk(nn), temp(nn), workspace(nn))
1368 :
1369 : ! Make a copy of the coefficients.
1370 76 : do j = 1_IK, nn
1371 76 : workspace(j) = coef(nn - j)
1372 : end do
1373 :
1374 : !write(*,*) "nn, workspace", nn, workspace
1375 :
1376 : ! Start the algorithm for one zero.
1377 : loopIter: do ! 30
1378 :
1379 44 : if (n <= 2_IK) then
1380 9 : if (n < 1_IK) return
1381 : ! Compute the final zero or pair of zeros.
1382 9 : if (n /= 2_IK) then
1383 6 : root(count)%re = -workspace(2) / workspace(1)
1384 6 : root(count)%im = 0._TKC
1385 6 : return
1386 : end if
1387 3 : call quad(workspace(1), workspace(2), workspace(3), root(count - 1), root(count))
1388 3 : return
1389 : end if
1390 :
1391 : ! Find largest and smallest moduli of coefficients.
1392 10 : max = 0._RKR
1393 10 : min = INF
1394 269 : do i = 1_IK, nn
1395 234 : x = abs(real(workspace(i), RKR))
1396 234 : if (x > max) max = x
1397 269 : if (x /= 0._RKR .and. x < min) min = x
1398 : end do
1399 :
1400 : ! Scale if there are large or very small coefficients.
1401 : ! Computes a scale factor to multiply the coefficients of the polynomial.
1402 : ! The scaling is done to avoid overflow and to avoid undetected underflow interfering with the convergence criterion.
1403 : ! The factor is a power of the `BASE`.
1404 35 : sc = LOW / min
1405 35 : if (sc <= 1._RKR) then
1406 35 : scalingEnabled = logical(max >= 10._RKR, LK)
1407 : else
1408 0 : scalingEnabled = logical(INF / sc >= max, LK)
1409 : end if
1410 35 : if (scalingEnabled) then
1411 8 : if (sc == 0._RKR) sc = TIN
1412 8 : factor = real(BASE, TKC) ** int(log(sc) / log(BASE) + .5_RKR, IK)
1413 8 : if (factor /= 1._TKC) then
1414 : do concurrent(i = 1_IK : nn)
1415 71 : workspace(i) = factor * workspace(i)
1416 : end do
1417 : end if
1418 : end if
1419 :
1420 : ! Compute lower bound on moduli of zeros.
1421 : do concurrent(i = 1_IK : nn)
1422 269 : pt(i) = real(abs(workspace(i)), RKR)
1423 : end do
1424 35 : pt(nn) = -pt(nn)
1425 :
1426 : ! Compute upper estimate of bound.
1427 35 : x = exp((log(-pt(nn)) - log(pt(1))) / n)
1428 35 : if (pt(n) /= 0._RKR) then
1429 : ! If newton step at the origin is better, use it.
1430 35 : xm = -pt(nn) / pt(n)
1431 35 : if (xm < x) x = xm
1432 : end if
1433 :
1434 : ! Chop the interval `(0, x)` until `ff <= 0`.
1435 0 : do
1436 25 : xm = x * .1_RKR
1437 10 : ff = pt(1)
1438 234 : do i = 2_IK, nn
1439 234 : ff = ff * xm + pt(i)
1440 : end do
1441 35 : if (ff <= 0._RKR) exit
1442 25 : x = xm
1443 : end do
1444 10 : dx = x
1445 :
1446 : ! Do newton iteration until `x` converges to two decimal places.
1447 118 : do
1448 153 : if (abs(dx / x) <= .005_RKR) exit
1449 33 : ff = pt(1)
1450 33 : df = ff
1451 652 : do i = 2_IK, n
1452 534 : ff = ff * x + pt(i)
1453 652 : df = df * x + ff
1454 : end do
1455 118 : ff = ff * x + pt(nn)
1456 118 : dx = ff / df
1457 118 : x = x - dx
1458 : end do
1459 10 : bnd = x
1460 :
1461 : ! Compute the derivative as the initial `k` polynomial and do `5` steps with no shift.
1462 35 : nm1 = n - 1_IK
1463 199 : do i = 2_IK, n
1464 199 : k(i) = (nn - i) * workspace(i) / n
1465 : end do
1466 35 : k(1) = workspace(1)
1467 35 : aa = workspace(nn)
1468 35 : bb = workspace(n)
1469 35 : zerok = logical(k(n) == 0._TKC, LK)
1470 210 : do jj = 1_IK, 5_IK
1471 175 : cc = k(n)
1472 210 : if (.not. zerok) then
1473 : ! Use scaled form of recurrence if value of k at 0 is nonzero
1474 175 : t = -aa / cc
1475 995 : do i = 1_IK, nm1
1476 820 : j = nn - i
1477 995 : k(j) = t * k(j - 1) + workspace(j)
1478 : end do
1479 175 : k(1) = workspace(1)
1480 175 : zerok = abs(k(n)) <= abs(bb) * EPS * 10._RKR
1481 : else
1482 : ! Use unscaled form of recurrence.
1483 0 : do i = 1_IK, nm1
1484 0 : j = nn - i
1485 0 : k(j) = k(j - 1)
1486 : end do
1487 0 : k(1) = 0._TKC
1488 0 : zerok = logical(k(n) == 0._TKC, LK)
1489 : end if
1490 : end do
1491 :
1492 : ! Save k for restarts with new shifts.
1493 234 : temp(1 : n) = k(1 : n)
1494 :
1495 : ! Loop to select the quadratic corresponding to each new shift.
1496 35 : do cnt = 1_IK, MAX_ITER
1497 : ! Quadratic corresponds to a double shift to a non-real point and niter complex conjugate.
1498 : ! The point has modulus bnd and amplitude rotated by 94 degrees from the previous shift.
1499 35 : xxx = COSR * xx - SINR * yy
1500 35 : yy = SINR * xx + COSR * yy
1501 10 : xx = xxx
1502 35 : s%re = bnd * xx
1503 35 : s%im = bnd * yy
1504 35 : u = -2._TKC * s%re
1505 35 : v = bnd
1506 : ! Second stage calculation, fixed quadratic.
1507 35 : call fxshfr(MAX_ITER * cnt, nz)
1508 35 : if (nz /= 0_IK) then
1509 : ! The second stage jumps directly to one of the third stage iterations and returns here if successful.
1510 : ! Deflate the polynomial, store the zero or zeros and return to the main algorithm.
1511 35 : j = count - n + 1_IK
1512 35 : root(j) = sz
1513 : !write(*,*) "j", j, count, n, nn, nz
1514 35 : nn = nn - nz
1515 35 : n = nn - 1_IK
1516 223 : workspace(1 : nn) = qp(1 : nn)
1517 35 : if (nz /= 1_IK) root(j + 1) = lz
1518 : cycle loopIter
1519 : end if
1520 : ! If the iteration is unsuccessful another quadratic is chosen after restoring `k`.
1521 0 : k(1 : nn) = temp(1 : nn)
1522 : end do
1523 : ! Return with failure if no convergence with `20` shifts.
1524 0 : count = count - n
1525 0 : return
1526 : end do loopIter
1527 :
1528 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1529 :
1530 : contains
1531 :
1532 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1533 :
1534 70 : subroutine fxshfr(l2, nz)
1535 : ! Compute up to `l2` fixed shift `k`-polynomials, testing for convergence in the linear or quadratic case.
1536 : ! Initiate one of the variable shift iterations and returns with the number of zeros found.
1537 : ! l2: limit of fixed shift steps
1538 : ! nz: number of zeros found
1539 : integer(IK) , intent(in) :: l2
1540 : integer(IK) , intent(out) :: nz
1541 : real(TKC) :: svu, svv, ui, vi, slocal
1542 : real(RKR) :: betas, betav, oss, ovv, ss, vv, ts, tv, ots, otv, tvv, tss
1543 : logical(LK) :: vpass, spass, vtry, stry
1544 : integer(IK) :: itype, j, iflag
1545 35 : nz = 0_IK
1546 10 : betav = .25_RKR
1547 10 : betas = .25_RKR
1548 35 : oss = real(s%re, RKR)
1549 35 : ovv = real(v, RKR)
1550 : ! Evaluate polynomial by synthetic division.
1551 : !write(*,*) "quadsd1, nn, u, v, workspace, qp, a, b", nn, u, v, workspace, qp, a, b
1552 35 : call quadsd(nn, u, v, workspace, qp, a, b)
1553 : !write(*,*) "quadsd2, nn, u, v, workspace, qp, a, b", nn, u, v, workspace, qp, a, b
1554 : !write(*,*) "calcsc1, itype", itype
1555 35 : call calcsc(itype)
1556 : !write(*,*) "calcsc2, itype", itype
1557 84 : do j = 1_IK, l2
1558 : ! Compute the next `k` polynomial and estimate `v`.
1559 84 : call nextk(itype)
1560 84 : call calcsc(itype)
1561 84 : call newest(itype, ui, vi)
1562 84 : vv = real(vi, RKR)
1563 : ! Estimate `s`.
1564 23 : ss = 0._RKR
1565 84 : if (k(n) /= 0._TKC) ss = real(-workspace(nn) / k(n), RKR)
1566 23 : tv = 1._RKR
1567 23 : ts = 1._RKR
1568 : !write(*,*) "BLOCK: ", j, itype
1569 84 : if (j /= 1_IK .and. itype /= 3_IK) then
1570 : ! Compute relative measures of convergence of `s` and `v` sequences.
1571 49 : if (vv /= 0._RKR) tv = abs((vv - ovv) / vv)
1572 49 : if (ss /= 0._RKR) ts = abs((ss - oss) / ss)
1573 : ! If decreasing, multiply two most recent convergence measures.
1574 13 : tvv = 1._RKR
1575 49 : if (tv < otv) tvv = tv * otv
1576 13 : tss = 1._RKR
1577 49 : if (ts < ots) tss = ts * ots
1578 : ! Compare with convergence criteria.
1579 49 : vpass = tvv < betav
1580 49 : spass = tss < betas
1581 49 : if (spass .or. vpass) then
1582 : ! At least one sequence has passed the convergence test. Store variables before iterating.
1583 37 : svu = u
1584 37 : svv = v
1585 242 : svk(1 : n) = k(1 : n)
1586 37 : slocal = ss
1587 : ! Choose iteration according to the fastest converging sequence.
1588 : vtry = .false._LK
1589 : stry = .false._LK
1590 37 : if (spass .and. ((.not. vpass) .or. tss < tvv)) goto 40
1591 : !write(*,*) "quadit1", ui, vi, nz
1592 19 : 20 call quadit(ui, vi, nz)
1593 : !write(*,*) "quadit2", ui, vi, nz
1594 43 : if (nz > 0_IK) return
1595 : ! Quadratic iteration has failed. flag that it has been tried and decrease the convergence criterion.
1596 : vtry = .true._LK
1597 8 : betav = betav * .25_RKR
1598 : ! Try linear iteration if it has not been tried and the `s` sequence is converging.
1599 8 : if (stry .or. (.not. spass)) goto 50
1600 48 : k(1 : n) = svk(1 : n)
1601 25 : 40 call realit(slocal, nz, iflag)
1602 25 : if (nz > 0_IK) return
1603 : ! Linear iteration has failed. flag that it has been tried and decrease the convergence criterion
1604 : stry = .true.
1605 1 : betas = betas * .25_RKR
1606 1 : if (iflag /= 0_IK) then
1607 : ! if linear iteration signals an almost double real zero attempt quadratic interation
1608 0 : ui = -(slocal + slocal)
1609 0 : vi = slocal * slocal
1610 0 : goto 20
1611 : end if
1612 : ! restore variables
1613 3 : 50 u = svu
1614 3 : v = svv
1615 12 : k(1 : n) = svk(1 : n)
1616 : ! try quadratic iteration if it has not been tried and the v sequence is converging
1617 3 : if (vpass .and. .not. vtry) goto 20
1618 : ! recompute qp and scalar values to continue the second stage
1619 2 : call quadsd(nn, u, v, workspace, qp, a, b)
1620 2 : call calcsc(itype)
1621 : end if
1622 : end if
1623 13 : ovv = vv
1624 13 : oss = ss
1625 13 : otv = tv
1626 97 : ots = ts
1627 : end do
1628 : end subroutine fxshfr
1629 :
1630 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1631 :
1632 19 : subroutine quadit(uu, vv, nz)
1633 : ! Variable-shift `k`-polynomial iteration for a quadratic factor, converges only if the zeros are equimodular or nearly so.
1634 : ! uu, vv: coefficients of starting quadratic
1635 : ! nz: number of zero found
1636 : real(TKC) , intent(in) :: uu
1637 : real(TKC) , intent(in) :: vv
1638 : integer(IK) , intent(out) :: nz
1639 : real(TKC) :: ui, vi
1640 : real(RKR) :: mp, omp, ee, relstp, t, zm
1641 : integer(IK) :: itype, i, j
1642 : logical(LK) :: tried
1643 19 : nz = 0_IK
1644 : tried = .false._LK
1645 19 : u = uu
1646 19 : v = vv
1647 : j = 0_IK
1648 : ! main loop
1649 43 : 10 call quad(1._TKC, u, v, sz, lz)
1650 : ! Return if roots of the quadratic are real and not close to multiple or nearly equal and of opposite sign.
1651 43 : if (abs(abs(sz%re) - abs(lz%re)) > .01_TKC * abs(lz%re)) return
1652 : ! Evaluate polynomial by quadratic synthetic division.
1653 35 : call quadsd(nn, u, v, workspace, qp, a, b)
1654 35 : mp = real(abs(a - sz%re * b) + abs(sz%im * b), RKR)
1655 : ! Compute a rigorous bound on the rounding error in evaluating P (`workspace`).
1656 35 : zm = sqrt(abs(real(v, RKR)))
1657 35 : ee = 2._RKR * abs(real(qp(1), RKR))
1658 10 : t = real(-sz%re * b, RKR)
1659 139 : do i = 2_IK, n
1660 139 : ee = ee * zm + abs(real(qp(i), RKR))
1661 : end do
1662 35 : ee = ee * zm + abs(real(a, RKR) + t)
1663 35 : ee = (5._RKR * MRE + 4._RKR * ARE) * ee - (5._RKR * MRE + 2._RKR * ARE) * (abs(real(a, RKR) + t) + abs(real(b, RKR)) * zm) + 2._RKR * ARE * abs(t)
1664 : ! Iteration has converged sufficiently if the polynomial value is less than 20 times this bound.
1665 35 : if (mp <= MAX_ITER * ee) then
1666 11 : nz = 2_IK
1667 11 : return
1668 : end if
1669 24 : j = j + 1_IK
1670 : ! Stop iteration after 20 steps.
1671 24 : if (j > MAX_ITER) return
1672 24 : if (j >= 2_IK) then
1673 14 : if (.not. (relstp > .01_RKR .or. mp < omp .or. tried)) then
1674 : ! A cluster appears to be stalling the convergence. Five fixed shift steps are taken with a `u, v` close to the cluster.
1675 0 : if (relstp < EPS) relstp = EPS
1676 0 : relstp = sqrt(relstp)
1677 0 : u = u - u * relstp
1678 0 : v = v + v * relstp
1679 0 : call quadsd(nn, u, v, workspace, qp, a, b)
1680 0 : do i = 1_IK, 5_IK
1681 0 : call calcsc(itype)
1682 0 : call nextk(itype)
1683 : end do
1684 : tried = .true._LK
1685 : j = 0_IK
1686 : end if
1687 : end if
1688 8 : omp = mp
1689 : ! Calculate next k polynomial and new u and v.
1690 24 : call calcsc(itype)
1691 24 : call nextk(itype)
1692 24 : call calcsc(itype)
1693 24 : call newest(itype, ui, vi)
1694 : ! If vi is zero the iteration is not converging.
1695 24 : if (vi == 0._TKC) return
1696 24 : relstp = real(abs((vi - v) / vi), RKR)
1697 24 : u = ui
1698 24 : v = vi
1699 24 : goto 10
1700 48 : end subroutine quadit
1701 :
1702 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1703 :
1704 25 : subroutine realit(sss,nz,iflag)
1705 : ! Variable-shift `h` polynomial iteration for a real zero.
1706 : ! sss : starting iterate.
1707 : ! nz : number of zero found.
1708 : ! iflag: flag to indicate a pair of zeros near real axis.
1709 : real(TKC) , intent(inout) :: sss
1710 : integer(IK) , intent(out) :: nz, iflag
1711 : real(TKC) :: pv, kv, t, s
1712 : real(RKR) :: ms, mp, omp, ee
1713 : integer(IK) :: i, j
1714 25 : nz = 0_IK
1715 25 : s = sss
1716 25 : iflag = 0_IK
1717 : j = 0_IK
1718 : ! main loop
1719 76 : do ! 10
1720 101 : pv = workspace(1)
1721 : ! Evaluate P (`workspace`) at s.
1722 101 : qp(1) = pv
1723 688 : do i = 2_IK, nn
1724 587 : pv = pv * s + workspace(i)
1725 688 : qp(i) = pv
1726 : end do
1727 101 : mp = real(abs(pv), RKR)
1728 : !write(*,*) "workspace", workspace
1729 : !write(*,*) "mp, pv, abs(pv)", mp, pv, abs(pv)
1730 : ! Compute a rigorous bound on the error in evaluating P (`workspace`).
1731 101 : ms = abs(real(s, RKR))
1732 101 : ee = (MRE / (ARE + MRE)) * abs(real(qp(1), RKR))
1733 688 : do i = 2_IK, nn
1734 688 : ee = ee * ms + abs(real(qp(i), RKR))
1735 : end do
1736 : ! Iteration has converged sufficiently if the polynomial value is less than 20 times this bound.
1737 101 : if (mp <= MAX_ITER * ((ARE + MRE) * ee - MRE * mp)) then
1738 24 : sz%re = s
1739 24 : sz%im = 0._TKC
1740 24 : nz = 1_IK
1741 24 : return
1742 : end if
1743 77 : j = j + 1_IK
1744 : ! Stop iteration after `10` steps.
1745 77 : if (j > 10_IK) return
1746 76 : if (j >= 2_IK) then
1747 52 : if (abs(t) <= .001_TKC * abs(s - t) .and. mp > omp) then
1748 : ! A cluster of zeros near the real axis has been encountered, return with iflag set to initiate a quadratic iteration.
1749 0 : iflag = 1_IK
1750 0 : sss = s
1751 0 : return
1752 : end if
1753 : end if
1754 : ! Return if the polynomial value has increased significantly.
1755 33 : omp = mp
1756 : ! Compute t, the next polynomial, and the new iterate.
1757 76 : kv = k(1)
1758 76 : qk(1) = kv
1759 428 : do i = 2_IK, n
1760 352 : kv = kv * s + k(i)
1761 428 : qk(i) = kv
1762 : end do
1763 76 : if (abs(kv) > abs(k(n)) * EPS * 10._TKC) then
1764 : ! Use the scaled form of the recurrence if the value of k at s is nonzero.
1765 76 : t = -pv / kv
1766 76 : k(1) = qp(1)
1767 428 : do i = 2_IK, n
1768 428 : k(i) = t * qk(i - 1) + qp(i)
1769 : end do
1770 : else
1771 : ! Use unscaled form.
1772 0 : k(1) = 0._TKC
1773 0 : do i = 2_IK, n
1774 0 : k(i) = qk(i - 1)
1775 : end do
1776 : end if
1777 76 : kv = k(1)
1778 428 : do i = 2_IK, n
1779 428 : kv = kv * s + k(i)
1780 : end do
1781 33 : t = 0._TKC
1782 76 : if (abs(kv) > abs(k(n)) * EPS * 10._TKC) t = -pv / kv
1783 76 : s = s + t
1784 : end do ! goto 10
1785 : end subroutine realit
1786 :
1787 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1788 :
1789 169 : subroutine calcsc(itype) ! checked
1790 : ! Compute scalar quantities need for the next `k` polynomial and new estimates of the quadratic coefficients.
1791 : ! itype: integer variable that indicates how the calculations are normalized to avoid overflow.
1792 : integer(IK), intent(out) :: itype
1793 : ! Synthetic division of k by the quadratic 1, u, v.
1794 169 : call quadsd(n, u, v, k, qk, c, d)
1795 169 : if (abs(c) <= abs(k(n)) * EPS * 100._TKC) then
1796 0 : if (abs(d) <= abs(k(n - 1)) * EPS * 100._TKC) then
1797 0 : itype = 3_IK ! Indicates the quadratic is almost a factor of `k`.
1798 0 : return
1799 : end if
1800 : end if
1801 169 : if (abs(d) >= abs(c)) then
1802 80 : itype = 2_IK ! Indicates that all formulas are divided by `d`.
1803 80 : e = a / d
1804 80 : f = c / d
1805 80 : g = u * b
1806 80 : h = v * b
1807 80 : a3 = (a + g) * e + h * (b / d)
1808 80 : a1 = b * f - a
1809 80 : a7 = (f + u) * a + h
1810 80 : return
1811 : end if
1812 89 : itype = 1_IK ! Indicates that all formulas are divided by `c`.
1813 89 : e = a / c
1814 89 : f = d / c
1815 89 : g = u * e
1816 89 : h = v * b
1817 89 : a3 = a * e + (h / c + g) * b
1818 89 : a1 = b - a * (d / c)
1819 89 : a7 = a + g * d + h * f
1820 89 : return
1821 : end subroutine calcsc
1822 :
1823 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1824 :
1825 108 : subroutine nextk(itype) ! checked
1826 : ! Compute the next `k` polynomials using scalars computed in `calcsc()`.
1827 : integer(IK) , intent(in) :: itype
1828 : real(TKC) :: temp
1829 : integer(IK) :: i
1830 108 : if (itype /= 3_IK) then
1831 108 : temp = a
1832 108 : if (itype == 1_IK) temp = b
1833 108 : if (abs(a1) <= abs(temp) * EPS * 10._TKC) then
1834 : ! If `a1` is nearly zero then use a special form of the recurrence.
1835 0 : k(1) = 0._TKC
1836 0 : k(2) = -a7 * qp(1)
1837 0 : do i = 3_IK, n
1838 0 : k(i) = a3 * qk(i - 2) - a7 * qp(i - 1)
1839 : end do
1840 : return
1841 : end if
1842 : ! Use scaled form of the recurrence.
1843 108 : a7 = a7 / a1
1844 108 : a3 = a3 / a1
1845 108 : k(1) = qp(1)
1846 108 : k(2) = qp(2) - a7 * qp(1)
1847 452 : do i = 3_IK, n
1848 452 : k(i) = a3 * qk(i - 2) - a7 * qp(i - 1) + qp(i)
1849 : end do
1850 : return
1851 : end if
1852 : ! Use unscaled form of the recurrence if itype is 3.
1853 0 : k(1) = 0._TKC
1854 0 : k(2) = 0._TKC
1855 0 : do i = 3_IK, n
1856 0 : k(i) = qk(i - 2)
1857 : end do
1858 : end subroutine nextk
1859 :
1860 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1861 :
1862 108 : subroutine newest(itype, uu, vv)
1863 : ! Compute new estimates of the quadratic coefficients using the scalars computed in calcsc.
1864 : integer(IK) , intent(in) :: itype
1865 : real(TKC) , intent(out) :: uu
1866 : real(TKC) , intent(out) :: vv
1867 : real(TKC) :: a4, a5, b1, b2, c1, c2, c3, c4, temp
1868 : ! Use formulas appropriate to setting of itype.
1869 108 : if (itype /= 3_IK) then
1870 108 : if (itype /= 2_IK) then
1871 54 : a4 = a + u * b + h * f
1872 54 : a5 = c + (u + v * f) * d
1873 : else
1874 54 : a4 = (a + g) * f + h
1875 54 : a5 = (f + u) * c + v * d
1876 : end if
1877 : ! Evaluate new quadratic coefficients.
1878 108 : b1 = -k(n) / workspace(nn)
1879 108 : b2 = -(k(n - 1) + b1 * workspace(n)) / workspace(nn)
1880 108 : c1 = v * b2 * a1
1881 108 : c2 = b1 * a7
1882 108 : c3 = b1 * b1 * a3
1883 108 : c4 = c1 - c2 - c3
1884 108 : temp = a5 + b1 * a4 - c4
1885 108 : if (temp /= 0._TKC) then
1886 108 : uu = u - (u * (c3 + c2) + v * (b1 * a1 + b2 * a7)) / temp
1887 108 : vv = v * (1._TKC + c4 / temp)
1888 108 : return
1889 : end if
1890 : end if
1891 : ! If itype = 3 the quadratic is zeroed.
1892 0 : uu = 0._TKC
1893 0 : vv = 0._TKC
1894 : end subroutine newest
1895 :
1896 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1897 :
1898 241 : pure subroutine quadsd(nn, u, v, workspace, q, a, b)
1899 : ! Divide P (`workspace`) by the quadratic `1, u, v` placing the quotient in `q` and the remainder in `a, b`.
1900 : integer(IK) , intent(in) :: nn
1901 : real(TKC) , intent(in) :: u, v, workspace(nn)
1902 : real(TKC) , intent(out) :: q(nn), a, b
1903 : real(TKC) :: c
1904 : integer(IK) :: i
1905 241 : b = workspace(1)
1906 241 : q(1) = b
1907 241 : a = workspace(2) - u * b
1908 241 : q(2) = a
1909 1036 : do i = 3_IK, nn
1910 795 : c = workspace(i) - u * a - v * b
1911 795 : q(i) = c
1912 795 : b = a
1913 1036 : a = c
1914 : end do
1915 : end subroutine quadsd
1916 :
1917 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1918 :
1919 46 : pure subroutine quad(a, b1, c, sroot, lroot)
1920 : ! Compute the zeros of the quadratic `a * z**2 + b1 * z + c`.
1921 : ! The quadratic formula, modified to avoid overflow, is used to find the larger zero if the zeros are real and both zeros are complex.
1922 : ! The smaller real zero is found directly from the product of the zeros `c / a`.
1923 : real(TKC) , intent(in) :: a, b1, c
1924 : complex(TKC), intent(out) :: sroot, lroot
1925 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC)
1926 : real(TKC) :: b, d, e
1927 46 : if (a == 0._TKC) then
1928 0 : lroot = ZERO
1929 0 : sroot = ZERO
1930 0 : if (b1 /= 0._TKC) sroot%re = -c / b1
1931 0 : return
1932 : end if
1933 46 : if (c == 0._TKC) then
1934 0 : lroot%re = -b1 / a
1935 0 : lroot%im = 0._TKC
1936 0 : sroot = ZERO
1937 0 : return
1938 : end if
1939 : ! Compute discriminant avoiding overflow.
1940 46 : b = b1 * .5_TKC
1941 46 : if (abs(b) >= abs(c)) then
1942 0 : e = 1._TKC - (a / b) * (c / b)
1943 0 : d = sqrt(abs(e)) * abs(b)
1944 : else
1945 13 : e = a
1946 46 : if (c < 0._TKC) e = -a
1947 46 : e = b * (b / abs(c)) - e
1948 46 : d = sqrt(abs(e)) * sqrt(abs(c))
1949 : end if
1950 46 : if (e >= 0._TKC) then
1951 : ! Real zeros.
1952 11 : if (b >= 0._TKC) d = -d
1953 11 : lroot%re = (-b + d) / a
1954 11 : lroot%im = 0._TKC
1955 11 : sroot = ZERO
1956 11 : if (lroot%re /= 0._TKC) sroot%re = (c / lroot%re) / a
1957 11 : return
1958 : end if
1959 : ! Complex conjugate zeros.
1960 35 : sroot%re = -b / a
1961 35 : lroot%re = sroot%re
1962 35 : sroot%im = abs(d / a)
1963 35 : lroot%im = -sroot%im
1964 : end subroutine
1965 :
1966 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1967 : #elif setPolyRoot_ENABLED && Lag_ENABLED
1968 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1969 :
1970 : integer(IK) :: ideg, jdeg, niter, degree
1971 : real(TKC), parameter :: EPS = 2 * epsilon(0._TKC) ! the estimated fractional roundoff error.
1972 23 : complex(TKC) :: temp, b, c, workspace(size(coef, 1, IK))
1973 23 : degree = size(coef, 1, IK) - 1_IK
1974 23 : CHECK_ASSERTION(__LINE__, 0_IK < degree, SK_"@setPolyRoot(): The condition `1 < size(coef)` must hold. size(coef) = "//getStr(size(coef, 1, IK)))
1975 23 : CHECK_ASSERTION(__LINE__, coef(degree) /= ZERO, SK_"@setPolyRoot(): The condition `coef(size(coef)) /= 0.` must hold. coef = "//getStr(coef))
1976 69 : CHECK_ASSERTION(__LINE__, size(root, 1, IK) == degree, SK_"@setPolyRoot(): The condition `size(root) == size(coef) - 1` must hold. size(root), size(coef) - 1 = "//getStr([size(root, 1, IK), degree + 1_IK]))
1977 23 : count = 0_IK
1978 231 : workspace = coef
1979 23 : if (degree < 1_IK) return
1980 208 : do ideg = degree, 1_IK, -1_IK
1981 185 : temp = ZERO
1982 185 : call setPolyRootPolished(temp, niter, workspace)
1983 208 : if (0_IK < niter) then
1984 185 : if (abs(temp%im) <= 2._TKC * EPS**2 * abs(temp%re)) temp%im = 0._TKC
1985 185 : count = count + 1_IK
1986 185 : root(count) = temp
1987 : ! deflate.
1988 185 : b = workspace(ideg + 1)
1989 1132 : do jdeg = ideg, 1_IK, -1_IK
1990 947 : c = workspace(jdeg)
1991 947 : workspace(jdeg) = b
1992 1132 : b = temp * b + c
1993 : end do
1994 : else
1995 : exit
1996 : end if
1997 : end do
1998 : ! Polish the roots using the un-deflated coefficients.
1999 : ideg = 0_IK
2000 : do
2001 208 : ideg = ideg + 1_IK
2002 208 : if (count < ideg) exit
2003 185 : call setPolyRootPolished(root(ideg), niter, coef)
2004 185 : if (0_IK < niter) cycle
2005 0 : count = count - 1_IK
2006 208 : root(ideg : count) = root(ideg + 1 : count + 1)
2007 : end do
2008 : ! Sort roots by their real parts by straight insertion.
2009 : !do jdeg = 2_IK, degree
2010 : ! x = root(jdeg)
2011 : ! do ideg = jdeg - 1, 1, -1
2012 : ! if(root(ideg)%re < x%re) exit
2013 : ! root(ideg + 1) = root(ideg)
2014 : ! end do
2015 : ! root(ideg + 1) = x
2016 : !end do
2017 :
2018 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2019 : #elif setPolyRootPolished_ENABLED && Lag_ENABLED && RK_RK_ENABLED
2020 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2021 :
2022 : complex(TKC) :: croot
2023 0 : croot = root
2024 0 : call setPolyRootPolished(croot, niter, coef)
2025 0 : root = croot%re
2026 :
2027 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2028 : #elif setPolyRootPolished_ENABLED && Lag_ENABLED && (CK_CK_ENABLED || RK_CK_ENABLED)
2029 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2030 :
2031 : ! Given the complex polynomial coefficients `coef(1 : degree + 1)` of the polynomial in increasing order,
2032 : ! a complex root `x`, this routine improves `x` by the Laguerre method until it converges,
2033 : ! within the achievable roundoff limit, to a root of the given polynomial.<br>
2034 : ! The number of iterations taken is returned as `niter`.<br>
2035 : ! Break (rare) limit cycles with different fractional values once every `nstep` steps, for `nitermax` total allowed iterations.<br>
2036 : real(TKC), parameter :: EPS = 10_IK * epsilon(0._TKC) ! the estimated fractional roundoff error.
2037 : real(TKC), parameter :: fraction(8) = [.5_TKC, .25_TKC, .75_TKC, .13_TKC, .38_TKC, .62_TKC, .88_TKC, 1._TKC] ! Fractions used to break a limit cycle.
2038 : integer(IK) , parameter :: nstep = 10, nitermax = size(fraction, 1, IK) * nstep * 10000 ! amir: Extra 10000 was added to allow convergence for extremely awkward coefficients.
2039 : complex(TKC) :: droot, rootnew, der0, der1, der2, der12, der12sq, hterm, dterm, denom1, denom2
2040 : real(TKC) :: absroot, absdenom1, absdenom2, relerr
2041 : integer(IK) :: degree, id, ifrac, counter
2042 376 : degree = size(coef, 1, IK) - 1
2043 376 : CHECK_ASSERTION(__LINE__, 0_IK < degree, SK_"@setPolyRootPolished(): The condition `1 < size(coef)` must hold. size(coef) = "//getStr(size(coef, 1, IK)))
2044 : counter = 0_IK
2045 : ifrac = 0_IK
2046 376 : niter = 0_IK
2047 376 : if (degree < 1_IK) return
2048 : do
2049 654938 : niter = niter + 1_IK
2050 654938 : der2 = coef(degree)
2051 654938 : relerr = abs(der2)
2052 654256 : der0 = ZERO
2053 654256 : der1 = ZERO
2054 654938 : absroot = abs(root)
2055 8507162 : do id = degree - 1, 0, -1
2056 7852224 : der0 = root * der0 + der1
2057 7852224 : der1 = root * der1 + der2
2058 7852224 : der2 = root * der2 + coef(id)
2059 8507162 : relerr = abs(der2) + absroot * relerr
2060 : end do
2061 654938 : relerr = EPS * relerr
2062 654938 : if(abs(der2) <= relerr) return
2063 : ! Use the Laguerre method.
2064 654562 : der12 = der1 / der2
2065 654562 : der12sq = der12 * der12
2066 654562 : hterm = der12sq - 2 * der0 / der2
2067 654562 : dterm = sqrt((degree - 1) * (degree * hterm - der12sq))
2068 654562 : denom1 = der12 + dterm
2069 654562 : denom2 = der12 - dterm
2070 654562 : absdenom1 = abs(denom1)
2071 654562 : absdenom2 = abs(denom2)
2072 654562 : if(absdenom1 < absdenom2) then
2073 215062 : absdenom1 = absdenom2
2074 326986 : denom1 = denom2
2075 : end if
2076 654562 : if (absdenom1 /= 0._TKC) then
2077 654546 : droot = degree / denom1
2078 : else
2079 16 : droot = exp(cmplx(log(1._TKC + absroot), real(niter, TKC), TKC))
2080 : endif
2081 654562 : rootnew = root - droot
2082 654562 : if(root == rootnew) return
2083 654562 : counter = counter + 1
2084 654562 : if (counter /= nstep) then
2085 589247 : root = rootnew
2086 : else
2087 : counter = 0_IK
2088 65315 : ifrac = ifrac + 1_IK
2089 65315 : root = root - droot * fraction(ifrac)
2090 65315 : if (ifrac == size(fraction)) ifrac = 0
2091 : end if
2092 654562 : if (niter == nitermax) exit
2093 : end do
2094 0 : niter = -niter
2095 :
2096 : #else
2097 : !%%%%%%%%%%%%%%%%%%%%%%%%
2098 : #error "Unrecognized interface."
2099 : !%%%%%%%%%%%%%%%%%%%%%%%%
2100 : #endif
2101 : #undef TYPE_KIND
2102 : #undef GET_RATIO
2103 : #undef GET_ABS
2104 : #undef COSD
2105 : #undef SIND
|
