Merge 707dc7a into dd63e02

ludvigak · Jan 23, 2021 · 1bfd173 · 1bfd173
2 parents dd63e02 + 707dc7a
commit 1bfd173
Show file tree

Hide file tree

Showing 2 changed files with 128 additions and 72 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "FINUFFT"
 uuid = "d8beea63-0952-562e-9c6a-8e8ef7364055"
 author = "Ludvig af Klinteberg <ludvigak@kth.se>"
-version = "1.1.0"
+version = "2.0.0"
 
 [deps]
 finufft_jll = "c41cd5a2-72a3-5203-9076-a500b088fc82"

diff --git a/src/FINUFFT.jl b/src/FINUFFT.jl
@@ -17,65 +17,105 @@ export nufft_c_opts # backward-compability
 ## External dependencies
 using finufft_jll
 
+const libfinufft = finufft_jll.libfinufft
+
 const BIGINT = Int64 # defined in src/finufft.h
 
 
 ## FINUFFT opts struct from src/finufft.h
 """
     mutable struct nufft_opts    
-        debug              :: Cint                
-        spread_debug       :: Cint         
-        spread_sort        :: Cint          
-        spread_kerevalmeth :: Cint   
-        spread_kerpad      :: Cint        
-        chkbnds            :: Cint              
-        fftw               :: Cint                 
-        modeord            :: Cint
-        upsampfac          :: Cdouble         
+        modeord             :: Cint
+        chkbnds             :: Cint              
+        # 
+        debug               :: Cint                
+        spread_debug        :: Cint         
+        showwarn            :: Cint
+        # 
+        nthreads            :: Cint
+        fftw                :: Cint                 
+        spread_sort         :: Cint          
+        spread_kerevalmeth  :: Cint   
+        spread_kerpad       :: Cint        
+        upsampfac           :: Cdouble         
+        spread_thread       :: Cint
+        maxbatchsize        :: Cint
+        spread_nthr_atomic  :: Cint
+        spread_max_sp_size  :: Cint
     end
 
 Options struct passed to the FINUFFT library.
 
 # Fields
 
-    debug :: Cint
-0: silent, 1: text basic timing output
-
-    spread_debug :: Cint
-passed to spread_opts, 0 (no text) 1 (some) or 2 (lots)
-
-    spread_sort :: Cint
-passed to spread_opts, 0 (don't sort) 1 (do) or 2 (heuristic)
-
-    spread_kerevalmeth :: Cint
-passed to spread_opts, 0: exp(sqrt()), 1: Horner ppval (faster)
-
-    spread_kerpad :: Cint
-passed to spread_opts, 0: don't pad to mult of 4, 1: do
+modeord :: Cint 
+    (type 1,2 only):    0: CMCL-style increasing mode order
+                        1: FFT-style mode order
+chkbnds :: Cint
+    0: don't check NU pts in [-3pi,3pi)
+    1: do (<few % slower)
+
+# diagnostic opts...
+debug :: Cint
+    0: silent
+    1: some timing/debug
+    2: more
+spread_debug :: Cint
+    0: silent
+    1: some timing/debug
+    2: tonnes
+showwarn :: Cint
+    0: don't print warnings to stderr
+    1: do
+
+# algorithm performance opts...
+nthreads :: Cint
+    number of threads to use, or 0 uses all available
+fftw :: Cint
+    plan flags to FFTW (FFTW_ESTIMATE=64, FFTW_MEASURE=0,...)
+spread_sort :: Cint
+    0: don't sort
+    1: do
+    2: heuristic choice
+spread_kerevalmeth :: Cint
+    0: exp(sqrt()) spreading kernel
+    1: Horner piecewise poly (faster)
+spread_kerpad :: Cint
+    option only for exp(sqrt())
+    0: don't pad kernel to 4n
+    1: do
+upsampfac :: Cdouble
+    upsampling ratio sigma: 2.0 std, 1.25 small FFT, 0.0 auto
+spread_thread :: Cint
+    (vectorized ntr>1 only):    0: auto, 1: seq multithreaded,
+                                2: parallel single-thread spread
+maxbatchsize :: Cint
+    option for vectorized ntr>1 only
+    max transform batch, 0 auto
+spread_nthr_atomic :: Cint
+    if >=0, threads above which spreader OMP critical goes atomic
+spread_max_sp_size :: Cint
+    if >0, overrides spreader (dir=1) max subproblem size
 
-    chkbnds :: Cint
-0: don't check if input NU pts in [-3pi,3pi], 1: do
-
-    fftw :: Cint
-0:`FFTW_ESTIMATE`, or 1:`FFTW_MEASURE` (slow plan but faster)
-
-    modeord :: Cint
-0: CMCL-style increasing mode ordering (neg to pos), or\\
-1: FFT-style mode ordering (affects type-1,2 only)
-
-    upsampfac::Cdouble
-upsampling ratio sigma, either 2.0 (standard) or 1.25 (small FFT)
 """
 mutable struct nufft_opts    
