SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.84k stars 224 forks source link

Failed while USING DifferentialEquations.jl in Julia 1.4.1 #623

Closed olaayeko closed 4 years ago

olaayeko commented 4 years ago

I am failing to use DifferentialEquations.jl (v6.14.0) in Julia 1.4.1 on Mac:

julia> Pkg.test("DifferentialEquations") Testing DifferentialEquations Status /private/var/folders/g9/yf_v1cnn01j5mplvzgxnlpnh0000gn/T/jl_w7zeiw/Manifest.toml [c3fe647b] AbstractAlgebra v0.9.2 [1520ce14] AbstractTrees v0.3.3 [79e6a3ab] Adapt v1.1.0 [ec485272] ArnoldiMethod v0.0.4 [7d9fca2a] Arpack v0.4.0 [68821587] Arpack_jll v3.5.0+3 [4fba245c] ArrayInterface v2.8.7 [4c555306] ArrayLayouts v0.3.4 [aae01518] BandedMatrices v0.15.11 [764a87c0] BoundaryValueDiffEq v2.5.0 [fa961155] CEnum v0.4.1 [a603d957] CanonicalTraits v0.2.1 [d360d2e6] ChainRulesCore v0.8.0 [861a8166] Combinatorics v1.0.2 [bbf7d656] CommonSubexpressions v0.2.0 [34da2185] Compat v3.12.0 [e66e0078] CompilerSupportLibraries_jll v0.3.3+0 [88cd18e8] ConsoleProgressMonitor v0.1.2 [187b0558] ConstructionBase v1.0.0 [adafc99b] CpuId v0.2.2 [9a962f9c] DataAPI v1.3.0 [864edb3b] DataStructures v0.17.18 [bcd4f6db] DelayDiffEq v5.24.1 [2b5f629d] DiffEqBase v6.38.3 [459566f4] DiffEqCallbacks v2.13.3 [5a0ffddc] DiffEqFinancial v2.3.0 [c894b116] DiffEqJump v6.9.1 [77a26b50] DiffEqNoiseProcess v4.2.0 [055956cb] DiffEqPhysics v3.5.0 [163ba53b] DiffResults v1.0.2 [b552c78f] DiffRules v1.0.1 [0c46a032] DifferentialEquations v6.14.0 [c619ae07] DimensionalPlotRecipes v1.2.0 [b4f34e82] Distances v0.9.0 [31c24e10] Distributions v0.23.4 [ffbed154] DocStringExtensions v0.8.2 [d4d017d3] ExponentialUtilities v1.6.0 [1a297f60] FillArrays v0.8.10 [6a86dc24] FiniteDiff v2.3.2 [59287772] Formatting v0.4.1 [f6369f11] ForwardDiff v0.10.10 [069b7b12] FunctionWrappers v1.1.1 [6b9d7cbe] GeneralizedGenerated v0.2.4 [01680d73] GenericSVD v0.3.0 [d25df0c9] Inflate v0.1.2 [42fd0dbc] IterativeSolvers v0.8.4 [82899510] IteratorInterfaceExtensions v1.0.0 [b14d175d] JuliaVariables v0.2.0 [b964fa9f] LaTeXStrings v1.1.0 [2ee39098] LabelledArrays v1.2.2 [23fbe1c1] Latexify v0.13.5 [1d6d02ad] LeftChildRightSiblingTrees v0.1.2 [093fc24a] LightGraphs v1.3.3 [d3d80556] LineSearches v7.0.1 [e6f89c97] LoggingExtras v0.4.1 [bdcacae8] LoopVectorization v0.8.5 [d00139f3] METIS_jll v5.1.0+4 [d8e11817] MLStyle v0.3.1 [1914dd2f] MacroTools v0.5.5 [e1d29d7a] Missings v0.4.3 [961ee093] ModelingToolkit v3.10.1 [46d2c3a1] MuladdMacro v0.2.2 [f9640e96] MultiScaleArrays v1.8.1 [d41bc354] NLSolversBase v7.6.1 [2774e3e8] NLsolve v4.4.0 [77ba4419] NaNMath v0.3.3 [71a1bf82] NameResolution v0.1.3 [6fe1bfb0] OffsetArrays v1.0.4 [4536629a] OpenBLAS_jll v0.3.9+4 [efe28fd5] OpenSpecFun_jll v0.5.3+3 [bac558e1] OrderedCollections v1.2.0 [1dea7af3] OrdinaryDiffEq v5.41.0 [90014a1f] PDMats v0.9.12 [65888b18] ParameterizedFunctions v5.3.0 [d96e819e] Parameters v0.12.1 [e409e4f3] PoissonRandom v0.4.0 [8162dcfd] PrettyPrint v0.1.0 [33c8b6b6] ProgressLogging v0.1.2 [92933f4c] ProgressMeter v1.3.1 [1fd47b50] QuadGK v2.3.1 [e6cf234a] RandomNumbers v1.4.0 [3cdcf5f2] RecipesBase v0.7.0 [731186ca] RecursiveArrayTools v2.4.4 [f2c3362d] RecursiveFactorization v0.1.2 [189a3867] Reexport v0.2.0 [ae029012] Requires v1.0.1 [ae5879a3] ResettableStacks v1.0.0 [79098fc4] Rmath v0.6.1 [f50d1b31] Rmath_jll v0.2.2+1 [f2b01f46] Roots v1.0.2 [21efa798] SIMDPirates v0.8.7 [476501e8] SLEEFPirates v0.5.1 [1bc83da4] SafeTestsets v0.0.1 [699a6c99] SimpleTraits v0.9.2 [a2af1166] SortingAlgorithms v0.3.1 [47a9eef4] SparseDiffTools v1.8.0 [276daf66] SpecialFunctions v0.10.3 [90137ffa] StaticArrays v0.12.3 [2913bbd2] StatsBase v0.32.2 [4c63d2b9] StatsFuns v0.9.5 [9672c7b4] SteadyStateDiffEq v1.5.1 [789caeaf] StochasticDiffEq v6.23.1 [bea87d4a] SuiteSparse_jll v5.4.0+8 [c3572dad] Sundials v4.2.3 [fb77eaff] Sundials_jll v5.2.0+0 [d1185830] SymbolicUtils v0.3.4 [3783bdb8] TableTraits v1.0.0 [5d786b92] TerminalLoggers v0.1.1 [a759f4b9] TimerOutputs v0.5.6 [a2a6695c] TreeViews v0.3.0 [3a884ed6] UnPack v1.0.1 [1986cc42] Unitful v1.2.1 [3d5dd08c] VectorizationBase v0.12.6 [19fa3120] VertexSafeGraphs v0.1.2 [700de1a5] ZygoteRules v0.2.0 [2a0f44e3] Base64 [ade2ca70] Dates [8bb1440f] DelimitedFiles [8ba89e20] Distributed [b77e0a4c] InteractiveUtils [76f85450] LibGit2 [8f399da3] Libdl [37e2e46d] LinearAlgebra [56ddb016] Logging [d6f4376e] Markdown [a63ad114] Mmap [44cfe95a] Pkg [de0858da] Printf [3fa0cd96] REPL [9a3f8284] Random [ea8e919c] SHA [9e88b42a] Serialization [1a1011a3] SharedArrays [6462fe0b] Sockets [2f01184e] SparseArrays [10745b16] Statistics [4607b0f0] SuiteSparse [8dfed614] Test [cf7118a7] UUIDs [4ec0a83e] Unicode Test Summary: | Pass Total Default Discrete Algorithm | 1 1 5.772186 seconds (16.05 M allocations: 868.707 MiB, 5.30% gc time) error: inline asm error: This value type register class is not natively supported! ERROR: Package DifferentialEquations errored during testing Stacktrace: [1] pkgerror(::String, ::Vararg{String,N} where N) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/Types.jl:53 [2] test(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}; coverage::Bool, julia_args::Cmd, test_args::Cmd, test_fn::Nothing) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/Operations.jl:1503 [3] test(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}; coverage::Bool, test_fn::Nothing, julia_args::Cmd, test_args::Cmd, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:316 [4] test(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:303 [5] #test#68 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:297 [inlined] [6] test at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:297 [inlined] [7] #test#67 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:296 [inlined] [8] test at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:296 [inlined] [9] test(::String; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:295 [10] test(::String) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/API.jl:295

olaayeko commented 4 years ago

Still getting the same error, is there a hard reset I need to complete first ?

chriselrod commented 4 years ago

Still getting the same error, is there a hard reset I need to complete first ?

Do you mean when testing DifferentialEquations, or trying to run vfmadd231(::Vec{8,Float64}...)?

If the latter, the point is: your CPU can do it, but LLVM doesn't think it can. So the change to VectorizationBase is meant to make some of DifferentialEquations.jl's dependencies use vfmadd231(::Vec{4,Float64}...) instead (basically, replace all the 8s with 4s, and it'll work).

olaayeko commented 4 years ago

When trying to run vfmadd231(::Vec{4,Float64}...), after changing the VectorizationBase settings

chriselrod commented 4 years ago

Are you sure that this does not work (replacing all the 8s with 4s)?

const Vec{W,T} = NTuple{W,Core.VecElement{T}}
@inline function vfmadd231(a::Vec{4,Float64}, b::Vec{4,Float64}, c::Vec{4,Float64})
    Base.llvmcall("""%res = call <4 x double> asm "vfmadd231pd \$3, \$2, \$1", "=v,0,v,v"(<4 x double> %2, <4 x double> %1, <4 x double> %0)
    ret <4 x double> %res""",
    Vec{4,Float64},
    Tuple{Vec{4,Float64},Vec{4,Float64},Vec{4,Float64}}, a, b, c)
end

x = ntuple(Val(4)) do _ Core.VecElement(rand()) end;
y = ntuple(Val(4)) do _ Core.VecElement(rand()) end;
z = ntuple(Val(4)) do _ Core.VecElement(rand()) end;
vfmadd231(x,y,z)
olaayeko commented 4 years ago

Sorry, yes that does work

maleadt commented 4 years ago

So, the problem is that the host's CPU features AVX512 support, but that LLVM thinks it is Haswell or something. Unfortunately, unsafe_string(LLVM.API.LLVMGetHostCPUFeatures()) returns -avx512f and all the other cpu feature flags.

So LLVM gets it wrong? https://github.com/llvm/llvm-project/commit/82921bf2baed96b700f90b090d5dc2530223d9c0 looks relevant then, and not included in our current LLVM (avx512f is set later based on HasLeaf7 && ((EBX >> 16) & 1) && HasAVX512Save).

olaayeko commented 4 years ago

Is there any action I can take to rectify this and any theory of why this bug appeared suddenly when it was working well before ?

chriselrod commented 4 years ago

So LLVM gets it wrong? llvm/llvm-project@82921bf looks relevant then, and not included in our current LLVM (avx512f is set later based on HasLeaf7 && ((EBX >> 16) & 1) && HasAVX512Save).

Yes, great. So it looks like the simplest approach to testing would be building Julia master with LLVM_VER=10.0.0 in the Make.user on that architecture. If that solves the problem, we could patch Julia's LLVM 9 for the 1.5 release?

EDIT: Or is there a way olaayeko could check if

bool HasAVX512Save = HasAVX && ((EAX & 0xe0) == 0xe0);

is currently set to false?

olaayeko commented 4 years ago

For the first one, how do I build that, sorry I am new to Julia. On the second one I am getting syntax/UndefVarError when I try to check that.

chriselrod commented 4 years ago

@maleadt Ha ha, I didn't realize that the feature string listed all the features, indicating with a + vs - whether or not they're present.

julia> using LLVM

julia> unsafe_string(LLVM.API.LLVMGetHostCPUFeatures())
"+sse2,+cx16,+sahf,-tbm,-avx512ifma,-sha,-gfni,-fma4,-vpclmulqdq,+prfchw,+bmi2,-cldemote,+fsgsbase,-ptwrite,+xsavec,+popcnt,+aes,-avx512bitalg,-movdiri,+xsaves,-avx512er,-avx512vnni,-avx512vpopcntdq,-pconfig,+clwb,+avx512f,-clzero,-pku,+mmx,-lwp,-rdpid,-xop,+rdseed,-waitpkg,-movdir64b,-sse4a,+avx512bw,+clflushopt,+xsave,-avx512vbmi2,+64bit,+avx512vl,+invpcid,+avx512cd,+avx,-vaes,+rtm,+fma,+bmi,+rdrnd,-mwaitx,+sse4.1,+sse4.2,+avx2,-wbnoinvd,+sse,+lzcnt,+pclmul,-prefetchwt1,+f16c,+ssse3,-sgx,-shstk,+cmov,-avx512vbmi,+movbe,+xsaveopt,+avx512dq,+adx,-avx512pf,+sse3"

julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i9-7900X CPU @ 3.30GHz

So the 7900X does correctly have +avx512f, but, for example, not -avx512er as it doesn't have exponentiation or accurate reciprocals. While olaayeko's CPU shows -avx512*, matching LLVM's behavior (although not the CPU's capabilities). So using this would let me generate code that agrees with LLVM. But it'd still be nice to stop underutilizing the host CPU on Apple.

Thanks @YingboMa for pointing out my misunderstanding.

So I think parsing this feature string is the right thing to do. I could drop the CpuId.jl dependency entirely if LLVM can somehow provide the number of physical cores as well as cache sizes? Is this possible at all?

I don't mind an LLVM.jl dependency too much, as I may consider using it in some cases instead of relying on Base.llvmcall. Any thoughts/concerns/recommendations? E.g., how do they compare on compile time? (llvmcall improved dramatically in Julia 1.4 versus early versions.) How robust is it versus changing underlying LLVM version? I've had a few issues with needing to checks based on Base.libllvm_version, but this is something I can easily handle myself. How stable/brittle is LLVM.jl with respect to this?

I can load the library and the feature string works with LLVM 10, but I get 1 test failure.

maleadt commented 4 years ago

I don't mind an LLVM.jl dependency too much, as I may consider using it in some cases instead of relying on Base.llvmcall. Any thoughts/concerns/recommendations?

LLVM.jl is a lightweight dependency, so OK for querying features. For replacing llvmcall, although it's a better (albeit sometimes verbose) way of emitting IR, the generated code is currently not precompilable because of how llvmcall with llvm::Function* arguments works (the package itself is perfectly precompilable though). Robustness is much better, since it uses the stable C ap vs. using the explicitly non-portable textual IR format. So it depends on your needs.

chriselrod commented 4 years ago

I'll wait until it is precompilable to switching to it more broadly. For now, I can use it for feature detection: https://github.com/chriselrod/VectorizationBase.jl/commit/1f1a22b3e1fbdcc67532320ad787a551545c7fcf

This also means I could add invariant.start and invariant.end methods, which I couldn't get to work with llvmcall (errors because it doesn't like {}*); these may be useful when working with arrays to which we have pointers but that are also static:

