Quad LAPACK

.github/workflows/CI.yml

@@ -58,7 +58,8 @@ jobs:
     #   run: fpm build --profile release
 
     - name: Run tests
-      run: fpm test --profile debug --flag -coverage
+      run: fpm test --profile debug --flag "-coverage -DUSE_STDLIB_LAPACK"
+      # run: fpm test --profile debug --flag "-coverage -lblas -llapack"
 
     - name: Create coverage report
       run: |

.gitignore

@@ -35,4 +35,5 @@
 .DS_Store
 /build
 /doc
-/tmp
+/tmp
+/*.png

README.md

@@ -103,13 +103,34 @@ Preprocessor flag | Kind  | Number of bytes
 `REAL64`  | `real(kind=real64)`  | 8
 `REAL128` | `real(kind=real128)` | 16
 
-For example, to build a single precision version of the library, use:
+<!-- For example, to build a single precision version of the library, use:
 
 ```
 fpm build --profile release --flag "-DREAL32"
 ```
 
-All routines, except for `polyroots` are available for any of the three real kinds. `polyroots` is not available for `real128` kinds since there is no corresponding LAPACK eigenvalue solver.
+All routines, except for `polyroots` are available for any of the three real kinds. `polyroots` is not available for `real128` kinds since there is no corresponding LAPACK eigenvalue solver. -->
+
+<!-- UPDATE: -->
+
+The library requires some LAPACK routines to be present. You can specify them to the linker:
+
+```
+fpm test --profile release --flag "-DREAL64 -lblas -llapack"
+```
+
+Or, on a Mac:
+```
+fpm test --profile release --flag "-DREAL64 -framework Accelerate"
+```
+
+Or, use the FPM package from [here](https://github.com/perazz/fortran-lapack) by defining the `USE_STDLIB_LAPACK` preprocessor flag:
+
+```
+fpm test --profile release --flag "-DREAL64 -DUSE_STDLIB_LAPACK"
+```
+
+Note that using `USE_STDLIB_LAPACK` is the only way to get a quad precision version of `polyroots` and `cpolyroots`, since the default LAPACK does not have a quad precision version.
 
 To run the unit tests:
 

fpm.toml

@@ -10,7 +10,10 @@ keywords = ["root"]
 auto-executables = false
 auto-examples = false
 auto-tests = true
-link = ["lapack", "blas"]
+# link = ["lapack", "blas"]
+
+[dependencies]
+    fortran-lapack = { git="https://github.com/perazz/fortran-lapack.git" }
 
 [dev-dependencies]
     mersenne-twister-fortran = { git="https://github.com/jacobwilliams/mersenne-twister-fortran.git", tag = "1.0.1" }

src/polyroots_lapack_module.F90

@@ -0,0 +1,61 @@
+!*****************************************************************************************
+!>
+!  LAPACK interface module
+!
+!### Author
+!  * Jacob Williams
+
+module polyroots_lapack_module
+
+#ifdef USE_STDLIB_LAPACK
+#ifdef REAL32
+    use stdlib_linalg_lapack, sgeev => stdlib_sgeev, cgeev => stdlib_cgeev
+#elif REAL128
+    use stdlib_linalg_lapack, qgeev => stdlib_qgeev, wgeev => stdlib_wgeev
+#else
+    use stdlib_linalg_lapack, dgeev => stdlib_dgeev, zgeev => stdlib_zgeev
+#endif
+#endif
+
+implicit none
+
+#ifndef USE_STDLIB_LAPACK
+
+#ifdef REAL32
+    interface
+        subroutine sgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
+            implicit none
+            character :: jobvl, jobvr
+            integer :: info, lda, ldvl, ldvr, lwork, n
+            real :: a(lda, *), vl(ldvl, *), vr(ldvr, *), wi(*), work(*), wr(*)
+        end subroutine sgeev
+        subroutine cgeev( jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info )
+            implicit none
+            character :: jobvl, jobvr
+            integer :: info, lda, ldvl, ldvr, lwork, n
+            real :: rwork( * )
+            complex :: a( lda, * ), vl( ldvl, * ), vr( ldvr, * ), w( * ), work( * )
+        end subroutine cgeev
+    end interface
+#elif REAL128
+    ! do not have a quad solver in lapack
+#else
+    interface
+        subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
+            implicit none
+            character :: jobvl, jobvr
+            integer :: info, lda, ldvl, ldvr, lwork, n
+            double precision :: a(lda, *), vl(ldvl, *), vr(ldvr, *), wi(*), work(*), wr(*)
+        end subroutine dgeev
+        subroutine zgeev( jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info )
+            implicit none
+            character :: jobvl, jobvr
+            integer :: info, lda, ldvl, ldvr, lwork, n
+            double precision :: rwork( * )
+            complex*16 :: a( lda, * ), vl( ldvl, * ), vr( ldvr, * ), w( * ), work( * )
+        end subroutine zgeev
+    end interface
+#endif
+#endif
+
+end module polyroots_lapack_module

src/polyroots_module.F90

@@ -2674,6 +2674,7 @@ end subroutine hqr
 !@note Works only for single and double precision.
 
 subroutine polyroots(n, p, wr, wi, info)
+    use polyroots_lapack_module
 
     implicit none
 
@@ -2691,28 +2692,6 @@ subroutine polyroots(n, p, wr, wi, info)
     real(wp), allocatable, dimension(:) :: work !! work array
     real(wp), dimension(1) :: vl, vr !! not used here
 
-#ifdef REAL32
-    interface
-        subroutine sgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
-            implicit none
-            character :: jobvl, jobvr
-            integer :: info, lda, ldvl, ldvr, lwork, n
-            real :: a(lda, *), vl(ldvl, *), vr(ldvr, *), wi(*), work(*), wr(*)
-        end subroutine sgeev
-    end interface
-#elif REAL128
-    ! do not have a quad solver in lapack
-#else
-    interface
-        subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
-            implicit none
-            character :: jobvl, jobvr
-            integer :: info, lda, ldvl, ldvr, lwork, n
-            double precision :: a(lda, *), vl(ldvl, *), vr(ldvr, *), wi(*), work(*), wr(*)
-        end subroutine dgeev
-    end interface
-#endif
-
     ! error check:
     if (abs(p(1)) == 0.0_wp) then
         info = -999
@@ -2737,7 +2716,12 @@ end subroutine dgeev
     ! single precision
     call sgeev('N', 'N', n, a, n, wr, wi, vl, 1, vr, 1, work, 3*n, info)
 #elif REAL128
+#ifdef USE_STDLIB_LAPACK
+    ! quad precision
+    call qgeev('N', 'N', n, a, n, wr, wi, vl, 1, vr, 1, work, 3*n, info)
+#else
     error stop 'do not have a quad solver in lapack'
+#endif
 #else
     ! by default, use double precision:
     call dgeev('N', 'N', n, a, n, wr, wi, vl, 1, vr, 1, work, 3*n, info)
@@ -2763,6 +2747,8 @@ end subroutine polyroots
 
 subroutine cpolyroots(n, p, w, info)
 
+    use polyroots_lapack_module
+
     implicit none
 
     integer, intent(in) :: n !! polynomial degree
@@ -2779,30 +2765,6 @@ subroutine cpolyroots(n, p, w, info)
     real(wp), allocatable, dimension(:) :: rwork !! rwork array (2*n)
     complex(wp), dimension(1) :: vl, vr !! not used here
 
-#ifdef REAL32
-    interface
-        subroutine cgeev( jobvl, jobvr, n, a, lda, w,      vl, ldvl, vr, ldvr, work, lwork, rwork, info )
-            implicit none
-            character :: jobvl, jobvr
-            integer :: info, lda, ldvl, ldvr, lwork, n
-            real :: rwork( * )
-            complex :: a( lda, * ), vl( ldvl, * ), vr( ldvr, * ), w( * ), work( * )
-        end subroutine cgeev
-    end interface
-#elif REAL128
-    ! do not have a quad solver in lapack
-#else
-    interface
-        subroutine zgeev( jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info )
-            implicit none
-            character :: jobvl, jobvr
-            integer :: info, lda, ldvl, ldvr, lwork, n
-            double precision :: rwork( * )
-            complex*16 :: a( lda, * ), vl( ldvl, * ), vr( ldvr, * ), w( * ), work( * )
-        end subroutine zgeev
-    end interface
-#endif
-
     ! error check:
     if (abs(p(1)) == 0.0_wp) then
         info = -999
@@ -2828,7 +2790,12 @@ end subroutine zgeev
     ! single precision
     call cgeev('N', 'N', n, a, n, w, vl, 1, vr, 1, work, 3*n, rwork, info)
 #elif REAL128
+#ifdef USE_STDLIB_LAPACK
+    ! quad precision
+    call wgeev('N', 'N', n, a, n, w, vl, 1, vr, 1, work, 3*n, rwork, info)
+#else
     error stop 'do not have a quad solver in lapack'
+#endif
 #else
     ! by default, use double precision:
     call zgeev('N', 'N', n, a, n, w, vl, 1, vr, 1, work, 3*n, rwork, info)
@@ -6117,6 +6084,12 @@ subroutine dpolz(ndeg,a,zr,zi,ierr)
         real(wp),parameter :: machep = eps        !! d1mach(4)
         integer,parameter :: base = radix(1.0_wp) !! i1mach(10)
         integer,parameter :: b2 = base*base
+#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
 
         ierr = 0
 
@@ -6306,7 +6279,7 @@ subroutine dpolz(ndeg,a,zr,zi,ierr)
                          endif
                       endif
                       en = enm2
-                   elseif ( its==30 ) then
+                   elseif ( its==maxiter ) then
                       ! ********** set error -- no convergence to an eigenvalue after 30 iterations **********
                       ierr = en
                       exit main

test/example.f90

@@ -17,15 +17,44 @@ program example
 real(wp),dimension(degree) :: zr !! real components of roots
 real(wp),dimension(degree) :: zi !! imaginary components of roots
 
-if (wp==real128) stop ! don't have a quad solver
+! if we have a quad solver:
+#ifdef USE_STDLIB_LAPACK
 
 call polyroots(degree, p, zr, zi, istatus)
 
 write(*,'(/A,1x,I3)') 'istatus: ', istatus
-write(*, '(*(a22,1x))') 'real part', 'imaginary part'
+write(*, '(*(a22,1x))') 'real part', 'imaginary part', 'root'
 do i = 1, degree
-    write(*,'(*(e22.15,1x))') zr(i), zi(i)
+    write(*,'(*(e22.15,1x))') zr(i), zi(i), abs(evaluate(p, cmplx(zr(i), zi(i), wp)))
 end do
 
+#endif
+
+contains
+
+    function evaluate(p, x) result(f)
+
+        !! 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
+
+        real(wp),dimension(:),intent(in) :: p !! coefficients
+        complex(wp), intent(in) :: x
+        complex(wp) :: f    !! function value at `x`
+
+        integer :: i !! counter
+
+        f = p(1)
+        do i = 2, size(p)-1 + 1
+            f = f*x + p(i)
+        end do
+
+    end function evaluate
+
 end program example
 !*****************************************************************************************

test/polyroots_test.f90

@@ -168,12 +168,15 @@ program polyroots_test
             call check_results('rroots_chebyshev_cubic', 0, zr, zi, degree)
         end if
 
-        if (wp /= real128) then
+        if (wp == real128) then
+#ifdef USE_STDLIB_LAPACK
+            ! if we have a quad solver:
             call polyroots(degree, p, zr, zi, istatus)
             call check_results('polyroots', istatus, zr, zi, degree)
 
             call cpolyroots(degree, cp, r, istatus)
             call check_results_complex('cpolyroots [complex coefficients]', istatus, real(r, wp), aimag(r), degree)
+#endif
         end if
 
         call rpoly(p, degree, zr, zi, istatus)