Update README and add NEWS

.gitignore

@@ -1 +1,10 @@
 *~
+*.jl.cov
+*.jl.*.cov
+*.jl.mem
+
+docs/build/
+docs/site/
+
+Manifest.toml
+

README.md

@@ -9,66 +9,40 @@ Type `QuasiUpperTriangular` stores a quasi upper triangular matrix in
 a square matrix. The various algorithms ignore the lower zero elements
 
 ## Functions
+Conventional `mul!` functions are defined to allow normal multiplication using `*`
 
-### Matrix - vector product
-- `A_mul_B!(a::QuasiUpperTriangular,b::AbstractVector,work::AbstractVector)`
-uses BLAS `trmv` and corrects the results for lower nonzero
-elements. For small matrices, (n < 27), uses regular matrix - vector
-product
-- `A_mul_B!(A::QuasiUpperTriangular, B::AbstractVector)` pure
-    Julia implementation. Doesn't need extra workspace.
-### Vector - matrix produc
-- `A_mul_B!(A::AbstractVector, B::QuasiUpperTriangular)` pure
-    Julia implementation. Doesn't need extra workspace.
-
-### Transposed matrix -vector product
-- `At_mul_B!(A::QuasiUpperTriangular, B::AbstractVector)` pure
-    Julia implementation. Doesn't need extra workspace.
-
-### Vector - transposed matrix product
-- `A_mul_Bt!(A::AbstractVector, B::QuasiUpperTriangular)`  pure
-Julia implementation. Doesn't need extra workspace.
-
-### Matrix - matrix product
-- `A_mul_B!(C::AbstractMatrix, A::QuasiUpperTriangular,
-      B::AbstractMatrix)`: `C = A*B`, `A` is quasi upper triangular
-- `A_mul_B!(C::AbstractMatrix, A::QuasiUpperTriangular,
-	    b::AbstractMatrix)`: `C = αA*B`, `A` is quasi upper triangular
-- `At_mul_B!(C::AbstractMatrix, A::QuasiUpperTriangular,
-          B::AbstractMatrix)`: `C = transpose(A)*B`, `A` is quasi upper triangular
-- `At_mul_B!(C::AbstractMatrix, A::QuasiUpperTriangular,
-	    b::AbstractMatrix)`: `C = α*transpose(A)*B`, `A` is quasi upper triangular
-- `At_mul_B!(A::QuasiUpperTriangular, B::AbstractMatrix)` `B = transpose(A)*B` in place computation,
-    `A` is quasi upper triangular.
-- `A_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::QuasiUpperTriangular)`:
-    `C = A*B`, `B` is quasi upper triangular.
-- `A_mul_B!(A::AbstractMatrix, B::QuasiUpperTriangular)`:
-    `C = A*B`, in place computation, `B` is quasi upper triangular.
-- `A_mul_Bt!(C::AbstractMatrix, A::AbstractMatrix, B::QuasiUpperTriangular)`:
-    `C = A*transpose(B)`, `B` is quasi upper triangular.
-- `A_mul_Bt!(A::AbstractMatrix, B::QuasiUpperTriangular)`:
-    `C = A*transpose(B)`, in place computation, `B` is quasi upper triangular.
-- `A_mul_B!(C:QuasiUpperTriangular, A::QuasiUpperTriangular, B::QuasiUpperTriangular)`,
-    `A`, `B` and `C` are quasi upper triangular.
+lets define $Q$ as a QuasiUpperTriangular matrix.
 
