UnitfulLinearAlgebra

function Base.Matrix(A::Union{AbstractUnitfulMatrix,AbstractUnitfulDimMatrix}) 

Expand A into array form
Useful for tests, display
pp. 193, Hart

DSVD <: Factorization

Matrix factorization type of the dimensioned singular value decomposition (DSVD) of a matrix A. This is the return type of dsvd(_), the corresponding matrix factorization function.

If F::DSVD is the factorization object, U, S, V and V⁻¹ can be obtained via F.U, F.S, F.V and F.V⁻¹, such that A = U * Diagonal(S) * V⁻¹. The singular values in S are sorted in descending order.

Iterating the decomposition produces the components U, S, and V.

Differences from SVD struct: Vt -> V⁻¹, U and V can have different types.

Functions available for DSVD: size, dsvdvals, inv. Function available for SVD that would be good to have to DSVD: ldiv!, transpose. ```

struct UnitfulMatrix

Extend DimArray to use dimensions for units, also add exact boolean flag

struct UnitfulMatrix{T,N,D<:Tuple,A<:AbstractArray{T,N}} <: AbstractUnitfulVecOrMat{T,N,D,A} data::A dims::D exact::Bool end

function UnitfulMatrix(A::AbstractMatrix)

Constructor to make inexact UnitfulMatrix. Satisfies algebraic interpretation of multipliable matrices.

function *(A::AbstractUnitfulType,b)

Matrix-scalar multiplication with units/dimensions.
Must account for change in the unitrange when the
 scalar has units.
Here, take product of dimension of the scalar and the unitrange.
Alternatively, divide the domain by the dimension of the scalar. 
Matrix-scalar multiplication is commutative.
Result is `exact` if input matrix is exact and scalar is dimensionless. 

function *(A,B)

Matrix-matrix multiplication with units/dimensions.
A*B represents two successive transformations.
Unitrange of B should equal domain of A in geometric interpretation.
Unitrange of B should be parallel to unitdomain of A in algebraic interpretation.

function +(A,B)

Matrix-matrix addition with units/dimensions.
A+B requires the two matrices to have dimensional similarity.

function -(A,B)

Matrix-matrix subtraction with units/dimensions.
A-B requires the two matrices to have dimensional similarity.

 function left divide

 Left divide of Multipliable Matrix.
 Reverse mapping from unitdomain to range.
 Is `exact` if input is exact.

function ldiv(F::LU{T,MultipliableMatrix{T},Vector{Int64}}, B::AbstractVector) where T<:Number

Perform matrix left divide on LU factorization object,
where LU object contains unit information.
Doesn't require LeftUniformMatrix.

 function inv

 Inverse of Multipliable Matrix.
 Only defined for nonsingular matrices.
 Inverse reverses mapping from unitdomain to range.
 Is `exact` if input is exact.

Hart, pp. 205.

Extend inv for Eigen factorizations of MultipliableMatrixs. Only defined for matrices with uniform units (pp. 101, Hart, 1995).

function vec(A::AbstractUnitfulType)

return a Vector{Quantity}
note ambiguity whether this function should return a Vector{Quantity} or an `AbstractUnitfulType` with one column

Arguments

• A::AbstractUnitfulType: input matrix

rebuild(A::UnitfulMatrix, data, [dims, exact]) => UnitfulMatrix
rebuild(A::UnitfulMatrix; kw...) => UnitfulMatrix

Rebuild a UnitfulMatrix with some field changes. All types that inherit from UnitfulMatrix must define this method if they have any additional fields or alternate field order.

This method can also be used with keyword arguments in place of regular arguments.

rebuild(A::UnitfulMatrix, data, dims, exact) => UnitfulMatrix
rebuild(A::UnitfulMatrix; kw...) => UnitfulMatrix

Rebuild a UnitfulMatrix with new fields. Handling partial field update is dealt with in rebuild for AbstractDimArray (still true?).

function cholesky(A::AbstractMultipliableMatrix)

Cholesky decomposition extended for matrices with units.
Requires unit (or dimensionally) symmetric matrix.
Functions available for LinearAlgebra.Cholesky objects: `size`, ``, `inv`, `det`, `logdet` and `isposdef`.
Functions available for UnitfulLinearAlgebra.Cholesky objects: `size`, `det`, and `isposdef`.

function det

Unitful matrix determinant.

function diag(A::AbstractUnitfulMatrix)

Diagonal elements of matrix with units.

Usual LinearAlgebra.diag function is not working due to different type elements on diagonal

function diagm(v::AbstractVector,r::Unitful.Units,d::Unitful.Units; exact = false)

Construct diagonal matrix with units where the diagonal has elements `v`.
If `v` has units, check that they conform with dimensional unit range `r`
 and dimensional unit domain `d`. Works for square or non-square matrices.

function eigen(A::T;permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T <: AbstractMultipliableMatrix

Thin wrapper for `UnitfulLinearAlgebra.eigen` with same keyword arguments as `LinearAlgebra.eigen`.
Only squarable matrices have eigenstructure (pp. 96, Hart, 1995).
There are multiple ways to distribute the units amongst the values and vectors.
If 𝐀 is endomorphic (i.e., the unit domain and range are the same), then the unit domain should
be taken as the units of the eigenvectors (pp. 205, Hart, 1995).  
In the general case, physical intuition and the equation 𝐀𝐱 = λ𝐱
dictate that the units of the eigenvectors are equal to the unit domain of 𝐀 (pp. 206, Hart, 1995).
Ideally the AbstractArray interface would automatically handle `eigen`,
but there is an unsolved issue with Unitful conversions.
The following functions are available for `Eigen` objects:  [`det`](@ref), [`inv`](@ref) and [`isposdef`](@ref). Some are restricted to uniform matrices.
`eigvals` of Eigen struct also available.

Extend isposdef for Eigen factorizations of MultipliableMatrixs. Should the units be stripped out of the function? Only defined for matrices with uniform units (pp. 101, Hart, 1995).

function lu(A::AbstractUnitfulVecOrMat{T})

Extend `lu` factorization to AbstractMultipliableMatrix.
Related to Gaussian elimination.
Store dimensional domain and range in "factors" attribute
even though this is not truly a MultipliableMatrix.
Returns `LU` type in analogy with `lu` for unitless matrices.
Based on LDU factorization, Hart, pp. 204.

svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD

Singular value decomposition (SVD) of `AbstractMultipliableMatrix`.
Only exists for uniform matrices (pp. 124, Hart, 1995).
Functions for `SVD{AbstractMultipliableMatrix}` object: `inv`, `size`, `adjoint`, `svdvals`.
Not implemented: `ldiv!`.

function convert_unitdomain(A, newdomain)

When using the geometric interpretation of matrices,
it is useful to convert the dimensional domain of the
matrix to match the expected vectors during multiplication.
Here we set the matrix to `exact=true` after this step.

function convert_unitrange(A, newrange)

When using the geometric interpretation of matrices,
it is useful to convert the dimensional range of the
matrix to match the desired output of multiplication.
Here we set the matrix to `exact=true` after this step.
Permits MatrixType to change.

function describe(A::UnitfulMatrix)

 Information regarding the type of multipliable matrix.

function dottable(a,b)

Are two quantities dimensionally compatible to take a dot product?

function dsvd(A::AbstractMultipliableMatrix,Prange::UnitSymmetricMatrix,Pdomain::UnitSymmetricMatrix;full=false,alg::LinearAlgebra.Algorithm = LinearAlgebra.default_svd_alg(A.numbers)) 

Dimensioned singular value decomposition (DSVD).
Appropriate version of SVD for non-uniform matrices.
`svd` can be computed for `Number`s, `Adjoint`s, `Transpose`s, and `Integers`; `dsvd` doesn't yet implement these.

Input

• A::AbstractMultipliableMatrix
• Pr::UnitSymmetricMatrix: square matrix defining norm of range
• Pd::UnitSymmetricMatrix: square matrix defining norm of domain
• full=false: optional argument
• alg: optional argument for algorithm

Output:

• F::DSVD: Dimensioned SVD object with units that can be deconstructed

function endomorphic(A)::Bool

Endomorphic matrices have a particular structure
 of the units/dimensions in the array.
It maps dimensioned vector space to itself.
Equivalent unit (dimensional) range and domain.

function exact(A)

• exact=true: geometric interpretation of unitdomain and unitrange
• exact=false: algebraic interpretation

function getindexqty

Get entry value of matrix including units.
Note: Calling B::UnitfulMatrix[i,j] doesn't currently return the units.

function identitymatrix(dimrange)

Input: dimensional (unit) range.
`A + I` only defined when `endomorphic(A)=true`
When accounting for units, there are many identity matrices.
This function returns a particular identity matrix
defined by its dimensional range.
Hart, pp. 200.

Maybe change the name to UnitfulI?

function left_uniform(A)

Definition: uniform unitrange of A
Left uniform matrix: output of matrix has uniform units

function multipliable(A)::Bool

Is an array multipliable?
It requires a particular structure of the units/dimensions in the array.

function parallel

Vector a is dimensionally parallel to vector b if they have the same length and a consistent dimensional change relates corresponding components. Guaranteed if two vectors are dimensionally similar. True for scalars in all cases.

pp. 188, Hart Note: Hart uses ≈, but this conflicts with an existing Julia function.

function right_uniform(A)

Does the unitdomain of A have uniform dimensions?
Right uniform matrix: input of matrix must have uniform units

function setindex!(A::AbstractUnitfulMatrix,v,i,j)

Set element (i,j) of a multipliable matrix.
Part of the AbstractArray interface.

#Input

• A::AbstractUnitfulMatrix
• v: new value
• i::Integer: row index
• j::Integer: column index

#Output

• Quantity: numerical value and units

function similarity(a,b)::Bool

Dimensional similarity of vectors, a binary relation Read "a has the same dimensional form as b" a and b may still have different units. A stronger condition than being parallel. pp. 184, Hart

function squarable(A::Matrix)

A squarable matrix is one where 𝐀² is defined.
Unit (dimensional) range and domain are parallel.
Key for solving difference and differential equations.
Have eigenstructure.

function square(A)

size(A)[1] == size(A)[2]

function trace(A)

Trace = sum of diagonal elements of a square matrix

function uniform(a)

Is the dimension of this quantity uniform?

There must be a way to inspect the Unitful type to answer this. Uniform matrix: All entries have the same units

function unit_symmetric(A::Matrix)

`UnitSymmetricMatrix`s have units that are symmetric about the main diagonal and define weighted norms. 
Definition: inverse dimensional range and dimensional domain are parallel.
Called "dimensionally symmetric" by Hart, 1995.

function unitdims(A::AbstractUnitfulVecOrMat) = dims(A)

Return tuple -> (unitrange, unitdomain)

for UnitfulMatrix, unit information overrides the `dims` function.

function unitdomain(A)

Find the dimensional (unit) domain of a matrix

function unitrange(A)

Find the dimensional (unit) range of a matrix