using LLVM, LLVM.Interop

struct Invariant{L,T}
    ivp::Ptr{Cvoid}
    ptr::Ptr{T}
end

@generated function invariant_start!(ptr::Ptr{T}, ::Val{L}) where {L,T}
    T_int = LLVM.IntType(sizeof(Int)*8, JuliaContext())
    paramtyps = [ T_int ]
    ret_typ = T_int # returning a Ptr{Cvoid}
    llvmf, _ = create_function(ret_typ, paramtyps)
    i8 = convert(LLVMType, Int8)
    i8_ptr = LLVM.PointerType(i8)

    emptystruct = LLVM.StructType(typeof(i8)[])
    emptystruct_ptr = LLVM.PointerType(emptystruct)

    mod = LLVM.parent(llvmf)
    intrinsic_typ = LLVM.FunctionType(emptystruct_ptr, [convert(LLVMType, Int), i8_ptr])
    intrinsic = LLVM.Function(mod, "llvm.invariant.start.p0i8", intrinsic_typ)

    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvmf, "entry", JuliaContext())
        position!(builder, entry)
        ptr = inttoptr!(builder, parameters(llvmf)[1], i8_ptr)
        ivt_ptr = call!(builder, intrinsic, [ConstantInt(T_int, sizeof(T)*L), ptr])
        ivt_int = ptrtoint!(builder, ivt_ptr, T_int)
        ret!(builder, ivt_int)
    end

    call_function(llvmf, Ptr{Cvoid}, Tuple{Ptr{T}}, :(ptr,))