-    debug              :: Cint                
-    spread_debug       :: Cint         
-    spread_sort        :: Cint          
-    spread_kerevalmeth :: Cint   
-    spread_kerpad      :: Cint        
-    chkbnds            :: Cint              
-    fftw               :: Cint                 
-    modeord            :: Cint
-    upsampfac          :: Cdouble         
+    modeord             :: Cint
+    chkbnds             :: Cint              
+    # 
+    debug               :: Cint                
+    spread_debug        :: Cint         
+    showwarn            :: Cint
+    # 
+    nthreads            :: Cint
+    fftw                :: Cint                 
+    spread_sort         :: Cint          
+    spread_kerevalmeth  :: Cint   
+    spread_kerpad       :: Cint        
+    upsampfac           :: Cdouble         
+    spread_thread       :: Cint
+    maxbatchsize        :: Cint
+    spread_nthr_atomic  :: Cint
+    spread_max_sp_size  :: Cint
 end
 
 const nufft_c_opts = nufft_opts # backward compability
@@ -87,25 +127,32 @@ Return a [`nufft_opts`](@ref) struct with the default FINUFFT settings.\\
 See: <https://finufft.readthedocs.io/en/latest/usage.html#options>
 """
 function finufft_default_opts()
-    opts = nufft_opts(0,0,0,0,0,0,0,0,0)
+    opts = nufft_opts(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
     ccall( (:finufft_default_opts, libfinufft),
            Nothing,
            (Ref{nufft_opts},),
            opts
            )
+    # default to number of julia threads
+    opts.nthreads = Threads.nthreads()
     return opts
 end
 
 ### Error handling
-const ERR_EPS_TOO_SMALL        = 1
-const ERR_MAXNALLOC            = 2
-const ERR_SPREAD_BOX_SMALL     = 3
-const ERR_SPREAD_PTS_OUT_RANGE = 4
-const ERR_SPREAD_ALLOC         = 5
-const ERR_SPREAD_DIR           = 6
-const ERR_UPSAMPFAC_TOO_SMALL  = 7
-const HORNER_WRONG_BETA        = 8
-const ERR_NDATA_NOTVALID       = 9
+const WARN_EPS_TOO_SMALL            = 1
+const ERR_MAXNALLOC                 = 2
+const ERR_SPREAD_BOX_SMALL          = 3
+const ERR_SPREAD_PTS_OUT_RANGE      = 4
+const ERR_SPREAD_ALLOC              = 5
+const ERR_SPREAD_DIR                = 6
+const ERR_UPSAMPFAC_TOO_SMALL       = 7
+const HORNER_WRONG_BETA             = 8
+const ERR_NDATA_NOTVALID            = 9
+const ERR_TYPE_NOTVALID             = 10
+# some generic internal allocation failure...
+const ERR_ALLOC                     = 11
+const ERR_DIM_NOTVALID              = 12
+const ERR_SPREAD_THREAD_NOTVALID    = 13
 
 struct FINUFFTError <: Exception
     errno::Cint
@@ -117,24 +164,33 @@ function check_ret(ret)
     # Check return value and output error messages
     if ret==0
         return
-    elseif ret==ERR_EPS_TOO_SMALL
-        msg = "requested tolerance epsilon too small"
+    elseif ret==WARN_EPS_TOO_SMALL
+        @warn "requested tolerance epsilon too small to achieve"
+        return
     elseif ret==ERR_MAXNALLOC
-        msg = "attemped to allocate internal arrays larger than MAX_NF (defined in common.h)"
+        msg = "attemped to allocate internal array larger than MAX_NF (defined in defs.h)"
     elseif ret==ERR_SPREAD_BOX_SMALL
-        msg = "spreader: fine grid too small"
+        msg = "spreader: fine grid too small compared to spread (kernel) width"
     elseif ret==ERR_SPREAD_PTS_OUT_RANGE
-        msg = "spreader: if chkbnds=1, a nonuniform point out of input range [-3pi,3pi]^d"
+        msg = "spreader: if chkbnds=1, a nonuniform point coordinate is out of input range [-3pi,3pi]^d"
     elseif ret==ERR_SPREAD_ALLOC
         msg = "spreader: array allocation error"
     elseif ret==ERR_SPREAD_DIR
         msg = "spreader: illegal direction (should be 1 or 2)"
     elseif ret==ERR_UPSAMPFAC_TOO_SMALL
-        msg = "upsampfac too small (should be >1)"
+        msg = "upsampfac too small (should be >1.0)"
     elseif ret==HORNER_WRONG_BETA
-        msg = "upsampfac not a value with known Horner eval: currently 2.0 or 1.25 only"
+        msg = "upsampfac not a value with known Horner poly eval rule (currently 2.0 or 1.25 only)"
     elseif ret==ERR_NDATA_NOTVALID
-        msg = "ndata not valid (should be >= 1)"
+        msg = "ntrans not valid in many (vectorized) or guru interface (should be >= 1)"
+    elseif ret==ERR_TYPE_NOTVALID
+        msg = "transform type invalid"
+    elseif ret==ERR_ALLOC
+        msg = "general allocation failure"
+    elseif ret==ERR_DIM_NOTVALID
+        msg = "dimension invalid"
+    elseif ret==ERR_SPREAD_THREAD_NOTVALID
+        msg = "spread_thread option invalid"
     else
         msg = "unknown error"
     end
@@ -422,7 +478,7 @@ function nufft1d1!(xj      :: Array{Float64},
                   Cdouble,
                   BIGINT,
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, cj, iflag, eps, ms, fk, opts
                  )
     check_ret(ret)
@@ -461,7 +517,7 @@ function nufft1d2!(xj      :: Array{Float64},
                   Cdouble,
                   BIGINT,
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, cj, iflag, eps, ms, fk, opts
                  )
     check_ret(ret)    
@@ -503,7 +559,7 @@ function nufft1d3!(xj      :: Array{Float64},
                   BIGINT,
                   Ref{Cdouble},            
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, cj, iflag, eps, nk, sk, fk, opts
                  )
     check_ret(ret)
@@ -549,7 +605,7 @@ function nufft2d1!(xj      :: Array{Float64},
                   BIGINT,
                   BIGINT,            
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                  )
     check_ret(ret)
@@ -593,7 +649,7 @@ function nufft2d2!(xj      :: Array{Float64},
                   BIGINT,
                   BIGINT,            
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                  )
     check_ret(ret)
@@ -642,7 +698,7 @@ function nufft2d3!(xj      :: Array{Float64},
                   Ref{Cdouble},
                   Ref{Cdouble},            
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, yj, cj, iflag, eps, nk, sk, tk, fk, opts
                  )
     check_ret(ret)
@@ -692,7 +748,7 @@ function nufft3d1!(xj      :: Array{Float64},
                   BIGINT,
                   BIGINT,
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                  )
     check_ret(ret)
@@ -740,7 +796,7 @@ function nufft3d2!(xj      :: Array{Float64},
                   BIGINT,
                   BIGINT,
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                  )
     check_ret(ret)
@@ -799,7 +855,7 @@ function nufft3d3!(xj      :: Array{Float64},
                   Ref{Cdouble},
                   Ref{Cdouble},                        
                   Ref{ComplexF64},
-                  nufft_opts),
+                  Ref{nufft_opts}),
                  nj, xj, yj, zj, cj, iflag, eps, nk, sk, tk, uk, fk, opts
                  )
     check_ret(ret)