This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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"), |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using DiffEqDevTools, BifurcationKit, Setfield, Plots | |
mutable struct TermCache | |
start | |
sign::Bool | |
count::Int | |
end | |
function stability_curve(tab) | |
function curve((x, ), (y, )) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#= | |
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
│ ─ %-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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) ) |
NewerOlder