jwscook

using ForwardDiff, StaticArrays

function jacobi(f, x)
  x_real = @SArray [real(x), imag(x)]
  function f_real(v)
    z = complex(v[1], v[2])
    result = f(z)
    return @SArray [real(result), imag(result)]
  end
  J = ForwardDiff.jacobian(f_real, x_real)

jwscook

using ForwardDiff
using SpecialFunctions

import SpecialFunctions.erfcx
function SpecialFunctions.erfcx(ϕ::Complex{<:ForwardDiff.Dual{TAG}}) where {TAG}
  # Split input into real and imaginary parts
  x, y = reim(ϕ)
  # This gives the 'finite' complex part of the input
  z = complex(ForwardDiff.value(x), ForwardDiff.value(y))

jwscook

using LinearAlgebra, SparseArrays

struct StaticCondensationMatrix{T, M<:AbstractMatrix{T}} <: AbstractMatrix{T}
  A::M
  indices::Vector{UnitRange{Int}}
  nlocalblocks::Int
  ncouplingblocks::Int
  reducedlocalindices::Vector{UnitRange{Int}}
  reducedcoupledindices::Vector{UnitRange{Int}}
  reducedlhs::M

jwscook

using LinearAlgebra, Test

# Initialize matrices with proper dimensions
A11 = rand(2, 2)
A12 = rand(2, 2)
A21 = rand(2, 2)  # Doesn't necessarily need symmetry
A22 = rand(2, 2)
A23 = rand(2, 2)
A32 = rand(2, 2)  # Doesn't necessarily need symmetry
A33 = rand(2, 2)

jwscook

using Base.Threads, LinearAlgebra, SparseArrays
using Random; Random.seed!(0)
using ChunkSplitters
using BlockArrays
using ThreadPinning
pinthreads(:cores)
using Polyester

lsolve!(A, L::AbstractSparseArray) = (A .= L \ A) # can't mutate L
lsolve!(A, L) = ldiv!(A, lu(L), A) # could use a work array here in lu!

jwscook

using LinearAlgebra
function qr_factorization(A)
  m, n = size(A)
  
  Q = zeros(float(eltype(A)), m, m)
  for i in 1:m Q[i, i] = 1 end  
  R = zeros(float(eltype(A)), m, n)
  R .= A
  G = zeros(float(eltype(A)), 2, 2)
  

jwscook

using Random, Distributed, Test
Random.seed!(0)

addprocs(4, exeflags="-t 2")

@everywhere begin
using LinearAlgebra, Random
using Distributed, DistributedArrays, SharedArrays

LinearAlgebra.BLAS.set_num_threads(Base.Threads.nthreads())

jwscook

using LinearAlgebra, Random, Base.Threads

Random.seed!(0)
alphafactor(x::Real) = -sign(x)
alphafactor(x::Complex) = -exp(im * angle(x))

function householder!(A::Matrix{T}) where T
  m, n = size(A)
  H = A
  α = zeros(T, min(m, n)) # Diagonal of R

jwscook

using LinearAlgebra, Test, Random
Random.seed!(0)
@testset "understand raw lapack calls" begin

@show Threads.nthreads()
@show BLAS.get_num_threads()

m = 5
n = 4
A = rand(m, n);

jwscook

Two stream instability

Electrostatic linear plasma dispersion relation

$$ 1 + \Sigma_s \frac{\Pi_s^2}{\omega}\int_{-\infty}^{\infty} \frac{\partial f_s(v_\parallel)}{\partial v_\parallel}\frac{v_\parallel}{\omega-k_\parallel v_\parallel} dv_\parallel=0 $$

where, in SI units: