From b49872e2e1259c69420ba90d4aa1dadf1d3ff543 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Mon, 27 Mar 2023 10:41:32 -0400 Subject: [PATCH 1/2] update Readme --- .gitignore | 9 +++++ README.md | 92 +++++++++++++++--------------------------- src/QuasiTriangular.jl | 16 +++----- 3 files changed, 47 insertions(+), 70 deletions(-) diff --git a/.gitignore b/.gitignore index b25c15b..2469202 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,10 @@ *~ +*.jl.cov +*.jl.*.cov +*.jl.mem + +docs/build/ +docs/site/ + +Manifest.toml + diff --git a/README.md b/README.md index d719039..0b89671 100644 --- a/README.md +++ b/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) diff --git a/src/QuasiTriangular.jl b/src/QuasiTriangular.jl index 45439dd..623a770 100644 --- a/src/QuasiTriangular.jl +++ b/src/QuasiTriangular.jl @@ -1,23 +1,17 @@ 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 @@ -25,8 +19,8 @@ struct QuasiUpperTriangular{T, S <: AbstractMatrix{T}} <: AbstractQuasiTriangula 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 From e1c6c850e5554fccbd8d40ba19fd70f4bb19ded1 Mon Sep 17 00:00:00 2001 From: Omar Elrefaei Date: Mon, 27 Mar 2023 10:52:52 -0400 Subject: [PATCH 2/2] add NEWS.md --- src/NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 src/NEWS.md diff --git a/src/NEWS.md b/src/NEWS.md new file mode 100644 index 0000000..8447752 --- /dev/null +++ b/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 \ No newline at end of file