end
@generated function invariant_end!(ivt_ptr::Ptr{Cvoid}, ptr::Ptr{T}, ::Val{L}) where {L,T}
    T_int = LLVM.IntType(sizeof(Int)*8, JuliaContext())
    voidtype = LLVM.VoidType(JuliaContext())

    paramtyps = [ T_int, T_int ]
    ret_typ = voidtype # returning a Ptr{Cvoid}
    llvmf, _ = create_function(ret_typ, paramtyps)
    i8 = convert(LLVMType, Int8)
    i8_ptr = LLVM.PointerType(i8)

    emptystruct = LLVM.StructType(typeof(i8)[])
    emptystruct_ptr = LLVM.PointerType(emptystruct)

    mod = LLVM.parent(llvmf)
    intrinsic_typ = LLVM.FunctionType(voidtype, [emptystruct_ptr, convert(LLVMType, Int), i8_ptr])
    intrinsic = LLVM.Function(mod, "llvm.invariant.end.p0i8", intrinsic_typ)

    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvmf, "entry", JuliaContext())
        position!(builder, entry)
        invt_ptr = inttoptr!(builder, parameters(llvmf)[1], emptystruct_ptr)
        ptr = inttoptr!(builder, parameters(llvmf)[2], i8_ptr)
        ret = call!(builder, intrinsic, [invt_ptr, ConstantInt(T_int, sizeof(T)*L), ptr])
        ret!(builder)
    end

    call_function(llvmf, Cvoid, Tuple{Ptr{Cvoid},Ptr{T}}, :(ivt_ptr, ptr))
