YingboMa

using LinearAlgebra, BandedMatrices

n = 3000
A = Matrix(SymTridiagonal(Symmetric(rand(n, n)))); A[end, 1] = A[1, end] = 1;
b = rand(size(A, 1))

function periodic_givens_rotation!(A::AbstractMatrix{T}) where T
    m, n = size(A)
    gs = Vector{LinearAlgebra.Givens{T}}(undef, 2n-2)
    jj = 0

YingboMa

using Plots, FastGaussQuadrature, FastTransforms
fejer1(N) = FastTransforms.fejernodes1(Float64, N), FastTransforms.fejerweights1(FastTransforms.FastTransforms.chebyshevmoments1(Float64, N))
fejer2(N) = FastTransforms.fejernodes2(Float64, N), FastTransforms.fejerweights2(FastTransforms.FastTransforms.chebyshevmoments2(Float64, N))
clenshawcurtis(N) = FastTransforms.clenshawcurtisnodes(Float64, N), FastTransforms.clenshawcurtisweights(FastTransforms.FastTransforms.chebyshevmoments1(Float64, N))
x = range(-1.2, stop = 1.2, length = 1000)
y = range(-0.4, stop = 0.4, length = 1000÷3)
r(z, x, w) = sum(k->w[k]/(z - x[k]), eachindex(x))
ϕ(z) = log((z+1)/(z-1))
quadrature_contour(xs, w; kwargs...) = contour(x, y, (x, y) -> abs(r(x + y*im, xs, w) - ϕ(x + y*im)); levels=map(i->1/10^i, 0:10), aspect_ratio=1, colorbar=false, xlimits=(-1.5, 1.5), ylimits=(-0.5, 0.5), kwargs...)
plts = [quadrature_contour(clenshawcurtis(64)...; title="Clenshaw-Curtis"),

YingboMa

using DiffEqDevTools, BifurcationKit, Setfield, Plots

mutable struct TermCache
    start
    sign::Bool
    count::Int
end

function stability_curve(tab)
    function curve((x, ), (y, ))

YingboMa

using ForwardDiff
goo((x, y, z),) = [x^2*z, x*y*z, abs(z)-y]
foo((x, y, z),) = [x^2*z, x*y*z, abs(z)-y]
function foo(u::Vector{ForwardDiff.Dual{T,V,P}}) where {T,V,P}
    # unpack: AoS -> SoA
    vs = ForwardDiff.value.(u)
    # you can play with the dimension here, sometimes it makes sense to transpose
    ps = mapreduce(ForwardDiff.partials, hcat, u)
    # get f(vs)
    val = foo(vs)

YingboMa

#=
Replace calcJ, calcW, linsolve, do_newJW to [u, p, t, newJ, newW, tol, error_estimate]
foo!(z, integrator, alg, nlsolver) -> z .= W\z
foo(z, integrator, alg, nlsolver) -> W\z

foo!(z, integrators, alg, nlsolver)

- User control inverse of W(_t)
- user defined factorization (object with ldiv! defined) (Info: u, p, t, newW,
  tol, error_estimate)

YingboMa

A = [11//45-7sqrt(6)/360         37//225-169sqrt(6)/1800   -2//225+sqrt(6)/75
     37//225+169sqrt(6)/1800     11//45+7sqrt(6)/360       -2//225-sqrt(6)/75
     4//9-sqrt(6)/36             4//9+sqrt(6)/36               1//9]

T = Float64
T11  = convert(T, 9.1232394870892942792e-02)
T12  = convert(T, -0.14125529502095420843e0)
T13  = convert(T, -3.0029194105147424492e-02)
T21  = convert(T, 0.24171793270710701896e0)
T22  = convert(T, 0.20412935229379993199e0)

YingboMa

│ ─ %-1  = invoke perform_step!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},Float64,DiffEqBase.NullParameters,Float
64,Float64,Float64,Array{Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},1},ODESolution{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},2,Array{
Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float
64,4},1},1},1},ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(lorenz),L
inearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple
{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFuncti

YingboMa

define void @"julia_perform_step!_19254"(%jl_value_t addrspace(10)* nonnull align 8 dereferenceable(304), %jl_value_t addrspace(10)* nonnull align 8 dereferenceable(544), i8) {
top:
  %3 = addrspacecast %jl_value_t addrspace(10)* %0 to %jl_value_t addrspace(11)*
  %4 = bitcast %jl_value_t addrspace(11)* %3 to i8 addrspace(11)*
  %5 = getelementptr inbounds i8, i8 addrspace(11)* %4, i64 32
  %6 = bitcast i8 addrspace(11)* %5 to double addrspace(11)*
  %7 = load double, double addrspace(11)* %6, align 8
  %8 = getelementptr inbounds i8, i8 addrspace(11)* %4, i64 48
  %9 = bitcast i8 addrspace(11)* %8 to %jl_value_t addrspace(10)* addrspace(11)*
  %10 = load %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)* addrspace(11)* %9, align 8

YingboMa

# by Mason Protter
using MacroTools: MacroTools, splitdef, combinedef, @capture
macro specialize_vararg(n::Int, fdef)
    d = splitdef(fdef)
    args = d[:args][end]
    @assert d[:args][end] isa Expr && d[:args][end].head == Symbol("...") 
    args_symbol = d[:args][end].args[]
    fdefs = Expr(:block)
    for i in 1:n-1
        di = deepcopy(d)

YingboMa

using ModelingToolkit
using ModelingToolkit: Differential, expand_derivatives, Expression, Operation, simplify_constants

function D(f)
    (args...) -> begin
        syms = map(_->Variable(gensym())(), args)
        ex = f(syms...)
        ds = map(s->Differential(s), syms)
        ops = [expand_derivatives(ds[i](ex)) for i in eachindex(syms)]
        rules = Dict( collect(zip(syms, args)) )