diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl index 1aa4b6f..d0ed264 100644 --- a/src/Characters/echelon_form.jl +++ b/src/Characters/echelon_form.jl @@ -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) @@ -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) @@ -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) @@ -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) @@ -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 @@ -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 diff --git a/src/actions.jl b/src/actions.jl index 2a9cb67..7bfdc70 100644 --- a/src/actions.jl +++ b/src/actions.jl @@ -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 diff --git a/src/matrix_projections.jl b/src/matrix_projections.jl index d8dd0fb..7816f3b 100644 --- a/src/matrix_projections.jl +++ b/src/matrix_projections.jl @@ -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 @@ -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 @@ -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 diff --git a/test/action_dihedral.jl b/test/action_dihedral.jl index 9314148..dab92dc 100644 --- a/test/action_dihedral.jl +++ b/test/action_dihedral.jl @@ -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, )