Template out storage types for Toeplitz, Hankel, etc. (#76)

* Refactor

see README.md for details.

* Update README.md

* Update README.md

* rollback changes in runtests.jl

* simplify boundscheck codes. complete interface

* add real and imag, bug fixes

* small fixes

* update

* move codes

* remove Infinities

* Update toeplitz.jl

* Update hankel.jl

* Update hankel.jl

* trying to resolve conflicts

* pass tests

* template hankel dims

* resolve conversations

* 78% coverage

* 90% coverage. forgot to import mul!🤣

* 94% coverage

* 95% coverage

* fix. reimplement size

* let Toeplitz(Matrix) make copies

* Small changes

---------

Co-authored-by: Sheehan Olver <[email protected]>

Project.toml

@@ -1,6 +1,6 @@
 name = "ToeplitzMatrices"
 uuid = "c751599d-da0a-543b-9d20-d0a503d91d24"
-version = "0.7.1"
+version = "0.8"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

README.md

@@ -8,7 +8,7 @@ ToeplitzMatrices.jl
 Fast matrix multiplication and division
 for Toeplitz, Hankel and circulant matrices in Julia
 
-# Note
+## Note
 
 Multiplication of large matrices and `sqrt`, `inv`, `LinearAlgebra.eigvals`,
 `LinearAlgebra.ldiv!`, and `LinearAlgebra.pinv` for circulant matrices
@@ -26,113 +26,137 @@ initialize the required arrays and plan the FFT only once. You can precompute
 the FFT factorization with `LinearAlgebra.factorize` and then use the factorization
 for the FFT-based computations.
 
-# Supported matrices
-
-## Toeplitz
+## Introduction
 
+### Toeplitz
 A Toeplitz matrix has constant diagonals. It can be constructed using
-
 ```julia
 Toeplitz(vc,vr)
+Toeplitz{T}(vc,vr)
 ```
-
 where `vc` are the entries in the first column and `vr` are the entries in the first row, where `vc[1]` must equal `vr[1]`. For example.
-
 ```julia
 Toeplitz(1:3, [1.,4.,5.])
 ```
-
 is a sparse representation of the matrix
-
 ```julia
 [ 1.0  4.0  5.0
   2.0  1.0  4.0
   3.0  2.0  1.0 ]
 ```
-
-## SymmetricToeplitz
-
-A symmetric Toeplitz matrix is a symmetric matrix with constant diagonals. It can be constructed with
-
+### Special toeplitz
+`SymmetricToeplitz`, `Circulant`, `UpperTriangularToeplitz` and `LowerTriangularToeplitz` only store one vector. By convention, `Circulant` stores the first column rather than the first row. They are constructed using `TYPE(v)` where `TYPE`∈{`SymmetricToeplitz`, `Circulant`, `UpperTriangularToeplitz`, `LowerTriangularToeplitz`}.
+
+### Hankel
+A Hankel matrix has constant anti-diagonals, for example,
 ```julia
-SymmetricToeplitz(vc)
+[ 1  2  3
+  2  3  4 ]
 ```
+There are a few ways to construct the above `Hankel{Int}`:
+- `Hankel([1,2,3,4], (2,3)) # Hankel(v, (h,w))`
+- `Hankel([1,2,3,4], 2,3) # Hankel(v, h, w)`
+- `Hankel([1,2], [2,3,4]) # Hankel(vc, vr)`
 
-where `vc` are the entries of the first column. For example,
+Note that the width is usually useless, since ideally, `w=length(v)-h+1`. It exists for infinite Hankel matrices. Its existence also means that `v` could be longer than necessary. `Hankel(v)`, where the size is not given, returns `Hankel(v, (l+1)÷2, (l+1)÷2)` where `l=length(v)`.
 
-```julia
-SymmetricToeplitz([1.0, 2.0, 3.0])
-```
+The `reverse` can transform between Hankel and Toeplitz. It is used to achieve fast Hankel algorithms.
 
-is a sparse representation of the matrix
+## Implemented interface
 
-```julia
-[ 1.0  2.0  3.0
-  2.0  1.0  2.0
-  3.0  2.0  1.0 ]
-```
+### Operations
+- ✓ implemented
+- ✗ error
+- _ fall back to `Matrix`
 
-## TriangularToeplitz
+||Toeplitz|Symmetric~|Circulant|UpperTriangular~|LowerTriangular~|Hankel|
+|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
+|getindex|✓|✓|✓|✓|✓|✓|
+|.vc|✓|✓|✓|✓|✓|✓|
+|.vr|✓|✓|✓|✓|✓|✓|
+|size|✓|✓|✓|✓|✓|✓|
+|copy|✓|✓|✓|✓|✓|✓|
+|similar|✓|✓|✓|✓|✓|✓|
+|zero|✓|✓|✓|✓|✓|✓|
+|real|✓|✓|✓|✓|✓|✓|
+|imag|✓|✓|✓|✓|✓|✓|
+|fill!|✓|✗|✗|✗|✗|✓|
+|conj|✓|✓|✓|✓|✓|✓|
+|transpose|✓|✓|✓|✓|✓|✓|
+|adjoint|✓|✓|✓|✓|✓|✓|
+|tril!|✓|✗|✗|✓|✓|✗|
+|triu!|✓|✗|✗|✓|✓|✗|
+|tril|✓|✓|✓|✓|✓|✗|
+|triu|✓|✓|✓|✓|✓|✗|
+|+|✓|✓|✓|✓|✓|✓|
+|-|✓|✓|✓|✓|✓|✓|
+|scalar<br>mult|✓|✓|✓|✓|✓|✓|
+|==|✓|✓|✓|✓|✓|✓|
+|issymmetric|||||||
+|istriu|||||||
+|istril|||||||
+|iszero|✓|✓|✓|✓|✓||
+|isone|||||||
+|copyto!|✓|✓|✓|✓|✓|✓|
+|reverse|✓|✓|✓|✓|✓|✓|
+|broadcast|||||||
+|broadcast!|||||||
 
-A triangular Toeplitz matrix can be constructed using
+Note that scalar multiplication, `conj`, `+` and `-` could be removed once `broadcast` is implemented.
 
-```julia
-TriangularToeplitz(ve,uplo)
-```
+`reverse(Hankel)` returns a `Toeplitz`, while `reverse(AbstractToeplitz)` returns a `Hankel`.
 
-where uplo is either `:L` or `:U` and `ve` are the rows or columns, respectively.  For example,
+### LinearAlgebra
 
-```julia
-TriangularToeplitz([1.,2.,3.],:L)
-```
+### Constructors and conversions
+||Toeplitz|Symmetric~|Circulant|UpperTriangular~|LowerTriangular~|Hankel|
+|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
+|from AbstractVector|✓|✓|✓|✓|✓|✓|
+|from AbstractMatrix|✓|✓|✓|✓|✓|✓|
+|from AbstractToeplitz|✓|✓|✓|✓|✓|✗|
+|to supertype|✓|✓|✓|✓|✓|✓|
+|to Toeplitz|-|✓|✓|✓|✓|✗|
+|to another eltype|✓|✓|✓|✓|✓|✓|
 
-is a sparse representation of the matrix
+When constructing `Toeplitz` from a matrix, the first row and the first column will be considered as `vr` and `vc`. Note that `vr` and `vc` are copied in construction to avoid the cases where they share memory. If you don't want copying, construct using vectors directly.
+
+When constructing `SymmetricToeplitz` or `Circulant` from `AbstractMatrix`, a second argument shall specify whether the first row or the first column is used. For example, for `A = [1 2; 3 4]`, 
+- `SymmetricToeplitz(A,:L)` gives `[1 3; 3 1]`, while
+- `SymmetricToeplitz(A,:U)` gives `[1 2; 2 1]`.
 
+For backward compatibility and consistency with `LinearAlgebra.Symmetric`,
 ```julia
-[ 1.0  0.0  0.0
-  2.0  1.0  0.0
-  3.0  2.0  1.0 ]
+SymmetricToeplitz(A) = SymmetricToeplitz(A, :U)
+Circulant(A) = Circulant(A, :L)
 ```
+`Hankel` constructor also accepts the second argument, `:L` denoting the first column and the last row while `:U` denoting the first row and the last column.
 
-## Hankel
-
-A Hankel matrix has constant anti-diagonals.  It can be constructed using
-
+`Symmetric`, `UpperTriangular` and `LowerTriangular` from `LinearAlgebra` are also overloaded for convenience.
 ```julia
-Hankel(vc,vr)
+Symmetric(T::Toeplitz) = SymmetricToeplitz(T)
+UpperTriangular(T::Toeplitz) = UpperTriangularToeplitz(T)
+LowerTriangular(T::Toeplitz) = LowerTriangularToeplitz(T)
 ```
 
-where `vc` are the entries in the first column and `vr` are the entries in the last row, where `vc[end]` must equal `vr[1]`.  For example.
-
+### TriangularToeplitz (obsolete)
+`TriangularToeplitz` is reserved for backward compatibility. 
 ```julia
-Hankel([1.,2.,3.], 3:5)
+TriangularToeplitz = Union{UpperTriangularToeplitz,LowerTriangularToeplitz}
 ```
-
-is a sparse representation of the matrix
-
+The old interface is implemented by
 ```julia
-[  1.0  2.0  3.0
-   2.0  3.0  4.0
-   3.0  4.0  5.0 ]
+getproperty(UpperTriangularToeplitz,:uplo) = :U
+getproperty(LowerTriangularToeplitz,:uplo) = :L
 ```
+This type is **obsolete** and will not be updated for features. Despite that, backward compatibility should be maintained. Codes that were using `TriangularToeplitz` should still work.
 
-## Circulant
+## Unexported interface
+Methods in this section are not exported.
 
-A circulant matrix is a special case of a Toeplitz matrix with periodic end conditions.
-It can be constructed using
+`_vr(A::AbstractMatrix)` returns the first row as a vector.
+`_vc(A::AbstractMatrix)` returns the first column as a vector.
+`_vr` and `_vc` are implemented for `AbstractToeplitz` as well. They are used to merge similar codes for `AbstractMatrix` and `AbstractToeplitz`.
 
-```julia
-Circulant(vc)
-```
-where `vc` is a vector with the entries for the first column.
-For example:
-```julia
-Circulant([1.0, 2.0, 3.0])
-```
-is a sparse representation of the matrix
+`_circulate(v::AbstractVector)` converts between the `vr` and `vc` of a `Circulant`.
 
-```julia
-[  1.0  3.0  2.0
-   2.0  1.0  3.0
-   3.0  2.0  1.0 ]
-```
+`isconcrete(A::Union{AbstractToeplitz,Hankel})` decides whether the stored vector(s) are concrete. It calls `Base.isconcretetype`.