jwscook

using LinearAlgebra

"Compute Givens rotation parameters to zero out b"
function givens_rotation(a, b)
  T = promote_type(typeof(a), typeof(b))
  if iszero(b)
    return one(T), zero(T)
  elseif abs(b) > abs(a)
    t = a / b
    s = 1 / sqrt(1 + t^2)

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);