Skip to content

Commit

Permalink
Merge pull request #46 from kalmarek/enh/image_basis_perf
Browse files Browse the repository at this point in the history
Enh/image basis perf
  • Loading branch information
Marek Kaluba authored Jan 19, 2022
2 parents fe363d2 + 3ff2008 commit 636964f
Show file tree
Hide file tree
Showing 4 changed files with 81 additions and 51 deletions.
69 changes: 30 additions & 39 deletions src/Characters/echelon_form.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ end
Base.@propagate_inbounds function _mul_row!(
A::AbstractMatrix,
val,
row_idx,
row_idx;
starting_at = 1,
)
@boundscheck @assert 1 starting_at size(A, 2)
Expand All @@ -28,7 +28,7 @@ end

Base.@propagate_inbounds function _find_pivot(
A::AbstractMatrix,
col_idx,
col_idx;
starting_at = 1,
)
k = findnext(!iszero, @view(A[:, col_idx]), starting_at)
Expand All @@ -39,7 +39,7 @@ end
Base.@propagate_inbounds function _reduce_column_by_pivot!(
A::AbstractMatrix,
row_idx,
col_idx,
col_idx;
starting_at = 1,
)
@boundscheck checkbounds(A, row_idx, col_idx)
Expand All @@ -56,57 +56,47 @@ Base.@propagate_inbounds function _reduce_column_by_pivot!(
return A
end

# version over AbstractFloat
_finalize_pivot_reduce!(A::AbstractMatrix, pivot) = A

Base.@propagate_inbounds function __findmax(
f,
A::AbstractArray{T};
atol=eps(real(eltype(A)))*max(size(A)...)
) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
@assert !isempty(A)

midx = first(eachindex(A))
mval = f(A[midx])

@inbounds for idx in eachindex(A)
iszero(A[idx]) && continue
v = f(A[idx])
if v > mval
midx, mval = idx, v
end
if v < atol
A[idx] = zero(T)
end
end
# version over AbstractFloat

return midx, mval
end
const FloatOrComplex = Union{AbstractFloat,Complex{<:AbstractFloat}}

Base.@propagate_inbounds function _find_pivot(
A::AbstractMatrix{T},
col_idx,
starting_at = 1;
col_idx;
starting_at = 1,
atol = eps(real(eltype(A))) * size(A, 1),
) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
) where {T<:FloatOrComplex}
isempty(starting_at:size(A, 1)) && return false, starting_at
@boundscheck @assert 1 starting_at
# find the largest entry in the column below the last pivot
@boundscheck @assert 1 col_idx size(A, 2)

midx, mval = __findmax(abs, @view(A[starting_at:end, col_idx]), atol=atol)
mval, midx = findmax(abs, @view A[starting_at:end, col_idx])
if mval < atol # the largest entry is below threshold so we zero everything in the column
@view(A[starting_at:end, col_idx]) .= zero(T)
return false, starting_at
end
return true, oftype(starting_at, starting_at + midx - 1)
end

function _finalize_pivot_reduce!(
A::AbstractSparseMatrix{T},
pivot,
) where {T<:FloatOrComplex}
m = T <: Complex ? 2abs(eps(T)) : eps(T)
droptol!(A, max(size(A)...) * m)
return A
end

Base.@propagate_inbounds function _reduce_column_by_pivot!(
A::AbstractMatrix{T},
row_idx,
col_idx,
col_idx;
starting_at = 1,
atol = eps(real(eltype(A))) * size(A, 1),
) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
) where {T<:FloatOrComplex}
@boundscheck checkbounds(A, row_idx, col_idx)
@boundscheck 1 starting_at size(A, 2)

