1module fastsum 2 3ending = ".so" 4 5if Sys.iswindows( ) 6 ending = ".dll" 7elseif Sys.isapple() 8 ending = ".dylib" 9end 10 11const lib_path = string( @__DIR__, "/libfastsumjulia", ending ) 12 13mutable struct fastsum_plan 14end 15 16mutable struct Plan 17 d::Integer # dimension 18 N::Integer # number of source nodes 19 M::Integer # number of target nodes 20 n::Integer # expansion degree 21 p::Integer # degree of smoothness 22 kernel::String # name of kernel 23 c::Vector{Float64} # kernel parameters 24 eps_I::Real # inner boundary 25 eps_B::Real # outer boundary 26 nn_x::Integer # oversampled nn in x 27 nn_y::Integer # oversampled nn in y 28 m_x::Integer # NFFT-cutoff in x 29 m_y::Integer # NFFT-cutoff in y 30 init_done::Bool # bool for plan init 31 finalized::Bool # bool for finalizer 32 flags::UInt32 # flags 33 x::Ref{Float64} # source nodes 34 y::Ref{Float64} # target nodes 35 alpha::Ref{ComplexF64} # source coefficients 36 f::Ref{ComplexF64} # target evaluations 37 plan::Ref{fastsum_plan} # plan (C pointer) 38 39 function Plan( d::Integer, N::Integer, M::Integer, n::Integer, p::Integer, kernel::String, c::Vector{<:Real}, eps_I::Real, eps_B::Real, nn_x::Integer, nn_y::Integer, m_x::Integer, m_y::Integer, flags::UInt32 ) 40 new( d, N, M, n, p, kernel, Vector{Float64}(c), eps_I, eps_B, nn_x, nn_y, m_x, m_y, false, false, flags ) 41 end 42end #struct fastsumplan 43 44function Plan( d::Integer, N::Integer, M::Integer, n::Integer, p::Integer, kernel::String, c::Real, eps_I::Real, eps_B::Real, nn::Integer, m::Integer ) 45 46 if N <= 0 47 error( "N has to be a positive Integer." ) 48 end 49 50 if M <= 0 51 error( "M has to be a positive Integer." ) 52 end 53 54 if n <= 0 55 error( "n has to be a positive Integer. ") 56 end 57 58 if m <= 0 59 error( "m has to be a positive Integer." ) 60 end 61 62 cv = Vector{Float64}( undef, 1 ) 63 cv[1] = Float64(c) 64 65 Plan( d, N, M, n, p, kernel, cv, eps_I, eps_B, nn, nn, m, m, UInt32(0) ) 66 67end #constructor 68 69function fastsum_init( p::Plan ) 70 71 ptr = ccall( ("jfastsum_alloc", lib_path), Ptr{fastsum_plan}, () ) 72 Core.setfield!( p, :plan, ptr ) 73 74 code = ccall( ("jfastsum_init",lib_path), Int64, (Ref{fastsum_plan}, Int32, Cstring, Ref{Float64}, UInt32, Int32, Int32, Float64, Float64, Int32, Int32, Int32, Int32, Int32, Int32), ptr, Int32(p.d), p.kernel, p.c, p.flags, Int32(p.n), Int32(p.p), Float64(p.eps_I), Float64(p.eps_B), Int32(p.N), Int32(p.M), Int32(p.nn_x), Int32(p.nn_y), Int32(p.m_x), Int32(p.m_y) ) 75 76 if code == 1 77 error( "Unkown kernel." ) 78 end 79 80 Core.setfield!( p, :init_done, true ) 81 finalizer( finalize_plan, p ) 82 83end #fastsum_init 84 85function finalize_plan( p::Plan ) 86 87 if !p.init_done 88 error( "Plan not initialized." ) 89 end 90 91 if !p.finalized 92 ccall( ("jfastsum_finalize",lib_path), Nothing, (Ref{fastsum_plan},), p.plan ) 93 Core.setfield!( p, :finalized, true ) 94 end 95 96end #finalize_plan 97 98function Base.setproperty!( p::Plan, v::Symbol, val ) 99 100 if !p.init_done 101 fastsum_init( p ) 102 end 103 104 if p.finalized 105 error( "Plan already finalized" ) 106 end 107 108 # edit source nodes 109 if v == :x 110 111 if p.d == 1 112 if typeof(val) != Vector{Float64} 113 error( "x has to be a Float64 vector." ) 114 end 115 if size(val)[1] != p.N 116 error( "x has to be a Float64 vector of length N." ) 117 end 118 else # => D >1 119 if typeof(val) != Array{Float64, 2} 120 error( "x has to be a Float64 matrix." ) 121 end 122 if size(val)[1] != p.N || size(val)[2] != p.d 123 error( "x has to be a Float64 matrix of size N." ) 124 end 125 end 126 127 ptr = ccall( ("jfastsum_set_x", lib_path), Ptr{Float64}, (Ref{fastsum_plan},Ref{Cdouble}), p.plan, val ) 128 Core.setfield!( p, v, ptr ) 129 130 # edit target nodes 131 elseif v == :y 132 133 if p.d == 1 134 if typeof(val) != Vector{Float64} 135 error( "y has to be a Float64 vector." ) 136 end 137 if size(val)[1] != p.M 138 error( "y has to be a Float64 vector of length M." ) 139 end 140 else # => D > 1 141 if typeof(val) != Array{Float64, 2} 142 error( "y has to be a Float64 matrix." ) 143 end 144 if size(val)[1] != p.M || size(val)[2] != p.d 145 error( "y has to be a Float64 matrix of size M." ) 146 end 147 end 148 149 ptr = ccall( ("jfastsum_set_y", lib_path), Ptr{Float64}, (Ref{fastsum_plan},Ref{Cdouble}), p.plan, val ) 150 Core.setfield!( p, v, ptr ) 151 152 # edit source coefficients 153 elseif v == :alpha 154 155 if typeof(val) != Vector{ComplexF64} 156 error( "alpha has to be a ComplexF64 vector." ) 157 end 158 if size(val)[1] != p.N 159 error( "alpha has to be a ComplexF64 vector of length N." ) 160 end 161 162 ptr = ccall( ("jfastsum_set_alpha", lib_path), Ptr{ComplexF64}, (Ref{fastsum_plan},Ref{ComplexF64}), p.plan, val ) 163 Core.setfield!( p, v, ptr ) 164 165 elseif v == :M 166 @warn("You can't modify the number of target nodes.") 167 elseif v == :N 168 @warn("You can't modify the number of source nodes.") 169 elseif v == :n 170 @warn("You can't modify the expansion degree.") 171 elseif v == :m 172 @warn("You can't modify the cut-off parameter.") 173 elseif v == :p 174 @warn("You can't modify the degree of smoothness.") 175 elseif v == :kernel 176 @warn("You can't modify the kernel.") 177 elseif v == :c 178 @warn("You can't modify the kernel parameters.") 179 elseif v == :eps_I 180 @warn("You can't modify the inner boundary.") 181 elseif v == :eps_B 182 @warn("You can't modify the outer boundary.") 183 elseif v == :plan 184 @warn("You can't modify the pointer to the fastsum plan.") 185 else 186 Core.setfield!( p, v, val ) 187 end 188 189end # Base.setproperty! 190 191# overwrite dot notation for plan struct in order to use C memory 192function Base.getproperty( p::Plan, v::Symbol ) 193 194 if v == :x 195 196 if !isdefined( p, :x ) 197 error("x is not set.") 198 end 199 200 ptr = Core.getfield( p, :x ) 201 202 if p.d == 1 203 return unsafe_wrap( Vector{Float64}, ptr, p.N ) # get source nodes from C memory and convert to Julia type 204 else 205 return unsafe_wrap( Matrix{Float64}, ptr, (p.d, p.N) ) # get source nodes from C memory and convert to Julia type 206 end 207 208 elseif v == :y 209 210 if !isdefined( p, :y ) 211 error("y is not set.") 212 end 213 214 ptr = Core.getfield( p, :y ) 215 216 if p.d == 1 217 return unsafe_wrap( Vector{Float64}, ptr, p.M ) 218 else 219 return unsafe_wrap( Matrix{Float64}, ptr, (p.d, p.M) ) 220 end 221 222 elseif v == :alpha 223 224 if !isdefined( p, :alpha ) 225 error("alpha is not set.") 226 end 227 228 ptr = Core.getfield( p, :alpha ) 229 return unsafe_wrap( Vector{ComplexF64}, ptr, p.N ) # get coefficients from C memory and convert to Julia type 230 231 elseif v == :f 232 233 if !isdefined( p, :f ) 234 error("f is not set.") 235 end 236 237 ptr = Core.getfield( p, :f ) 238 return unsafe_wrap( Vector{ComplexF64}, ptr, p.M ) # get function values from C memory and convert to Julia type 239 240 else 241 return Core.getfield( p, v ) 242 end 243end # Base.getproperty 244 245function trafo( p::Plan ) 246 247 if p.finalized 248 error( "Plan already finalized." ) 249 end 250 251 if !isdefined( p, :x ) 252 error( "x has not been set." ) 253 end 254 255 if !isdefined( p, :y ) 256 error( "y has not been set." ) 257 end 258 259 if !isdefined( p, :alpha ) 260 error( "alpha has not been set." ) 261 end 262 263 ptr = ccall( ( "jfastsum_trafo", lib_path), Ptr{ComplexF64}, (Ref{fastsum_plan},), p.plan ) 264 Core.setfield!( p, :f, ptr ) 265 266end #trafo 267 268function trafo_exact( p::Plan ) 269 270 if p.finalized 271 error( "Plan already finalized." ) 272 end 273 274 if !isdefined( p, :x ) 275 error( "x has not been set." ) 276 end 277 278 if !isdefined( p, :y ) 279 error( "y has not been set." ) 280 end 281 282 if !isdefined( p, :alpha ) 283 error( "alpha has not been set." ) 284 end 285 286 ptr = ccall( ("jfastsum_exact", lib_path), Ptr{ComplexF64}, (Ref{fastsum_plan},), p.plan ) 287 Core.setfield!( p, :f, ptr ) 288 289end #trafo 290 291end #module 292