-### Linear problem solver
-- `A_ldiv_B!(A::QuasiUpperTriangular, B::AbstractMatrix)` solves `A*X = B`,
-    where `A` is quasi upper triangular. Solves by back substitution. Lower off-diagonal elements
-    make 2 * 2 problems that are solved explicitely.
-- `A_rdiv_B!(A::AbstractMatrix, B::QuasiUpperTriangular)` solves `X*B = A`,
-    where `B` is quasi upper triangular. Solves by back substitution. Lower off-diagonal elements
-    make 2 * 2 problems that are solved explicitely.
-- `A_rdiv_Bt!(A::AbstractMatrix, B::QuasiUpperTriangular)` solves `X*transpose(B) = A`,
-    where `B` is quasi upper triangular. Solves by back substitution. Lower off-diagonal elements
-    make 2 * 2 problems that are solved explicitely.
-- `I_plus_rA_ldiv_B!(r::Float64,a::QuasiUpperTriangular, b::AbstractVector)` solves
-    `(I + r*A)*x = b` where `A` is quasi upper triangular.
+### Matrix - vector product
+- $Q * \vec{v}$
+- $Q^T * \vec{v}$
+
+### Matrix - Matrix product
+- $Q * {A}$
+- $Q^T * A$
+- $A * Q$
+- $A * Q^T$
+
+
+### Linear problem solvers
+- `ldiv!(Q, A)` solves $Q*X = A$,
+    Solves by back substitution. Lower off-diagonal elements make 2 * 2 problems that are solved explicitly.
+- `rdiv!(A, Q)` solves $X*Q = A$,
+    Solves by back substitution. Lower off-diagonal elements make 2 * 2 problems that are solved explicitly.
+- `rdiv!(A, Q')` solves $X*Q^T = A$,
+    Solves by back substitution. Lower off-diagonal elements make 2 * 2 problems that are solved explicitly.
+
+lets define $r$ and $s$ as floats
+
+> Note: these functions break conventions and mutate their last argument
+- `I_plus_rA_ldiv_B!(r, Q, b)` solves $(I + rQ)*\vec{x} = \vec{b}$
+- `I_plus_rA_ldiv_B!(r, Q, B)` solves $(I + rQ)*X = B$ 
+- `I_plus_rA_plus_sB_ldiv_C!(r, s, Q1, Q2, c)` solves $(I + rQ_1 + sQ_2)*\vec{x} = \vec{c}$
+- `I_plus_rA_plus_sB_ldiv_C!(r, s, Q1, Q2, C)` solves $(I + rQ_1 + sQ_2)*X = C$
 
 ## TODO
-- replace function names for product by `mul!`
-- introduce lazy transpose evaluation
-- handle quasi lower triangular matrices
-- benchmark cases that have two different implementations
-
-
-
+- [x] replace function names for product by `mul!`
+- [x] introduce lazy transpose evaluation
+- [ ] assert that sub-diagonal does not contain consecutive non-zero elements 
+- [ ] handle quasi lower triangular matrices
+- [ ] profile, benchmark, and reintroduce BLAS based implementations if needed (for specific strided-matrix element-types)

src/NEWS.md

@@ -0,0 +1,5 @@
+0.2.0
+======
+- move most functions that operate on submatrices to KroneckerTools package
+- rename and refactor all `mul!` functions to be consistent with Julia's post v1.0 conventions
+- delete all implementations that use BLAS to keep the core package generic

src/QuasiTriangular.jl

@@ -1,32 +1,26 @@
 module QuasiTriangular
-## QuasiUpperTriangular Matrix of real numbers
 
 using LinearAlgebra
+
+import Base.copy
 import Base.size
 import Base.similar
 import Base.getindex
-import Base.require_one_based_indexing
 import Base.setindex!
-import Base.copy
 import LinearAlgebra.mul!
 import LinearAlgebra.ldiv!
 import LinearAlgebra.rdiv!
-import LinearAlgebra.checksquare
-import LinearAlgebra.BlasInt
-import LinearAlgebra.BLAS.@blasfunc
-import LinearAlgebra.BLAS.libblas
 
-export QuasiUpperTriangular, I_plus_rA_ldiv_B!, I_plus_rA_plus_sB_ldiv_C!,
- rdiv!, mul!, ldiv!, rdiv!
+export QuasiUpperTriangular, mul!, ldiv!, rdiv!, I_plus_rA_ldiv_B!, I_plus_rA_plus_sB_ldiv_C!
 
 abstract type AbstractQuasiTriangular{T, S <: AbstractMatrix} <: AbstractMatrix{T} end
 
 struct QuasiUpperTriangular{T, S <: AbstractMatrix{T}} <: AbstractQuasiTriangular{T,S}
     data::S
 
     function QuasiUpperTriangular{T,S}(data) where {T,S<:AbstractMatrix{T}}
-        require_one_based_indexing(data)
-        checksquare(data)
+        Base.require_one_based_indexing(data)
+        LinearAlgebra.checksquare(data)
         new(data)
     end
 end