Expand All @@ -117,9 +107,8 @@ Base.@propagate_inbounds function _reduce_column_by_pivot!(
@inbounds A[ridx, col_idx] = zero(T)
continue
end
for cidx in starting_at:size(A, 2)
@inbounds A[ridx, cidx] -=
cidx == col_idx ? A[ridx, cidx] : v * A[row_idx, cidx]
@inbounds for cidx in starting_at:size(A, 2)
A[ridx, cidx] -= v * A[row_idx, cidx]
end
end
return A
Expand All @@ -129,24 +118,26 @@ Base.@propagate_inbounds function row_echelon_form!(A::AbstractMatrix)
pivots = Int[]
row_idx = 0 # also: current_rank
@inbounds for col_idx in 1:size(A, 2)
found, j = _find_pivot(A, col_idx, row_idx + 1)
found, j = _find_pivot(A, col_idx, starting_at = row_idx + 1)
found || continue
row_idx += 1
# @info col_idx
push!(pivots, col_idx)

# swap the rows so that A[row_idx, :] is the row with leading nz in
# i-th column
_swap_rows!(A, row_idx, j)
A = _swap_rows!(A, row_idx, j)

w = inv(A[row_idx, col_idx])

# multiply A[row_idx, :] by w
_mul_row!(A, w, row_idx, col_idx)
A = _mul_row!(A, w, row_idx, starting_at = col_idx)
# A[current_rank, col_idx] is 1 now

# zero the whole col_idx-th column above and below pivot:
# to the left of col_idx everything is zero
_reduce_column_by_pivot!(A, row_idx, col_idx, col_idx)
A = _reduce_column_by_pivot!(A, row_idx, col_idx, starting_at = col_idx)
A = _finalize_pivot_reduce!(A, col_idx)
end
return A, pivots
end
Expand Down
11 changes: 9 additions & 2 deletions src/actions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -137,8 +137,15 @@ function CachedExtensionHomomorphism(
if precompute
# to make sure that the future access to chom is read-only, i.e. thread-safe
# one may choose to precompute the images
for g in G
induce(chom, g)

lck = Threads.SpinLock()
@sync for g in G
Threads.@spawn begin
val = induce(action, chom.ehom, g)
lock(lck) do
chom.cache[g] = val
end
end
end
end
return chom
Expand Down
50 changes: 41 additions & 9 deletions src/matrix_projections.jl
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,6 @@ end

_mproj_fitsT!(args...) = _mproj_outsT!(args...)


function _mproj_outsT!(mproj::AbstractMatrix{T}, χ::Character) where {T}
mproj .= sum(
c .* matrix_projection_irr(Character(table(χ), i)) for
Expand Down Expand Up @@ -337,8 +336,11 @@ For characters or algebra elements return a basis of the row-space of
`matrix_projection([hom, ]χ)` or `matrix_representation([hom, ]α)`.
Two methods are employed to achieve the goal.
* By default (symbolic, exact) row echelon form is produced, and therefore there is no guarantee on the orthogonality of the returned basis vectors (rows).
* If `eltype(A) <: AbstractFloat` this is followed by a call to `svd` and the appropriate rows of its `Vt` are returned, thus the basis is (numerically) orthonormal.
* By default (symbolic, exact) row echelon form is produced, and therefore there
is no guarantee on the orthogonality of the returned basis vectors (rows).
* If `eltype(A) <: AbstractFloat` this is followed by a call to `svd` (or
`qr` in the sparse case) and the appropriate rows of its orthogonal factor are
returned, thus the basis is (numerically) orthonormal.
# Examples:
```julia
Expand Down Expand Up @@ -396,12 +398,42 @@ function image_basis!(A::AbstractSparseMatrix)
return A, p
end

function image_basis!(A::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Complex}}
_eps(::Type{T}) where {T} = T <: Complex ? 2eps(real(T)) : eps(T)

function _orth!(M::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Complex}}
F = svd!(convert(Matrix{T}, M))
M_rank = count(>(maximum(size(M)) * _eps(T)), F.S)
return F.Vt, M_rank
end

function _orth!(
M::AbstractSparseMatrix{T},
) where {T<:Union{AbstractFloat,Complex}}
F = qr(M)
M_rank = rank(F)
result = let tmp = F.Q * Matrix(I, size(F.Q, 2), M_rank)
sparse(tmp[invperm(F.prow), :]')
end
result = droptol!(result, maximum(size(result)) * _eps(T))
return result, M_rank
end

function image_basis!(
A::AbstractMatrix{T},
) where {T<:Union{AbstractFloat,Complex}}
A, p = row_echelon_form!(A)
fact = svd!(convert(Matrix{T}, @view(A[1:length(p), :])))
m = T <: Complex ? 2eps(real(T)) : eps(T)
A_rank = sum(fact.S .> maximum(size(A)) * m)
@assert A_rank == length(p)
return fact.Vt, 1:length(p)
A_orth, A_rank = _orth!(@view A[1:length(p), :])
@assert A_rank == length(p) "_orth rank doesn't agree with rref rank!"
return A_orth, 1:A_rank
end

function image_basis!(
A::AbstractSparseMatrix{T},
) where {T<:Union{AbstractFloat,Complex}}
N = LinearAlgebra.checksquare(A)
droptol!(A, N * _eps(T))
# A, p = row_echelon_form!(A)
A_orth, A_rank = _orth!(A)
# @assert A_rank == length(p) "_orth rank doesn't agree with rref rank!"
return A_orth, 1:A_rank
end
2 changes: 1 addition & 1 deletion test/action_dihedral.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ const OPTIMIZER = optimizer_with_attributes(
SCS.Optimizer,
"acceleration_lookback" => 10,
"max_iters" => 3_000,
"alpha" => 1.2,
"alpha" => 1.0,
"eps" => 1e-6,
"linear_solver" => SCS.DirectSolver,
)
Expand Down

0 comments on commit 636964f

Please sign in to comment.