end
@inline function freeze!(A_ptr::Ptr{T}, ::Val{L}) where {T,L}
    Invariant{L,T}(invariant_start!(A_ptr, Val(L)), A_ptr)
end
@inline function melt!(ivt::Invariant{L,T}) where {L,T}
    invariant_end!(ivt.ivp, ivt.ptr, Val(L))
end

The new release of VectorizationBase should be registered soon. It can only use LLVMGetHostCPUFeatures with LLVM is >= 8, so this wont be fixed on Julia versions <= 1.3. I hope not many people are still using these old versions.

olaayeko commented 4 years ago

I made the changes suggested here https://github.com/JuliaLang/julia/issues/36329, DifferentialEquations is passing Pkg.test(). In Jupyter notebook the DifferentialEquations module is still crashing the kernel. Any suggestions on next steps for debugging ? Thank you for the help so far as well.

chriselrod commented 4 years ago

When is DifferentialEquations crashing the kernel? Also, could you try running RecursiveFactorization's tests?

olaayeko commented 4 years ago

The kernel crashes when I attempt to run "solve(prob, Rodas4())". The RecursiveFactorization tests passed

olaayeko commented 4 years ago

and fails when I let it choose the solving method by emitting Rodas4()

olaayeko commented 4 years ago

What I find especially bizarre is that it was working perfectly fine earlier this month. Of note was I downloaded a packaged "FourierTransform", coincidentally(or not) it stopped working after this

chriselrod commented 4 years ago

And solve(prob, Rodas4()) works outside of Jupyter? Is versioninfo() the same inside and outside of Jupyter?

olaayeko commented 4 years ago

No its 1.5 in Julia and 1.4.2 in Jupyter, and it is working fine in terminal strangely

chriselrod commented 4 years ago

Could you try either of

  1. using Pkg; Pkg.update() in Jupyter, to make sure you're on the latest versions.
  2. using IJulia; installkernel() in the terminal, and then using the newly installed kernel in Jupyter so that you're using the same Julia version in both

and then let me know if you're still seeing the problem?

olaayeko commented 4 years ago

I am using the 1.5 kernel in Jupyter, its sorted now thank you for the help!

chriselrod commented 4 years ago

Great! Sounds like we can close this?

YingboMa commented 4 years ago

Thanks, everyone!