Assorted fixes/improvements

src/SmithNormalForm.jl

@@ -1,48 +1,88 @@
 module SmithNormalForm
 
-# import Base: show, getindex
 using LinearAlgebra
 using SparseArrays
 using Base.CoreLogging
 
-import Base: show
-import LinearAlgebra: diagm
+import Base: show, summary
+import LinearAlgebra: diagm, diag
 
-export snf, smith
+export snf, smith, Smith
 
 include("bezout.jl")
 include("snf.jl")
 
-struct Smith{P,Q<:AbstractMatrix} <: Factorization{P}
+# ---------------------------------------------------------------------------------------- #
+
+struct Smith{P,Q<:AbstractMatrix{P},V<:AbstractVector{P}} <: Factorization{P}
     S::Q
     Sinv::Q
     T::Q
     Tinv::Q
-    SNF::AbstractVector{P}
-    Smith{P,Q}(S::AbstractMatrix{P}, Sinv::AbstractMatrix{P},
-               T::AbstractMatrix{P}, Tinv::AbstractMatrix{P},
-               A::AbstractVector{P}) where {P,Q} = new(S, Sinv, T, Tinv, A)
+    SNF::V
+    Smith{P,Q,V}(S::AbstractMatrix{P}, Sinv::AbstractMatrix{P},
+                 T::AbstractMatrix{P}, Tinv::AbstractMatrix{P},
+                 D::AbstractVector{P}) where {P,Q,V} = new(S, Sinv, T, Tinv, D)
 end
-Smith(S::AbstractMatrix{P}, T::AbstractMatrix{P}, A::AbstractVector{P}) where {P} =
-    Smith{P,typeof(A)}(S, similar(S, 0, 0), T, similar(T, 0, 0), A)
+Smith(S::AbstractMatrix{P}, T::AbstractMatrix{P}, SNF::AbstractVector{P}) where {P} =
+    Smith{P,typeof(S),typeof(SNF)}(S, similar(S, 0, 0), T, similar(T, 0, 0), SNF)
+
+# ---------------------------------------------------------------------------------------- #
+
+"""
+    smith(X::AbstractMatrix{P}; inverse::Bool=true) --> Smith{P,Q,V}
+
+Return a Smith normal form of an integer matrix `X` as a `Smith` structure (of element type
+`P`, matrix type `Q`, and invariant factor type `V`).
+
+The Smith normal form is well-defined for any matrix ``m×n`` matrix `X` with elements in a
+principal domain (PID; e.g., integers) and provides a decomposition of `X` into ``m×m``,
+``m×n``, and `S`, `Λ`, and `T` as `X = SΛT`, where `Λ` is a diagonal matrix with entries 
+("invariant factors") Λᵢ ≥ Λᵢ₊₁ ≥ 0 with nonzero entries divisible in the sense Λᵢ | Λᵢ₊₁.
+The invariant factors can be obtained from [`diag(::Smith)`](@ref).
 
-function smith(X::AbstractMatrix{P}; inverse=true) where {P}
-    S, T, A, Sinv, Tinv = snf(X, inverse=inverse)
-    return Smith{P, typeof(X)}(S, Sinv, T, Tinv, diag(A))
+`S` and `T` are invertible matrices; if the keyword argument `inverse` is true (default),
+the inverse matrices are computed and returned as part of the `Smith` factorization.
+"""
+function smith(X::AbstractMatrix{P}; inverse::Bool=true) where {P}
+    S, T, D, Sinv, Tinv = snf(X, inverse=inverse)
+    SNF = diag(D)
+    return Smith{P, typeof(X), typeof(SNF)}(S, Sinv, T, Tinv, SNF)
 end
 
-"""Retrive SNF matrix from the factorization"""
-function diagm(F::Smith{P,Q}) where {P,Q}
-    base = issparse(F.SNF) ? spzeros(P, size(F.S,1), size(F.T,1)) : zeros(P, size(F.S,1), size(F.T,1))
-    for i in 1:length(F.SNF)
-        base[i,i] = F.SNF[i]
+# ---------------------------------------------------------------------------------------- #
+
+"""
+    diagm(F::Smith) --> AbstractMatrix
+
+Return the Smith normal form diagonal matrix `Λ` from a factorization `F`.
+"""
+function diagm(F::Smith{P}) where P
+    rows = size(F.S, 1)
+    cols = size(F.T, 1)
+    D    = issparse(F.SNF) ? spzeros(P, rows, cols) : zeros(P, rows, cols)
+    for (i,Λᵢ) in enumerate(diag(F))
+        D[i,i] = Λᵢ
     end
-    return base
+    return D
 end
 
-function show(io::IO, F::Smith)
-    println(io, "Smith normal form:")
-    show(io, diagm(F))
+"""
+    diag(F::Smith) --> AbstractVector
+
+Return the invariant factors (or, equivalently, the elementary divisors) of a Smith normal
+form `F`.
+"""
+diag(F::Smith) = F.SNF
+
+summary(io::IO, F::Smith{P,Q,V}) where {P,Q,V} = print(io, "Smith{", P, ",", Q, ",", V, "}")
+function show(io::IO, ::MIME"text/plain", F::Smith)
+    summary(io, F)
+    println(io, " with:")
+    Base.print_matrix(io, diagm(F), " Λ = "); println(io)
+    Base.print_matrix(io, F.S, " S = "); println(io)
+    Base.print_matrix(io, F.T, " T = ")
+    return nothing
 end
 
 end # module

src/snf.jl

@@ -1,6 +1,6 @@
 # Smith Normal Form
 
-function divisable(y::R, x::R ) where {R}
+function divisable(y::R, x::R) where {R}
   x == zero(R) && return y == zero(R)
   return div(y,x)*x == y
 end
@@ -115,11 +115,11 @@ function colpivot(V::AbstractArray{R,2},
     end
 end
 
-function smithpivot(U::AbstractArray{R,2},
-                    Uinv::AbstractArray{R,2},
-                    V::AbstractArray{R,2},
-                    Vinv::AbstractArray{R,2},
-                    D::AbstractArray{R,2},
+function smithpivot(U::AbstractMatrix{R},
+                    Uinv::AbstractMatrix{R},
+                    V::AbstractMatrix{R},
+                    Vinv::AbstractMatrix{R},
+                    D::AbstractMatrix{R},
                     i, j; inverse=true) where {R}
 
     pivot = D[i,j]
@@ -240,11 +240,11 @@ function snf(M::AbstractMatrix{R}; inverse=true) where {R}
         j > cols && break
         Λⱼ = D[j,j]
         if Λⱼ < 0
-            @views V[j,:]    .*= -1 # T′   = sign(Λ)*T    [rows]
+            @views V[j,:] .*= -1        # T′   = sign(Λ)*T    [rows]
             if inverse
                 @views Vinv[:,j] .*= -1 # T⁻¹′ = T⁻¹*sign(Λ)  [columns]
             end
-            D[j,j] = abs(Λⱼ)        # Λ′ = Λ*sign(Λ)
+            D[j,j] = abs(Λⱼ)            # Λ′ = Λ*sign(Λ)
         end
     end
     @logmsg (Base.CoreLogging.Debug-1) "Factorization" D=formatmtx(D) U=formatmtx(U) V=formatmtx(V) U⁻¹=formatmtx(Uinv) V⁻¹=formatmtx(Vinv)