From 759491156d937f569fd983dd6644b7e25bcd5d18 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 18 Sep 2022 02:12:22 -0500 Subject: [PATCH 1/3] some interface cleanups --- src/polyroots_module.F90 | 198 ++++++++++++++++++++------------------- test/rpoly_test.f90 | 7 +- 2 files changed, 102 insertions(+), 103 deletions(-) diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index 2cce654..089beb0 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -20,44 +20,44 @@ module polyroots_module -use iso_fortran_env + use iso_fortran_env -implicit none + implicit none -private + private #ifdef REAL32 -integer, parameter, public :: polyroots_module_rk = real32 !! real kind used by this module [4 bytes] + integer, parameter, public :: polyroots_module_rk = real32 !! real kind used by this module [4 bytes] #elif REAL64 -integer, parameter, public :: polyroots_module_rk = real64 !! real kind used by this module [8 bytes] + integer, parameter, public :: polyroots_module_rk = real64 !! real kind used by this module [8 bytes] #elif REAL128 -integer, parameter, public :: polyroots_module_rk = real128 !! real kind used by this module [16 bytes] + integer, parameter, public :: polyroots_module_rk = real128 !! real kind used by this module [16 bytes] #else -integer, parameter, public :: polyroots_module_rk = real64 !! real kind used by this module [8 bytes] + integer, parameter, public :: polyroots_module_rk = real64 !! real kind used by this module [8 bytes] #endif -integer, parameter :: wp = polyroots_module_rk !! local copy of `polyroots_module_rk` with a shorter name + integer, parameter :: wp = polyroots_module_rk !! local copy of `polyroots_module_rk` with a shorter name -real(wp), parameter :: eps = epsilon(1.0_wp) !! machine epsilon -real(wp), parameter :: pi = acos(-1.0_wp) + real(wp), parameter :: eps = epsilon(1.0_wp) !! machine epsilon + real(wp), parameter :: pi = acos(-1.0_wp) -! general polynomial routines: -public :: polyroots -public :: rpoly -public :: cpzero -public :: rpzero -public :: rpqr79 -public :: cpqr79 -public :: qr_algeq_solver + ! general polynomial routines: + public :: polyroots + public :: rpoly + public :: cpzero + public :: rpzero + public :: rpqr79 + public :: cpqr79 + public :: qr_algeq_solver -! special polynomial routines: -public :: dcbcrt -public :: dqdcrt + ! special polynomial routines: + public :: dcbcrt + public :: dqdcrt -! utility routines: -public :: dcbrt -public :: cpevl -public :: newton_root_polish + ! utility routines: + public :: dcbrt + public :: cpevl + public :: newton_root_polish contains !***************************************************************************************** @@ -75,10 +75,10 @@ subroutine rpoly(op, degree, zeror, zeroi, istat) implicit none - real(wp), intent(in) :: op(:) !! vector of coefficients in order of decreasing powers integer, intent(in) :: degree !! degree of polynomial - real(wp), intent(out) :: zeror(:) !! output vector of real parts of the zeros - real(wp), intent(out) :: zeroi(:) !! output vector of imaginary parts of the zeros + real(wp), dimension(degree+1), intent(in) :: op !! vector of coefficients in order of decreasing powers + real(wp), dimension(degree), intent(out) :: zeror !! output vector of real parts of the zeros + real(wp), dimension(degree), intent(out) :: zeroi !! output vector of imaginary parts of the zeros integer, intent(out) :: istat !! status output: !! !! * 0 : success @@ -86,17 +86,12 @@ subroutine rpoly(op, degree, zeror, zeroi, istat) !! * -2 : no roots found !! * >0 : the number of zeros found - ! these were formerly in a common block: - real(wp), allocatable :: p(:), qp(:), k(:), qk(:), svk(:) + real(wp), dimension(:), allocatable :: p, qp, k, qk, svk, temp, pt real(wp) :: sr, si, u, v, a, b, c, d, a1, a3, & - a7, e, f, g, h, szr, szi, lzr, lzi - integer :: n, nn - !-------------------------------------------------------- - - real(wp), allocatable :: temp(:) - real(wp), allocatable :: pt(:) - real(wp) :: t, aa, bb, cc, factor, mx, mn, xx, yy, xxx, x, sc, bnd, xm, ff, df, dx - integer :: cnt, nz, i, j, jj, l, nm1 + a7, e, f, g, h, szr, szi, lzr, lzi, & + t, aa, bb, cc, factor, mx, mn, xx, yy, & + xxx, x, sc, bnd, xm, ff, df, dx + integer :: cnt, nz, i, j, jj, l, nm1, n, nn logical :: zerok real(wp), parameter :: deg2rad = pi/180.0_wp @@ -1253,7 +1248,7 @@ end subroutine dqdcrt ! * [`/opt/companion.tgz`](https://netlib.org/opt/companion.tgz) on Netlib ! from [Edelman & Murakami (1995)](https://www.ams.org/journals/mcom/1995-64-210/S0025-5718-1995-1262279-2/S0025-5718-1995-1262279-2.pdf), -subroutine qr_algeq_solver(n, c, zr, zi, detil, istatus) +subroutine qr_algeq_solver(n, c, zr, zi, istatus, detil) implicit none @@ -1261,13 +1256,13 @@ subroutine qr_algeq_solver(n, c, zr, zi, detil, istatus) real(wp), intent(in) :: c(n + 1) !! coefficients of the polynomial. in order of decreasing powers. real(wp), intent(out) :: zr(n) !! real part of output roots real(wp), intent(out) :: zi(n) !! imaginary part of output roots - real(wp), intent(out) :: detil !! accuracy hint. integer, intent(out) :: istatus !! return code: !! !! * -1 : degree <= 0 !! * -2 : leading coefficient `c(1)` is zero !! * 0 : success !! * otherwise, the return code from `hqr_eigen_hessenberg` + real(wp), intent(out), optional :: detil !! accuracy hint. real(wp), allocatable :: a(:, :) !! work matrix integer, allocatable :: cnt(:) !! work area for counting the qr-iterations @@ -1307,14 +1302,18 @@ subroutine qr_algeq_solver(n, c, zr, zi, detil, istatus) return end if - ! count the total qr iteration. - iter = 0 - do i = 1, n - if (cnt(i) > 0) iter = iter + cnt(i) - end do + if (present(detil)) then - ! calculate the accuracy hint. - detil = eps*n*iter*afnorm + ! count the total qr iteration. + iter = 0 + do i = 1, n + if (cnt(i) > 0) iter = iter + cnt(i) + end do + + ! calculate the accuracy hint. + detil = eps*n*iter*afnorm + + end if contains @@ -1678,13 +1677,13 @@ subroutine cpevl(n, m, a, z, c, b, kbd) !! !! if M > N+1 function and all N derivatives will be calculated. complex(wp), intent(in) :: a(*) !! vector containing the N+1 coefficients of polynomial. - !! A(I) = coefficient of Z**(N+1-I) + !! A(I) = coefficient of Z**(N+1-I) complex(wp), intent(in) :: z !! point at which the evaluation is to take place complex(wp), intent(out) :: c(*) !! Array of 2(M+1) words: C(I+1) contains the complex value of the I-th !! derivative at Z, I=0,...,M complex(wp), intent(out) :: b(*) !! Array of 2(M+1) words: B(I) contains the bounds on the real and imaginary parts - !! of C(I) if they were requested. only needed if bounds are to be calculated. - !! It is not used otherwise. + !! of C(I) if they were requested. only needed if bounds are to be calculated. + !! It is not used otherwise. logical, intent(in) :: kbd !! A logical variable, e.g. .TRUE. or .FALSE. which is !! to be set .TRUE. if bounds are to be computed. @@ -1731,16 +1730,15 @@ end subroutine cpevl ! * 891214 Prologue converted to Version 4.0 format. (BAB) ! * Jacob Williams, 9/13/2022 : modernized this routine -subroutine cpzero(in, a, r, t, iflg, s) +subroutine cpzero(in, a, r, iflg, s) implicit none - integer, intent(in) :: in !! degree of P(Z) - complex(wp), intent(in) :: a(*) !! complex vector containing coefficients of P(Z), - !! A(I) = coefficient of Z**(N+1-i) - complex(wp), intent(inout) :: r(*) !! N word complex vector. On input: containing initial estimates for zeros - !! if these are known. On output: Ith zero - complex(wp) :: t(*) !! `4(N+1)` word array used for temporary storage + integer, intent(in) :: in !! `N`, the degree of `P(Z)` + complex(wp), dimension(in+1), intent(in) :: a !! complex vector containing coefficients of `P(Z)`, + !! `A(I)` = coefficient of `Z**(N+1-i)` + complex(wp), dimension(in), intent(inout) :: r !! `N` word complex vector. On input: containing initial + !! estimates for zeros if these are known. On output: Ith zero integer, intent(inout) :: iflg !!### On Input: !! !! flag to indicate if initial estimates of zeros are input: @@ -1760,16 +1758,19 @@ subroutine cpzero(in, a, r, t, iflg, s) !! * If IFLG == 2 on return, the program failed to converge !! after 25*N iterations. Best current estimates of the !! zeros are in R(I). Error bounds are not calculated. - real(wp), intent(out) :: s(*) !! an `N` word array. S(I) = bound for R(I) + real(wp), intent(out) :: s(in) !! an `N` word array. `S(I)` = bound for `R(I)` integer :: i, imax, j, n, n1, nit, nmax, nr real(wp) :: u, v, x complex(wp) :: pn, temp complex(wp) :: ctmp(1), btmp(1) + complex(wp), dimension(:), allocatable :: t !! `4(N+1)` word array used for temporary storage if (in <= 0 .or. abs(a(1)) == 0.0_wp) then iflg = 1 else + ! work array: + allocate(t(4*(in+1))) ! check for easily obtained zeros n = in n1 = n + 1 @@ -1880,16 +1881,15 @@ end subroutine cpzero ! !@note This is just a wrapper to [[cpzero]] -subroutine rpzero(n, a, r, t, iflg, s) +subroutine rpzero(n, a, r, iflg, s) implicit none - integer, intent(in) :: n !! degree of P(X) - real(wp), intent(in) :: a(*) !! real vector containing coefficients of P(X), - !! A(I) = coefficient of X**(N+1-I) - complex(wp), intent(inout) :: r(*) !! N word complex vector. On Input: containing initial estimates for zeros - !! if these are known. On output: ith zero. - complex(wp) :: t(*) !! `6(N+1)` word array used for temporary storage + integer, intent(in) :: n !! degree of `P(X)` + real(wp), dimension(n+1), intent(in) :: a !! real vector containing coefficients of `P(X)`, + !! `A(I)` = coefficient of `X**(N+1-I)` + complex(wp), dimension(n), intent(inout) :: r !! N word complex vector. On Input: containing initial estimates for zeros + !! if these are known. On output: ith zero. integer, intent(inout) :: iflg !!### On Input: !! !! flag to indicate if initial estimates of zeros are input: @@ -1909,15 +1909,17 @@ subroutine rpzero(n, a, r, t, iflg, s) !! * If IFLG == 2 on return, the program failed to converge !! after 25*N iterations. Best current estimates of the !! zeros are in R(I). Error bounds are not calculated. - real(wp), intent(out) :: s(*) !! an `N` word array. bound for R(I). + real(wp), dimension(n), intent(out) :: s !! an `N` word array. bound for R(I). + + integer :: i + complex(wp), dimension(:), allocatable :: p !! complex coefficients - integer :: i, n1 + allocate(p(n+1)) - n1 = n + 1 - do i = 1, n1 - t(i) = cmplx(a(i), 0.0_wp, wp) + do i = 1, n + 1 + p(i) = cmplx(a(i), 0.0_wp, wp) end do - call cpzero(n, t, r, t(n + 2), iflg, s) + call cpzero(n, p, r, iflg, s) end subroutine rpzero !***************************************************************************************** @@ -1942,9 +1944,9 @@ subroutine rpqr79(ndeg, coeff, root, ierr) implicit none integer, intent(in) :: ndeg !! degree of polynomial - real(wp), intent(in) :: coeff(*) !! `NDEG+1` coefficients in descending order. i.e., - !! `P(Z) = COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)` - complex(wp), intent(out) :: root(*) !! `NDEG` vector of roots + real(wp), dimension(ndeg+1), intent(in) :: coeff !! `NDEG+1` coefficients in descending order. i.e., + !! `P(Z) = COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)` + complex(wp), dimension(ndeg), intent(out) :: root !! `NDEG` vector of roots integer, intent(out) :: ierr !! Output Error Code !! !!### Normal Code: @@ -2285,16 +2287,16 @@ subroutine polyroots(n, p, wr, wi, info) implicit none - integer, intent(in) :: n !! polynomial degree - real(wp), intent(in) :: p(n + 1) !! polynomial coefficients array, in order of decreasing powers - real(wp), intent(out) :: wr(n) !! real parts of roots - real(wp), intent(out) :: wi(n) !! imaginary parts of roots - integer, intent(out) :: info !! output from the lapack solver. - !! if `info=0` the routine converged. - !! if `info=-999`, then the leading coefficient is zero. + integer, intent(in) :: n !! polynomial degree + real(wp), dimension(n+1), intent(in) :: p !! polynomial coefficients array, in order of decreasing powers + real(wp), dimension(n), intent(out) :: wr !! real parts of roots + real(wp), dimension(n), intent(out) :: wi !! imaginary parts of roots + integer, intent(out) :: info !! output from the lapack solver. + !! if `info=0` the routine converged. + !! if `info=-999`, then the leading coefficient is zero. integer :: i !! counter - real(wp), allocatable, dimension(:, :) :: a !! companion matrix + real(wp), allocatable, dimension(:,:) :: a !! companion matrix real(wp), allocatable, dimension(:) :: work !! work array real(wp), dimension(1) :: vl, vr !! not used here @@ -2373,9 +2375,9 @@ subroutine cpqr79(ndeg, coeff, root, ierr) implicit none integer, intent(in) :: ndeg !! degree of polynomial - complex(wp), intent(in) :: coeff(*) !! `(NDEG+1)` coefficients in descending order. i.e., - !! `P(Z)= COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)` - complex(wp), intent(out) :: root(*) !! `(NDEG)` vector of roots + complex(wp), dimension(ndeg+1), intent(in) :: coeff !! `(NDEG+1)` coefficients in descending order. i.e., + !! `P(Z)= COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)` + complex(wp), dimension(ndeg), intent(out) :: root !! `(NDEG)` vector of roots integer, intent(out) :: ierr !! Output Error Code. !! !!### Normal Code: @@ -2808,19 +2810,19 @@ subroutine newton_root_polish(n, p, zr, zi, ftol, ztol, maxiter, istat) implicit none - integer, intent(in) :: n !! degree of polynomial - real(wp), intent(in) :: p(n + 1) !! vector of coefficients in order of decreasing powers - real(wp), intent(inout) :: zr !! output vector of real parts of the zeros - real(wp), intent(inout) :: zi !! output vector of imaginary parts of the zeros - real(wp), intent(in) :: ftol !! convergence tolerance for the root - real(wp), intent(in) :: ztol !! convergence tolerance for `x` - integer, intent(in) :: maxiter !! maximum number of iterations - integer, intent(out) :: istat !! status flag: - !! - !! * 0 = converged in `f` - !! * 1 = converged in `x` - !! * -1 = singular - !! * -2 = max iterations reached + integer, intent(in) :: n !! degree of polynomial + real(wp), dimension(n+1), intent(in) :: p !! vector of coefficients in order of decreasing powers + real(wp), intent(inout) :: zr !! real part of the zero + real(wp), intent(inout) :: zi !! imaginary part of the zero + real(wp), intent(in) :: ftol !! convergence tolerance for the root + real(wp), intent(in) :: ztol !! convergence tolerance for `x` + integer, intent(in) :: maxiter !! maximum number of iterations + integer, intent(out) :: istat !! status flag: + !! + !! * 0 = converged in `f` + !! * 1 = converged in `x` + !! * -1 = singular + !! * -2 = max iterations reached complex(wp) :: z, f, z_prev, z_best, f_best, dfdx integer :: i !! counter diff --git a/test/rpoly_test.f90 b/test/rpoly_test.f90 index 044809d..9a595fd 100644 --- a/test/rpoly_test.f90 +++ b/test/rpoly_test.f90 @@ -14,7 +14,6 @@ program rpoly_test real(wp),dimension(:),allocatable :: p, zr, zi, s complex(wp),dimension(:),allocatable :: r, cp - complex(wp),dimension(:),allocatable :: t !! work array for [[rpzero]] integer :: degree, i, istatus, icase, n integer,dimension(:),allocatable :: seed real(wp) :: detil @@ -103,7 +102,7 @@ program rpoly_test write(*, '(/A,1x,i3)') 'rpzero' write(*, '(a)') ' real part imaginary part root' istatus = 0 ! no estimates input - call rpzero(degree,p,r,t,istatus,s) + call rpzero(degree,p,r,istatus,s) if (istatus/=0) error stop ' ** failure in rpzero **' call check_results(real(r,wp), aimag(r), degree) @@ -125,7 +124,7 @@ program rpoly_test write(*, '(/A,1x,i3)') 'qr_algeq_solver' write(*, '(a)') ' real part imaginary part root' - call qr_algeq_solver(degree,p,zr,zi,detil,istatus) + call qr_algeq_solver(degree,p,zr,zi,istatus,detil=detil) if (istatus/=0) error stop ' ** failure in qr_algeq_solver **' call check_results(zr, zi, degree) @@ -154,7 +153,6 @@ subroutine allocate_arrays(d) if (allocated(s)) deallocate(s) if (allocated(r)) deallocate(r) if (allocated(cp)) deallocate(cp) - if (allocated(t)) deallocate(t) allocate(p(max_degree+1)) allocate(zr(max_degree+1)) @@ -162,7 +160,6 @@ subroutine allocate_arrays(d) allocate(s(max_degree+1)) allocate(r(max_degree+1)) allocate(cp(max_degree+1)) - allocate(t(6*(max_degree+1))) end subroutine allocate_arrays !******************************************************************** From 51ec610b92edcb72b5ed1276a2f4be126968d6ad Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 18 Sep 2022 10:44:17 -0500 Subject: [PATCH 2/3] some minor cleanups test updates --- src/polyroots_module.F90 | 50 +++++++++++++++++----------------- test/dcbcrt_test.f90 | 4 +-- test/rpoly_test.f90 | 59 +++++++++++++++++++++++++++++++++++----- 3 files changed, 79 insertions(+), 34 deletions(-) diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index 089beb0..ece01ae 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -97,10 +97,10 @@ subroutine rpoly(op, degree, zeror, zeroi, istat) real(wp), parameter :: deg2rad = pi/180.0_wp real(wp), parameter :: cosr = cos(94.0_wp*deg2rad) real(wp), parameter :: sinr = sin(86.0_wp*deg2rad) - real(wp), parameter :: base = radix(0.0_wp) + real(wp), parameter :: base = radix(1.0_wp) real(wp), parameter :: eta = eps - real(wp), parameter :: infin = huge(0.0_wp) - real(wp), parameter :: smalno = tiny(0.0_wp) + real(wp), parameter :: infin = huge(1.0_wp) + real(wp), parameter :: smalno = tiny(1.0_wp) real(wp), parameter :: sqrthalf = sqrt(0.5_wp) real(wp), parameter :: are = eta !! unit error in + real(wp), parameter :: mre = eta !! unit error in * @@ -642,17 +642,17 @@ subroutine calcsc(type) a3 = (a + g)*e + h*(b/d) a1 = b*f - a a7 = (f + u)*a + h - return + else + type = 1 + ! type=1 indicates that all formulas are divided by c + e = a/c + f = d/c + g = u*e + h = v*b + a3 = a*e + (h/c + g)*b + a1 = b - a*(d/c) + a7 = a + g*d + h*f end if - type = 1 - ! type=1 indicates that all formulas are divided by c - e = a/c - f = d/c - g = u*e - h = v*b - a3 = a*e + (h/c + g)*b - a1 = b - a*(d/c) - a7 = a + g*d + h*f end subroutine calcsc @@ -686,15 +686,15 @@ subroutine nextk(type) do i = 3, n k(i) = a3*qk(i - 2) - a7*qp(i - 1) + qp(i) end do - return - end if - ! use unscaled form of the recurrence if type is 3 - k(1) = 0.0_wp - k(2) = 0.0_wp - do i = 3, n - k(i) = qk(i - 2) - end do + else + ! use unscaled form of the recurrence if type is 3 + k(1) = 0.0_wp + k(2) = 0.0_wp + do i = 3, n + k(i) = qk(i - 2) + end do + end if end subroutine nextk @@ -1289,21 +1289,21 @@ subroutine qr_algeq_solver(n, c, zr, zi, istatus, detil) ! balancing the a itself. call balance_companion(n, a) - ! compute the frobenius norm of the balanced companion matrix a. - afnorm = frobenius_norm_companion(n, a) - ! qr iterations from a. call hqr_eigen_hessenberg(n, a, zr, zi, cnt, istatus) if (istatus /= 0) then write (*, '(A,1X,I4)') 'abnormal return from hqr_eigen_hessenberg. istatus=', istatus if (istatus == 1) write (*, '(A)') 'matrix is completely zero.' - if (istatus == 2) write (*, '(A)') 'qr iteration does not converged.' + if (istatus == 2) write (*, '(A)') 'qr iteration did not converge.' if (istatus > 3) write (*, '(A)') 'arguments violate the condition.' return end if if (present(detil)) then + ! compute the frobenius norm of the balanced companion matrix a. + afnorm = frobenius_norm_companion(n, a) + ! count the total qr iteration. iter = 0 do i = 1, n diff --git a/test/dcbcrt_test.f90 b/test/dcbcrt_test.f90 index ffbd816..7dc5c1e 100644 --- a/test/dcbcrt_test.f90 +++ b/test/dcbcrt_test.f90 @@ -23,7 +23,7 @@ program dcbcrt_test do i = 1, 3 z = cmplx(zr(i), zi(i), wp) root = a(1) + a(2)*z + a(3)*z**2 + a(4)*z**3 - write(*,*) 'root is: ', real(root,wp) + write(*,*) 'root is: ', zr(i), zi(i), abs(root) if (abs(root) > 100*epsilon(1.0_wp)) error stop 'Error: insufficient accuracy' end do @@ -32,7 +32,7 @@ program dcbcrt_test do i = 1, 2 z = cmplx(zr(i), zi(i), wp) root = a(1) + a(2)*z + a(3)*z**2 - write(*,*) 'root is: ', real(root,wp) + write(*,*) 'root is: ', zr(i), zi(i), abs(root) if (abs(root) > 100*epsilon(1.0_wp)) error stop 'Error: insufficient accuracy' end do diff --git a/test/rpoly_test.f90 b/test/rpoly_test.f90 index 9a595fd..8210e3c 100644 --- a/test/rpoly_test.f90 +++ b/test/rpoly_test.f90 @@ -9,7 +9,7 @@ program rpoly_test implicit none - integer,parameter :: max_degree = 10 + integer,parameter :: max_degree = 10 !! max degree polynomials to test for random cases integer,parameter :: n_cases = 30 !! number of cases to run real(wp),dimension(:),allocatable :: p, zr, zi, s @@ -79,6 +79,13 @@ program rpoly_test call allocate_arrays(3) p = [ -8.0e18_wp,3.0e12_wp,5.0e6_wp,1.0_wp] + case(13) + call allocate_arrays(3) + p = [4.0_wp, 3.0_wp, 2.0_wp, 1.0_wp] + case(14) + call allocate_arrays(2) + p = [3.0_wp, 2.0_wp, 1.0_wp] + case default ! random coefficients call allocate_arrays(get_random_integer_number(3,max_degree)) @@ -90,6 +97,22 @@ program rpoly_test write(*,'(A,1X,I3)') ' Degree: ', degree write(*,'(A,1X/,*(g23.15/))') ' Coefficients: ', p(1:degree+1) + if (degree==2) then + ! also test this one (only for quadratic equations): + write(*, '(A,1x,i3)') 'dqdcrt' + write(*, '(a)') ' real part imaginary part root' + call dqdcrt(reverse(p), zr, zi) + call check_results(zr, zi, degree) + end if + + if (degree==3) then + ! also test this one (only for cubic equations): + write(*, '(A,1x,i3)') 'dcbcrt' + write(*, '(a)') ' real part imaginary part root' + call dcbcrt(reverse(p), zr, zi) + call check_results(zr, zi, degree) + end if + write(*, '(A,1x,i3)') 'rpoly' write(*, '(a)') ' real part imaginary part root' call rpoly(p, degree, zr, zi, istatus) @@ -140,6 +163,28 @@ program rpoly_test contains + !******************************************************************** + pure function reverse(x) result(y) + + !! reverse a `real(wp)` vector + + implicit none + + real(wp), dimension(:), intent(in) :: x + real(wp), dimension(size(x)) :: y + + integer :: i !! counter + integer :: n !! size of `x` + + n = size(x) + + do i = 1, n + y(i) = x(n-i+1) + end do + + end function reverse + !******************************************************************** + !******************************************************************** subroutine allocate_arrays(d) @@ -154,12 +199,12 @@ subroutine allocate_arrays(d) if (allocated(r)) deallocate(r) if (allocated(cp)) deallocate(cp) - allocate(p(max_degree+1)) - allocate(zr(max_degree+1)) - allocate(zi(max_degree+1)) - allocate(s(max_degree+1)) - allocate(r(max_degree+1)) - allocate(cp(max_degree+1)) + allocate(p(degree+1)) + allocate(zr(degree+1)) + allocate(zi(degree+1)) + allocate(s(degree+1)) + allocate(r(degree+1)) + allocate(cp(degree+1)) end subroutine allocate_arrays !******************************************************************** From 8f4d15b8a0f3ca038e436745134526675d0937c6 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 18 Sep 2022 14:44:34 -0500 Subject: [PATCH 3/3] added cmplx_roots_gen (Skowron & Gould) fixed an issue in qr_algeq_solver for quad precision. public only procedures for ford documentation. --- LICENSE.md | 185 +++++++++- README.md | 5 +- ford.md | 9 +- src/polyroots_module.F90 | 719 ++++++++++++++++++++++++++++++++++++++- test/rpoly_test.f90 | 39 ++- 5 files changed, 939 insertions(+), 18 deletions(-) diff --git a/LICENSE.md b/LICENSE.md index 5f0f282..7690b03 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -26,4 +26,187 @@ ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +---------------------------------------------------------------------------------- + +"General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses" +Copyright 2012 Jan Skowron & Andrew Gould + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + diff --git a/README.md b/README.md index b828559..c256d2d 100644 --- a/README.md +++ b/README.md @@ -21,11 +21,12 @@ Method name | Polynomial type | Coefficients | Roots | Reference [`rpoly`](https://jacobwilliams.github.io/polyroots-fortran/proc/rpoly.html) | General | real | complex | [Jenkins & Traub](https://dl.acm.org/doi/10.1145/355637.355643) [`dcbcrt`](https://jacobwilliams.github.io/polyroots-fortran/proc/dcbcrt.html) | Cubic | real | complex | [NSWC Library](https://github.com/jacobwilliams/nswc) [`dqdcrt`](https://jacobwilliams.github.io/polyroots-fortran/proc/dqdcrt.html) | Quadratic | real | complex | [NSWC Library](https://github.com/jacobwilliams/nswc) -[`qr_algeq_solver`](https://jacobwilliams.github.io/polyroots-fortran/proc/qr_algeq_solver.html) | General | real | complex | [Edelman & Murakami (1995)](https://www.ams.org/journals/mcom/1995-64-210/S0025-5718-1995-1262279-2/S0025-5718-1995-1262279-2.pdf) [`rpzero`](https://jacobwilliams.github.io/polyroots-fortran/proc/rpzero.html) | General | real | complex | [SLATEC](https://netlib.org/slatec/src/rpzero.f) [`cpzero`](https://jacobwilliams.github.io/polyroots-fortran/proc/cpzero.html) | General | complex | complex | [SLATEC](https://netlib.org/slatec/src/cpzero.f) [`rpqr79`](https://jacobwilliams.github.io/polyroots-fortran/proc/rpqr79.html) | General | real | complex | [SLATEC](https://netlib.org/slatec/src/rpqr79.f) [`cpqr79`](https://jacobwilliams.github.io/polyroots-fortran/proc/cpqr79.html) | General | complex | complex | [SLATEC](https://netlib.org/slatec/src/cpqr79.f) +[`qr_algeq_solver`](https://jacobwilliams.github.io/polyroots-fortran/proc/qr_algeq_solver.html) | General | real | complex | [Edelman & Murakami (1995)](https://www.ams.org/journals/mcom/1995-64-210/S0025-5718-1995-1262279-2/S0025-5718-1995-1262279-2.pdf) +[`cmplx_roots_gen`](https://jacobwilliams.github.io/polyroots-fortran/proc/cmplx_roots_gen.html) | General | complex | complex | [Skowron & Gould (2012)](http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/) [`polyroots`](https://jacobwilliams.github.io/polyroots-fortran/proc/polyroots.html) | General | real | complex | -- ## Compiling @@ -109,7 +110,7 @@ The polyroots-fortran source code and related files and documentation are distri * [PA03](https://www.hsl.rl.ac.uk/archive/specs/pa03.pdf) HSL Archive code for computing all the roots of a cubic polynomial * [PA05](https://www.hsl.rl.ac.uk/archive/specs/pa05.pdf) HSL Archive code for computing all the roots of a quartic polynomial * [PA16](https://www.hsl.rl.ac.uk/catalogue/pa16.html), [PA17](https://www.hsl.rl.ac.uk/catalogue/pa17.html) HSL codes for computing zeros of polynomials using method of Madsen and Reid -* `POLZEROS` - Polynomial zeros using Aberth's method, [_Numer Algor_ **13**, 179–200 (1996).](https://doi.org/10.1007/BF02207694) Fortran source code at Netlib: [`numeralgo/na10`](https://netlib.org/numeralgo/na10), several versions provided also by [Alan Miller](https://jblevins.org/mirror/amiller/) +* `POLZEROS` - Polynomial zeros using Aberth's method, [_Numer Algor_ **13**, 179–200 (1996).](https://doi.org/10.1007/BF02207694) Fortran source code at Netlib: [`numeralgo/na10`](https://netlib.org/numeralgo/na10), several versions provided also by [Alan Miller](https://jblevins.org/mirror/amiller/) ## See also diff --git a/ford.md b/ford.md index 3aeea23..128d850 100644 --- a/ford.md +++ b/ford.md @@ -11,15 +11,10 @@ predocmark: < docmark_alt: docmark: ! display: public - private source: true graph: false -search: false +search: true preprocessor: gfortran -E -exclude: pyplot_module.f90 - face.F90 -exclude_dir: ./src/tests -extra_mods: pyplot_module:https://github.com/jacobwilliams/pyplot-fortran - iso_fortran_env:https://gcc.gnu.org/onlinedocs/gfortran/ISO_005fFORTRAN_005fENV.html +extra_mods: iso_fortran_env:https://gcc.gnu.org/onlinedocs/gfortran/ISO_005fFORTRAN_005fENV.html {!README.md!} \ No newline at end of file diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index ece01ae..e8824d0 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -49,14 +49,13 @@ module polyroots_module public :: rpqr79 public :: cpqr79 public :: qr_algeq_solver + public :: cmplx_roots_gen ! special polynomial routines: public :: dcbcrt public :: dqdcrt ! utility routines: - public :: dcbrt - public :: cpevl public :: newton_root_polish contains @@ -1346,7 +1345,7 @@ subroutine balance_companion(n, a) !! blancing the unsymmetric matrix `a`. !! !! this fortran code is based on the algol code "balance" from paper: - !! "balancing a matrixfor calculation of eigenvalues and eigenvectors" + !! "balancing a matrix for calculation of eigenvalues and eigenvectors" !! by b.n.parlett and c.reinsch, numer. math. 13, 293-304(1969). !! !! note: the only non-zero elements of the companion matrix are touched. @@ -1468,7 +1467,7 @@ subroutine hqr_eigen_hessenberg(n0, h, wr, wi, cnt, istatus) !! parts in the array wr(1:n0) and the imaginary parts in the !! array wi(1:n0). !! the procedure fails if any eigenvalue takes more than - !! 30 iterations. + !! `maxiter` iterations. implicit none @@ -1483,6 +1482,13 @@ subroutine hqr_eigen_hessenberg(n0, h, wr, wi, cnt, istatus) real(wp) :: p, q, r, s, t, w, x, y, z logical :: notlast, found +#if REAL128 + integer, parameter :: maxiter = 100 !! max iterations. It seems we need more than 30 + !! for quad precision (see test case 11) +#else + integer, parameter :: maxiter = 30 !! max iterations +#endif + ! note: n is changing in this subroutine. n = n0 istatus = 0 @@ -1544,7 +1550,7 @@ subroutine hqr_eigen_hessenberg(n0, h, wr, wi, cnt, istatus) n = n - 2 cycle main else - if (its == 30) then + if (its == maxiter) then ! 30 for the orignial double precision code istatus = 1 return end if @@ -2909,6 +2915,709 @@ end subroutine func end subroutine newton_root_polish !***************************************************************************************** +!***************************************************************************************** +!> +! This subroutine finds roots of a complex polynomial. +! It uses a new dynamic root finding algorithm (see the Paper). +! +! It can use Laguerre's method (subroutine [[cmplx_laguerre]]) +! or Laguerre->SG->Newton method (subroutine +! [[cmplx_laguerre2newton]] - this is default choice) to find +! roots. It divides polynomial one by one by found roots. At the +! end it finds last root from Viete's formula for quadratic +! equation. Finally, it polishes all found roots using a full +! polynomial and Newton's or Laguerre's method (default is +! Laguerre's - subroutine [[cmplx_laguerre]]). +! You can change default choices by commenting out and uncommenting +! certain lines in the code below. +! +!### Reference +! * J. Skowron & A. Gould, +! "[General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses](https://arxiv.org/pdf/1203.1034.pdf)" +! (2012) +! +!### History +! * Original code here (Apache license): http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/ +! * Jacob Williams, 9/18/2022 : refactored this code a bit +! +!### Notes: +! +! * we solve for the last root with Viete's formula rather +! than doing full Laguerre step (which is time consuming +! and unnecessary) +! * we do not introduce any preference to real roots +! * in Laguerre implementation we omit unneccesarry calculation of +! absolute values of denominator +! * we do not sort roots. + +subroutine cmplx_roots_gen(degree, poly, roots, polish_roots_after, use_roots_as_starting_points) + + implicit none + + integer, intent(in) :: degree !! degree of the polynomial and size of 'roots' array + complex(kind=wp), dimension(degree+1), intent(in) :: poly !! coeffs of the polynomial, in order of increasing powers. + complex(kind=wp), dimension(degree), intent(inout) :: roots !! array which will hold all roots that had been found. + !! If the flag 'use_roots_as_starting_points' is set to + !! .true., then instead of point (0,0) we use value from + !! this array as starting point for cmplx_laguerre + logical, intent(in), optional :: polish_roots_after !! after all roots have been found by dividing + !! original polynomial by each root found, + !! you can opt in to polish all roots using full + !! polynomial. [default is false] + logical, intent(in), optional :: use_roots_as_starting_points !! usually we start Laguerre's + !! method from point (0,0), but you can decide to use the + !! values of 'roots' array as starting point for each new + !! root that is searched for. This is useful if you have + !! very rough idea where some of the roots can be. + !! [default is false] + + complex(kind=wp), dimension(:), allocatable :: poly2 !! `degree+1` array + integer :: i, n, iter + logical :: success + complex(kind=wp) :: coef, prev + + integer, parameter :: MAX_ITERS=50 + ! constants needed to break cycles in the scheme + integer, parameter :: FRAC_JUMP_EVERY=10 + integer, parameter :: FRAC_JUMP_LEN=10 + real(kind=wp), dimension(FRAC_JUMP_LEN), parameter :: FRAC_JUMPS=& + [0.64109297_wp, & + 0.91577881_wp, 0.25921289_wp, 0.50487203_wp, & + 0.08177045_wp, 0.13653241_wp, 0.306162_wp , & + 0.37794326_wp, 0.04618805_wp, 0.75132137_wp] !! some random numbers + real(kind=wp), parameter :: FRAC_ERR = 10.0_wp * eps !! fractional error (see. Adams 1967 Eqs 9 and 10) [2.0d-15 in original code] + complex(kind=wp), parameter :: zero = cmplx(0.0_wp,0.0_wp,wp) + complex(kind=wp), parameter :: c_one=cmplx(1.0_wp,0.0_wp,wp) + + ! initialize starting points + if (present(use_roots_as_starting_points)) then + if (.not.use_roots_as_starting_points) roots = zero + else + roots = zero + end if + + ! skip small degree polynomials from doing Laguerre's method + if (degree<=1) then + if (degree==1) roots(1)=-poly(1)/poly(2) + return + endif + + allocate(poly2(degree+1)) + poly2=poly + + do n=degree, 3, -1 + + ! find root with Laguerre's method + !call cmplx_laguerre(poly2, n, roots(n), iter, success) + ! or + ! find root with (Laguerre's method -> SG method -> Newton's method) + call cmplx_laguerre2newton(poly2, n, roots(n), iter, success, 2) + if (.not.success) then + roots(n)=zero + call cmplx_laguerre(poly2, n, roots(n), iter, success) + endif + + ! divide the polynomial by this root + coef=poly2(n+1) + do i=n,1,-1 + prev=poly2(i) + poly2(i)=coef + coef=prev+roots(n)*coef + enddo + ! variable coef now holds a remainder - should be close to 0 + + enddo + + ! find all but last root with Laguerre's method + !call cmplx_laguerre(poly2, 2, roots(2), iter, success) + ! or + call cmplx_laguerre2newton(poly2, 2, roots(2), iter, success, 2) + if (.not.success) then + call solve_quadratic_eq(roots(2),roots(1),poly2) + else + ! calculate last root from Viete's formula + roots(1)=-(roots(2)+poly2(2)/poly2(3)) + endif + + if (present(polish_roots_after)) then + if (polish_roots_after) then + do n=1, degree ! polish roots one-by-one with a full polynomial + call cmplx_laguerre(poly, degree, roots(n), iter, success) + !call cmplx_newton_spec(poly, degree, roots(n), iter, success) + enddo + endif + end if + + contains + + recursive subroutine cmplx_laguerre(poly, degree, root, iter, success) + + ! Subroutine finds one root of a complex polynomial using + ! Laguerre's method. In every loop it calculates simplified + ! Adams' stopping criterion for the value of the polynomial. + ! + ! For a summary of the method go to: + ! http://en.wikipedia.org/wiki/Laguerre's_method + + implicit none + + integer, intent(in) :: degree !! a degree of the polynomial + complex(kind=wp), dimension(degree+1), intent(in) :: poly !! an array of polynomial cooefs + !! length = degree+1, poly(1) is constant + !!``` + !! 1 2 3 + !! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... + !!``` + integer, intent(out) :: iter !! number of iterations performed (the number of polynomial + !! evaluations and stopping criterion evaluation) + complex(kind=wp), intent(inout) :: root !! input: guess for the value of a root + !! output: a root of the polynomial + !! + !! Uses 'root' value as a starting point (!!!!!) + !! Remember to initialize 'root' to some initial guess or to + !! point (0,0) if you have no prior knowledge. + + logical, intent(out) :: success !! is false if routine reaches maximum number of iterations + + real(kind=wp) :: faq !! jump length + complex(kind=wp) :: p !! value of polynomial + complex(kind=wp) :: dp !! value of 1st derivative + complex(kind=wp) :: d2p_half !! value of 2nd derivative + integer :: i, k + logical :: good_to_go + complex(kind=wp) :: denom, denom_sqrt, dx, newroot + real(kind=wp) :: ek, absroot, abs2p + complex(kind=wp) :: fac_netwon, fac_extra, F_half, c_one_nth + real(kind=wp) :: one_nth, n_1_nth, two_n_div_n_1 + real(kind=wp) :: stopping_crit2 + + iter=0 + success=.true. + + ! next if-endif block is an EXTREME failsafe, not usually needed, and thus turned off in this version. + !if (.false.) then ! change false-->true if you would like to use caution about having first coefficient == 0 + if (degree<0) then + write(*,*) 'Error: cmplx_laguerre: degree<0' + return + endif + if (poly(degree+1)==zero) then + if (degree==0) return + call cmplx_laguerre(poly, degree-1, root, iter, success) + return + endif + if (degree<=1) then + if (degree==0) then ! we know from previous check than poly(1) not equal zero + success=.false. + write(*,*) 'Warning: cmplx_laguerre: degree=0 and poly(1)/=0, no roots' + return + else + root=-poly(1)/poly(2) + return + endif + endif + !endif + ! end EXTREME failsafe + + good_to_go=.false. + one_nth=1.0_wp/degree + n_1_nth=(degree-1.0_wp)*one_nth + two_n_div_n_1=2.0_wp/n_1_nth + c_one_nth=cmplx(one_nth,0.0_wp,wp) + + do i=1,MAX_ITERS + ! prepare stoping criterion + ek=abs(poly(degree+1)) + absroot=abs(root) + ! calculate value of polynomial and its first two derivatives + p =poly(degree+1) + dp =zero + d2p_half=zero + do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives + d2p_half=dp + d2p_half*root + dp =p + dp*root + p =poly(k)+p*root ! b_k + ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", + ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 + ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf + ! Eq 8. + ek=absroot*ek+abs(p) + enddo + iter=iter+1 + + abs2p=real(conjg(p)*p) + if (abs2p==0.0_wp) return + stopping_crit2=(FRAC_ERR*ek)**2 + if (abs2p=0.0_wp) then + ! real part of a square root is positive for probably all compilers. You can + ! test this on your compiler and if so, you can omit this check + denom=c_one_nth+n_1_nth*denom_sqrt + else + denom=c_one_nth-n_1_nth*denom_sqrt + endif + endif + if (denom==zero) then !test if demoninators are > 0.0 not to divide by zero + dx=(absroot+1.0_wp)*exp(cmplx(0.0_wp,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi,wp)) ! make some random jump + else + dx=fac_netwon/denom + !dx=degree/denom + endif + + newroot=root-dx + if (newroot==root) return ! nothing changes -> return + if (good_to_go) then ! this was jump already after stopping criterion was met + root=newroot + return + endif + + if (mod(i,FRAC_JUMP_EVERY)==0) then ! decide whether to do a jump of modified length (to break cycles) + faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) + newroot=root-faq*dx ! do jump of some semi-random length (0=0.0_wp ) then ! scalar product to decide the sign yielding bigger magnitude + x0=-0.5_wp*(b+delta) + else + x0=-0.5_wp*(b-delta) + endif + if (x0==cmplx(0.0_wp,0.0_wp,wp)) then + x1=cmplx(0.0_wp,0.0_wp,wp) + else ! Viete's formula + x1=c/x0 + x0=x0/a + endif + + end subroutine solve_quadratic_eq + + recursive subroutine cmplx_laguerre2newton(poly, degree, root, iter, success, starting_mode) + + ! Subroutine finds one root of a complex polynomial using + ! Laguerre's method, Second-order General method and Newton's + ! method - depending on the value of function F, which is a + ! combination of second derivative, first derivative and + ! value of polynomial [F=-(p"*p)/(p'p')]. + ! + ! Subroutine has 3 modes of operation. It starts with mode=2 + ! which is the Laguerre's method, and continues until F + ! becames F<0.50, at which point, it switches to mode=1, + ! i.e., SG method (see paper). While in the first two + ! modes, routine calculates stopping criterion once per every + ! iteration. Switch to the last mode, Newton's method, (mode=0) + ! happens when becomes F<0.05. In this mode, routine calculates + ! stopping criterion only once, at the beginning, under an + ! assumption that we are already very close to the root. + ! If there are more than 10 iterations in Newton's mode, + ! it means that in fact we were far from the root, and + ! routine goes back to Laguerre's method (mode=2). + ! + ! For a summary of the method see the paper: Skowron & Gould (2012) + + implicit none + + integer, intent(in) :: degree !! a degree of the polynomial + complex(kind=wp), dimension(degree+1), intent(in) :: poly !! is an array of polynomial cooefs + !! length = degree+1, poly(1) is constant + !!``` + !! 1 2 3 + !! poly(1) x^0 + poly(2) x^1 + poly(3) x^2 + ... + !!``` + complex(kind=wp), intent(inout) :: root !! input: guess for the value of a root + !! output: a root of the polynomial + !! + !! Uses 'root' value as a starting point (!!!!!) + !! Remember to initialize 'root' to some initial guess or to + !! point (0,0) if you have no prior knowledge. + integer, intent(in) :: starting_mode !! this should be by default = 2. However if you + !! choose to start with SG method put 1 instead. + !! Zero will cause the routine to + !! start with Newton for first 10 iterations, and + !! then go back to mode 2. + integer, intent(out) :: iter !! number of iterations performed (the number of polynomial + !! evaluations and stopping criterion evaluation) + logical, intent(out) :: success !! is false if routine reaches maximum number of iterations + + real(kind=wp) :: faq ! jump length + complex(kind=wp) :: p ! value of polynomial + complex(kind=wp) :: dp ! value of 1st derivative + complex(kind=wp) :: d2p_half ! value of 2nd derivative + integer :: i, j, k + logical :: good_to_go + complex(kind=wp) :: denom, denom_sqrt, dx, newroot + real(kind=wp) :: ek, absroot, abs2p, abs2_F_half + complex(kind=wp) :: fac_netwon, fac_extra, F_half, c_one_nth + real(kind=wp) :: one_nth, n_1_nth, two_n_div_n_1 + integer :: mode + real(kind=wp) :: stopping_crit2 + + iter=0 + success=.true. + stopping_crit2 = 0.0_wp ! value not important, will be initialized anyway on the first loop (because mod(1,10)==1) + + ! next if-endif block is an EXTREME failsafe, not usually needed, and thus turned off in this version. + !if(.false.)then ! change false-->true if you would like to use caution about having first coefficient == 0 + if (degree<0) then + write(*,*) 'Error: cmplx_laguerre2newton: degree<0' + return + endif + if (poly(degree+1)==zero) then + if (degree==0) return + call cmplx_laguerre2newton(poly, degree-1, root, iter, success, starting_mode) + return + endif + if (degree<=1) then + if (degree==0) then ! we know from previous check than poly(1) not equal zero + success=.false. + write(*,*) 'Warning: cmplx_laguerre2newton: degree=0 and poly(1)/=0, no roots' + return + else + root=-poly(1)/poly(2) + return + endif + endif + !endif + ! end EXTREME failsafe + + j=1 + good_to_go=.false. + mode=starting_mode ! mode=2 full laguerre, mode=1 SG, mode=0 newton + + do ! infinite loop, just to be able to come back from newton, if more than 10 iteration there + + !------------------------------------------------------------- mode 2 + if (mode>=2) then ! LAGUERRE'S METHOD + one_nth=1.0_wp/degree + n_1_nth=(degree-1.0_wp)*one_nth + two_n_div_n_1=2.0_wp/n_1_nth + c_one_nth=cmplx(one_nth,0.0_wp,wp) + + do i=1,MAX_ITERS ! + faq=1.0_wp + + ! prepare stoping criterion + ek=abs(poly(degree+1)) + absroot=abs(root) + ! calculate value of polynomial and its first two derivatives + p =poly(degree+1) + dp =zero + d2p_half=zero + do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives + d2p_half=dp + d2p_half*root + dp =p + dp*root + p =poly(k)+p*root ! b_k + ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", + ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 + ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf + ! Eq 8. + ek=absroot*ek+abs(p) + enddo + abs2p=real(conjg(p)*p, wp) !abs(p) + iter=iter+1 + if (abs2p==0.0_wp) return + + stopping_crit2=(FRAC_ERR*ek)**2 + if (abs2p=0.0_wp) then + ! real part of a square root is positive for probably all compilers. You can + ! test this on your compiler and if so, you can omit this check + denom=c_one_nth+n_1_nth*denom_sqrt + else + denom=c_one_nth-n_1_nth*denom_sqrt + endif + endif + if (denom==zero) then !test if demoninators are > 0.0 not to divide by zero + dx=(abs(root)+1.0_wp)*exp(cmplx(0.0_wp,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi,wp)) ! make some random jump + else + dx=fac_netwon/denom + endif + + newroot=root-dx + if (newroot==root) return ! nothing changes -> return + if (good_to_go) then ! this was jump already after stopping criterion was met + root=newroot + return + endif + + if (mode/=2) then + root=newroot + j=i+1 ! remember iteration index + exit ! go to Newton's or SG + endif + + if (mod(i,FRAC_JUMP_EVERY)==0) then ! decide whether to do a jump of modified length (to break cycles) + faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) + newroot=root-faq*dx ! do jump of some semi-random length (0=MAX_ITERS) then + success=.false. + return + endif + + endif ! if mode 2 + + !------------------------------------------------------------- mode 1 + if (mode==1) then ! SECOND-ORDER GENERAL METHOD (SG) + + do i=j,MAX_ITERS ! + faq=1.0_wp + + ! calculate value of polynomial and its first two derivatives + p =poly(degree+1) + dp =zero + d2p_half=zero + if (mod(i-j,10)==0) then + ! prepare stoping criterion + ek=abs(poly(degree+1)) + absroot=abs(root) + do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives + d2p_half=dp + d2p_half*root + dp =p + dp*root + p =poly(k)+p*root ! b_k + ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", + ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 + ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf + ! Eq 8. + ek=absroot*ek+abs(p) + enddo + stopping_crit2=(FRAC_ERR*ek)**2 + else + do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives + d2p_half=dp + d2p_half*root + dp =p + dp*root + p =poly(k)+p*root ! b_k + enddo + endif + + abs2p=real(conjg(p)*p, wp) !abs(p)**2 + iter=iter+1 + if (abs2p==0.0_wp) return + + if (abs2p 0.0 not to divide by zero + dx=(abs(root)+1.0_wp)*exp(cmplx(0.0_wp,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi,wp)) ! make some random jump + else + fac_netwon=p/dp + fac_extra=d2p_half/dp + F_half=fac_netwon*fac_extra + + abs2_F_half=real(conjg(F_half)*F_half, wp) + if (abs2_F_half<=0.000625_wp) then ! F<0.05, F/2<0.025 + mode=0 ! set Newton's, go there after jump + endif + + dx=fac_netwon*(c_one+F_half) ! SG + endif + + newroot=root-dx + if (newroot==root) return ! nothing changes -> return + if (good_to_go) then ! this was jump already after stopping criterion was met + root=newroot + return + endif + + if (mode/=1) then + root=newroot + j=i+1 ! remember iteration number + exit ! go to Newton's + endif + + if (mod(i,FRAC_JUMP_EVERY)==0) then ! decide whether to do a jump of modified length (to break cycles) + faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) + newroot=root-faq*dx ! do jump of some semi-random length (0=MAX_ITERS) then + success=.false. + return + endif + + endif ! if mode 1 + + !------------------------------------------------------------- mode 0 + if (mode==0) then ! NEWTON'S METHOD + + do i=j,j+10 ! do only 10 iterations the most, then go back to full Laguerre's + faq=1.0_wp + + ! calculate value of polynomial and its first two derivatives + p =poly(degree+1) + dp =zero + if (i==j) then ! calculate stopping crit only once at the begining + ! prepare stoping criterion + ek=abs(poly(degree+1)) + absroot=abs(root) + do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives + dp =p + dp*root + p =poly(k)+p*root ! b_k + ! Adams, Duane A., 1967, "A stopping criterion for polynomial root finding", + ! Communications of the ACM, Volume 10 Issue 10, Oct. 1967, p. 655 + ! ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/55/CS-TR-67-55.pdf + ! Eq 8. + ek=absroot*ek+abs(p) + enddo + stopping_crit2=(FRAC_ERR*ek)**2 + else ! + do k=degree,1,-1 ! Horner Scheme, see for eg. Numerical Recipes Sec. 5.3 how to evaluate polynomials and derivatives + dp =p + dp*root + p =poly(k)+p*root ! b_k + enddo + endif + abs2p=real(conjg(p)*p, wp) !abs(p)**2 + iter=iter+1 + if (abs2p==0.0_wp) return + + if (abs2p 0.0 not to divide by zero + dx=(abs(root)+1.0_wp)*exp(cmplx(0.0_wp,FRAC_JUMPS(mod(i,FRAC_JUMP_LEN)+1)*2*pi,wp)) ! make some random jump + else + dx=p/dp + endif + + newroot=root-dx + if (newroot==root) return ! nothing changes -> return + if (good_to_go) then + root=newroot + return + endif + + ! this loop is done only 10 times. So skip this check + !if (mod(i,FRAC_JUMP_EVERY)==0) then ! decide whether to do a jump of modified length (to break cycles) + ! faq=FRAC_JUMPS(mod(i/FRAC_JUMP_EVERY-1,FRAC_JUMP_LEN)+1) + ! newroot=root-faq*dx ! do jump of some semi-random length (0=MAX_ITERS) then + ! too many iterations here + success=.false. + return + endif + mode=2 ! go back to Laguerre's. This happens when we were unable to converge in 10 iterations with Newton's + + endif ! if mode 0 + + enddo ! end of infinite loop + + success=.false. + + end subroutine cmplx_laguerre2newton + +end subroutine cmplx_roots_gen +!***************************************************************************************** + !***************************************************************************************** end module polyroots_module !***************************************************************************************** \ No newline at end of file diff --git a/test/rpoly_test.f90 b/test/rpoly_test.f90 index 8210e3c..d79427b 100644 --- a/test/rpoly_test.f90 +++ b/test/rpoly_test.f90 @@ -12,7 +12,7 @@ program rpoly_test integer,parameter :: max_degree = 10 !! max degree polynomials to test for random cases integer,parameter :: n_cases = 30 !! number of cases to run - real(wp),dimension(:),allocatable :: p, zr, zi, s + real(wp),dimension(:),allocatable :: p, zr, zi, s, q complex(wp),dimension(:),allocatable :: r, cp integer :: degree, i, istatus, icase, n integer,dimension(:),allocatable :: seed @@ -97,11 +97,13 @@ program rpoly_test write(*,'(A,1X,I3)') ' Degree: ', degree write(*,'(A,1X/,*(g23.15/))') ' Coefficients: ', p(1:degree+1) + q = reverse(p) ! the following two accept the coefficients in reverse order + if (degree==2) then ! also test this one (only for quadratic equations): write(*, '(A,1x,i3)') 'dqdcrt' write(*, '(a)') ' real part imaginary part root' - call dqdcrt(reverse(p), zr, zi) + call dqdcrt(q, zr, zi) call check_results(zr, zi, degree) end if @@ -109,7 +111,7 @@ program rpoly_test ! also test this one (only for cubic equations): write(*, '(A,1x,i3)') 'dcbcrt' write(*, '(a)') ' real part imaginary part root' - call dcbcrt(reverse(p), zr, zi) + call dcbcrt(q, zr, zi) call check_results(zr, zi, degree) end if @@ -149,6 +151,13 @@ program rpoly_test write(*, '(a)') ' real part imaginary part root' call qr_algeq_solver(degree,p,zr,zi,istatus,detil=detil) if (istatus/=0) error stop ' ** failure in qr_algeq_solver **' + call check_results(real(r,wp), aimag(r), degree) + + write(*, '(/A,1x,i3)') 'cmplx_roots_gen' + write(*, '(a)') ' real part imaginary part root' + cp = reversez(cp) + call cmplx_roots_gen(degree, cp, r) + if (istatus/=0) error stop ' ** failure in cmplx_roots_gen **' call check_results(zr, zi, degree) if (wp /= REAL128) then @@ -185,6 +194,28 @@ pure function reverse(x) result(y) end function reverse !******************************************************************** + !******************************************************************** + pure function reversez(x) result(y) + + !! reverse a `complex(wp)` vector + + implicit none + + complex(wp), dimension(:), intent(in) :: x + complex(wp), dimension(size(x)) :: y + + integer :: i !! counter + integer :: n !! size of `x` + + n = size(x) + + do i = 1, n + y(i) = x(n-i+1) + end do + + end function reversez + !******************************************************************** + !******************************************************************** subroutine allocate_arrays(d) @@ -198,8 +229,10 @@ subroutine allocate_arrays(d) if (allocated(s)) deallocate(s) if (allocated(r)) deallocate(r) if (allocated(cp)) deallocate(cp) + if (allocated(q)) deallocate(q) allocate(p(degree+1)) + allocate(q(degree+1)) allocate(zr(degree+1)) allocate(zi(degree+1)) allocate(s(degree+1))