From 19161f987ab7ea2e4e30d0dd1b0936dcd879d574 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sat, 24 Aug 2019 17:21:27 -0400 Subject: [PATCH 01/15] Fix index in BlockMap adjoint/transpose multiply --- src/blockmap.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/blockmap.jl b/src/blockmap.jl index bbe5e889..8e9deda4 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -244,19 +244,19 @@ function At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!")) maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges mapind = 0 - # first block row (rowind = 1), fill all of y + # first block row (rowind = 1) of A, meaning first block column of A', fill all of y @views @inbounds begin xcol = x[xinds[1]] for colind in 1:rows[1] mapind +=1 - A_mul_B!(y[yinds[colind]], transpose(maps[mapind]), xcol) + A_mul_B!(y[yinds[mapind]], transpose(maps[mapind]), xcol) end - # subsequent block rows, add results to corresponding parts of y + # subsequent block rows of A, add results to corresponding parts of y for rowind in 2:length(rows) xcol = x[xinds[rowind]] for colind in 1:rows[rowind] mapind +=1 - mul!(y[yinds[colind]], transpose(maps[mapind]), xcol, true, true) + mul!(y[yinds[mapind]], transpose(maps[mapind]), xcol, true, true) end end end @@ -268,19 +268,19 @@ function Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!")) maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges mapind = 0 - # first block row (rowind = 1), fill all of y + # first block row (rowind = 1) of A, fill all of y @views @inbounds begin xcol = x[xinds[1]] for colind in 1:rows[1] mapind +=1 - A_mul_B!(y[yinds[colind]], adjoint(maps[mapind]), xcol) + A_mul_B!(y[yinds[mapind]], adjoint(maps[mapind]), xcol) end - # subsequent block rows, add results to corresponding parts of y + # subsequent block rows of A, add results to corresponding parts of y for rowind in 2:length(rows) xcol = x[xinds[rowind]] for colind in 1:rows[rowind] mapind +=1 - mul!(y[yinds[colind]], adjoint(maps[mapind]), xcol, true, true) + mul!(y[yinds[mapind]], adjoint(maps[mapind]), xcol, true, true) end end end From 65738fa40ac94cfdb509e642abb9ed52248322a8 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 25 Aug 2019 17:18:20 -0400 Subject: [PATCH 02/15] Add hvcat test for "misaligned" blockmap --- test/runtests.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 3289b921..9936a22f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -509,6 +509,7 @@ end @test_throws DimensionMismatch vcat(LinearMap(A11), LinearMap(A21)) end end + @testset "hvcat" begin for elty in (Float32, Float64, ComplexF64) A11 = rand(elty, 10, 10) @@ -546,8 +547,19 @@ end @test_throws DimensionMismatch A = [I A21; A12 I] @test_throws DimensionMismatch A = [A12 A12; A21 A21] @test_throws DimensionMismatch A = [A12 A21; A12 A21] + + # basic test of "misaligned" blocks + M = ones(elty, 3, 2) # non-square + A = LinearMap(M) + B = [I A; A I] + C = [I M; M I] + @test B isa LinearMaps.BlockMap{elty} + @test Matrix(B) == C + @test Matrix(transpose(B)) == C' + @test Matrix(adjoint(B)) == C' end end + @testset "adjoint/transpose" begin for elty in (Float32, Float64, ComplexF64), transform in (transpose, adjoint) A12 = rand(elty, 10, 10) From 1f2501b888a9bacefc290f93f0335d0d11a0248d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 26 Aug 2019 09:09:13 +0200 Subject: [PATCH 03/15] minor code style changes --- test/runtests.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9936a22f..ceb40f84 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -509,7 +509,6 @@ end @test_throws DimensionMismatch vcat(LinearMap(A11), LinearMap(A21)) end end - @testset "hvcat" begin for elty in (Float32, Float64, ComplexF64) A11 = rand(elty, 10, 10) @@ -549,17 +548,16 @@ end @test_throws DimensionMismatch A = [A12 A21; A12 A21] # basic test of "misaligned" blocks - M = ones(elty, 3, 2) # non-square + M = ones(elty, 3, 2) # non-square A = LinearMap(M) B = [I A; A I] C = [I M; M I] @test B isa LinearMaps.BlockMap{elty} @test Matrix(B) == C - @test Matrix(transpose(B)) == C' + @test Matrix(transpose(B)) == transpose(C) @test Matrix(adjoint(B)) == C' end end - @testset "adjoint/transpose" begin for elty in (Float32, Float64, ComplexF64), transform in (transpose, adjoint) A12 = rand(elty, 10, 10) From eb1d1e373f8491168a4a8de2e9d02f6431c6c1da Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sat, 31 Aug 2019 14:35:16 -0400 Subject: [PATCH 04/15] Initial version of IndexMap --- src/LinearMaps.jl | 1 + src/indexmap.jl | 138 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 src/indexmap.jl diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index c3c6f813..d556223f 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -112,6 +112,7 @@ include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby al include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I include("functionmap.jl") # using a function as linear map include("blockmap.jl") # block linear maps +include("indexmap.jl") # single block map defined by row/col index """ LinearMap(A; kwargs...) diff --git a/src/indexmap.jl b/src/indexmap.jl new file mode 100644 index 00000000..00fa38db --- /dev/null +++ b/src/indexmap.jl @@ -0,0 +1,138 @@ +#= +indexmap.jl +2019-08-29 Jeff Fessler, University of Michigan +=# + +export IndexMap, hcat_new + +#using LinearMaps +#import LinearMaps: A_mul_B! +#import LinearAlgebra + + +""" +`check_index(row::AbstractVector{Int}, dims::Dims{2})` + +Do not allow repeated index values because that would complicate the adjoint. +And insist on monotone increasing too because it is the easiest way to +ensure no repeats and probably better for cache too. +""" +function check_index(index::AbstractVector{Int}, dimA::Int, dimB::Int) + length(index) != dimA && throw("invalid index length") + minimum(index) < 0 && throw("negative index") + maximum(index) > dimB && + throw("index $(maximum(index)) > dimB=$dimB") + !(index isa UnitRange{Int64} || + ((index isa StepRange{Int64,Int64}) && (index.step > 0)) || + all(diff(index) .> 0)) && throw("non-monotone index") + nothing +end + + +""" +`IndexMap` + +The original motivation was to support construction of linear operator +of the form `B = [0 0 0; 0 A 0; 0 0 0]` +from a given linear operator `A` +where `size(B) == dims` + +This version is more general than the form above +because it allows the rows and colums of `A` +to be interspersed among the rows and columns of `B` arbitrarily. +Essentially like the above form but with some permutation-like matrices +before and after. +""" +struct IndexMap{T, As <: LinearMap, Rs <: AbstractVector{Int}, + Cs <: AbstractVector{Int}} <: LinearMap{T} + map::As + dims::Dims{2} + rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1) + cols::Cs # typically j1:j2 + function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, rows::Rs, cols::Cs) where + {T, As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} + + check_index(rows, size(map,1), dims[1]) + check_index(cols, size(map,2), dims[2]) + + return new{T,As,Rs,Cs}(map, dims, rows, cols) + end +end + +IndexMap{T}(map::As, dims::Dims{2} ; offset::Dims{2}) where {T, As <: LinearMap} = + IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, + offset[1] .+ (1:size(map,1)), + offset[2] .+ (1:size(map,2))) + +IndexMap(map::LinearMap, dims::Dims{2} ; offset::Dims{2}) = + IndexMap{eltype(map)}(map, dims ; offset=offset) + +Base.size(A::IndexMap) = A.dims + + +############ +# basic methods +############ + +# the following rules are sufficient but perhaps not necessary +LinearAlgebra.issymmetric(A::IndexMap) = issymmetric(A.map) && + (A.dims[1] == A.dims[2]) && (A.rows == A.cols) + +LinearAlgebra.ishermitian(A::IndexMap) = ishermitian(A.map) && + (A.dims[1] == A.dims[2]) && (A.rows == A.cols) + +LinearAlgebra.ishermitian(A::IndexMap{<:Real}) = issymmetric(A) + + +# comparison of IndexMap objects +Base.:(==)(A::IndexMap, B::IndexMap) = (eltype(A) == eltype(B)) && + (A.map == B.map) && (A.dims == B.dims) && + (A.rows == B.rows) && (A.cols == B.cols) + +# special transposition behavior +LinearAlgebra.transpose(A::IndexMap) = TransposeMap(A) +LinearAlgebra.adjoint(A::IndexMap) = AdjointMap(A) + + +############ +# multiplication with vectors +############ + +function A_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) + m, n = size(A) + @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) + y[:] .= 0 + @views @inbounds mul!(y[A.rows], A.map, x[A.cols], true, true) + return y +end + +function At_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) + m, n = size(A) + @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!")) + y[:] .= 0 + @views @inbounds mul!(y[A.cols], transpose(A.map), x[A.rows], true, true) + return y +end + +function Ac_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) + m, n = size(A) + y[:] .= 0 + @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("Ac_mul_B!")) + @views @inbounds mul!(y[A.cols], adjoint(A.map), x[A.rows], true, true) + return y +end + +#= +=# +"simplified version of hcat" +function hcat_new(As::LinearMap...) + rows = collect(map(A -> size(A,1), As)) + any(rows .!= rows[1]) && throw("row mismatch") + cols = collect(map(A -> size(A,2), As)) + + nmap = length(As) + col_offsets = [0; cumsum(cols)[1:nmap-1]] + dims = (rows[1], sum(cols)) + return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii])) + for ii in 1:nmap]...) +end From a87f2314827e0ae99f5b2da272cdca37f5935765 Mon Sep 17 00:00:00 2001 From: "fessler@umich.edu" Date: Tue, 3 Sep 2019 14:48:59 -0400 Subject: [PATCH 05/15] Add block_diag --- src/indexmap.jl | 90 +++++++++++++++++++++++++++++-------------------- 1 file changed, 54 insertions(+), 36 deletions(-) diff --git a/src/indexmap.jl b/src/indexmap.jl index 00fa38db..3e006f95 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -3,12 +3,7 @@ indexmap.jl 2019-08-29 Jeff Fessler, University of Michigan =# -export IndexMap, hcat_new - -#using LinearMaps -#import LinearMaps: A_mul_B! -#import LinearAlgebra - +export IndexMap, hcat_new, block_diag """ `check_index(row::AbstractVector{Int}, dims::Dims{2})` @@ -47,25 +42,35 @@ struct IndexMap{T, As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} <: LinearMap{T} map::As dims::Dims{2} + set0::Bool # set to zero the output array? default true rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1) cols::Cs # typically j1:j2 - function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, rows::Rs, cols::Cs) where - {T, As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} + function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, set0::Bool, + rows::Rs, cols::Cs) where + {T, As <: LinearMap, + Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} check_index(rows, size(map,1), dims[1]) check_index(cols, size(map,2), dims[2]) - return new{T,As,Rs,Cs}(map, dims, rows, cols) + return new{T,As,Rs,Cs}(map, dims, set0, rows, cols) end end -IndexMap{T}(map::As, dims::Dims{2} ; offset::Dims{2}) where {T, As <: LinearMap} = - IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, +IndexMap{T}(map::As, dims::Dims{2} ; set0::Bool = true, offset::Dims{2}) where + {T, As <: LinearMap} = + IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, set0, offset[1] .+ (1:size(map,1)), offset[2] .+ (1:size(map,2))) -IndexMap(map::LinearMap, dims::Dims{2} ; offset::Dims{2}) = - IndexMap{eltype(map)}(map, dims ; offset=offset) +IndexMap(map::LinearMap, dims::Dims{2} ; set0::Bool = true, offset::Dims{2}) = + IndexMap{eltype(map)}(map, dims ; offset=offset, set0=set0) + +IndexMap(map::LinearMap, dims::Dims{2}, + rows::AbstractVector{Int}, cols::AbstractVector{Int}, + ; set0::Bool = true) = + IndexMap{eltype(map), typeof(map), typeof(rows), typeof(cols)}( + map, dims, set0, rows, cols) Base.size(A::IndexMap) = A.dims @@ -87,11 +92,14 @@ LinearAlgebra.ishermitian(A::IndexMap{<:Real}) = issymmetric(A) # comparison of IndexMap objects Base.:(==)(A::IndexMap, B::IndexMap) = (eltype(A) == eltype(B)) && (A.map == B.map) && (A.dims == B.dims) && - (A.rows == B.rows) && (A.cols == B.cols) + (A.rows == B.rows) && (A.cols == B.cols) && + (A.set0 == B.set0) # special transposition behavior -LinearAlgebra.transpose(A::IndexMap) = TransposeMap(A) -LinearAlgebra.adjoint(A::IndexMap) = AdjointMap(A) +LinearAlgebra.transpose(A::IndexMap) = + IndexMap(transpose(A.map), reverse(A.dims), A.cols, A.rows ; set0=A.set0) +LinearAlgebra.adjoint(A::IndexMap) = + IndexMap(adjoint(A.map), reverse(A.dims), A.cols, A.rows ; set0=A.set0) ############ @@ -101,30 +109,16 @@ LinearAlgebra.adjoint(A::IndexMap) = AdjointMap(A) function A_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) m, n = size(A) @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) - y[:] .= 0 +# @show A.set0, A.cols # todo: debug + if A.set0 + fill!(y, zero(eltype(y))) + end @views @inbounds mul!(y[A.rows], A.map, x[A.cols], true, true) return y end -function At_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) - m, n = size(A) - @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!")) - y[:] .= 0 - @views @inbounds mul!(y[A.cols], transpose(A.map), x[A.rows], true, true) - return y -end - -function Ac_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) - m, n = size(A) - y[:] .= 0 - @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("Ac_mul_B!")) - @views @inbounds mul!(y[A.cols], adjoint(A.map), x[A.rows], true, true) - return y -end -#= -=# -"simplified version of hcat" +"(possibly) simplified version of hcat - initial draft" function hcat_new(As::LinearMap...) rows = collect(map(A -> size(A,1), As)) any(rows .!= rows[1]) && throw("row mismatch") @@ -133,6 +127,30 @@ function hcat_new(As::LinearMap...) nmap = length(As) col_offsets = [0; cumsum(cols)[1:nmap-1]] dims = (rows[1], sum(cols)) - return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii])) + return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii]), + set0=true) # todo: for now always set to 0 to see timing penalty +# todo set0=(ii==1)) +# this was foiled by reversed order in +# https://github.com/Jutho/LinearMaps.jl/blob/master/src/linearcombination.jl#L33 + for ii in 1:nmap]...) +end + + +""" +first attempt at block diagonal. +it works, but has the overhead of zeroing out every block's output vector +""" +function block_diag(As::LinearMap...) + rows = collect(map(A -> size(A,1), As)) + cols = collect(map(A -> size(A,2), As)) + dims = (sum(rows), sum(cols)) + nmap = length(As) + col_offsets = [0; cumsum(cols)[1:nmap-1]] + row_offsets = [0; cumsum(rows)[1:nmap-1]] + return +([IndexMap(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii]), + set0=true) # todo: for now always set to 0 to see timing penalty +# todo set0=(ii==1)) for ii in 1:nmap]...) end + +# todo: instead of relying on LinearCombination, make a special sum type? From 6395500ceb8b544d82a89fa907ca073ec1eaa32e Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 29 Sep 2019 21:52:57 -0400 Subject: [PATCH 06/15] Flag possible errors in mul! for LinearCombination --- src/linearcombination.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 55e2e6d3..9ab93325 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -98,7 +98,7 @@ function LinearAlgebra.mul!(y::AbstractVector, A::LinearCombination{T,As}, x::Ab length(y) == size(A, 1) || throw(DimensionMismatch("mul!")) if isone(α) iszero(β) && (A_mul_B!(y, A, x); return y) - !isone(β) && rmul!(y, β) + !isone(β) && rmul!(y, β) # todo: this seems incorrect, cf LinearMaps.jl elseif iszero(α) iszero(β) && (fill!(y, zero(eltype(y))); return y) isone(β) && return y @@ -111,7 +111,7 @@ function LinearAlgebra.mul!(y::AbstractVector, A::LinearCombination{T,As}, x::Ab rmul!(y, α) return y elseif !isone(β) - rmul!(y, β) + rmul!(y, β) # todo: also seems incorrect, missing A * x ? end # β-cases end # α-cases From 4bdc059765e5fcac0dfaaf9319ee39eb636048b7 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 29 Sep 2019 22:12:33 -0400 Subject: [PATCH 07/15] Add indexmap tests --- test/indexmap.jl | 77 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ 2 files changed, 79 insertions(+) create mode 100644 test/indexmap.jl diff --git a/test/indexmap.jl b/test/indexmap.jl new file mode 100644 index 00000000..d7263b0a --- /dev/null +++ b/test/indexmap.jl @@ -0,0 +1,77 @@ +#= +test/indexmap.jl +2019-09-29 Jeff Fessler, University of Michigan +=# + +using LinearAlgebra: mul! +using LinearMaps: LinearMap +using Random: seed! +using Test: @test, @testset + +@testset "indexmap" begin + m = 5; n = 6 + M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) + M2 = randn(m,n+2); L2 = LinearMap(M2) + M3 = randn(m,n+3); L3 = LinearMap(M3) + Mh = [M1 M2 M3] + Lh = [L1 L2 L3] + offset = (3,4) + BM1 = [zeros(offset...) zeros(offset[1],size(M1,2)); + zeros(size(M1,1),offset[2]) M1] + BL1 = @inferred LinearMap(L1, size(BM1) ; offset=offset) + (s1,s2) = size(BM1) + + seed!(0) + x = randn(s2) + y = BM1 * x + + @test @inferred BL1 * x ≈ BM1 * x + @test @inferred BL1' * y ≈ BM1' * y + + # todo: add @test_throws + + # test all flavors of 5-arg mul! for IndexMap + for α in (0, 0.3, 1) + for β in (0, 0.4, 1) + # @show α, β + yl1 = randn(s1) + ym1 = copy(yl1) + # mul!(ym1, BM1, x, α, β) # fails in 1.2 + ym1 = α * (BM1 * x) .+ β * ym1 + @test @inferred mul!(yl1, BL1, x, α, β) ≈ ym1 + end + end + + @test @inferred Matrix(BL1) == BM1 + @test @inferred Matrix(BL1') == Matrix(BL1)' + +#end + + +#@testset "hcat_new" begin + + x = randn(size(Mh,2)) + Bh = hcat_new(L1, L2, L3) +# Bh = @inferred hcat_new(L1, L2, L3) # todo: type instability + @test Mh == Matrix(Lh) + @test Mh == Matrix(Bh) + @test Matrix(Bh') == Matrix(Bh)' + @test Mh * x ≈ Bh * x + +end + +#= +todo +@testset "blockdiag" begin + + using SparseArrays: blockdiag, sparse + Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) + # Md = diag(M1, M2, M3) + x = randn(size(Md,2)) + Bd = blockdiag(L1, L2, L3, L2, L1) + @test @inferred Bd * x ≈ Md * x + @test @inferred Matrix(Bd) == Md + @test @inferred Matrix(Bd') == Matrix(Bd)' + +end +=# diff --git a/test/runtests.jl b/test/runtests.jl index 8ed8ba4b..ded36544 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,7 @@ using Test, LinearMaps +include("indexmap.jl") + include("linearmaps.jl") include("transpose.jl") From f2e6e9a9a51d971236b8b9a1de9281baefdcc8fb Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 29 Sep 2019 22:13:00 -0400 Subject: [PATCH 08/15] Add Random for tests --- Manifest.toml | 19 +++++++++++++++++++ Project.toml | 1 + 2 files changed, 20 insertions(+) create mode 100644 Manifest.toml diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 00000000..00f61e2c --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,19 @@ +# This file is machine-generated - editing it directly is not advised + +[[Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[LinearAlgebra]] +deps = ["Libdl"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[Random]] +deps = ["Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/Project.toml b/Project.toml index 5fdb8be1..e9afbb53 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "2.5.2" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] From 7e8b4476febe5fc590920e6389eea0e595b3ca02 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 29 Sep 2019 22:13:38 -0400 Subject: [PATCH 09/15] Update for 5-arg mul! --- src/indexmap.jl | 184 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 129 insertions(+), 55 deletions(-) diff --git a/src/indexmap.jl b/src/indexmap.jl index 3e006f95..079707b2 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -3,14 +3,15 @@ indexmap.jl 2019-08-29 Jeff Fessler, University of Michigan =# -export IndexMap, hcat_new, block_diag +export hcat_new, blockdiag + """ `check_index(row::AbstractVector{Int}, dims::Dims{2})` Do not allow repeated index values because that would complicate the adjoint. -And insist on monotone increasing too because it is the easiest way to -ensure no repeats and probably better for cache too. +Insist on monotone increasing too because it is the easiest way to +ensure no repeats and is probably better for cache too. """ function check_index(index::AbstractVector{Int}, dimA::Int, dimB::Int) length(index) != dimA && throw("invalid index length") @@ -27,7 +28,7 @@ end """ `IndexMap` -The original motivation was to support construction of linear operator +The original motivation was to support construction of linear operators of the form `B = [0 0 0; 0 A 0; 0 0 0]` from a given linear operator `A` where `size(B) == dims` @@ -36,16 +37,16 @@ This version is more general than the form above because it allows the rows and colums of `A` to be interspersed among the rows and columns of `B` arbitrarily. Essentially like the above form but with some permutation-like matrices -before and after. +before and after, implicitly. """ struct IndexMap{T, As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} <: LinearMap{T} map::As dims::Dims{2} - set0::Bool # set to zero the output array? default true rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1) - cols::Cs # typically j1:j2 - function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, set0::Bool, + cols::Cs # typically j1:j2 with 1 <= j1 <= j2 <= size(map,2) + + function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, rows::Rs, cols::Cs) where {T, As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} @@ -53,33 +54,47 @@ struct IndexMap{T, As <: LinearMap, Rs <: AbstractVector{Int}, check_index(rows, size(map,1), dims[1]) check_index(cols, size(map,2), dims[2]) - return new{T,As,Rs,Cs}(map, dims, set0, rows, cols) + return new{T,As,Rs,Cs}(map, dims, rows, cols) end end -IndexMap{T}(map::As, dims::Dims{2} ; set0::Bool = true, offset::Dims{2}) where +IndexMap{T}(map::As, dims::Dims{2} ; offset::Dims{2}) where {T, As <: LinearMap} = - IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, set0, + IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, offset[1] .+ (1:size(map,1)), offset[2] .+ (1:size(map,2))) -IndexMap(map::LinearMap, dims::Dims{2} ; set0::Bool = true, offset::Dims{2}) = - IndexMap{eltype(map)}(map, dims ; offset=offset, set0=set0) +IndexMap(map::LinearMap, dims::Dims{2} ; offset::Dims{2}) = + IndexMap{eltype(map)}(map, dims, offset=offset) IndexMap(map::LinearMap, dims::Dims{2}, - rows::AbstractVector{Int}, cols::AbstractVector{Int}, - ; set0::Bool = true) = + rows::AbstractVector{Int}, cols::AbstractVector{Int}) = IndexMap{eltype(map), typeof(map), typeof(rows), typeof(cols)}( - map, dims, set0, rows, cols) + map, dims, rows, cols) Base.size(A::IndexMap) = A.dims -############ -# basic methods -############ +# Provide constructors via LinearMap + +LinearMap(A::LinearMap, dims::Dims{2}, + index::Tuple{AbstractVector{Int}, AbstractVector{Int}}) = + IndexMap(A, dims, index[1], index[2]) # given (rows,cols) tuple + +LinearMap(A::LinearMap, dims::Dims{2} ; offset::Dims{2}) = + IndexMap(A, dims ; offset=offset) # given offset + +#= +LinearMap(A::LinearMap ; offset::Dims{2}, dims::Dims{2}=size(A) .+ offset) = + IndexMap(A, dims, offset) # possible alternative offset syntax +=# + + +#= +basic methods +=# -# the following rules are sufficient but perhaps not necessary +# the following symmetry rules are sufficient but perhaps not necessary LinearAlgebra.issymmetric(A::IndexMap) = issymmetric(A.map) && (A.dims[1] == A.dims[2]) && (A.rows == A.cols) @@ -89,68 +104,127 @@ LinearAlgebra.ishermitian(A::IndexMap) = ishermitian(A.map) && LinearAlgebra.ishermitian(A::IndexMap{<:Real}) = issymmetric(A) -# comparison of IndexMap objects +# compare IndexMap objects Base.:(==)(A::IndexMap, B::IndexMap) = (eltype(A) == eltype(B)) && (A.map == B.map) && (A.dims == B.dims) && - (A.rows == B.rows) && (A.cols == B.cols) && - (A.set0 == B.set0) + (A.rows == B.rows) && (A.cols == B.cols) -# special transposition behavior +# transpose/adjoint by simple swap LinearAlgebra.transpose(A::IndexMap) = - IndexMap(transpose(A.map), reverse(A.dims), A.cols, A.rows ; set0=A.set0) + IndexMap(transpose(A.map), reverse(A.dims), A.cols, A.rows) LinearAlgebra.adjoint(A::IndexMap) = - IndexMap(adjoint(A.map), reverse(A.dims), A.cols, A.rows ; set0=A.set0) - + IndexMap(adjoint(A.map), reverse(A.dims), A.cols, A.rows) -############ -# multiplication with vectors -############ - -function A_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) - m, n = size(A) - @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) -# @show A.set0, A.cols # todo: debug - if A.set0 - fill!(y, zero(eltype(y))) - end - @views @inbounds mul!(y[A.rows], A.map, x[A.cols], true, true) - return y -end +#= +combinations of IndexMap objects +=# -"(possibly) simplified version of hcat - initial draft" +""" +(possibly) simplified version of hcat - initial draft +Requires 5-arg mul! to be efficient. +""" function hcat_new(As::LinearMap...) rows = collect(map(A -> size(A,1), As)) any(rows .!= rows[1]) && throw("row mismatch") cols = collect(map(A -> size(A,2), As)) - nmap = length(As) col_offsets = [0; cumsum(cols)[1:nmap-1]] dims = (rows[1], sum(cols)) - return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii]), - set0=true) # todo: for now always set to 0 to see timing penalty -# todo set0=(ii==1)) -# this was foiled by reversed order in -# https://github.com/Jutho/LinearMaps.jl/blob/master/src/linearcombination.jl#L33 + return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii])) for ii in 1:nmap]...) end """ -first attempt at block diagonal. -it works, but has the overhead of zeroing out every block's output vector +`B = blockdiag(As::LinearMap...)` + +Block diagonal `LinearMap`, akin to `SparseArrays.blockdiag()` +Requires 5-arg mul! to be efficient. """ -function block_diag(As::LinearMap...) +function blockdiag(As::LinearMap...) rows = collect(map(A -> size(A,1), As)) cols = collect(map(A -> size(A,2), As)) dims = (sum(rows), sum(cols)) nmap = length(As) col_offsets = [0; cumsum(cols)[1:nmap-1]] row_offsets = [0; cumsum(rows)[1:nmap-1]] - return +([IndexMap(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii]), - set0=true) # todo: for now always set to 0 to see timing penalty -# todo set0=(ii==1)) + return +([IndexMap(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii])) for ii in 1:nmap]...) end -# todo: instead of relying on LinearCombination, make a special sum type? + +#= +multiplication with vectors + +When IndexMaps are put in a LinearCombination, hopefully A_mul_B! +is called only on the "first" block, so that the overhead of `fill` +happens just once. +=# + +#if VERSION < v"1.3.0-alpha.115" + +""" +Multiply a single IndexMap with a vector. +Basic *(A,x) in LinearMaps.jl#L22 uses similar() so must zero `y` here. +In principle we could zero out just the complement of `A.rows`. +""" +function A_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) + m, n = size(A) +# @show "todoIM-AB:", size(A) + @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) + fill!(y, zero(eltype(y))) + @views @inbounds mul!(y[A.rows], A.map, x[A.cols], true, true) + return y +end + + +""" +This 5-arg mul! is based on LinearAlgebra.mul! in LinearMaps.jl +with modifications to support the indexing used by an IndexMap. +y = α B x + β y +""" +function LinearAlgebra.mul!(y_in::AbstractVector, B::IndexMap, + x_in::AbstractVector, α::Number=true, β::Number=false) + length(y_in) == size(B, 1) || throw(DimensionMismatch("mul!")) + A = B.map +# @show "todoIM:", size(B), α, β + xv = @views @inbounds x_in[B.cols] + yv = @views @inbounds y_in[B.rows] + if isone(α) # α = 1, so B x + β y + if iszero(β) # β = 0: basic y = B*x case + fill!(y_in, zero(eltype(y_in))) + A_mul_B!(yv, A, xv) + elseif isone(β) # β = 1 + yv .+= A * xv + else # β != 0, 1 + rmul!(y_in, β) + yv .+= A * xv + end # β-cases + elseif iszero(α) # α = 0, so β y + if iszero(β) # β = 0 + fill!(y_in, zero(eltype(y_in))) + elseif !isone(β) # β != 1 + rmul!(y_in, β) # β != 0, 1 + end # β-cases + else # α != 0, 1, so α B x + β y + if iszero(β) # β = 0, so α B x + fill!(y_in, zero(eltype(y_in))) + A_mul_B!(yv, A, xv) + rmul!(yv, α) + else + if !isone(β) # β != 0, 1 + rmul!(y_in, β) + end + yv .+= rmul!(A * xv, α) + end # β-cases + end # α-cases + + return y_in +end + +#else # 5-arg mul! is available for matrices + + # throw("todo: not done") + +#end # VERSION From 905ab6d7d368581d6bbf590d84e299a05be4cecd Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 29 Sep 2019 23:06:59 -0400 Subject: [PATCH 10/15] Rename block_diag --- src/indexmap.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/indexmap.jl b/src/indexmap.jl index 079707b2..bec369d5 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -3,7 +3,8 @@ indexmap.jl 2019-08-29 Jeff Fessler, University of Michigan =# -export hcat_new, blockdiag +export hcat_new +export block_diag # use _ to avoid odd conflict with SparseArrays.blockdiag """ @@ -14,6 +15,7 @@ Insist on monotone increasing too because it is the easiest way to ensure no repeats and is probably better for cache too. """ function check_index(index::AbstractVector{Int}, dimA::Int, dimB::Int) +#@show index, dimA, dimB length(index) != dimA && throw("invalid index length") minimum(index) < 0 && throw("negative index") maximum(index) > dimB && @@ -137,12 +139,12 @@ end """ -`B = blockdiag(As::LinearMap...)` +`B = block_diag(As::LinearMap...)` Block diagonal `LinearMap`, akin to `SparseArrays.blockdiag()` Requires 5-arg mul! to be efficient. """ -function blockdiag(As::LinearMap...) +function block_diag(As::LinearMap...) rows = collect(map(A -> size(A,1), As)) cols = collect(map(A -> size(A,2), As)) dims = (sum(rows), sum(cols)) From 9cd78571db9db60fd6f12f12f4ff95a7b652ae29 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Sun, 29 Sep 2019 23:07:19 -0400 Subject: [PATCH 11/15] Update tests --- test/block_diag.jl | 24 ++++++++++++++++++++++++ test/indexmap.jl | 24 +++++++----------------- test/runtests.jl | 2 ++ 3 files changed, 33 insertions(+), 17 deletions(-) create mode 100644 test/block_diag.jl diff --git a/test/block_diag.jl b/test/block_diag.jl new file mode 100644 index 00000000..0faeacde --- /dev/null +++ b/test/block_diag.jl @@ -0,0 +1,24 @@ +#= +test/block_diag.jl +2019-09-29 Jeff Fessler, University of Michigan +=# + +using LinearMaps: LinearMap +using SparseArrays: sparse, blockdiag +using Random: seed! +using Test: @test, @testset + +@testset "block_diag" begin + m = 5; n = 6 + M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) + M2 = randn(m,n+2); L2 = LinearMap(M2) + M3 = randn(m,n+3); L3 = LinearMap(M3) + + # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: + Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) + x = randn(size(Md,2)) + Bd = block_diag(L1, L2, L3, L2, L1) + @test @inferred Bd * x ≈ Md * x + @test @inferred Matrix(Bd) == Md + @test @inferred Matrix(Bd') == Matrix(Bd)' +end diff --git a/test/indexmap.jl b/test/indexmap.jl index d7263b0a..771246d8 100644 --- a/test/indexmap.jl +++ b/test/indexmap.jl @@ -3,7 +3,7 @@ test/indexmap.jl 2019-09-29 Jeff Fessler, University of Michigan =# -using LinearAlgebra: mul! +using LinearAlgebra: mul!, issymmetric, ishermitian using LinearMaps: LinearMap using Random: seed! using Test: @test, @testset @@ -21,6 +21,11 @@ using Test: @test, @testset BL1 = @inferred LinearMap(L1, size(BM1) ; offset=offset) (s1,s2) = size(BM1) + @test @inferred !issymmetric(BL1) + @test @inferred !ishermitian(BL1) + @test @inferred LinearMap(L1, size(BM1), + (offset[1] .+ (1:m), offset[2] .+ (1:(n+1)))) == BL1 # test tuple and == + seed!(0) x = randn(s2) y = BM1 * x @@ -44,6 +49,7 @@ using Test: @test, @testset @test @inferred Matrix(BL1) == BM1 @test @inferred Matrix(BL1') == Matrix(BL1)' + @test @inferred transpose(BL1) * y ≈ transpose(BM1) * y #end @@ -59,19 +65,3 @@ using Test: @test, @testset @test Mh * x ≈ Bh * x end - -#= -todo -@testset "blockdiag" begin - - using SparseArrays: blockdiag, sparse - Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) - # Md = diag(M1, M2, M3) - x = randn(size(Md,2)) - Bd = blockdiag(L1, L2, L3, L2, L1) - @test @inferred Bd * x ≈ Md * x - @test @inferred Matrix(Bd) == Md - @test @inferred Matrix(Bd') == Matrix(Bd)' - -end -=# diff --git a/test/runtests.jl b/test/runtests.jl index ded36544..53210315 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,8 @@ using Test, LinearMaps include("indexmap.jl") +include("block_diag.jl") + include("linearmaps.jl") include("transpose.jl") From 5ff80ff90aa67736ad8e14d7004bda21f605a537 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Mon, 30 Sep 2019 06:16:29 -0400 Subject: [PATCH 12/15] Tab to 4-space --- src/indexmap.jl | 202 +++++++++++++++++++++++------------------------ test/indexmap.jl | 84 ++++++++++---------- 2 files changed, 143 insertions(+), 143 deletions(-) diff --git a/src/indexmap.jl b/src/indexmap.jl index bec369d5..0739ebde 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -16,14 +16,14 @@ ensure no repeats and is probably better for cache too. """ function check_index(index::AbstractVector{Int}, dimA::Int, dimB::Int) #@show index, dimA, dimB - length(index) != dimA && throw("invalid index length") - minimum(index) < 0 && throw("negative index") - maximum(index) > dimB && - throw("index $(maximum(index)) > dimB=$dimB") - !(index isa UnitRange{Int64} || - ((index isa StepRange{Int64,Int64}) && (index.step > 0)) || - all(diff(index) .> 0)) && throw("non-monotone index") - nothing + length(index) != dimA && throw("invalid index length") + minimum(index) < 0 && throw("negative index") + maximum(index) > dimB && + throw("index $(maximum(index)) > dimB=$dimB") + !(index isa UnitRange{Int64} || + ((index isa StepRange{Int64,Int64}) && (index.step > 0)) || + all(diff(index) .> 0)) && throw("non-monotone index") + nothing end @@ -42,37 +42,37 @@ Essentially like the above form but with some permutation-like matrices before and after, implicitly. """ struct IndexMap{T, As <: LinearMap, Rs <: AbstractVector{Int}, - Cs <: AbstractVector{Int}} <: LinearMap{T} - map::As - dims::Dims{2} - rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1) - cols::Cs # typically j1:j2 with 1 <= j1 <= j2 <= size(map,2) - - function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, - rows::Rs, cols::Cs) where - {T, As <: LinearMap, - Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} - - check_index(rows, size(map,1), dims[1]) - check_index(cols, size(map,2), dims[2]) - - return new{T,As,Rs,Cs}(map, dims, rows, cols) - end + Cs <: AbstractVector{Int}} <: LinearMap{T} + map::As + dims::Dims{2} + rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1) + cols::Cs # typically j1:j2 with 1 <= j1 <= j2 <= size(map,2) + + function IndexMap{T,As,Rs,Cs}(map::As, dims::Dims{2}, + rows::Rs, cols::Cs) where + {T, As <: LinearMap, + Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}} + + check_index(rows, size(map,1), dims[1]) + check_index(cols, size(map,2), dims[2]) + + return new{T,As,Rs,Cs}(map, dims, rows, cols) + end end IndexMap{T}(map::As, dims::Dims{2} ; offset::Dims{2}) where - {T, As <: LinearMap} = - IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, - offset[1] .+ (1:size(map,1)), - offset[2] .+ (1:size(map,2))) + {T, As <: LinearMap} = + IndexMap{T, As, UnitRange{Int64}, UnitRange{Int64}}(map, dims, + offset[1] .+ (1:size(map,1)), + offset[2] .+ (1:size(map,2))) IndexMap(map::LinearMap, dims::Dims{2} ; offset::Dims{2}) = - IndexMap{eltype(map)}(map, dims, offset=offset) + IndexMap{eltype(map)}(map, dims, offset=offset) IndexMap(map::LinearMap, dims::Dims{2}, - rows::AbstractVector{Int}, cols::AbstractVector{Int}) = - IndexMap{eltype(map), typeof(map), typeof(rows), typeof(cols)}( - map, dims, rows, cols) + rows::AbstractVector{Int}, cols::AbstractVector{Int}) = + IndexMap{eltype(map), typeof(map), typeof(rows), typeof(cols)}( + map, dims, rows, cols) Base.size(A::IndexMap) = A.dims @@ -80,15 +80,15 @@ Base.size(A::IndexMap) = A.dims # Provide constructors via LinearMap LinearMap(A::LinearMap, dims::Dims{2}, - index::Tuple{AbstractVector{Int}, AbstractVector{Int}}) = - IndexMap(A, dims, index[1], index[2]) # given (rows,cols) tuple + index::Tuple{AbstractVector{Int}, AbstractVector{Int}}) = + IndexMap(A, dims, index[1], index[2]) # given (rows,cols) tuple LinearMap(A::LinearMap, dims::Dims{2} ; offset::Dims{2}) = - IndexMap(A, dims ; offset=offset) # given offset + IndexMap(A, dims ; offset=offset) # given offset #= LinearMap(A::LinearMap ; offset::Dims{2}, dims::Dims{2}=size(A) .+ offset) = - IndexMap(A, dims, offset) # possible alternative offset syntax + IndexMap(A, dims, offset) # possible alternative offset syntax =# @@ -98,24 +98,24 @@ basic methods # the following symmetry rules are sufficient but perhaps not necessary LinearAlgebra.issymmetric(A::IndexMap) = issymmetric(A.map) && - (A.dims[1] == A.dims[2]) && (A.rows == A.cols) + (A.dims[1] == A.dims[2]) && (A.rows == A.cols) LinearAlgebra.ishermitian(A::IndexMap) = ishermitian(A.map) && - (A.dims[1] == A.dims[2]) && (A.rows == A.cols) + (A.dims[1] == A.dims[2]) && (A.rows == A.cols) LinearAlgebra.ishermitian(A::IndexMap{<:Real}) = issymmetric(A) # compare IndexMap objects Base.:(==)(A::IndexMap, B::IndexMap) = (eltype(A) == eltype(B)) && - (A.map == B.map) && (A.dims == B.dims) && - (A.rows == B.rows) && (A.cols == B.cols) + (A.map == B.map) && (A.dims == B.dims) && + (A.rows == B.rows) && (A.cols == B.cols) # transpose/adjoint by simple swap LinearAlgebra.transpose(A::IndexMap) = - IndexMap(transpose(A.map), reverse(A.dims), A.cols, A.rows) + IndexMap(transpose(A.map), reverse(A.dims), A.cols, A.rows) LinearAlgebra.adjoint(A::IndexMap) = - IndexMap(adjoint(A.map), reverse(A.dims), A.cols, A.rows) + IndexMap(adjoint(A.map), reverse(A.dims), A.cols, A.rows) #= @@ -127,14 +127,14 @@ combinations of IndexMap objects Requires 5-arg mul! to be efficient. """ function hcat_new(As::LinearMap...) - rows = collect(map(A -> size(A,1), As)) - any(rows .!= rows[1]) && throw("row mismatch") - cols = collect(map(A -> size(A,2), As)) - nmap = length(As) - col_offsets = [0; cumsum(cols)[1:nmap-1]] - dims = (rows[1], sum(cols)) - return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii])) - for ii in 1:nmap]...) + rows = collect(map(A -> size(A,1), As)) + any(rows .!= rows[1]) && throw("row mismatch") + cols = collect(map(A -> size(A,2), As)) + nmap = length(As) + col_offsets = [0; cumsum(cols)[1:nmap-1]] + dims = (rows[1], sum(cols)) + return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii])) + for ii in 1:nmap]...) end @@ -145,14 +145,14 @@ Block diagonal `LinearMap`, akin to `SparseArrays.blockdiag()` Requires 5-arg mul! to be efficient. """ function block_diag(As::LinearMap...) - rows = collect(map(A -> size(A,1), As)) - cols = collect(map(A -> size(A,2), As)) - dims = (sum(rows), sum(cols)) - nmap = length(As) - col_offsets = [0; cumsum(cols)[1:nmap-1]] - row_offsets = [0; cumsum(rows)[1:nmap-1]] - return +([IndexMap(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii])) - for ii in 1:nmap]...) + rows = collect(map(A -> size(A,1), As)) + cols = collect(map(A -> size(A,2), As)) + dims = (sum(rows), sum(cols)) + nmap = length(As) + col_offsets = [0; cumsum(cols)[1:nmap-1]] + row_offsets = [0; cumsum(rows)[1:nmap-1]] + return +([IndexMap(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii])) + for ii in 1:nmap]...) end @@ -172,12 +172,12 @@ Basic *(A,x) in LinearMaps.jl#L22 uses similar() so must zero `y` here. In principle we could zero out just the complement of `A.rows`. """ function A_mul_B!(y::AbstractVector, A::IndexMap, x::AbstractVector) - m, n = size(A) -# @show "todoIM-AB:", size(A) - @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) - fill!(y, zero(eltype(y))) - @views @inbounds mul!(y[A.rows], A.map, x[A.cols], true, true) - return y + m, n = size(A) +# @show "todoIM-AB:", size(A) + @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) + fill!(y, zero(eltype(y))) + @views @inbounds mul!(y[A.rows], A.map, x[A.cols], true, true) + return y end @@ -187,46 +187,46 @@ with modifications to support the indexing used by an IndexMap. y = α B x + β y """ function LinearAlgebra.mul!(y_in::AbstractVector, B::IndexMap, - x_in::AbstractVector, α::Number=true, β::Number=false) - length(y_in) == size(B, 1) || throw(DimensionMismatch("mul!")) - A = B.map -# @show "todoIM:", size(B), α, β - xv = @views @inbounds x_in[B.cols] - yv = @views @inbounds y_in[B.rows] - if isone(α) # α = 1, so B x + β y - if iszero(β) # β = 0: basic y = B*x case - fill!(y_in, zero(eltype(y_in))) - A_mul_B!(yv, A, xv) - elseif isone(β) # β = 1 - yv .+= A * xv - else # β != 0, 1 - rmul!(y_in, β) - yv .+= A * xv - end # β-cases - elseif iszero(α) # α = 0, so β y - if iszero(β) # β = 0 - fill!(y_in, zero(eltype(y_in))) - elseif !isone(β) # β != 1 - rmul!(y_in, β) # β != 0, 1 - end # β-cases - else # α != 0, 1, so α B x + β y - if iszero(β) # β = 0, so α B x - fill!(y_in, zero(eltype(y_in))) - A_mul_B!(yv, A, xv) - rmul!(yv, α) - else - if !isone(β) # β != 0, 1 - rmul!(y_in, β) - end - yv .+= rmul!(A * xv, α) - end # β-cases - end # α-cases - - return y_in + x_in::AbstractVector, α::Number=true, β::Number=false) + length(y_in) == size(B, 1) || throw(DimensionMismatch("mul!")) + A = B.map +# @show "todoIM:", size(B), α, β + xv = @views @inbounds x_in[B.cols] + yv = @views @inbounds y_in[B.rows] + if isone(α) # α = 1, so B x + β y + if iszero(β) # β = 0: basic y = B*x case + fill!(y_in, zero(eltype(y_in))) + A_mul_B!(yv, A, xv) + elseif isone(β) # β = 1 + yv .+= A * xv + else # β != 0, 1 + rmul!(y_in, β) + yv .+= A * xv + end # β-cases + elseif iszero(α) # α = 0, so β y + if iszero(β) # β = 0 + fill!(y_in, zero(eltype(y_in))) + elseif !isone(β) # β != 1 + rmul!(y_in, β) # β != 0, 1 + end # β-cases + else # α != 0, 1, so α B x + β y + if iszero(β) # β = 0, so α B x + fill!(y_in, zero(eltype(y_in))) + A_mul_B!(yv, A, xv) + rmul!(yv, α) + else + if !isone(β) # β != 0, 1 + rmul!(y_in, β) + end + yv .+= rmul!(A * xv, α) + end # β-cases + end # α-cases + + return y_in end #else # 5-arg mul! is available for matrices - # throw("todo: not done") + # throw("todo: not done") #end # VERSION diff --git a/test/indexmap.jl b/test/indexmap.jl index 771246d8..840ad966 100644 --- a/test/indexmap.jl +++ b/test/indexmap.jl @@ -9,59 +9,59 @@ using Random: seed! using Test: @test, @testset @testset "indexmap" begin - m = 5; n = 6 - M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) - M2 = randn(m,n+2); L2 = LinearMap(M2) - M3 = randn(m,n+3); L3 = LinearMap(M3) - Mh = [M1 M2 M3] - Lh = [L1 L2 L3] - offset = (3,4) - BM1 = [zeros(offset...) zeros(offset[1],size(M1,2)); - zeros(size(M1,1),offset[2]) M1] - BL1 = @inferred LinearMap(L1, size(BM1) ; offset=offset) - (s1,s2) = size(BM1) + m = 5; n = 6 + M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) + M2 = randn(m,n+2); L2 = LinearMap(M2) + M3 = randn(m,n+3); L3 = LinearMap(M3) + Mh = [M1 M2 M3] + Lh = [L1 L2 L3] + offset = (3,4) + BM1 = [zeros(offset...) zeros(offset[1],size(M1,2)); + zeros(size(M1,1),offset[2]) M1] + BL1 = @inferred LinearMap(L1, size(BM1) ; offset=offset) + (s1,s2) = size(BM1) - @test @inferred !issymmetric(BL1) - @test @inferred !ishermitian(BL1) - @test @inferred LinearMap(L1, size(BM1), - (offset[1] .+ (1:m), offset[2] .+ (1:(n+1)))) == BL1 # test tuple and == + @test @inferred !ishermitian(1im * BL1) + @test @inferred !issymmetric(BL1) + @test @inferred LinearMap(L1, size(BM1), + (offset[1] .+ (1:m), offset[2] .+ (1:(n+1)))) == BL1 # test tuple and == - seed!(0) - x = randn(s2) - y = BM1 * x + seed!(0) + x = randn(s2) + y = BM1 * x - @test @inferred BL1 * x ≈ BM1 * x - @test @inferred BL1' * y ≈ BM1' * y + @test @inferred BL1 * x ≈ BM1 * x + @test @inferred BL1' * y ≈ BM1' * y - # todo: add @test_throws + # todo: add @test_throws - # test all flavors of 5-arg mul! for IndexMap - for α in (0, 0.3, 1) - for β in (0, 0.4, 1) - # @show α, β - yl1 = randn(s1) - ym1 = copy(yl1) - # mul!(ym1, BM1, x, α, β) # fails in 1.2 - ym1 = α * (BM1 * x) .+ β * ym1 - @test @inferred mul!(yl1, BL1, x, α, β) ≈ ym1 - end - end + # test all flavors of 5-arg mul! for IndexMap + for α in (0, 0.3, 1) + for β in (0, 0.4, 1) + # @show α, β + yl1 = randn(s1) + ym1 = copy(yl1) + # mul!(ym1, BM1, x, α, β) # fails in 1.2 + ym1 = α * (BM1 * x) .+ β * ym1 + @test @inferred mul!(yl1, BL1, x, α, β) ≈ ym1 + end + end - @test @inferred Matrix(BL1) == BM1 - @test @inferred Matrix(BL1') == Matrix(BL1)' - @test @inferred transpose(BL1) * y ≈ transpose(BM1) * y + @test @inferred Matrix(BL1) == BM1 + @test @inferred Matrix(BL1') == Matrix(BL1)' + @test @inferred transpose(BL1) * y ≈ transpose(BM1) * y #end #@testset "hcat_new" begin - x = randn(size(Mh,2)) - Bh = hcat_new(L1, L2, L3) -# Bh = @inferred hcat_new(L1, L2, L3) # todo: type instability - @test Mh == Matrix(Lh) - @test Mh == Matrix(Bh) - @test Matrix(Bh') == Matrix(Bh)' - @test Mh * x ≈ Bh * x + x = randn(size(Mh,2)) + Bh = hcat_new(L1, L2, L3) +# Bh = @inferred hcat_new(L1, L2, L3) # todo: type instability + @test Mh == Matrix(Lh) + @test Mh == Matrix(Bh) + @test Matrix(Bh') == Matrix(Bh)' + @test Mh * x ≈ Bh * x end From fd8a30defcca363eb87333f674be9845d882953f Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Mon, 30 Sep 2019 06:46:32 -0400 Subject: [PATCH 13/15] Try to fix type instability hcat_new block_diag --- src/indexmap.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/indexmap.jl b/src/indexmap.jl index 0739ebde..a3d527df 100644 --- a/src/indexmap.jl +++ b/src/indexmap.jl @@ -86,6 +86,9 @@ LinearMap(A::LinearMap, dims::Dims{2}, LinearMap(A::LinearMap, dims::Dims{2} ; offset::Dims{2}) = IndexMap(A, dims ; offset=offset) # given offset +#LinearMap{T}(A::LinearMap, dims::Dims{2} ; offset::Dims{2}) where {T} = +# IndexMap{T}(A, dims ; offset=offset) # given offset + #= LinearMap(A::LinearMap ; offset::Dims{2}, dims::Dims{2}=size(A) .+ offset) = IndexMap(A, dims, offset) # possible alternative offset syntax @@ -133,8 +136,11 @@ function hcat_new(As::LinearMap...) nmap = length(As) col_offsets = [0; cumsum(cols)[1:nmap-1]] dims = (rows[1], sum(cols)) - return +([IndexMap(As[ii], dims ; offset=(0, col_offsets[ii])) - for ii in 1:nmap]...) + T = promote_type(map(eltype, As)...) +# return +([LinearMap{T}(As[ii], dims ; offset=(0, col_offsets[ii])) +# for ii in 1:nmap]...) + return +([IndexMap{T}(As[ii], dims ; offset=(0, col_offsets[ii])) + for ii in 1:nmap]...) end @@ -151,7 +157,8 @@ function block_diag(As::LinearMap...) nmap = length(As) col_offsets = [0; cumsum(cols)[1:nmap-1]] row_offsets = [0; cumsum(rows)[1:nmap-1]] - return +([IndexMap(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii])) + T = promote_type(map(eltype, As)...) + return +([IndexMap{T}(As[ii], dims ; offset=(row_offsets[ii], col_offsets[ii])) for ii in 1:nmap]...) end From d767d0711f90138016f260f09e5ccdcdf4fcf6bc Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Mon, 30 Sep 2019 06:47:56 -0400 Subject: [PATCH 14/15] Update tests --- test/block_diag.jl | 24 +++++++++++++----------- test/indexmap.jl | 5 ++--- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/test/block_diag.jl b/test/block_diag.jl index 0faeacde..5fd34d84 100644 --- a/test/block_diag.jl +++ b/test/block_diag.jl @@ -9,16 +9,18 @@ using Random: seed! using Test: @test, @testset @testset "block_diag" begin - m = 5; n = 6 - M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) - M2 = randn(m,n+2); L2 = LinearMap(M2) - M3 = randn(m,n+3); L3 = LinearMap(M3) + m = 5; n = 6 + M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) + M2 = randn(m,n+2); L2 = LinearMap(M2) + M3 = randn(m,n+3); L3 = LinearMap(M3) - # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: - Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) - x = randn(size(Md,2)) - Bd = block_diag(L1, L2, L3, L2, L1) - @test @inferred Bd * x ≈ Md * x - @test @inferred Matrix(Bd) == Md - @test @inferred Matrix(Bd') == Matrix(Bd)' + # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: + Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) + x = randn(size(Md,2)) +# Bd = @inferred block_diag(L1, L2, L3, L2, L1) # todo: type instability + Bd = block_diag(L1, L2, L3, L2, L1) + @test Bd isa LinearMaps.LinearCombination + @test @inferred Bd * x ≈ Md * x + @test @inferred Matrix(Bd) == Md + @test @inferred Matrix(Bd') == Matrix(Bd)' end diff --git a/test/indexmap.jl b/test/indexmap.jl index 840ad966..bf8f4c95 100644 --- a/test/indexmap.jl +++ b/test/indexmap.jl @@ -33,8 +33,6 @@ using Test: @test, @testset @test @inferred BL1 * x ≈ BM1 * x @test @inferred BL1' * y ≈ BM1' * y - # todo: add @test_throws - # test all flavors of 5-arg mul! for IndexMap for α in (0, 0.3, 1) for β in (0, 0.4, 1) @@ -57,8 +55,9 @@ using Test: @test, @testset #@testset "hcat_new" begin x = randn(size(Mh,2)) +# Bh = @inferred hcat_new(L1, L2, L3) # todo: type instability Bh = hcat_new(L1, L2, L3) -# Bh = @inferred hcat_new(L1, L2, L3) # todo: type instability + @test Bh isa LinearMaps.LinearCombination @test Mh == Matrix(Lh) @test Mh == Matrix(Bh) @test Matrix(Bh') == Matrix(Bh)' From d1dc812e68e0fd5bd2a7e7f5167f472ac3d6daa1 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Mon, 30 Sep 2019 07:10:46 -0400 Subject: [PATCH 15/15] Test complex case --- test/indexmap.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/indexmap.jl b/test/indexmap.jl index bf8f4c95..657b012a 100644 --- a/test/indexmap.jl +++ b/test/indexmap.jl @@ -21,10 +21,13 @@ using Test: @test, @testset BL1 = @inferred LinearMap(L1, size(BM1) ; offset=offset) (s1,s2) = size(BM1) - @test @inferred !ishermitian(1im * BL1) + @test @inferred !ishermitian(BL1) @test @inferred !issymmetric(BL1) @test @inferred LinearMap(L1, size(BM1), (offset[1] .+ (1:m), offset[2] .+ (1:(n+1)))) == BL1 # test tuple and == + Wc = @inferred LinearMap([2 im; -im 0]) + Bc = @inferred LinearMap(Wc, (4,4) ; offset=(2,2)) + @test @inferred ishermitian(Bc) seed!(0) x = randn(s2)