From 435700f4e497d28d33b9d2664a674dea91c445b9 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Thu, 6 Oct 2022 16:01:55 -0500 Subject: [PATCH 1/5] initial cpolz refactor from MATH77 --- README.md | 1 + src/polyroots_module.F90 | 438 +++++++++++++++++++++++++++++++++++++++ test/polyroots_test.f90 | 5 +- 3 files changed, 443 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 0a602ff..b6598e8 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,7 @@ Method name | Polynomial type | Coefficients | Roots | Reference [`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) [`dpolz`](https://jacobwilliams.github.io/polyroots-fortran/proc/dpolz.html) | General | real | complex | [MATH77 Library](https://netlib.org/math/) +[`cpolz`](https://jacobwilliams.github.io/polyroots-fortran/proc/cpolz.html) | General | complex | complex | [MATH77 Library](https://netlib.org/math/) [`polyroots`](https://jacobwilliams.github.io/polyroots-fortran/proc/polyroots.html) | General | real | complex | [LAPACK](https://netlib.org/lapack/explore-html/index.html) [`cpolyroots`](https://jacobwilliams.github.io/polyroots-fortran/proc/cpolyroots.html) | General | complex | complex | [LAPACK](https://netlib.org/lapack/explore-html/index.html) [`rroots_chebyshev_cubic`](https://jacobwilliams.github.io/polyroots-fortran/proc/rroots_chebyshev_cubic.html) | Cubic | real | complex | [Lebedev (1991)](https://doi.org/10.1515/rnam.1991.6.4.315) diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index 12dea7f..36bfae6 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -56,6 +56,7 @@ module polyroots_module public :: polzeros public :: fpml public :: dpolz + public :: cpolz ! special polynomial routines: public :: dcbcrt @@ -5934,6 +5935,443 @@ subroutine dpolz(ndeg,a,zr,zi,ierr) end subroutine dpolz !***************************************************************************************** +!***************************************************************************************** +!> +! In the discussion below, the notation A([*,],k} should be interpreted +! as the complex number A(k) if A is declared complex, and should be +! interpreted as the complex number A(1,k) + i * A(2,k) if A is not +! declared to be of type complex. Similar statements apply for Z(k). +! +! Given complex coefficients A([*,[1),...,A([*,]NDEG+1) this +! subr computes the NDEG roots of the polynomial +! A([*,]1)*X**NDEG + ... + A([*,]NDEG+1) +! storing the roots as complex numbers in the array Z( ). +! Require NDEG >= 1 and A([*,]1) /= 0. +! +!### Reference +! * Original code from [JPL MATH77 Library](https://netlib.org/math/) +! +!### License +! Copyright (c) 1996 California Institute of Technology, Pasadena, CA. +! ALL RIGHTS RESERVED. +! Based on Government Sponsored Research NAS7-03001. +! +!### History +! * C.L.Lawson & S.Y.Chan, JPL, June 3,1986. +! * 1987-02-25 CPOLZ Lawson Initial code. +! * 1989-10-20 CLL Delcared all variables. +! * 1992-05-11 CLL IERR was not being set when N = 0 or 1. Fixed this. +! * 1995-01-18 CPOLZ Krogh More M77CON for conversion to C. +! * 1995-11-17 CPOLZ Krogh Added M77CON statements for conversion to C +! * 1995-11-17 CPOLZ Krogh Converted SFTRAN to Fortran 77. +! * 1996-03-30 CPOLZ Krogh Added external statement. +! * 1996-04-27 CPOLZ Krogh Changes to use .C. and C%%. +! * 2001-05-25 CPOLZ Krogh Minor change for making .f90 version. +! * 2022-10-06, Jacob Williams modernized this routine + + subroutine cpolz(a,ndeg,z,ierr) + + integer,intent(in) :: ndeg !! degree of the polynomial + complex(wp),intent(in) :: a(ndeg+1) !! contains the complex coefficients of a polynomial + !! high order coefficient first, with a([*,]1)/=0. the + !! real and imaginary parts of the jth coefficient must + !! be provided in a([*],j). the contents of this array will + !! not be modified by the subroutine. + complex(wp),intent(out) :: z(ndeg) !! contains the polynomial roots stored as complex + !! numbers. the real and imaginary parts of the jth root + !! will be stored in z([*,]j). + integer,intent(out) :: ierr !! error flag. set by the subroutine to 0 on normal + !! termination. set to -1 if a([*,]1)=0. set to -2 if ndeg + !! <= 0. set to j > 0 if the iteration count limit + !! has been exceeded and roots 1 through j have not been + !! determined. + + complex(wp) :: temp + integer :: i, j, n + real(wp) :: c, f, g, r, s + logical :: more, first + real(wp) :: h(ndeg,ndeg,2) !! array of work space + + real(wp),parameter :: zero = 0.0_wp + real(wp),parameter :: one = 1.0_wp + real(wp),parameter :: c95 = 0.95_wp + integer,parameter :: base = radix(1.0_wp) !! i1mach(10) + integer,parameter :: b2 = base * base + + if (ndeg <= 0) then + ierr = -2 + write(*,*) 'ndeg <= 0' + return + end if + + if (a(1) == cmplx(zero, zero, wp)) then + ierr = -1 + write(*,*) 'a(*,1) == zero' + return + end if + + n = ndeg + ierr = 0 + + ! build first row of companion matrix. + + do i = 2,n+1 + temp = -(a(i)/a(1)) + h(1,i-1,1) = real(temp,wp) + h(1,i-1,2) = aimag(temp) + end do + + ! extract any exact zero roots and set n = degree of + ! remaining polynomial. + + do j = ndeg,1,-1 + if (h(1,j,1)/=zero .or. h(1,j,2)/=zero) exit + z(j) = zero + n = n - 1 + end do + + ! special for n = 0 or 1. + + if (n == 0) return + if (n == 1) then + z(1) = cmplx(h(1,1,1),h(1,1,2),wp) + return + end if + + ! build rows 2 thru n of the companion matrix. + + do i = 2,n + do j = 1,n + if (j == i-1) then + h(i,j,1) = one + h(i,j,2) = zero + else + h(i,j,1) = zero + h(i,j,2) = zero + end if + end do + end do + + ! ***************** BALANCE THE MATRIX *********************** + ! + ! This is an adaption of the EISPACK subroutine BALANC to + ! the special case of a complex companion matrix. 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). + + ! ********** ITERATIVE LOOP FOR NORM REDUCTION ********** + do + MORE = .FALSE. + DO I = 1, N + ! Compute R = sum of magnitudes in row I skipping diagonal. + ! C = sum of magnitudes in col I skipping diagonal. + IF (I == 1) THEN + R = ABS(H(1,2,1)) + ABS(H(1,2,2)) + DO J = 3,N + R = R + ABS(H(1,J,1)) + ABS(H(1,J,2)) + end do + C = ABS(H(2,1,1)) + ABS(H(2,1,2)) + ELSE + R = ABS(H(I,I-1,1)) + ABS(H(I,I-1,2)) + C = ABS(H(1,I,1)) + ABS(H(1,I,2)) + IF (I /= N) THEN + C = C + ABS(H(I+1,I,1)) + ABS(H(I+1,I,2)) + END IF + END IF + + ! Determine column scale factor, F. + + G = R / BASE + F = ONE + S = C + R + + do + IF (C >= G) exit + F = F * BASE + C = C * B2 + end do + G = R * BASE + do + IF (C < G) exit + F = F / BASE + C = C / B2 + end do + ! Will the factor F have a significant effect ? + + IF ((C + R) / F < C95 * S) THEN + + ! Yes, so do the scaling. + + G = ONE / F + MORE = .TRUE. + + ! Scale Row I + + IF (I == 1) THEN + DO J = 1,N + H(1,J,1) = H(1,J,1)*G + H(1,J,2) = H(1,J,2)*G + end do + ELSE + H(I,I-1,1) = H(I,I-1,1)*G + H(I,I-1,2) = H(I,I-1,2)*G + END IF + + ! Scale Column I + + H(1,I,1) = H(1,I,1) * F + H(1,I,2) = H(1,I,2) * F + IF (I /= N) THEN + H(I+1,I,1) = H(I+1,I,1) * F + H(I+1,I,2) = H(I+1,I,2) * F + END IF + + END IF + end do + IF (.not. MORE) exit + end do + + call scomqr(ndeg,n,1,n,h(1,1,1),h(1,1,2),z,ierr) + + if (ierr /= 0) write(*,*) 'Convergence failure in cpolz' + + end subroutine cpolz + +! Copyright (c) 1996 California Institute of Technology, Pasadena, CA. +! ALL RIGHTS RESERVED. +! Based on Government Sponsored Research NAS7-03001. +!>> 2001-01-24 SCOMQR Krogh CSQRT -> CSQRTX to avoid C lib. conflicts. +!>> 1996-04-27 SCOMQR Krogh Changes to use .C. and C%%. +!>> 1996-03-30 SCOMQR Krogh Added external statement. +!>> 1996-01-18 SCOMQR Krogh Added M77CON statements for conv. to C. +!>> 1995-01-03 SCOMQR WVS Added EXTERNAL CQUO, CSQRT so VAX won't gripe +!>> 1992-03-13 SCOMQR FTK Removed implicit statements. +!>> 1987-11-12 SCOMQR Lawson Initial code. +! ------------------------------------------------------------------ +! Version of the Eispack subr, COMQR, for use in the JPL MATH77 +! library. C. L. Lawson, JPL, 1987 Feb 17. +! ------------------------------------------------------------------ +! +! THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE +! ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN +! AND WILKINSON. +! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). +! THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS +! (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. +! +! THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX +! UPPER HESSENBERG MATRIX BY THE QR METHOD. +! +! ON INPUT- +! +! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL +! ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM +! DIMENSION STATEMENT, +! +! N IS THE ORDER OF THE MATRIX, +! +! LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING +! SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, +! SET LOW=1, IGH=N, +! +! HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, +! RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. +! THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN +! INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN +! THE REDUCTION BY CORTH, IF PERFORMED. +! +! ON OUTPUT- +! +! THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN +! DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE +! CALLING COMQR IF SUBSEQUENT CALCULATION OF +! EIGENVECTORS IS TO BE PERFORMED, +! +! WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, +! RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR +! EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT +! FOR INDICES IERR+1,...,N, +! +! IERR IS SET TO +! ZERO FOR NORMAL RETURN, +! J IF THE J-TH EIGENVALUE HAS NOT BEEN +! DETERMINED AFTER 30 ITERATIONS. +! +! ARITHMETIC IS REAL EXCEPT FOR THE REPLACEMENT OF THE ALGOL +! PROCEDURE CDIV BY COMPLEX DIVISION AND USE OF THE SUBROUTINES +! SQRT AND CMPLX IN COMPUTING COMPLEX SQUARE ROOTS. +! +! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, +! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY + + subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) + + integer :: en,enm1,i,ierr,igh,its,j,l,ll,low,lp1,n,nm + real(wp) :: hi(nm,n),hr(nm,n) + real(wp) :: norm,si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr + complex(wp) :: z(n), z3 + + ierr = 0 + if (low /= igh) then + ! ********** create real subdiagonal elements ********** + l = low + 1 + + do i = l, igh + ll = min(i+1,igh) + if (hi(i,i-1) == 0.0_wp) cycle + norm = abs(cmplx(hr(i,i-1),hi(i,i-1),wp)) + yr = hr(i,i-1) / norm + yi = hi(i,i-1) / norm + hr(i,i-1) = norm + hi(i,i-1) = 0.0_wp + + do j = i, igh + si = yr * hi(i,j) - yi * hr(i,j) + hr(i,j) = yr * hr(i,j) + yi * hi(i,j) + hi(i,j) = si + end do + + do j = low, ll + si = yr * hi(j,i) + yi * hr(j,i) + hr(j,i) = yr * hr(j,i) - yi * hi(j,i) + hi(j,i) = si + end do + end do + end if + ! ********** store roots isolated by cbal ********** + do i = 1, n + if (i >= low .and. i <= igh) cycle + z(i) = cmplx(hr(i,i),hi(i,i),wp) + end do + + en = igh + tr = 0.0_wp + ti = 0.0_wp + + do + ! ********** search for next eigenvalue ********** + if (en < low) return + its = 0 + enm1 = en - 1 + ! ********** look for single small sub-diagonal element + ! for l=en step -1 until low -- ********** + 240 do ll = low, en + l = en + low - ll + if (l == low) exit + if (abs(hr(l,l-1)) <= & + eps * (abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) & + + abs(hr(l,l)) +abs(hi(l,l)))) exit + end do + ! ********** form shift ********** + if (l == en) go to 660 + if (its == 30) go to 1000 + if (its == 10 .or. its == 20) go to 320 + sr = hr(en,en) + si = hi(en,en) + xr = hr(enm1,en) * hr(en,enm1) + xi = hi(enm1,en) * hr(en,enm1) + if (xr == 0.0_wp .and. xi == 0.0_wp) go to 340 + yr = (hr(enm1,enm1) - sr) / 2.0_wp + yi = (hi(enm1,enm1) - si) / 2.0_wp + z3 = sqrt(cmplx(yr**2-yi**2+xr,2.0_wp*yr*yi+xi,wp)) + zzr = real(z3,wp) + zzi = aimag(z3) + if (yr * zzr + yi * zzi < 0.0_wp) then + zzr = -zzr + zzi = -zzi + end if + z3 = cmplx(xr,xi,wp) / cmplx(yr+zzr,yi+zzi,wp) + sr = sr - real(z3,wp) + si = si - aimag(z3) + go to 340 + + ! ********** form exceptional shift ********** + 320 sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2)) + si = 0.0_wp + + 340 do i = low, en + hr(i,i) = hr(i,i) - sr + hi(i,i) = hi(i,i) - si + end do + + tr = tr + sr + ti = ti + si + its = its + 1 + ! ********** reduce to triangle (rows) ********** + lp1 = l + 1 + + do i = lp1, en + sr = hr(i,i-1) + hr(i,i-1) = 0.0_wp + norm = sqrt(hr(i-1,i-1)*hr(i-1,i-1)+hi(i-1,i-1)*hi(i-1,i-1)+sr*sr) + xr = hr(i-1,i-1) / norm + xi = hi(i-1,i-1) / norm + z(i-1) = cmplx(xr,xi,wp) + hr(i-1,i-1) = norm + hi(i-1,i-1) = 0.0_wp + hi(i,i-1) = sr / norm + do j = i, en + yr = hr(i-1,j) + yi = hi(i-1,j) + zzr = hr(i,j) + zzi = hi(i,j) + hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr + hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi + hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr + hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi + end do + end do + + si = hi(en,en) + if (si /= 0.0_wp) then + norm = abs(cmplx(hr(en,en),si,wp)) + sr = hr(en,en) / norm + si = si / norm + hr(en,en) = norm + hi(en,en) = 0.0_wp + end if + ! ********** inverse operation (columns) ********** + do j = lp1, en + xr = real(z(j-1),wp) + xi = aimag(z(j-1)) + do i = l, j + yr = hr(i,j-1) + yi = 0.0 + zzr = hr(i,j) + zzi = hi(i,j) + if (i /= j) then + yi = hi(i,j-1) + hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi + end if + hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi + end do + end do + + if (si == 0.0) go to 240 + + do i = l, en + yr = hr(i,en) + yi = hi(i,en) + hr(i,en) = sr * yr - si * yi + hi(i,en) = sr * yi + si * yr + end do + go to 240 + + ! ********** a root found ********** + 660 continue + z(en) = cmplx(hr(en,en)+tr,hi(en,en)+ti,wp) + en = enm1 + end do + + ! ********** set error -- no convergence to an + ! eigenvalue after 30 iterations ********** + 1000 ierr = en + + end subroutine scomqr + !***************************************************************************************** end module polyroots_module !***************************************************************************************** \ No newline at end of file diff --git a/test/polyroots_test.f90 b/test/polyroots_test.f90 index af4ab51..29f685f 100644 --- a/test/polyroots_test.f90 +++ b/test/polyroots_test.f90 @@ -150,7 +150,10 @@ program polyroots_test call dpolz(degree,p,zr,zi,istatus) call check_results('dpolz', istatus, zr, zi, degree) - ! for now, just test the following two with the real coefficients only: + ! for now, just test the following with the real coefficients only: + + call cpolz(cp,degree,r,istatus) + call check_results('cpolz', istatus, real(r,wp), aimag(r), degree) q = 0.0_wp istatus = 0 From ed9e3b82462cd1f8d8e9e41bc267cd18911d9ca4 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Thu, 6 Oct 2022 23:04:08 -0500 Subject: [PATCH 2/5] some cpolz refactoring expanded unit tests to include non-zero complex coefficients. added newton polish solver for complex coefficients --- src/polyroots_module.F90 | 325 ++++++++++++++++++++++++++------------- test/example.f90 | 2 +- test/polyroots_test.f90 | 141 +++++++++++++---- 3 files changed, 332 insertions(+), 136 deletions(-) diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index 36bfae6..fc0a862 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -67,6 +67,11 @@ module polyroots_module public :: newton_root_polish public :: sort_roots + interface newton_root_polish + module procedure :: newton_root_polish_real, & + newton_root_polish_complex + end interface + contains !***************************************************************************************** @@ -2911,7 +2916,7 @@ end function pythag !### History ! * Jacob Williams, 9/15/2023, created this routine. -subroutine newton_root_polish(n, p, zr, zi, ftol, ztol, maxiter, istat) +subroutine newton_root_polish_real(n, p, zr, zi, ftol, ztol, maxiter, istat) implicit none @@ -3016,7 +3021,123 @@ subroutine func(x, f, dfdx) end subroutine func -end subroutine newton_root_polish +end subroutine newton_root_polish_real +!***************************************************************************************** + +!***************************************************************************************** +!> +! "Polish" a root using a complex Newton Raphson method. +! This routine will perform a Newton iteration and update the roots only if they improve, +! otherwise, they are left as is. +! +!@note Same as [[newton_root_polish_real]] except that `p` is `complex(wp)` + +subroutine newton_root_polish_complex(n, p, zr, zi, ftol, ztol, maxiter, istat) + + implicit none + + integer, intent(in) :: n !! degree of polynomial + complex(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 !! complex number for `(zr,zi)` + complex(wp) :: f !! function evaluation + complex(wp) :: z_prev !! previous value of `z` + complex(wp) :: z_best !! best `z` so far + complex(wp) :: f_best !! best `f` so far + complex(wp) :: dfdx !! derivative evaluation + integer :: i !! counter + + real(wp), parameter :: alpha = 1.0_wp !! newton step size + + ! first evaluate initial point: + z = cmplx(zr, zi, wp) + call func(z, f, dfdx) + + ! initialize: + istat = 0 + z_prev = z + z_best = z + f_best = f + + main: do i = 1, maxiter + + if (i > 1) call func(z, f, dfdx) + if (abs(f) < abs(f_best)) then + ! best so far + zr = real(z, wp) + zi = aimag(z) + z_best = z + f_best = f + end if + if (abs(f) <= ftol) exit main + + if (i == maxiter) then ! max iterations reached + istat = -2 + exit main + end if + + if (dfdx == 0.0_wp) then ! can't proceed + istat = -1 + exit main + end if + + ! Newton correction for next step: + z = z - alpha*(f/dfdx) + + if (abs(z - z_prev) <= ztol) then + ! convergence in x. try this point and see if there is any improvement + istat = 1 + call func(z, f, dfdx) + if (abs(f) < abs(f_best)) then ! best so far + zr = real(z, wp) + zi = aimag(z) + end if + exit main + end if + z_prev = z + + end do main + +contains + + subroutine func(x, f, dfdx) + + !! evaluate the polynomial: + !! `f = p(1)*x**n + p(2)*x**(n-1) + ... + p(n)*x + p(n+1)` + !! and its derivative using Horner's Rule. + !! + !! See: "Roundoff in polynomial evaluation", W. Kahan, 1986 + !! https://rosettacode.org/wiki/Horner%27s_rule_for_polynomial_evaluation#Fortran + + implicit none + + complex(wp), intent(in) :: x + complex(wp), intent(out) :: f !! function value at `x` + complex(wp), intent(out) :: dfdx !! function derivative at `x` + + integer :: i !! counter + + f = p(1) + dfdx = 0.0_wp + do i = 2, n + 1 + dfdx = dfdx*x + f + f = f*x + p(i) + end do + + end subroutine func + +end subroutine newton_root_polish_complex !***************************************************************************************** !***************************************************************************************** @@ -4523,7 +4644,7 @@ subroutine newton(n, poly, apoly, apolyr, z, small, radius, corr, again) az = abs(z) ! if |z|<=1 then apply ruffini-horner's rule for p(z)/p'(z) ! and for the computation of the inclusion radius - if (az <= 1) then + if (az <= 1.0_wp) then p = poly(n + 1) ap = apoly(n + 1) p1 = p @@ -4539,7 +4660,6 @@ subroutine newton(n, poly, apoly, apolyr, z, small, radius, corr, again) ap = ap again = (absp > (small + ap)) if (.not. again) radius = n*(absp + ap)/abs(p1) - return else ! if |z|>1 then apply ruffini-horner's rule to the reversed polynomial ! and use formula (2) for p(z)/p'(z). analogously do for the inclusion @@ -6137,80 +6257,74 @@ subroutine cpolz(a,ndeg,z,ierr) if (ierr /= 0) write(*,*) 'Convergence failure in cpolz' end subroutine cpolz +!***************************************************************************************** -! Copyright (c) 1996 California Institute of Technology, Pasadena, CA. -! ALL RIGHTS RESERVED. -! Based on Government Sponsored Research NAS7-03001. -!>> 2001-01-24 SCOMQR Krogh CSQRT -> CSQRTX to avoid C lib. conflicts. -!>> 1996-04-27 SCOMQR Krogh Changes to use .C. and C%%. -!>> 1996-03-30 SCOMQR Krogh Added external statement. -!>> 1996-01-18 SCOMQR Krogh Added M77CON statements for conv. to C. -!>> 1995-01-03 SCOMQR WVS Added EXTERNAL CQUO, CSQRT so VAX won't gripe -!>> 1992-03-13 SCOMQR FTK Removed implicit statements. -!>> 1987-11-12 SCOMQR Lawson Initial code. -! ------------------------------------------------------------------ -! Version of the Eispack subr, COMQR, for use in the JPL MATH77 -! library. C. L. Lawson, JPL, 1987 Feb 17. -! ------------------------------------------------------------------ -! -! THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE -! ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN -! AND WILKINSON. -! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). -! THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS -! (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. -! -! THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX -! UPPER HESSENBERG MATRIX BY THE QR METHOD. -! -! ON INPUT- -! -! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL -! ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM -! DIMENSION STATEMENT, -! -! N IS THE ORDER OF THE MATRIX, -! -! LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING -! SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, -! SET LOW=1, IGH=N, -! -! HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, -! RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. -! THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN -! INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN -! THE REDUCTION BY CORTH, IF PERFORMED. -! -! ON OUTPUT- -! -! THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN -! DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE -! CALLING COMQR IF SUBSEQUENT CALCULATION OF -! EIGENVECTORS IS TO BE PERFORMED, +!***************************************************************************************** +!> +! This subroutine finds the eigenvalues of a complex +! upper hessenberg matrix by the qr method. ! -! WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, -! RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR -! EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT -! FOR INDICES IERR+1,...,N, +! This subroutine is a translation of a unitary analogue of the +! algol procedure comlr, num. math. 12, 369-376(1968) by martin +! and wilkinson. +! handbook for auto. comp., vol.ii-linear algebra, 396-403(1971). +! the unitary analogue substitutes the qr algorithm of francis +! (comp. jour. 4, 332-345(1962)) for the lr algorithm. ! -! IERR IS SET TO -! ZERO FOR NORMAL RETURN, -! J IF THE J-TH EIGENVALUE HAS NOT BEEN -! DETERMINED AFTER 30 ITERATIONS. +!### Reference +! * Original code from [JPL MATH77 Library](https://netlib.org/math/) ! -! ARITHMETIC IS REAL EXCEPT FOR THE REPLACEMENT OF THE ALGOL -! PROCEDURE CDIV BY COMPLEX DIVISION AND USE OF THE SUBROUTINES -! SQRT AND CMPLX IN COMPUTING COMPLEX SQUARE ROOTS. +!### License +! Copyright (c) 1996 California Institute of Technology, Pasadena, CA. +! ALL RIGHTS RESERVED. +! Based on Government Sponsored Research NAS7-03001. ! -! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, -! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY +!### History +! * 1987-11-12 SCOMQR Lawson Initial code. +! * 1992-03-13 SCOMQR FTK Removed implicit statements. +! * 1995-01-03 SCOMQR WVS Added EXTERNAL CQUO, CSQRT so VAX won't gripe +! * 1996-01-18 SCOMQR Krogh Added M77CON statements for conv. to C. +! * 1996-03-30 SCOMQR Krogh Added external statement. +! * 1996-04-27 SCOMQR Krogh Changes to use .C. and C%%. +! * 2001-01-24 SCOMQR Krogh CSQRT -> CSQRTX to avoid C lib. conflicts. +! * 2022-10-06, Jacob Williams modernized this routine subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) - integer :: en,enm1,i,ierr,igh,its,j,l,ll,low,lp1,n,nm - real(wp) :: hi(nm,n),hr(nm,n) + integer,intent(in) :: nm !! the row dimension of two-dimensional array + !! parameters as declared in the calling program + !! dimension statement + integer,intent(in) :: n !! the order of the matrix + integer,intent(in) :: low !! low and igh are integers determined by the balancing + !! subroutine cbal. if cbal has not been used, + !! set low=1, igh=n + integer,intent(in) :: igh !! low and igh are integers determined by the balancing + !! subroutine cbal. if cbal has not been used, + !! set low=1, igh=n + real(wp),intent(inout) :: hi(nm,n) !! Input: hr and hi contain the real and imaginary parts, + !! respectively, of the complex upper hessenberg matrix. + !! their lower triangles below the subdiagonal contain + !! information about the unitary transformations used in + !! the reduction by corth, if performed. + !! + !! Output: the upper hessenberg portions of hr and hi have been + !! destroyed. therefore, they must be saved before + !! calling comqr if subsequent calculation of + !! eigenvectors is to be performed, + real(wp),intent(inout) :: hr(nm,n) !! see `hi` description + complex(wp),intent(out) :: z(n) !! the real and imaginary parts, + !! respectively, of the eigenvalues. if an error + !! exit is made, the eigenvalues should be correct + !! for indices ierr+1,...,n, + integer,intent(out) :: ierr !! is set to: + !! + !! * zero -- for normal return, + !! * j -- if the j-th eigenvalue has not been + !! determined after 30 iterations. + + integer :: en,enm1,i,its,j,l,ll,lp1 real(wp) :: norm,si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr - complex(wp) :: z(n), z3 + complex(wp) :: z3 ierr = 0 if (low /= igh) then @@ -6239,6 +6353,7 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) end do end do end if + ! ********** store roots isolated by cbal ********** do i = 1, n if (i >= low .and. i <= igh) cycle @@ -6249,11 +6364,12 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) tr = 0.0_wp ti = 0.0_wp - do + main : do ! ********** search for next eigenvalue ********** if (en < low) return its = 0 enm1 = en - 1 + ! ********** look for single small sub-diagonal element ! for l=en step -1 until low -- ********** 240 do ll = low, en @@ -6264,33 +6380,39 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) + abs(hr(l,l)) +abs(hi(l,l)))) exit end do ! ********** form shift ********** - if (l == en) go to 660 - if (its == 30) go to 1000 - if (its == 10 .or. its == 20) go to 320 - sr = hr(en,en) - si = hi(en,en) - xr = hr(enm1,en) * hr(en,enm1) - xi = hi(enm1,en) * hr(en,enm1) - if (xr == 0.0_wp .and. xi == 0.0_wp) go to 340 - yr = (hr(enm1,enm1) - sr) / 2.0_wp - yi = (hi(enm1,enm1) - si) / 2.0_wp - z3 = sqrt(cmplx(yr**2-yi**2+xr,2.0_wp*yr*yi+xi,wp)) - zzr = real(z3,wp) - zzi = aimag(z3) - if (yr * zzr + yi * zzi < 0.0_wp) then - zzr = -zzr - zzi = -zzi + if (l == en) then + ! ********** a root found ********** + z(en) = cmplx(hr(en,en)+tr,hi(en,en)+ti,wp) + en = enm1 + cycle main + end if + if (its == 30) exit main + if (its == 10 .or. its == 20) then + ! ********** form exceptional shift ********** + sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2)) + si = 0.0_wp + else + sr = hr(en,en) + si = hi(en,en) + xr = hr(enm1,en) * hr(en,enm1) + xi = hi(enm1,en) * hr(en,enm1) + if (xr /= 0.0_wp .or. xi /= 0.0_wp) then + yr = (hr(enm1,enm1) - sr) / 2.0_wp + yi = (hi(enm1,enm1) - si) / 2.0_wp + z3 = sqrt(cmplx(yr**2-yi**2+xr,2.0_wp*yr*yi+xi,wp)) + zzr = real(z3,wp) + zzi = aimag(z3) + if (yr * zzr + yi * zzi < 0.0_wp) then + zzr = -zzr + zzi = -zzi + end if + z3 = cmplx(xr,xi,wp) / cmplx(yr+zzr,yi+zzi,wp) + sr = sr - real(z3,wp) + si = si - aimag(z3) + end if end if - z3 = cmplx(xr,xi,wp) / cmplx(yr+zzr,yi+zzi,wp) - sr = sr - real(z3,wp) - si = si - aimag(z3) - go to 340 - - ! ********** form exceptional shift ********** - 320 sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2)) - si = 0.0_wp - 340 do i = low, en + do i = low, en hr(i,i) = hr(i,i) - sr hi(i,i) = hi(i,i) - si end do @@ -6358,17 +6480,14 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) hr(i,en) = sr * yr - si * yi hi(i,en) = sr * yi + si * yr end do + go to 240 - ! ********** a root found ********** - 660 continue - z(en) = cmplx(hr(en,en)+tr,hi(en,en)+ti,wp) - en = enm1 - end do + end do main - ! ********** set error -- no convergence to an - ! eigenvalue after 30 iterations ********** - 1000 ierr = en + ! set error -- no convergence to an + ! eigenvalue after 30 iterations + ierr = en end subroutine scomqr diff --git a/test/example.f90 b/test/example.f90 index a6a5807..7b1a0ae 100644 --- a/test/example.f90 +++ b/test/example.f90 @@ -19,7 +19,7 @@ program example call polyroots(degree, p, zr, zi, istatus) -write(*,'(A,1x,I3)') 'istatus: ', istatus +write(*,'(/A,1x,I3)') 'istatus: ', istatus write(*, '(*(a22,1x))') 'real part', 'imaginary part' do i = 1, degree write(*,'(*(e22.15,1x))') zr(i), zi(i) diff --git a/test/polyroots_test.f90 b/test/polyroots_test.f90 index 29f685f..f6a8024 100644 --- a/test/polyroots_test.f90 +++ b/test/polyroots_test.f90 @@ -11,9 +11,9 @@ program polyroots_test integer,parameter :: n_cases = 14 + 110 !! number of cases to run - real(wp),dimension(:),allocatable :: p, zr, zi, s, q, radius,rr,rc, berr,cond + real(wp),dimension(:),allocatable :: p, pi, zr, zi, s, q, radius, berr,cond integer,dimension(:),allocatable :: conv - complex(wp),dimension(:),allocatable :: r, cp + complex(wp),dimension(:),allocatable :: r, cp, cp_ integer :: degree, i, istatus, icase, n integer,dimension(:),allocatable :: seed real(wp) :: detil @@ -51,46 +51,60 @@ program polyroots_test 12753576._wp, & -10628640._wp, & 3628800._wp ] + pi = 0.0_wp case(2) call allocate_arrays(4) p = [1,-3,20,44,54] + pi = 0.0_wp case(3) call allocate_arrays(6) p = [1,-2,2,1,6,-6,8] + pi = 0.0_wp case(4) call allocate_arrays(5) p = [1,1,-8,-16,7,15] + pi = 0.0_wp case(5) call allocate_arrays(5) p = [1,7,5,6,3,2] + pi = 0.0_wp case(6) call allocate_arrays(5) p = [2,3,6,5,7,1] + pi = 0.0_wp case(7) call allocate_arrays(6) p = [1,0,-14,0,49,0,-36] + pi = 0.0_wp case(8) call allocate_arrays(8) p = [1,0,-30,0,273,0,-820,0,576] + pi = 0.0_wp case(9) call allocate_arrays(4) p = [1,0,0,0,-16] + pi = 0.0_wp case(10) call allocate_arrays(6) p = [1,-2,2,1,6,-6,8] + pi = 0.0_wp case(11) ! a case where 1 is an obvious root call allocate_arrays(5) + pi = 0.0_wp p = [8,-8,16,-16,8,-8] case(12) call allocate_arrays(3) p = [ -8.0e18_wp,3.0e12_wp,5.0e6_wp,1.0_wp] + pi = 0.0_wp case(13) call allocate_arrays(3) p = [4.0_wp, 3.0_wp, 2.0_wp, 1.0_wp] + pi = 0.0_wp case(14) call allocate_arrays(2) p = [3.0_wp, 2.0_wp, 1.0_wp] + pi = 0.0_wp case default ! test a set of random coefficients for each degree: @@ -103,17 +117,18 @@ program polyroots_test do i = 1, degree+1 p(i) = get_random_number(-1000.0_wp,1000.0_wp) + pi(i) = get_random_number(-10000.0_wp,10000.0_wp) end do end select + do i = 1, degree+1 + cp(i) = cmplx(p(i), pi(i), wp) ! put in a complex number + end do + q = reverse(p) ! + cp_ = reversez(cp) ! for the ones that require reverse order 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 - do i = 1, degree+1 - cp(i) = cmplx(p(i), 0.0_wp, wp) ! put in a complex number - end do - if (degree==2) then ! also test this one (only for quadratic equations): call dqdcrt(q, zr, zi) @@ -134,7 +149,7 @@ program polyroots_test call check_results('polyroots', istatus, zr, zi, degree) call cpolyroots(degree, cp, r, istatus) - call check_results('cpolyroots', istatus, real(r, wp), aimag(r), degree) + call check_results_complex('cpolyroots [complex coefficients]', istatus, real(r, wp), aimag(r), degree) end if call rpoly(p, degree, zr, zi, istatus) @@ -153,37 +168,30 @@ program polyroots_test ! for now, just test the following with the real coefficients only: call cpolz(cp,degree,r,istatus) - call check_results('cpolz', istatus, real(r,wp), aimag(r), degree) + call check_results_complex('cpolz [complex coefficients]', istatus, real(r,wp), aimag(r), degree) - q = 0.0_wp istatus = 0 - call cpoly(p,q,degree,zr,zi,fail) + call cpoly(real(cp, wp),aimag(cp),degree,zr,zi,fail) if (fail) istatus = -1 - call check_results('cpoly', istatus, zr, zi, degree) + call check_results_complex('cpoly [complex coefficients]', istatus, zr, zi, degree) call cpqr79(degree,cp,r,istatus) - call check_results('cpqr79', istatus, real(r,wp), aimag(r), degree) + call check_results_complex('cpqr79 [complex coefficients]', istatus, real(r,wp), aimag(r), degree) call qr_algeq_solver(degree,p,zr,zi,istatus,detil=detil) - call check_results('qr_algeq_solver', istatus, real(r,wp), aimag(r), degree) + call check_results('qr_algeq_solver', istatus, zr,zi, degree) - !..... these accept the complex coefficients in reverse order - cp = reversez(cp) - call cmplx_roots_gen(degree, cp, r) ! no status flag - rr = real(r, wp) - rc = aimag(r) - call check_results('cmplx_roots_gen', 0, rr, rc, degree) + !.... + ! these accept the complex coefficients in reverse order + call cmplx_roots_gen(degree, cp_, r) ! no status flag + call check_results_complex('cmplx_roots_gen [complex coefficients]', 0, real(r, wp), aimag(r), degree) - istatus = 0 - call polzeros(degree, cp, 100, r, radius, err) - if (any(err)) istatus = -1 - rr = real(r, wp) - rc = aimag(r) - call check_results('polzeros', istatus, rr, rc, degree) - - call fpml(cp, degree, r, berr, cond, conv, itmax=100) - call check_results('fpml', 0, real(r, wp), aimag(r), degree) + call polzeros(degree, cp_, 100, r, radius, err); istatus = 0; if (any(err)) istatus = -1 + call check_results_complex('polzeros [complex coefficients]', istatus, real(r, wp), aimag(r), degree) + call fpml(cp_, degree, r, berr, cond, conv, itmax=100) + call check_results_complex('fpml [complex coefficients]', 0, real(r, wp), aimag(r), degree) + !.... end do if (failure) error stop 'At least one test failed' @@ -244,8 +252,10 @@ subroutine allocate_arrays(d) degree = d p = [(0, i=1,degree+1)] + pi = [(0, i=1,degree+1)] q = [(0, i=1,degree+1)] cp = [(0, i=1,degree+1)] + cp_ = [(0, i=1,degree+1)] berr = [(0, i=1,degree+1)] zr = [(0, i=1,degree)] @@ -254,8 +264,6 @@ subroutine allocate_arrays(d) r = [(0, i=1,degree)] radius = [(0, i=1,degree)] err = [(.false., i=1,degree)] - rr = [(0, i=1,degree)] - rc = [(0, i=1,degree)] cond = [(0, i=1,degree)] conv = [(0, i=1,degree)] @@ -304,7 +312,7 @@ subroutine check_results(name, istatus, zr, zi, degree) z = cmplx(re(j), im(j), wp) root = p(1) do i = 2, degree+1 - root = root * z + p(i) ! horner's rule + root = root * z + p(i) ! horner's rule !>..need to test ones with complex coefficients end do write(*, '(3(2g23.15,1x))') re(j), im(j), abs(root) if (polish .and. abs(root) > ftol) then @@ -331,6 +339,75 @@ subroutine check_results(name, istatus, zr, zi, degree) end subroutine check_results !***************************************************************************************** + !***************************************************************************************** + subroutine check_results_complex(name, istatus, zr, zi, degree) + + !! check the results (for complex coefficients). + !! if any are not within the tolerance, + !! then also try to polish them using the newton method. + + character(len=*),intent(in) :: name !! name of method + integer,intent(in) :: istatus !! status flag (0 = success) + real(wp),dimension(:),intent(in) :: zr, zi + integer,intent(in) :: degree + + real(wp) :: zr_, zi_ ! copy of inputs for polishing + real(wp),dimension(size(zr)) :: re, im ! copy of inputs for sorting + complex(wp) :: z, root + integer :: i,j !! counter + integer :: istat + + real(wp),parameter :: tol = 1.0e-2_wp !! acceptable root tolerance for tests + real(wp),parameter :: ftol = 1.0e-8_wp !! desired root tolerance + real(wp),parameter :: ztol = 10*epsilon(1.0_wp) !! newton tol for x + logical,parameter :: polish = .true. + + write(*, '(/A,1x,i3)') trim(name) + + if (istatus /= 0) then + failure = .true. + write(*,'(A,1x,i3)') 'Error: method did not converge. istatus = ', istatus + return + end if + + ! sort them in increasing order: + re = zr + im = zi + call sort_roots(re, im) + + write(*, '(a)') ' real part imaginary part root' + + do j = 1, degree + z = cmplx(re(j), im(j), wp) + root = cp(1) + do i = 2, degree+1 + root = root * z + cp(i) ! horner's rule !>..need to test ones with complex coefficients + end do + write(*, '(3(2g23.15,1x))') re(j), im(j), abs(root) + if (polish .and. abs(root) > ftol) then + ! attempt to polish the root: + zr_ = re(j) + zi_ = im(j) + call newton_root_polish(degree, cp, zr_, zi_, & + ftol=ftol, ztol=ztol, maxiter=10, & + istat=istat) + z = cmplx(zr_, zi_, wp) ! recompute root with possibly updated values + root = cp(1) + do i = 2, degree+1 + root = root * z + cp(i) ! horner's rule + end do + write(*, '(3(2g23.15,1x),1X,A)') zr_, zi_, abs(root), 'POLISHED' + if (abs(root) > tol) then + failure = .true. + write(*,'(A)') 'Error: insufficient accuracy *******' + !error stop 'Error: insufficient accuracy' + end if + end if + end do + + end subroutine check_results_complex + !***************************************************************************************** + !***************************************************************************************** !> author: Jacob Williams ! From 218cd2c13a087a13aea802506bebe3a63b63ebd3 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Thu, 6 Oct 2022 23:08:36 -0500 Subject: [PATCH 3/5] cpoly refactoring --- src/polyroots_module.F90 | 331 ++++++++++++++++++++------------------- 1 file changed, 166 insertions(+), 165 deletions(-) diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index fc0a862..fd2c32e 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -6172,84 +6172,84 @@ subroutine cpolz(a,ndeg,z,ierr) end do end do - ! ***************** BALANCE THE MATRIX *********************** + ! ***************** balance the matrix *********************** ! - ! This is an adaption of the EISPACK subroutine BALANC to - ! the special case of a complex companion matrix. 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). + ! this is an adaption of the eispack subroutine balanc to + ! the special case of a complex companion matrix. 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). - ! ********** ITERATIVE LOOP FOR NORM REDUCTION ********** + ! ********** iterative loop for norm reduction ********** do - MORE = .FALSE. - DO I = 1, N - ! Compute R = sum of magnitudes in row I skipping diagonal. - ! C = sum of magnitudes in col I skipping diagonal. - IF (I == 1) THEN - R = ABS(H(1,2,1)) + ABS(H(1,2,2)) - DO J = 3,N - R = R + ABS(H(1,J,1)) + ABS(H(1,J,2)) + more = .false. + do i = 1, n + ! compute r = sum of magnitudes in row i skipping diagonal. + ! c = sum of magnitudes in col i skipping diagonal. + if (i == 1) then + r = abs(h(1,2,1)) + abs(h(1,2,2)) + do j = 3,n + r = r + abs(h(1,j,1)) + abs(h(1,j,2)) end do - C = ABS(H(2,1,1)) + ABS(H(2,1,2)) - ELSE - R = ABS(H(I,I-1,1)) + ABS(H(I,I-1,2)) - C = ABS(H(1,I,1)) + ABS(H(1,I,2)) - IF (I /= N) THEN - C = C + ABS(H(I+1,I,1)) + ABS(H(I+1,I,2)) - END IF - END IF + c = abs(h(2,1,1)) + abs(h(2,1,2)) + else + r = abs(h(i,i-1,1)) + abs(h(i,i-1,2)) + c = abs(h(1,i,1)) + abs(h(1,i,2)) + if (i /= n) then + c = c + abs(h(i+1,i,1)) + abs(h(i+1,i,2)) + end if + end if - ! Determine column scale factor, F. + ! determine column scale factor, f. - G = R / BASE - F = ONE - S = C + R + g = r / base + f = one + s = c + r do - IF (C >= G) exit - F = F * BASE - C = C * B2 + if (c >= g) exit + f = f * base + c = c * b2 end do - G = R * BASE + g = r * base do - IF (C < G) exit - F = F / BASE - C = C / B2 + if (c < g) exit + f = f / base + c = c / b2 end do - ! Will the factor F have a significant effect ? + ! will the factor f have a significant effect ? - IF ((C + R) / F < C95 * S) THEN + if ((c + r) / f < c95 * s) then - ! Yes, so do the scaling. + ! yes, so do the scaling. - G = ONE / F - MORE = .TRUE. + g = one / f + more = .true. - ! Scale Row I + ! scale row i - IF (I == 1) THEN - DO J = 1,N - H(1,J,1) = H(1,J,1)*G - H(1,J,2) = H(1,J,2)*G + if (i == 1) then + do j = 1,n + h(1,j,1) = h(1,j,1)*g + h(1,j,2) = h(1,j,2)*g end do - ELSE - H(I,I-1,1) = H(I,I-1,1)*G - H(I,I-1,2) = H(I,I-1,2)*G - END IF + else + h(i,i-1,1) = h(i,i-1,1)*g + h(i,i-1,2) = h(i,i-1,2)*g + end if - ! Scale Column I + ! scale column i - H(1,I,1) = H(1,I,1) * F - H(1,I,2) = H(1,I,2) * F - IF (I /= N) THEN - H(I+1,I,1) = H(I+1,I,1) * F - H(I+1,I,2) = H(I+1,I,2) * F - END IF + h(1,i,1) = h(1,i,1) * f + h(1,i,2) = h(1,i,2) * f + if (i /= n) then + h(i+1,i,1) = h(i+1,i,1) * f + h(i+1,i,2) = h(i+1,i,2) * f + end if - END IF + end if end do - IF (.not. MORE) exit + if (.not. more) exit end do call scomqr(ndeg,n,1,n,h(1,1,1),h(1,1,2),z,ierr) @@ -6328,7 +6328,7 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) ierr = 0 if (low /= igh) then - ! ********** create real subdiagonal elements ********** + ! create real subdiagonal elements l = low + 1 do i = l, igh @@ -6354,7 +6354,7 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) end do end if - ! ********** store roots isolated by cbal ********** + ! store roots isolated by cbal do i = 1, n if (i >= low .and. i <= igh) cycle z(i) = cmplx(hr(i,i),hi(i,i),wp) @@ -6365,124 +6365,125 @@ subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr) ti = 0.0_wp main : do - ! ********** search for next eigenvalue ********** + ! search for next eigenvalue if (en < low) return its = 0 enm1 = en - 1 - ! ********** look for single small sub-diagonal element - ! for l=en step -1 until low -- ********** - 240 do ll = low, en - l = en + low - ll - if (l == low) exit - if (abs(hr(l,l-1)) <= & - eps * (abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) & - + abs(hr(l,l)) +abs(hi(l,l)))) exit - end do - ! ********** form shift ********** - if (l == en) then - ! ********** a root found ********** - z(en) = cmplx(hr(en,en)+tr,hi(en,en)+ti,wp) - en = enm1 - cycle main - end if - if (its == 30) exit main - if (its == 10 .or. its == 20) then - ! ********** form exceptional shift ********** - sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2)) - si = 0.0_wp - else - sr = hr(en,en) - si = hi(en,en) - xr = hr(enm1,en) * hr(en,enm1) - xi = hi(enm1,en) * hr(en,enm1) - if (xr /= 0.0_wp .or. xi /= 0.0_wp) then - yr = (hr(enm1,enm1) - sr) / 2.0_wp - yi = (hi(enm1,enm1) - si) / 2.0_wp - z3 = sqrt(cmplx(yr**2-yi**2+xr,2.0_wp*yr*yi+xi,wp)) - zzr = real(z3,wp) - zzi = aimag(z3) - if (yr * zzr + yi * zzi < 0.0_wp) then - zzr = -zzr - zzi = -zzi + do + ! look for single small sub-diagonal element + ! for l=en step -1 until low + do ll = low, en + l = en + low - ll + if (l == low) exit + if (abs(hr(l,l-1)) <= & + eps * (abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) & + + abs(hr(l,l)) +abs(hi(l,l)))) exit + end do + ! form shift + if (l == en) then + ! a root found + z(en) = cmplx(hr(en,en)+tr,hi(en,en)+ti,wp) + en = enm1 + cycle main + end if + if (its == 30) exit main + if (its == 10 .or. its == 20) then + ! form exceptional shift + sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2)) + si = 0.0_wp + else + sr = hr(en,en) + si = hi(en,en) + xr = hr(enm1,en) * hr(en,enm1) + xi = hi(enm1,en) * hr(en,enm1) + if (xr /= 0.0_wp .or. xi /= 0.0_wp) then + yr = (hr(enm1,enm1) - sr) / 2.0_wp + yi = (hi(enm1,enm1) - si) / 2.0_wp + z3 = sqrt(cmplx(yr**2-yi**2+xr,2.0_wp*yr*yi+xi,wp)) + zzr = real(z3,wp) + zzi = aimag(z3) + if (yr * zzr + yi * zzi < 0.0_wp) then + zzr = -zzr + zzi = -zzi + end if + z3 = cmplx(xr,xi,wp) / cmplx(yr+zzr,yi+zzi,wp) + sr = sr - real(z3,wp) + si = si - aimag(z3) end if - z3 = cmplx(xr,xi,wp) / cmplx(yr+zzr,yi+zzi,wp) - sr = sr - real(z3,wp) - si = si - aimag(z3) - end if - end if + end if - do i = low, en - hr(i,i) = hr(i,i) - sr - hi(i,i) = hi(i,i) - si - end do + do i = low, en + hr(i,i) = hr(i,i) - sr + hi(i,i) = hi(i,i) - si + end do - tr = tr + sr - ti = ti + si - its = its + 1 - ! ********** reduce to triangle (rows) ********** - lp1 = l + 1 - - do i = lp1, en - sr = hr(i,i-1) - hr(i,i-1) = 0.0_wp - norm = sqrt(hr(i-1,i-1)*hr(i-1,i-1)+hi(i-1,i-1)*hi(i-1,i-1)+sr*sr) - xr = hr(i-1,i-1) / norm - xi = hi(i-1,i-1) / norm - z(i-1) = cmplx(xr,xi,wp) - hr(i-1,i-1) = norm - hi(i-1,i-1) = 0.0_wp - hi(i,i-1) = sr / norm - do j = i, en - yr = hr(i-1,j) - yi = hi(i-1,j) - zzr = hr(i,j) - zzi = hi(i,j) - hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr - hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi - hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr - hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi - end do - end do + tr = tr + sr + ti = ti + si + its = its + 1 + ! reduce to triangle (rows) + lp1 = l + 1 + + do i = lp1, en + sr = hr(i,i-1) + hr(i,i-1) = 0.0_wp + norm = sqrt(hr(i-1,i-1)*hr(i-1,i-1)+hi(i-1,i-1)*hi(i-1,i-1)+sr*sr) + xr = hr(i-1,i-1) / norm + xi = hi(i-1,i-1) / norm + z(i-1) = cmplx(xr,xi,wp) + hr(i-1,i-1) = norm + hi(i-1,i-1) = 0.0_wp + hi(i,i-1) = sr / norm + do j = i, en + yr = hr(i-1,j) + yi = hi(i-1,j) + zzr = hr(i,j) + zzi = hi(i,j) + hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr + hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi + hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr + hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi + end do + end do - si = hi(en,en) - if (si /= 0.0_wp) then - norm = abs(cmplx(hr(en,en),si,wp)) - sr = hr(en,en) / norm - si = si / norm - hr(en,en) = norm - hi(en,en) = 0.0_wp - end if - ! ********** inverse operation (columns) ********** - do j = lp1, en - xr = real(z(j-1),wp) - xi = aimag(z(j-1)) - do i = l, j - yr = hr(i,j-1) - yi = 0.0 - zzr = hr(i,j) - zzi = hi(i,j) - if (i /= j) then - yi = hi(i,j-1) - hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi - end if - hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr - hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr - hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi - end do - end do + si = hi(en,en) + if (si /= 0.0_wp) then + norm = abs(cmplx(hr(en,en),si,wp)) + sr = hr(en,en) / norm + si = si / norm + hr(en,en) = norm + hi(en,en) = 0.0_wp + end if + ! inverse operation (columns) + do j = lp1, en + xr = real(z(j-1),wp) + xi = aimag(z(j-1)) + do i = l, j + yr = hr(i,j-1) + yi = 0.0 + zzr = hr(i,j) + zzi = hi(i,j) + if (i /= j) then + yi = hi(i,j-1) + hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi + end if + hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi + end do + end do - if (si == 0.0) go to 240 + if (si /= 0.0_wp) then + do i = l, en + yr = hr(i,en) + yi = hi(i,en) + hr(i,en) = sr * yr - si * yi + hi(i,en) = sr * yi + si * yr + end do + end if - do i = l, en - yr = hr(i,en) - yi = hi(i,en) - hr(i,en) = sr * yr - si * yi - hi(i,en) = sr * yi + si * yr end do - go to 240 - end do main ! set error -- no convergence to an From f84c9dabe6e4adeaee248b77cd76e495bc2d9338 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Thu, 6 Oct 2022 23:09:54 -0500 Subject: [PATCH 4/5] commenting --- test/polyroots_test.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/polyroots_test.f90 b/test/polyroots_test.f90 index f6a8024..b844515 100644 --- a/test/polyroots_test.f90 +++ b/test/polyroots_test.f90 @@ -312,7 +312,7 @@ subroutine check_results(name, istatus, zr, zi, degree) z = cmplx(re(j), im(j), wp) root = p(1) do i = 2, degree+1 - root = root * z + p(i) ! horner's rule !>..need to test ones with complex coefficients + root = root * z + p(i) ! horner's rule end do write(*, '(3(2g23.15,1x))') re(j), im(j), abs(root) if (polish .and. abs(root) > ftol) then @@ -381,7 +381,7 @@ subroutine check_results_complex(name, istatus, zr, zi, degree) z = cmplx(re(j), im(j), wp) root = cp(1) do i = 2, degree+1 - root = root * z + cp(i) ! horner's rule !>..need to test ones with complex coefficients + root = root * z + cp(i) ! horner's rule end do write(*, '(3(2g23.15,1x))') re(j), im(j), abs(root) if (polish .and. abs(root) > ftol) then From 009f9e6bdc12ad120cc8bd429d44f8d15677a8d2 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Fri, 7 Oct 2022 20:55:52 -0500 Subject: [PATCH 5/5] removed nonstandard isnan --- src/polyroots_module.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/polyroots_module.F90 b/src/polyroots_module.F90 index fd2c32e..c4c51d7 100644 --- a/src/polyroots_module.F90 +++ b/src/polyroots_module.F90 @@ -21,6 +21,7 @@ module polyroots_module use iso_fortran_env + use ieee_arithmetic implicit none @@ -5402,7 +5403,8 @@ function check_nan_inf(a) result(res) ! check for nan and inf re_a = real(a,wp) im_a = aimag(a) - res = isnan(re_a) .or. isnan(im_a) .or. (abs(re_a)>big) .or. (abs(im_a)>big) + res = ieee_is_nan(re_a) .or. ieee_is_nan(im_a) .or. & + (abs(re_a)>big) .or. (abs(im_a)>big) end function check_nan_inf