From ca52a1e76ccf2c31e739f2dcd62f92727c78a50f Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Sun, 12 Oct 2025 06:44:09 +0530 Subject: [PATCH 1/4] Tutorial update --- docs/src/tutorial.md | 69 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 68 insertions(+), 1 deletion(-) diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 876cd0a0..2f825a88 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -1 +1,68 @@ -You can check an [Introduction to LinearOperators.jl](https://jso.dev/tutorials/introduction-to-linear-operators/) on our site, [jso.dev](https://jso.dev/). +# Tutorial + +This page follows the Documenter.jl tutorial layout. It includes runnable examples for common usage patterns in LinearOperators.jl. + +```@contents +Pages = ["tutorial.md"] +``` + +## Getting started + +Calling into LinearOperators is simple: + +```@docs +using LinearOperators +``` + +## Simple examples + +Construct an operator from a matrix and use it like a matrix: + +```@example ex_matrix +using LinearOperators + +A = [1.0 2.0; 3.0 4.0] +op = LinearOperator(A) +y = op * [1.0, 1.0] +M = Matrix(op) +println("y = ", y) +println("Matrix(op) = \n", M) +``` + +Create function-based operators (showing the 5-arg mul! signature): + +```@example ex_fun +n = 4 +function mymul!(res, v, α, β) + if β == 0 + res .= α .* v .* (1:n) + else + res .= α .* v .* (1:n) .+ β .* res + end +end +opfun = LinearOperator(Float64, n, n, false, false, mymul!) +println(opfun * ones(n)) +``` + +Diagonal operators (use `opDiagonal`): + +```@example ex_diag +d = [2.0, 3.0, 4.0] +D = opDiagonal(d) +println(D * ones(3)) +``` + +Composing operators with vertical concatenation: + +```@example ex_cat +A = rand(3,3); B = rand(3,3) +opA = LinearOperator(A); opB = LinearOperator(B) +opcat = [opA; opB] +println(size(opcat)) +``` + +## Return values and tips + +- The `LinearOperator` type implements `*` and can be converted to a `Matrix` when necessary. +- Prefer function-based operators when you want to avoid materializing large matrices. +- See the full introduction: [Introduction to LinearOperators.jl](https://jso.dev/tutorials/introduction-to-linear-operators/) From ce862061403122afc72542cddd9b678952ebdef7 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Sun, 19 Oct 2025 02:53:45 +0530 Subject: [PATCH 2/4] Storage-type promotion --- src/cat.jl | 22 ++++++++++++++++++---- src/operations.jl | 22 ++++++++++++++++++---- 2 files changed, 36 insertions(+), 8 deletions(-) diff --git a/src/cat.jl b/src/cat.jl index 764c7893..7fc4bd5d 100644 --- a/src/cat.jl +++ b/src/cat.jl @@ -46,8 +46,15 @@ function hcat(A::AbstractLinearOperator, B::AbstractLinearOperator) hcat_ctprod!(res, adjoint(A), adjoint(B), Ancol, Ancol + Bncol, w, α, β) args5 = (has_args5(A) && has_args5(B)) S = promote_type(storage_type(A), storage_type(B)) - isconcretetype(S) || - throw(LinearOperatorException("storage types cannot be promoted to a concrete type")) + if !isconcretetype(S) + if isconcretetype(storage_type(A)) + S = storage_type(A) + elseif isconcretetype(storage_type(B)) + S = storage_type(B) + else + S = Vector{T} + end + end CompositeLinearOperator(T, nrow, ncol, false, false, prod!, tprod!, ctprod!, args5, S = S) end @@ -105,8 +112,15 @@ function vcat(A::AbstractLinearOperator, B::AbstractLinearOperator) vcat_ctprod!(res, adjoint(A), adjoint(B), Anrow, Anrow + Bnrow, w, α, β) args5 = (has_args5(A) && has_args5(B)) S = promote_type(storage_type(A), storage_type(B)) - isconcretetype(S) || - throw(LinearOperatorException("storage types cannot be promoted to a concrete type")) + if !isconcretetype(S) + if isconcretetype(storage_type(A)) + S = storage_type(A) + elseif isconcretetype(storage_type(B)) + S = storage_type(B) + else + S = Vector{T} + end + end CompositeLinearOperator(T, nrow, ncol, false, false, prod!, tprod!, ctprod!, args5, S = S) end diff --git a/src/operations.jl b/src/operations.jl index 18259d40..a65e3407 100644 --- a/src/operations.jl +++ b/src/operations.jl @@ -140,8 +140,15 @@ function *(op1::AbstractLinearOperator, op2::AbstractLinearOperator) throw(LinearOperatorException("shape mismatch")) end S = promote_type(storage_type(op1), storage_type(op2)) - isconcretetype(S) || - throw(LinearOperatorException("storage types cannot be promoted to a concrete type")) + if !isconcretetype(S) + if isconcretetype(storage_type(op1)) + S = storage_type(op1) + elseif isconcretetype(storage_type(op2)) + S = storage_type(op2) + else + S = Vector{T} + end + end #tmp vector for products vtmp = fill!(S(undef, m2), zero(T)) utmp = fill!(S(undef, n1), zero(T)) @@ -211,8 +218,15 @@ function +(op1::AbstractLinearOperator, op2::AbstractLinearOperator) herm = (ishermitian(op1) && ishermitian(op2)) args5 = (has_args5(op1) && has_args5(op2)) S = promote_type(storage_type(op1), storage_type(op2)) - isconcretetype(S) || - throw(LinearOperatorException("storage types cannot be promoted to a concrete type")) + if !isconcretetype(S) + if isconcretetype(storage_type(op1)) + S = storage_type(op1) + elseif isconcretetype(storage_type(op2)) + S = storage_type(op2) + else + S = Vector{T} + end + end return CompositeLinearOperator(T, m1, n1, symm, herm, prod!, tprod!, ctprod!, args5, S = S) end From a4f8f70f15f676649a544cf7162b624d6ba8c7a6 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Sun, 19 Oct 2025 02:57:02 +0530 Subject: [PATCH 3/4] Update tutorial.md --- docs/src/tutorial.md | 69 +------------------------------------------- 1 file changed, 1 insertion(+), 68 deletions(-) diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 2f825a88..876cd0a0 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -1,68 +1 @@ -# Tutorial - -This page follows the Documenter.jl tutorial layout. It includes runnable examples for common usage patterns in LinearOperators.jl. - -```@contents -Pages = ["tutorial.md"] -``` - -## Getting started - -Calling into LinearOperators is simple: - -```@docs -using LinearOperators -``` - -## Simple examples - -Construct an operator from a matrix and use it like a matrix: - -```@example ex_matrix -using LinearOperators - -A = [1.0 2.0; 3.0 4.0] -op = LinearOperator(A) -y = op * [1.0, 1.0] -M = Matrix(op) -println("y = ", y) -println("Matrix(op) = \n", M) -``` - -Create function-based operators (showing the 5-arg mul! signature): - -```@example ex_fun -n = 4 -function mymul!(res, v, α, β) - if β == 0 - res .= α .* v .* (1:n) - else - res .= α .* v .* (1:n) .+ β .* res - end -end -opfun = LinearOperator(Float64, n, n, false, false, mymul!) -println(opfun * ones(n)) -``` - -Diagonal operators (use `opDiagonal`): - -```@example ex_diag -d = [2.0, 3.0, 4.0] -D = opDiagonal(d) -println(D * ones(3)) -``` - -Composing operators with vertical concatenation: - -```@example ex_cat -A = rand(3,3); B = rand(3,3) -opA = LinearOperator(A); opB = LinearOperator(B) -opcat = [opA; opB] -println(size(opcat)) -``` - -## Return values and tips - -- The `LinearOperator` type implements `*` and can be converted to a `Matrix` when necessary. -- Prefer function-based operators when you want to avoid materializing large matrices. -- See the full introduction: [Introduction to LinearOperators.jl](https://jso.dev/tutorials/introduction-to-linear-operators/) +You can check an [Introduction to LinearOperators.jl](https://jso.dev/tutorials/introduction-to-linear-operators/) on our site, [jso.dev](https://jso.dev/). From 3dfa5ff415aa8d745f972ed3e5a43b2c3de06e73 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Wed, 22 Oct 2025 04:48:59 +0530 Subject: [PATCH 4/4] Centralized storage selection is now used by both operations.jl and cat.jl consistently (tests still pass functionally; only the known two tiny allocation assertions remain). --- src/abstract.jl | 17 +++++++++++++++++ src/cat.jl | 22 ++-------------------- src/operations.jl | 22 ++-------------------- 3 files changed, 21 insertions(+), 40 deletions(-) diff --git a/src/abstract.jl b/src/abstract.jl index c3940a33..2f33b0f7 100644 --- a/src/abstract.jl +++ b/src/abstract.jl @@ -179,6 +179,23 @@ storage_type(op::Adjoint) = storage_type(parent(op)) storage_type(op::Transpose) = storage_type(parent(op)) storage_type(op::Diagonal) = typeof(parent(op)) +@inline function _select_storage_type( + op1::AbstractLinearOperator, + op2::AbstractLinearOperator, + T::Type, +) + S = promote_type(storage_type(op1), storage_type(op2)) + if isconcretetype(S) + return S + elseif isconcretetype(storage_type(op1)) + return storage_type(op1) + elseif isconcretetype(storage_type(op2)) + return storage_type(op2) + else + return Vector{T} + end +end + """ reset!(op) diff --git a/src/cat.jl b/src/cat.jl index 7fc4bd5d..b20e1be4 100644 --- a/src/cat.jl +++ b/src/cat.jl @@ -45,16 +45,7 @@ function hcat(A::AbstractLinearOperator, B::AbstractLinearOperator) ctprod! = @closure (res, w, α, β) -> hcat_ctprod!(res, adjoint(A), adjoint(B), Ancol, Ancol + Bncol, w, α, β) args5 = (has_args5(A) && has_args5(B)) - S = promote_type(storage_type(A), storage_type(B)) - if !isconcretetype(S) - if isconcretetype(storage_type(A)) - S = storage_type(A) - elseif isconcretetype(storage_type(B)) - S = storage_type(B) - else - S = Vector{T} - end - end + S = _select_storage_type(A, B, T) CompositeLinearOperator(T, nrow, ncol, false, false, prod!, tprod!, ctprod!, args5, S = S) end @@ -111,16 +102,7 @@ function vcat(A::AbstractLinearOperator, B::AbstractLinearOperator) ctprod! = @closure (res, w, α, β) -> vcat_ctprod!(res, adjoint(A), adjoint(B), Anrow, Anrow + Bnrow, w, α, β) args5 = (has_args5(A) && has_args5(B)) - S = promote_type(storage_type(A), storage_type(B)) - if !isconcretetype(S) - if isconcretetype(storage_type(A)) - S = storage_type(A) - elseif isconcretetype(storage_type(B)) - S = storage_type(B) - else - S = Vector{T} - end - end + S = _select_storage_type(A, B, T) CompositeLinearOperator(T, nrow, ncol, false, false, prod!, tprod!, ctprod!, args5, S = S) end diff --git a/src/operations.jl b/src/operations.jl index a65e3407..b02803d3 100644 --- a/src/operations.jl +++ b/src/operations.jl @@ -139,16 +139,7 @@ function *(op1::AbstractLinearOperator, op2::AbstractLinearOperator) if m2 != n1 throw(LinearOperatorException("shape mismatch")) end - S = promote_type(storage_type(op1), storage_type(op2)) - if !isconcretetype(S) - if isconcretetype(storage_type(op1)) - S = storage_type(op1) - elseif isconcretetype(storage_type(op2)) - S = storage_type(op2) - else - S = Vector{T} - end - end + S = _select_storage_type(op1, op2, T) #tmp vector for products vtmp = fill!(S(undef, m2), zero(T)) utmp = fill!(S(undef, n1), zero(T)) @@ -217,16 +208,7 @@ function +(op1::AbstractLinearOperator, op2::AbstractLinearOperator) symm = (issymmetric(op1) && issymmetric(op2)) herm = (ishermitian(op1) && ishermitian(op2)) args5 = (has_args5(op1) && has_args5(op2)) - S = promote_type(storage_type(op1), storage_type(op2)) - if !isconcretetype(S) - if isconcretetype(storage_type(op1)) - S = storage_type(op1) - elseif isconcretetype(storage_type(op2)) - S = storage_type(op2) - else - S = Vector{T} - end - end + S = _select_storage_type(op1, op2, T) return CompositeLinearOperator(T, m1, n1, symm, herm, prod!, tprod!, ctprod!, args5, S = S) end