Merge pull request #15 from aeberspaecher/amos

Add AMOS sources and convenience wrappers

src/CMakeLists.txt

@@ -5,6 +5,7 @@ set(SRC
     ppm.f90
     optimize.f90
     special.f90
+    amos.f90
 )
 
 if(WITH_LAPACK)
@@ -18,4 +19,8 @@ if(WITH_HDF5)
     include_directories($ENV{SPKG_LOCAL}/include)
 endif()
 
+
+add_definitions(-Wno-implicit-interface)
+add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/legacy/amos)
+
 add_library(fortran_utils ${SRC})

src/amos.f90

@@ -0,0 +1,95 @@
+module amos
+  implicit none
+
+  ! This module contains interfaces for the amos routines (in unmodified form
+  ! found on http://netlib.org/amos).
+  ! The slightly modified sources taken from SciPy's "legacy/amos/special"
+  ! module are linked against.
+
+  ! The original code was written by Donald E. Amos at Sandia National Labratories
+  ! See also "Portable package for Bessel functions of a complex argument and
+  ! non-negative order" by D. E. Amos in Trans. Math. Software (1986)
+
+  ! Note that we made all routines pure routines. We do this to be able to
+  ! use these routines in OpenMP loops and forall statements. The F77 routines
+  ! are, of course, not declared as pure. However, they behave as such - the do
+  ! not have side effects. Also, the intent of every argument is clear.
+  ! Technically, there is one exception: the ?1mach.f90 routines (used in
+  ! finding machine constants) each contain a write statement in case the query
+  ! fails. We ignore this possibility (as we trust the queries made by amos).
+  ! The author of the wrapper uses this without any problems in production code.
+
+  interface
+
+     ! Bessel J function:
+     pure subroutine zbesj(re, im, order, scaling, length, zOut_r, zOut_i, &
+          underflows, ierr)
+       implicit none
+       double precision, intent(in) :: re, im, order
+       integer, intent(in) :: scaling, length
+       double precision, intent(out) :: zOut_r(length), zOut_i(length)
+       integer, intent(out) :: underflows, ierr
+     end subroutine zbesj
+
+     ! Bessel Y function of second kind:
+     pure subroutine zbesy(re, im, order, scaling, length, zOut_r, zOut_i, &
+          underflows, workr, worki, ierr)
+       implicit none
+       double precision, intent(in) :: re, im, order
+       integer, intent(in) :: scaling, length
+       double precision, intent(out) :: workr(length), worki(length), &
+            zOut_r(length), zOut_i(length)
+       integer, intent(out) :: underflows, ierr
+     end subroutine zbesy
+
+     ! modified Bessel I function of first kind:
+     pure subroutine zbesi(re, im, order, scaling, length, zOut_r, zOut_i, &
+          underflows, ierr)
+       implicit none
+       double precision, intent(in) :: re, im, order
+       integer, intent(in) :: scaling, length
+       double precision, intent(out) :: zOut_r(length), zOut_i(length)
+       integer, intent(out) :: underflows, ierr
+     end subroutine zbesi
+
+     ! modified Bessel K function of second kind:
+     pure subroutine zbesk(re, im, order, scaling, length, zOut_r, zOut_i, &
+          underflows, ierr)
+       implicit none
+       double precision, intent(in) :: re, im, order
+       integer, intent(in) :: scaling, length
+       double precision, intent(out) :: zOut_r(length), zOut_i(length)
+       integer, intent(out) :: underflows, ierr
+     end subroutine zbesk
+
+     ! Hankel H functions:
+     pure subroutine zbesh(re, im, order, scaling, hankelkind, length, zOut_r, &
+          zOut_i, underflows, ierr)
+       implicit none
+       double precision, intent(in) :: re, im, order
+       integer, intent(in) :: scaling, hankelkind, length
+       double precision, intent(out) :: zOut_r(length), zOut_i(length)
+       integer, intent(out) :: underflows, ierr
+     end subroutine zbesh
+
+     ! Airy function Ai (or its derivative):
+     pure subroutine zairy(re, im, deriv, scaling, zOut_r, zOut_i, underflow, ierr)
+       implicit none
+       double precision, intent(in) :: re, im
+       integer, intent(in) :: deriv, scaling
+       double precision, intent(out) :: zOut_r, zOut_i
+       integer, intent(out) :: underflow, ierr
+     end subroutine zairy
+
+     ! Airy function Bi (or its derivative)
+     pure subroutine zbiry(re, im, deriv, scaling, zOut_r, zOut_i, ierr)
+       implicit none
+       double precision, intent(in) :: re, im
+       integer, intent(in) :: deriv, scaling
+       double precision, intent(out) :: zOut_r, zOut_i
+       integer, intent(out) :: ierr
+     end subroutine zbiry
+
+  end interface
+
+end module amos

src/legacy/amos/CMakeLists.txt

@@ -0,0 +1,49 @@
+
+set(SRC
+    ${SRC}
+    ${CMAKE_CURRENT_SOURCE_DIR}/d1mach.f90
+    ${CMAKE_CURRENT_SOURCE_DIR}/i1mach.f90
+    ${CMAKE_CURRENT_SOURCE_DIR}/r1mach.f90
+    ${CMAKE_CURRENT_SOURCE_DIR}/dgamln.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/dsclmr.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/fdump.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zabs.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zacai.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zacon.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zairy.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zasyi.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbesh.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbesi.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbesj.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbesk.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbesy.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbinu.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbiry.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbknu.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbuni.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zbunk.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zdiv.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zexp.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zkscl.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zlog.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zmlri.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zmlt.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zrati.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zs1s2.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zseri.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zshch.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zsqrt.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zuchk.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zunhj.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zuni1.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zuni2.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zunik.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zunk1.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zunk2.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zuoik.f
+    ${CMAKE_CURRENT_SOURCE_DIR}/zwrsk.f
+    PARENT_SCOPE
+)
+
+add_definitions(-Wno-implicit-interface)
+# TODO: the line above is useless

src/legacy/amos/d1mach.f90

@@ -0,0 +1,77 @@
+!DECK D1MACH
+      DOUBLE PRECISION FUNCTION D1MACH (I)
+      IMPLICIT NONE
+      INTEGER :: I
+      DOUBLE PRECISION :: B, X
+!***BEGIN PROLOGUE  D1MACH
+!***PURPOSE  Return floating point machine dependent constants.
+!***LIBRARY   SLATEC
+!***CATEGORY  R1
+!***TYPE      SINGLE PRECISION (D1MACH-S, D1MACH-D)
+!***KEYWORDS  MACHINE CONSTANTS
+!***AUTHOR  Fox, P. A., (Bell Labs)
+!           Hall, A. D., (Bell Labs)
+!           Schryer, N. L., (Bell Labs)
+!***DESCRIPTION
+!
+!   D1MACH can be used to obtain machine-dependent parameters for the
+!   local machine environment.  It is a function subprogram with one
+!   (input) argument, and can be referenced as follows:
+!
+!        A = D1MACH(I)
+!
+!   where I=1,...,5.  The (output) value of A above is determined by
+!   the (input) value of I.  The results for various values of I are
+!   discussed below.
+!
+!   D1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
+!   D1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
+!   D1MACH(3) = B**(-T), the smallest relative spacing.
+!   D1MACH(4) = B**(1-T), the largest relative spacing.
+!   D1MACH(5) = LOG10(B)
+!
+!   Assume single precision numbers are represented in the T-digit,
+!   base-B form
+!
+!              sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
+!
+!   where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
+!   EMIN .LE. E .LE. EMAX.
+!
+!   The values of B, T, EMIN and EMAX are provided in I1MACH as
+!   follows:
+!   I1MACH(10) = B, the base.
+!   I1MACH(11) = T, the number of base-B digits.
+!   I1MACH(12) = EMIN, the smallest exponent E.
+!   I1MACH(13) = EMAX, the largest exponent E.
+!
+!
+!***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
+!                 a portable library, ACM Transactions on Mathematical
+!                 Software 4, 2 (June 1978), pp. 177-188.
+!***ROUTINES CALLED  XERMSG
+!***REVISION HISTORY  (YYMMDD)
+!   790101  DATE WRITTEN
+!   960329  Modified for Fortran 90 (BE after suggestions by EHG)      
+!***END PROLOGUE  D1MACH
+!      
+      X = 1.0D0
+      B = RADIX(X)
+      SELECT CASE (I)
+        CASE (1)
+          D1MACH = B**(MINEXPONENT(X)-1) ! the smallest positive magnitude.
+        CASE (2)
+          D1MACH = HUGE(X)               ! the largest magnitude.
+        CASE (3)
+          D1MACH = B**(-DIGITS(X))       ! the smallest relative spacing.
+        CASE (4)
+          D1MACH = B**(1-DIGITS(X))      ! the largest relative spacing.
+        CASE (5)
+          D1MACH = LOG10(B)
+        CASE DEFAULT
+          WRITE (*, FMT = 9000)
+ 9000     FORMAT ('1ERROR    1 IN D1MACH - I OUT OF BOUNDS')
+          STOP
+      END SELECT
+      RETURN
+      END