From d562ce1c3fff9dfdb8546bc9c76fa3668064e93a Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 11:07:22 +0200 Subject: [PATCH 01/17] Added algebra interfaces --- src/Algebra/Algebra.jl | 19 +++- src/Algebra/AlgebraInterfaces.jl | 108 ++++++++++++++++++++ src/Algebra/CompressedSparseMatrices.jl | 34 ------ src/Algebra/MethodsAbstractMatrices.jl | 3 - src/Algebra/SparseMatrices.jl | 62 ++--------- src/Algebra/SparseMatrixCSC.jl | 2 +- src/Arrays/Arrays.jl | 2 - src/Arrays/Interface.jl | 21 ---- src/MultiField/MultiFieldArrays.jl | 2 +- test/AlgebraTests/AlgebraInterfacesTests.jl | 56 ++++++++++ test/AlgebraTests/runtests.jl | 2 + 11 files changed, 197 insertions(+), 114 deletions(-) create mode 100644 src/Algebra/AlgebraInterfaces.jl delete mode 100644 src/Algebra/MethodsAbstractMatrices.jl create mode 100644 test/AlgebraTests/AlgebraInterfacesTests.jl diff --git a/src/Algebra/Algebra.jl b/src/Algebra/Algebra.jl index 4b478601b..4e1627f68 100644 --- a/src/Algebra/Algebra.jl +++ b/src/Algebra/Algebra.jl @@ -13,6 +13,17 @@ using Test using Gridap.Helpers +import Base: convert, size, getindex, show, count, * +import LinearAlgebra: mul! +import SparseArrays: nnz, nonzeros, nzrange, findnz, rowvals + +export allocate_vector +export allocate_in_domain +export allocate_in_range +export add_entries! +export scale_entries! +export muladd! + export push_coo! export is_entry_stored export finalize_coo! @@ -71,12 +82,18 @@ export get_vector export SparseMatrixCSR export SymSparseMatrixCSR -include("MethodsAbstractMatrices.jl") +include("AlgebraInterfaces.jl") include("SparseMatrices.jl") include("CompressedSparseMatrices.jl") +include("SparseMatrixCSC.jl") + +include("SparseMatrixCSR.jl") + +include("SymSparseMatrixCSR.jl") + include("NonlinearOperators.jl") include("NonlinearSolvers.jl") diff --git a/src/Algebra/AlgebraInterfaces.jl b/src/Algebra/AlgebraInterfaces.jl new file mode 100644 index 000000000..f9755a000 --- /dev/null +++ b/src/Algebra/AlgebraInterfaces.jl @@ -0,0 +1,108 @@ + +""" + allocate_vector(::Type{V},indices) where V + +Allocate a vector of type `V` indexable at the incides `indices` +""" +function allocate_vector(::Type{V},indices) where V + T = eltype(V) + n = length(indices) + zeros(T,n) +end + +""" + allocate_in_range(::Type{V},matrix) where V + +Allocate a vector of type `V` in the range of matrix `matrix`. +""" +function allocate_in_range(::Type{V},matrix) where V + indices = 1:size(matrix,1) + allocate_vector(V,indices) +end + +""" + allocate_in_domain(::Type{V},matrix) where V + +Allocate a vector of type `V` in the domain of matrix `matrix`. +""" +function allocate_in_domain(::Type{V},matrix) where V + indices = 1:size(matrix,2) + allocate_vector(V,indices) +end + +""" + fill_entries!(a,v) + +Fill the entries of array `a` with the value `v`. Returns `a`. +""" +function fill_entries!(a,v) + fill!(a,v) + a +end + +""" + copy_entries!(a,b) + +Copy the entries of array `b` into array `a`. Returns `a`. +""" +function copy_entries!(a,b) + if a !== b + copyto!(a,b) + end + a +end + +""" + add_entries!(a,b,combine=+) + +Perform the operation `combine` element-wise in the entries of arrays `a` and `b` +and store the result in array `a`. Returns `a`. +""" +function add_entries!(a,b,combine=+) + @assert length(a) == length(b) + @inbounds for i in eachindex(a) + a[i] = combine(a[i],b[i]) + end + a +end + +""" + scale_entries!(a,v) + +Scale the entries of array `a` with the value `v`. Returns `a`. +""" +function scale_entries!(a,b) + @inbounds for i in eachindex(a) + a[i] = b*a[i] + end + a +end + +# Base.mul! + +""" + muladd!(c,a,b) + +Matrix multiply a*b and add to result to c. Returns c. +""" +function muladd!(c,a,b) + _muladd!(c,a,b) + c +end + +@static if VERSION >= v"1.3" + function _muladd!(c,a,b) + mul!(c,a,b,1,1) + end +else + function _muladd!(c,a,b) + @assert length(c) == size(a,1) + @assert length(b) == size(a,2) + @inbounds for j in 1:size(a,2) + for i in 1:size(a,1) + c[i] += a[i,j]*b[j] + end + end + end +end + diff --git a/src/Algebra/CompressedSparseMatrices.jl b/src/Algebra/CompressedSparseMatrices.jl index 1b7412a3e..4daf79b7d 100644 --- a/src/Algebra/CompressedSparseMatrices.jl +++ b/src/Algebra/CompressedSparseMatrices.jl @@ -1,34 +1,4 @@ -# CompressedSparseMatrices implementation contains: -# -# - Data types: -# + `SparseMatrixCSR`:Compressed Sparse Row (CSR) sparse matrix implementation with Bi-based indexing. -# + `SymSparseMatrixCSR`: Symmetric Compressed Sparse Row sparse matrix implementation with Bi-based indexing. -# -# - Procedures: -# + `sparsecsr`: Analogous method to [`sparse`](@ref) function to build a SparseMatrixCSR. -# + `symsparsecsr`: Analogous method to [`sparse`](@ref) function to build a SymSparseMatrixCSR. -# + `colvals`: Analogous method to [`rowvals`](@ref) to return `colvals` array. -# + `hasrowmajororder`: Return `true` if matrix values are ordered by row. -# + `hascolmajororder`: Return `true` if matrix values are ordered by column. -# + `getptr`: Return the pointer array of a SparseMatrix (`rowptr` or `colptr` depending on the SparseMatrix type) -# + `getindices`: Return the indices array of a SparseMatrix (`rowval` or `colval` depending on the SparseMatrix type) -# -# - Overloaded procedures: -# + `*`: SparseMatrix-Vector product. -# + `mul!`: SparseMatrix-Vector product. -# + `nnz`: Return the number of stored (filled) elements in a sparse array. -# + `nonzeros`: Return `nzval` array. -# + `nzrange`: Return the range of indices for a particular row or column (Depending on the SparseMatrix type) -# + `findnz`: Return a tuple (I, J, V) where I and J are the row and column indices. -# + `rowvals`: Return row indices or raises an error (Depending on the SparseMatrix type) -# + `convert`: Type conversion -# + `count`: Count the number of elements in nonzeros(S) for which predicate pred returns true. - -import Base: convert, size, getindex, show, count, * -import LinearAlgebra: mul! -import SparseArrays: nnz, nonzeros, nzrange, findnz, rowvals - """ hasrowmajororder(::Type{AbstractSparseMatrix}) @@ -157,7 +127,3 @@ function fill_entries!(A::AbstractSparseMatrix{Tv,Ti},v::Number) where {Tv,Ti} nonzeros(A) .= convert(Tv,v) end - -include("SparseMatrixCSC.jl") -include("SparseMatrixCSR.jl") -include("SymSparseMatrixCSR.jl") diff --git a/src/Algebra/MethodsAbstractMatrices.jl b/src/Algebra/MethodsAbstractMatrices.jl deleted file mode 100644 index 16f526c57..000000000 --- a/src/Algebra/MethodsAbstractMatrices.jl +++ /dev/null @@ -1,3 +0,0 @@ -function fill_entries!(J::AbstractArray,v::Number) - J .= convert(eltype(J),v) -end diff --git a/src/Algebra/SparseMatrices.jl b/src/Algebra/SparseMatrices.jl index 2c34d5378..9704c139a 100644 --- a/src/Algebra/SparseMatrices.jl +++ b/src/Algebra/SparseMatrices.jl @@ -1,43 +1,28 @@ -# SparseMatrices implementation contains: -# -# - Data types: -# + `SparseMatrixCSR`:Compressed Sparse Row (CSR) sparse matrix implementation with Bi-based indexing. -# + `SymSparseMatrixCSR`: Symmetric Compressed Sparse Row sparse matrix implementation with Bi-based indexing. -# -# - Procedures: -# + `push_coo!`: Helper function to build COO arrays for further building a SparseMatrix -# + `finalize_coo!`: Finalization of COO arrays building. -# + `sparse_from_coo`: Return a SparseMatrix from COO data given. -# + `add_entry!`: Add an entry given its position and the operation to perform. -# + `fill_entries!`: Fills an existing sparse matrix with a given value. - # Extended Sparse matrix interface """ - sparse_from_coo(::Type{T} where T,args...) - -`args...` are the same as for function `sparse` + sparse_from_coo(::Type{T} where T,I,J,V,m,n) """ -function sparse_from_coo(::Type{T} where T,args...) +function sparse_from_coo(::Type{T} where T,I,J,V,m,n) @abstractmethod end """ - add_entry!(A,v::Number,i::Integer,j::Integer,combine::Function=+) + add_entry!(A,v,i,j,combine=+) Add an entry given its position and the operation to perform. """ -function add_entry!(A,v::Number,i::Integer,j::Integer,combine::Function=+) +function add_entry!(A,v,i,j,combine=+) @abstractmethod end """ - fill_entry!(A,v::Number) + fill_entry!(A::AbstractSparseMatrix,v) -Fills an existing sparse matrix A with a given value v. +Fills the non-zero entries in the sparse matrix A with a given value v. """ -function fill_entries!(A,v::Number) +function fill_entries!(A::AbstractSparseMatrix,v) @abstractmethod end @@ -46,16 +31,16 @@ end Inserts entries in COO vectors for further building a sparse matrix of type T. """ -function push_coo!(::Type{T} where T,I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number) +function push_coo!(::Type{T} where T,I,J,V,ik,jk,vk) @abstractmethod end """ - is_entry_stored(::Type{T} where T,i::Integer,j::Integer) -> Bool + is_entry_stored(::Type{T} where T,i,j) -> Bool Tells if the entry with coordinates `[i,j]` will be stored when calling function push_coo! """ -function is_entry_stored(::Type{T} where T,i::Integer,j::Integer) +function is_entry_stored(::Type{T} where T,i,j) @abstractmethod end @@ -65,20 +50,10 @@ end Check and insert diagonal entries in COO vectors if needed. """ -function finalize_coo!(::Type{T} where T,I::Vector,J::Vector,V::Vector,m::Integer,n::Integer) +function finalize_coo!(::Type{T} where T,I,J,V,m,n) @abstractmethod end -""" -""" -function create_coo_vectors(::Type{M}) where {M} - return (Int[], Int[], Float64[]) -end - -function create_coo_vectors(::Type{M}) where {Tv,Ti,M<:AbstractSparseMatrix{Tv,Ti}} - return (Ti[], Ti[], Tv[]) -end - """ """ function allocate_coo_vectors(::Type{M},n::Integer) where M @@ -89,18 +64,3 @@ function allocate_coo_vectors(::Type{M},n::Integer) where {Tv,Ti,M<:AbstractSpar return (zeros(Ti,n), zeros(Ti,n), zeros(Tv,n)) end -""" -""" -function copy_entries!(a::AbstractMatrix,b::AbstractMatrix) - if a !== b - _copy!(a,b) - end -end - -# We define this, since its not in 1.0 -function _copy!(a,b) - @assert size(a) == size(b) "Array dimension mismatch when copying" - for i in eachindex(a) - a[i] = b[i] - end -end diff --git a/src/Algebra/SparseMatrixCSC.jl b/src/Algebra/SparseMatrixCSC.jl index f60897ecd..08825e8d2 100644 --- a/src/Algebra/SparseMatrixCSC.jl +++ b/src/Algebra/SparseMatrixCSC.jl @@ -61,7 +61,7 @@ function _copy_entries_sparse!(a,b) na = nonzeros(a) nb = nonzeros(b) if na !== nb - _copy!(na,nb) + copyto!(na,nb) end end diff --git a/src/Arrays/Arrays.jl b/src/Arrays/Arrays.jl index 8f18113f9..9074604d9 100644 --- a/src/Arrays/Arrays.jl +++ b/src/Arrays/Arrays.jl @@ -84,8 +84,6 @@ export SubVector export pair_arrays export unpair_arrays -export matvec_muladd! - export lazy_append export lazy_split export AppendedArray diff --git a/src/Arrays/Interface.jl b/src/Arrays/Interface.jl index 64ac7ab76..4be55eb72 100644 --- a/src/Arrays/Interface.jl +++ b/src/Arrays/Interface.jl @@ -335,26 +335,5 @@ function add_to_array!(a::AbstractArray,b::Number,combine=+) end end -""" -""" -function matvec_muladd!(c::AbstractVector,a::AbstractMatrix,b::AbstractVector) - _matvec_muladd!(c,a,b) -end - -@static if VERSION >= v"1.3" - function _matvec_muladd!(c,a,b) - mul!(c,a,b,1,1) - end -else - function _matvec_muladd!(c,a,b) - @assert length(c) == size(a,1) - @assert length(b) == size(a,2) - @inbounds for j in 1:size(a,2) - for i in 1:size(a,1) - c[i] += a[i,j]*b[j] - end - end - end -end diff --git a/src/MultiField/MultiFieldArrays.jl b/src/MultiField/MultiFieldArrays.jl index f7bbd5c1a..181ce738d 100644 --- a/src/MultiField/MultiFieldArrays.jl +++ b/src/MultiField/MultiFieldArrays.jl @@ -184,7 +184,7 @@ function mul!(c::MultiFieldArray{Tc,1},a::MultiFieldArray{Ta,2},b::MultiFieldArr bk = b.blocks[p] q = c.ptrs[ci] ck = c.blocks[q] - matvec_muladd!(ck,ak,bk) + muladd!(ck,ak,bk) end end diff --git a/test/AlgebraTests/AlgebraInterfacesTests.jl b/test/AlgebraTests/AlgebraInterfacesTests.jl new file mode 100644 index 000000000..d4108c9b1 --- /dev/null +++ b/test/AlgebraTests/AlgebraInterfacesTests.jl @@ -0,0 +1,56 @@ +module AlgebraInterfacesTests + +using Gridap.Algebra +using Test + +a = allocate_vector(Vector{Int},1:10) +@test isa(a,Vector{Int}) +@test length(a) == 10 + +a = zeros(4,6) + +b = allocate_in_range(Vector{Int},a) +@test length(b) == 4 + +b = allocate_in_domain(Vector{Int},a) +@test length(b) == 6 + +fill_entries!(b,4) +@test all( b .== 4 ) + +a = rand(6) +b = rand(6) + +copy_entries!(a,b) +@test all( a .== b ) + +a = rand(6) +c = copy(a) +b = rand(6) + +add_entries!(a,b) +@test all( a .== ( c .+ b) ) + +a = rand(6) +c = copy(a) +b = rand(6) + +add_entries!(a,b,-) +@test all( a .== ( c .- b) ) + +a = rand(6) +c = copy(a) +scale_entries!(a,10) +@test all( a .== 10*c) + + +a = rand(4,6) +c = rand(4) +d = copy(c) +b = rand(6) + +muladd!(c,a,b) + +@test all( c .≈ (d .+ a*b ) ) + +end # module diff --git a/test/AlgebraTests/runtests.jl b/test/AlgebraTests/runtests.jl index 32457aa65..90bca8289 100644 --- a/test/AlgebraTests/runtests.jl +++ b/test/AlgebraTests/runtests.jl @@ -2,6 +2,8 @@ module AlgebraTests using Test +@testset "AlgebraInterfaces" begin include("AlgebraInterfacesTests.jl") end + @testset "SparseMatrixCSC" begin include("SparseMatrixCSC.jl") end @testset "SparseMatrixCSR" begin include("SparseMatrixCSR.jl") end From 0b06cf3a8441036029242a0bf68e32d69c9b88e1 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 12:07:03 +0200 Subject: [PATCH 02/17] Simplifying zero_initial_guess --- src/Algebra/LinearSolvers.jl | 8 ++++---- src/Algebra/NonlinearOperators.jl | 22 +++++++--------------- src/FESpaces/FEOperators.jl | 4 +++- 3 files changed, 14 insertions(+), 20 deletions(-) diff --git a/src/Algebra/LinearSolvers.jl b/src/Algebra/LinearSolvers.jl index dbb3a1aa8..2ab8088fe 100644 --- a/src/Algebra/LinearSolvers.jl +++ b/src/Algebra/LinearSolvers.jl @@ -33,13 +33,13 @@ function jacobian(op::AffineOperator,x::AbstractVector) op.matrix end -function zero_initial_guess(::Type{T},op::AffineOperator) where T - n = size(op.matrix,2) - zeros(T,n) +function zero_initial_guess(op::AffineOperator) + x = allocate_in_domain(typeof(op.vector),op.matrix) + fill_entries!(x,zero(eltype(x))) end function allocate_residual(op::AffineOperator,x::AbstractVector) - similar(x,eltype(x),length(op.vector)) + allocate_in_range(typeof(op.vector),op.matrix) end function allocate_jacobian(op::AffineOperator,x::AbstractVector) diff --git a/src/Algebra/NonlinearOperators.jl b/src/Algebra/NonlinearOperators.jl index 31df78930..320b44bd9 100644 --- a/src/Algebra/NonlinearOperators.jl +++ b/src/Algebra/NonlinearOperators.jl @@ -3,7 +3,7 @@ - [`residual!(b::AbstractVector,op::NonlinearOperator,x::AbstractVector)`](@ref) - [`jacobian!(A::AbstractMatrix,op::NonlinearOperator,x::AbstractVector)`](@ref) -- [`zero_initial_guess(::Type{T},op::NonlinearOperator) where T`](@ref) +- [`zero_initial_guess(op::NonlinearOperator)`](@ref) - [`allocate_residual(op::NonlinearOperator,x::AbstractVector)`](@ref) - [`allocate_jacobian(op::NonlinearOperator,x::AbstractVector)`](@ref) @@ -63,18 +63,11 @@ function residual_and_jacobian(op::NonlinearOperator,x::AbstractVector) (b, A) end -""" - zero_initial_guess(::Type{T},op::NonlinearOperator) where T -""" -function zero_initial_guess(::Type{T},op::NonlinearOperator) where T - @abstractmethod -end - """ zero_initial_guess(op::NonlinearOperator) """ function zero_initial_guess(op::NonlinearOperator) - zero_initial_guess(Float64,op) + @abstractmethod end """ @@ -120,8 +113,6 @@ function test_nonlinear_operator( @test pred(b,b1) x0 = zero_initial_guess(op) - x0 = zero_initial_guess(Int,op) - @assert eltype(x0) == Int if jac != nothing nrows, ncols = size(jac) @@ -156,15 +147,16 @@ function jacobian!(A::AbstractMatrix,::NonlinearOperatorMock,x::AbstractVector) A[2,2] = 1 end -function zero_initial_guess(::Type{T},op::NonlinearOperatorMock) where T - x = T[] - allocate_residual(op,x) +function zero_initial_guess(op::NonlinearOperatorMock) + x = allocate_residual(op,Float64[]) + fill!(x,zero(eltype(x))) + x end function allocate_residual(op::NonlinearOperatorMock,x::AbstractVector) T = eltype(x) n = 2 - zeros(T,n) + similar(x,T,n) end function allocate_jacobian(op::NonlinearOperatorMock,x::AbstractVector) diff --git a/src/FESpaces/FEOperators.jl b/src/FESpaces/FEOperators.jl index 617b45f90..47f443672 100644 --- a/src/FESpaces/FEOperators.jl +++ b/src/FESpaces/FEOperators.jl @@ -160,7 +160,9 @@ function residual_and_jacobian(op::AlgebraicOpFromFEOp,x::AbstractVector) residual_and_jacobian(op.feop,u) end -function zero_initial_guess(::Type{T},op::AlgebraicOpFromFEOp) where T +function zero_initial_guess(op::AlgebraicOpFromFEOp) + # TODO + T = Float64 trial = get_trial(op.feop) x = zero_free_values(T,trial) end From b6844065194e21f05ab74b752a2cc1636196c52d Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 12:23:42 +0200 Subject: [PATCH 03/17] Writing operators and solvers in terms of the new algebra interface --- src/Algebra/LinearSolvers.jl | 17 ++++++++++------- src/Algebra/NLSolvers.jl | 6 +++--- src/Algebra/NonlinearSolvers.jl | 3 +-- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/Algebra/LinearSolvers.jl b/src/Algebra/LinearSolvers.jl index 2ab8088fe..1dd42e142 100644 --- a/src/Algebra/LinearSolvers.jl +++ b/src/Algebra/LinearSolvers.jl @@ -20,7 +20,7 @@ get_vector(op::AffineOperator) = op.vector function residual!(b::AbstractVector,op::AffineOperator,x::AbstractVector) mul!(b,op.matrix,x) - b .-= op.vector + add_entries!(b,op.vector,-) b end @@ -36,6 +36,7 @@ end function zero_initial_guess(op::AffineOperator) x = allocate_in_domain(typeof(op.vector),op.matrix) fill_entries!(x,zero(eltype(x))) + x end function allocate_residual(op::AffineOperator,x::AbstractVector) @@ -151,12 +152,12 @@ function solve!(x::AbstractVector, ls::LinearSolver, op::NonlinearOperator, cache::Nothing) - x .= zero(eltype(x)) + fill_entries!(x,zero(eltype(x))) b = residual(op, x) A = jacobian(op, x) ss = symbolic_setup(ls, A) ns = numerical_setup(ss,A) - broadcast!(*,b,b,-1) + scale_entries!(b,-1) solve!(x,ns,b) LinearSolverCache(A,b,ns) end @@ -165,13 +166,13 @@ function solve!(x::AbstractVector, ls::LinearSolver, op::NonlinearOperator, cache) - x .= zero(eltype(x)) + fill_entries!(x,zero(eltype(x))) b = cache.b A = cache.A ns = cache.ns residual!(b, op, x) numerical_setup!(ns,A) - broadcast!(*,b,b,-1) + scale_entries!(b,-1) solve!(x,ns,b) cache end @@ -238,7 +239,8 @@ end function solve!( x::AbstractVector,ns::LUNumericalSetup,b::AbstractVector) y = ns.factors\b # the allocation of y can be avoided - x .= y + copy_entries!(x,y) + x end """ @@ -265,5 +267,6 @@ end function solve!( x::AbstractVector,ns::BackslashNumericalSetup,b::AbstractVector) - copyto!(x, ns.A\b) + copy_entries!(x, ns.A\b) end + diff --git a/src/Algebra/NLSolvers.jl b/src/Algebra/NLSolvers.jl index b8e417365..6e56d475c 100644 --- a/src/Algebra/NLSolvers.jl +++ b/src/Algebra/NLSolvers.jl @@ -54,9 +54,9 @@ function _solve_nr!(x,A,b,dx,ns,nls,op) for nliter in 1:nls.max_nliters # Solve linearized problem - broadcast!(*,b,b,-1) + scale_entries!(b,-1) solve!(dx,ns,b) - broadcast!(+,x,x,dx) + add_entries!(x,dx) # Check convergence for the current residual residual!(b, op, x) @@ -157,7 +157,7 @@ function _nlsolve_with_updated_cache!(x,nls,op,cache) end r = nlsolve(df,x;linsolve=linsolve!,kwargs...) cache.result = r - x[:] .= r.zero + copy_entries!(x,r.zero) end function _new_nlsolve_cache(x0,nls,op) diff --git a/src/Algebra/NonlinearSolvers.jl b/src/Algebra/NonlinearSolvers.jl index c2abf087d..d831ccc93 100644 --- a/src/Algebra/NonlinearSolvers.jl +++ b/src/Algebra/NonlinearSolvers.jl @@ -25,8 +25,7 @@ end """ solve!(x::AbstractVector,nls::NonlinearSolver,op::NonlinearOperator,cache) -Solve using the cache object from a previous solve. If cache === nothing, a new -cache is created a returned. +Solve using the cache object from a previous solve. """ function solve!(x::AbstractVector,nls::NonlinearSolver,op::NonlinearOperator,cache) @abstractmethod From c61b1b5ad765a1a59375f7b1aa536c7d6e946ab5 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 12:42:33 +0200 Subject: [PATCH 04/17] Simplified API of zero_free_values --- src/FESpaces/CLagrangianFESpaces.jl | 4 ++-- src/FESpaces/DirichletFESpaces.jl | 7 +++---- src/FESpaces/ExtendedFESpaces.jl | 4 ++-- src/FESpaces/FEOperators.jl | 4 +--- src/FESpaces/FESpacesInterfaces.jl | 23 +++++++++++----------- src/FESpaces/FESpacesWithLastDofRemoved.jl | 4 ++-- src/FESpaces/TrialFESpaces.jl | 2 +- src/FESpaces/UnconstrainedFESpaces.jl | 7 +++---- src/FESpaces/ZeroMeanFESpaces.jl | 2 +- src/MultiField/MultiFieldFESpaces.jl | 8 ++++---- 10 files changed, 30 insertions(+), 35 deletions(-) diff --git a/src/FESpaces/CLagrangianFESpaces.jl b/src/FESpaces/CLagrangianFESpaces.jl index 810c0e4b5..012571a42 100644 --- a/src/FESpaces/CLagrangianFESpaces.jl +++ b/src/FESpaces/CLagrangianFESpaces.jl @@ -37,8 +37,8 @@ function get_cell_basis(f::CLagrangianFESpace) get_cell_basis(f.space) end -function zero_free_values(::Type{T},f::CLagrangianFESpace) where T - zero_free_values(T,f.space) +function zero_free_values(f::CLagrangianFESpace) + zero_free_values(f.space) end # SingleFieldFESpace interface diff --git a/src/FESpaces/DirichletFESpaces.jl b/src/FESpaces/DirichletFESpaces.jl index 5e879160a..dea555685 100644 --- a/src/FESpaces/DirichletFESpaces.jl +++ b/src/FESpaces/DirichletFESpaces.jl @@ -32,8 +32,8 @@ function num_free_dofs(f::DirichletFESpace) num_dirichlet_dofs(f.space) end -function zero_free_values(::Type{T},f::DirichletFESpace) where T - zeros(T,num_free_dofs(f)) +function zero_free_values(f::DirichletFESpace) + zero_dirichlet_values(f.space) end function get_cell_dofs(f::DirichletFESpace) @@ -45,8 +45,7 @@ function num_dirichlet_dofs(f::DirichletFESpace) end function zero_dirichlet_values(f::DirichletFESpace) - T = Float64 - zero_free_values(T,f.space) + zero_free_values(f.space) end function num_dirichlet_tags(f::DirichletFESpace) diff --git a/src/FESpaces/ExtendedFESpaces.jl b/src/FESpaces/ExtendedFESpaces.jl index f9eede12c..fcd32b126 100644 --- a/src/FESpaces/ExtendedFESpaces.jl +++ b/src/FESpaces/ExtendedFESpaces.jl @@ -236,8 +236,8 @@ function num_free_dofs(f::ExtendedFESpace) num_free_dofs(f.space) end -function zero_free_values(::Type{T},f::ExtendedFESpace) where T - zeros(T,num_free_dofs(f)) +function zero_free_values(f::ExtendedFESpace) + zero_free_values(f.space) end function num_dirichlet_dofs(f::ExtendedFESpace) diff --git a/src/FESpaces/FEOperators.jl b/src/FESpaces/FEOperators.jl index 47f443672..cc225d6b4 100644 --- a/src/FESpaces/FEOperators.jl +++ b/src/FESpaces/FEOperators.jl @@ -161,8 +161,6 @@ function residual_and_jacobian(op::AlgebraicOpFromFEOp,x::AbstractVector) end function zero_initial_guess(op::AlgebraicOpFromFEOp) - # TODO - T = Float64 trial = get_trial(op.feop) - x = zero_free_values(T,trial) + x = zero_free_values(trial) end diff --git a/src/FESpaces/FESpacesInterfaces.jl b/src/FESpaces/FESpacesInterfaces.jl index 68fad9a7e..f42a8e9b3 100644 --- a/src/FESpaces/FESpacesInterfaces.jl +++ b/src/FESpaces/FESpacesInterfaces.jl @@ -20,15 +20,11 @@ end """ abstract type FESpace <: GridapType end -""" -""" -function num_free_dofs(f::FESpace) - @abstractmethod -end +# Minimal FE interface (used by FEOperator) """ """ -function get_cell_basis(f::FESpace) +function num_free_dofs(f::FESpace) @abstractmethod end @@ -44,12 +40,8 @@ end """ """ -function zero_free_values(::Type{T},fs::FESpace) where T - @abstractmethod -end - function zero_free_values(fs::FESpace) - zero_free_values(Float64,fs) + @abstractmethod end """ @@ -59,6 +51,14 @@ function Base.zero(f::FESpace) FEFunction(f,free_values) end +# Extended FEInterface used by FEOperatorFromTerms and Assemblers + +""" +""" +function get_cell_basis(f::FESpace) + @abstractmethod +end + """ """ function constraint_style(::Type{<:FESpace}) @@ -89,7 +89,6 @@ end """ function test_fe_space(f::FESpace) free_values = zero_free_values(f) - @test eltype(zero_free_values(Int,f)) == Int fe_function = FEFunction(f,free_values) test_fe_function(fe_function) fe_basis = get_cell_basis(f) diff --git a/src/FESpaces/FESpacesWithLastDofRemoved.jl b/src/FESpaces/FESpacesWithLastDofRemoved.jl index 9711574cf..fb4d3f7c5 100644 --- a/src/FESpaces/FESpacesWithLastDofRemoved.jl +++ b/src/FESpaces/FESpacesWithLastDofRemoved.jl @@ -25,8 +25,8 @@ function num_free_dofs(f::FESpaceWithLastDofRemoved) num_free_dofs(f.space) - 1 end -function zero_free_values(::Type{T},f::FESpaceWithLastDofRemoved) where T - zeros(T,num_free_dofs(f)) +function zero_free_values(f::FESpaceWithLastDofRemoved) + zeros(num_free_dofs(f)) end function get_cell_dofs(f::FESpaceWithLastDofRemoved) diff --git a/src/FESpaces/TrialFESpaces.jl b/src/FESpaces/TrialFESpaces.jl index 97c5f18ff..2e468ac53 100644 --- a/src/FESpaces/TrialFESpaces.jl +++ b/src/FESpaces/TrialFESpaces.jl @@ -78,7 +78,7 @@ get_cell_dof_basis(f::TrialFESpace) = get_cell_dof_basis(f.space) num_free_dofs(f::TrialFESpace) = num_free_dofs(f.space) -zero_free_values(::Type{T},f::TrialFESpace) where T = zero_free_values(T,f.space) +zero_free_values(f::TrialFESpace) = zero_free_values(f.space) get_cell_dofs(f::TrialFESpace) = get_cell_dofs(f.space) diff --git a/src/FESpaces/UnconstrainedFESpaces.jl b/src/FESpaces/UnconstrainedFESpaces.jl index 8cdf8b046..cd83635af 100644 --- a/src/FESpaces/UnconstrainedFESpaces.jl +++ b/src/FESpaces/UnconstrainedFESpaces.jl @@ -54,8 +54,8 @@ function get_cell_basis(f::UnconstrainedFESpace) f.cell_basis end -function zero_free_values(::Type{T},f::UnconstrainedFESpace) where T - zeros(T,num_free_dofs(f)) +function zero_free_values(f::UnconstrainedFESpace) + zeros(num_free_dofs(f)) end # SingleFieldFESpace interface @@ -77,8 +77,7 @@ function num_dirichlet_tags(f::UnconstrainedFESpace) end function zero_dirichlet_values(f::UnconstrainedFESpace) - T = Float64 # TODO - zeros(T,num_dirichlet_dofs(f)) + zeros(num_dirichlet_dofs(f)) end function get_dirichlet_dof_tag(f::UnconstrainedFESpace) diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl index 8aed5ade8..3ef825d28 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -95,7 +95,7 @@ get_cell_dof_basis(f::ZeroMeanFESpace) = get_cell_dof_basis(f.space) num_free_dofs(f::ZeroMeanFESpace) = num_free_dofs(f.space) -zero_free_values(::Type{T},f::ZeroMeanFESpace) where T = zero_free_values(T,f.space) +zero_free_values(f::ZeroMeanFESpace) = zero_free_values(f.space) get_cell_dofs(f::ZeroMeanFESpace) = get_cell_dofs(f.space) diff --git a/src/MultiField/MultiFieldFESpaces.jl b/src/MultiField/MultiFieldFESpaces.jl index 2f8dbb639..695d55906 100644 --- a/src/MultiField/MultiFieldFESpaces.jl +++ b/src/MultiField/MultiFieldFESpaces.jl @@ -94,13 +94,13 @@ function EvaluationFunction(spaces::Vector{<:SingleFieldFESpace}, free_values) EvaluationFunction(f,free_values) end -function zero_free_values(::Type{T},fs::MultiFieldFESpace) where T - zeros(T,num_free_dofs(fs)) +function zero_free_values(fs::MultiFieldFESpace) + zeros(num_free_dofs(fs)) end -function zero_free_values(::Type{T},spaces::Vector{<:SingleFieldFESpace}) where T +function zero_free_values(spaces::Vector{<:SingleFieldFESpace}) f = MultiFieldFESpace(spaces) - zero_free_values(T,f) + zero_free_values(f) end # API for multi field case From e21242eada6c95a1443189f27e38d054d44340ae Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 12:59:16 +0200 Subject: [PATCH 05/17] Moving to v0.10.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index fca9ac916..cdff7e492 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gridap" uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" authors = ["Santiago Badia ", "Francesc Verdugo "] -version = "0.9.2" +version = "0.10.0" [deps] BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" From def8485c6077c9ccbee7211186bdf971a6ad624e Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 14:46:36 +0200 Subject: [PATCH 06/17] moving allocate_{vector,matrix,matrix_and_vector} to Gridap.Algebra --- src/Algebra/Algebra.jl | 2 ++ src/Algebra/AlgebraInterfaces.jl | 3 +++ src/FESpaces/FESpaces.jl | 9 ++++++--- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/Algebra/Algebra.jl b/src/Algebra/Algebra.jl index 4e1627f68..f3a4d4ba4 100644 --- a/src/Algebra/Algebra.jl +++ b/src/Algebra/Algebra.jl @@ -18,6 +18,8 @@ import LinearAlgebra: mul! import SparseArrays: nnz, nonzeros, nzrange, findnz, rowvals export allocate_vector +export allocate_matrix +export allocate_matrix_and_vector export allocate_in_domain export allocate_in_range export add_entries! diff --git a/src/Algebra/AlgebraInterfaces.jl b/src/Algebra/AlgebraInterfaces.jl index f9755a000..3535109a3 100644 --- a/src/Algebra/AlgebraInterfaces.jl +++ b/src/Algebra/AlgebraInterfaces.jl @@ -1,4 +1,7 @@ +function allocate_matrix end +function allocate_matrix_and_vector end + """ allocate_vector(::Type{V},indices) where V diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index 3ca28154c..267e094ea 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -74,6 +74,9 @@ import Gridap.Algebra: get_matrix import Gridap.Algebra: get_vector import Gridap.Algebra: solve! import Gridap.Algebra: solve +import Gridap.Algebra: allocate_vector +import Gridap.Algebra: allocate_matrix +import Gridap.Algebra: allocate_matrix_and_vector export FEFunctionStyle export is_a_fe_function @@ -102,16 +105,16 @@ export test_fe_space export Assembler export get_test export get_trial -export allocate_matrix export assemble_matrix! export assemble_matrix_add! export assemble_matrix -export allocate_vector export assemble_vector! export assemble_vector -export allocate_matrix_and_vector export assemble_matrix_and_vector! export assemble_matrix_and_vector +export allocate_vector +export allocate_matrix +export allocate_matrix_and_vector export test_assembler export SingleFieldFESpace From a505c58695990f0b0bed750c3b8090dfbdc559f7 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 5 May 2020 14:57:55 +0200 Subject: [PATCH 07/17] More work in Algebra interface --- src/Algebra/AlgebraInterfaces.jl | 34 ++++++++++++++++++--- src/Algebra/SparseMatrices.jl | 11 +------ src/Algebra/SparseMatrixCSC.jl | 4 +-- src/Algebra/SparseMatrixCSR.jl | 4 +-- src/Algebra/SymSparseMatrixCSR.jl | 4 +-- test/AlgebraTests/AlgebraInterfacesTests.jl | 4 +++ 6 files changed, 40 insertions(+), 21 deletions(-) diff --git a/src/Algebra/AlgebraInterfaces.jl b/src/Algebra/AlgebraInterfaces.jl index 3535109a3..7adad261f 100644 --- a/src/Algebra/AlgebraInterfaces.jl +++ b/src/Algebra/AlgebraInterfaces.jl @@ -8,8 +8,12 @@ function allocate_matrix_and_vector end Allocate a vector of type `V` indexable at the incides `indices` """ function allocate_vector(::Type{V},indices) where V - T = eltype(V) n = length(indices) + allocate_vector(V,n) +end + +function allocate_vector(::Type{V},n::Integer) where V + T = eltype(V) zeros(T,n) end @@ -19,8 +23,8 @@ end Allocate a vector of type `V` in the range of matrix `matrix`. """ function allocate_in_range(::Type{V},matrix) where V - indices = 1:size(matrix,1) - allocate_vector(V,indices) + n = size(matrix,1) + allocate_vector(V,n) end """ @@ -29,8 +33,8 @@ end Allocate a vector of type `V` in the domain of matrix `matrix`. """ function allocate_in_domain(::Type{V},matrix) where V - indices = 1:size(matrix,2) - allocate_vector(V,indices) + n = size(matrix,2) + allocate_vector(V,n) end """ @@ -69,6 +73,26 @@ function add_entries!(a,b,combine=+) a end +""" + add_entry!(A,v,i,j,combine=+) + +Add an entry given its position and the operation to perform. +""" +function add_entry!(A,v,i::Integer,j::Integer,combine=+) + aij = A[i,j] + A[i,j] = combine(aij,v) +end + +""" + add_entry!(A,v,i,combine=+) + +Add an entry given its position and the operation to perform. +""" +function add_entry!(A,v,i::Integer,combine=+) + ai = A[i] + A[i] = combine(ai,v) +end + """ scale_entries!(a,v) diff --git a/src/Algebra/SparseMatrices.jl b/src/Algebra/SparseMatrices.jl index 9704c139a..b44eb1706 100644 --- a/src/Algebra/SparseMatrices.jl +++ b/src/Algebra/SparseMatrices.jl @@ -9,16 +9,7 @@ function sparse_from_coo(::Type{T} where T,I,J,V,m,n) end """ - add_entry!(A,v,i,j,combine=+) - -Add an entry given its position and the operation to perform. -""" -function add_entry!(A,v,i,j,combine=+) - @abstractmethod -end - -""" - fill_entry!(A::AbstractSparseMatrix,v) + fill_entries!(A::AbstractSparseMatrix,v) Fills the non-zero entries in the sparse matrix A with a given value v. """ diff --git a/src/Algebra/SparseMatrixCSC.jl b/src/Algebra/SparseMatrixCSC.jl index 08825e8d2..33ef38bbb 100644 --- a/src/Algebra/SparseMatrixCSC.jl +++ b/src/Algebra/SparseMatrixCSC.jl @@ -1,6 +1,6 @@ -function sparse_from_coo(::Type{<:SparseMatrixCSC}, args...) - sparse(args...) +function sparse_from_coo(::Type{<:SparseMatrixCSC}, I,J,V,m,n) + sparse(I,J,V,m,n) end @inline function push_coo!( diff --git a/src/Algebra/SparseMatrixCSR.jl b/src/Algebra/SparseMatrixCSR.jl index 74a165325..421019e92 100644 --- a/src/Algebra/SparseMatrixCSR.jl +++ b/src/Algebra/SparseMatrixCSR.jl @@ -163,8 +163,8 @@ function findnz(S::SparseMatrixCSR{Bi,Tv,Ti}) where {Bi,Tv,Ti} return (I, J, V) end -function sparse_from_coo(M::Type{<:SparseMatrixCSR}, args...) - sparsecsr(M, args...) +function sparse_from_coo(M::Type{<:SparseMatrixCSR}, I,J,V,m,n) + sparsecsr(M, I,J,V,m,n) end function push_coo!(::Type{<:SparseMatrixCSR}, diff --git a/src/Algebra/SymSparseMatrixCSR.jl b/src/Algebra/SymSparseMatrixCSR.jl index a28885f50..be5b1f58f 100644 --- a/src/Algebra/SymSparseMatrixCSR.jl +++ b/src/Algebra/SymSparseMatrixCSR.jl @@ -11,8 +11,8 @@ struct SymSparseMatrixCSR{Bi,T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti} uppertrian :: SparseMatrixCSR{Bi,T,Ti} end -function sparse_from_coo(M::Type{<:SymSparseMatrixCSR}, args...) - symsparsecsr(M, args...) +function sparse_from_coo(M::Type{<:SymSparseMatrixCSR}, I,J,V,m,n) + symsparsecsr(M, I,J,V,m,n) end # SparseMatrix interface implementation diff --git a/test/AlgebraTests/AlgebraInterfacesTests.jl b/test/AlgebraTests/AlgebraInterfacesTests.jl index d4108c9b1..907289798 100644 --- a/test/AlgebraTests/AlgebraInterfacesTests.jl +++ b/test/AlgebraTests/AlgebraInterfacesTests.jl @@ -3,6 +3,10 @@ module AlgebraInterfacesTests using Gridap.Algebra using Test +a = allocate_vector(Vector{Int},10) +@test isa(a,Vector{Int}) +@test length(a) == 10 + a = allocate_vector(Vector{Int},1:10) @test isa(a,Vector{Int}) @test length(a) == 10 From d0c748c2b678a812e9c6ad498fe531e5f846fbfc Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 6 May 2020 10:24:33 +0200 Subject: [PATCH 08/17] Added AssemblyStrategy and changed test/trial in single field assembler --- src/Algebra/NonlinearOperators.jl | 2 + src/FESpaces/AffineFEOperators.jl | 11 +- src/FESpaces/Assemblers.jl | 46 +++ src/FESpaces/FEOperatorsFromTerms.jl | 9 +- src/FESpaces/SparseMatrixAssemblers.jl | 331 ++++++++++-------- src/FESpaces/ZeroMeanFESpaces.jl | 2 +- test/FESpacesTests/AffineFEOperatorsTests.jl | 2 +- test/FESpacesTests/FESolversTests.jl | 2 +- test/FESpacesTests/FETermsTests.jl | 2 +- .../SparseMatrixAssemblersTests.jl | 2 +- 10 files changed, 253 insertions(+), 156 deletions(-) diff --git a/src/Algebra/NonlinearOperators.jl b/src/Algebra/NonlinearOperators.jl index 320b44bd9..a91625998 100644 --- a/src/Algebra/NonlinearOperators.jl +++ b/src/Algebra/NonlinearOperators.jl @@ -138,6 +138,7 @@ struct NonlinearOperatorMock <: NonlinearOperator end function residual!(b::AbstractVector,::NonlinearOperatorMock,x::AbstractVector) b[1] = (x[1]-2)*(x[1]-1) b[2] = x[2]-3 + b end function jacobian!(A::AbstractMatrix,::NonlinearOperatorMock,x::AbstractVector) @@ -145,6 +146,7 @@ function jacobian!(A::AbstractMatrix,::NonlinearOperatorMock,x::AbstractVector) A[1,2] = 0 A[2,1] = 0 A[2,2] = 1 + A end function zero_initial_guess(op::NonlinearOperatorMock) diff --git a/src/FESpaces/AffineFEOperators.jl b/src/FESpaces/AffineFEOperators.jl index f27964182..e6ca9f6e9 100644 --- a/src/FESpaces/AffineFEOperators.jl +++ b/src/FESpaces/AffineFEOperators.jl @@ -40,15 +40,22 @@ function AffineFEOperator(trial::FESpace,test::FESpace,assem::Assembler,terms::A end function AffineFEOperator(trial::FESpace,test::FESpace,terms::AffineFETerm...) - assem = SparseMatrixAssembler(test,trial) + assem = SparseMatrixAssembler(trial,test) AffineFEOperator(trial,test,assem,terms...) end -function AffineFEOperator(mat::Type{<:AbstractSparseMatrix},trial::FESpace,test::FESpace,terms::AffineFETerm...) +function AffineFEOperator( + mat::Type{<:AbstractSparseMatrix},trial::FESpace,test::FESpace,terms::AffineFETerm...) assem = SparseMatrixAssembler(mat,trial,test) AffineFEOperator(trial,test,assem,terms...) end +function AffineFEOperator( + mat::Type{<:AbstractSparseMatrix},vec::Type{<:AbstractVector},trial::FESpace,test::FESpace,terms::AffineFETerm...) + assem = SparseMatrixAssembler(mat,vec,trial,test) + AffineFEOperator(trial,test,assem,terms...) +end + get_test(feop::AffineFEOperator) = feop.test get_trial(feop::AffineFEOperator) = feop.trial diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index 26434740a..cb5188862 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -1,4 +1,42 @@ +""" +""" +abstract type AssemblyStrategy end + +""" +""" +function row_map(a::AssemblyStrategy,row) + @abstractmethod +end + +""" +""" +function col_map(a::AssemblyStrategy,col) + @abstractmethod +end + +""" +""" +function row_mask(a::AssemblyStrategy,row) + @abstractmethod +end + +""" +""" +function col_mask(a::AssemblyStrategy,col) + @abstractmethod +end + +struct DefaultAssemblyStrategy <: AssemblyStrategy end + +row_map(a::DefaultAssemblyStrategy,row) = row + +col_map(a::DefaultAssemblyStrategy,col) = col + +row_mask(a::DefaultAssemblyStrategy,row) = true + +col_mask(a::DefaultAssemblyStrategy,col) = true + """ """ abstract type Assembler <: GridapType end @@ -15,6 +53,12 @@ function get_trial(a::Assembler) @abstractmethod end +""" +""" +function get_assembly_strategy(a::Assembler) + @abstractmethod +end + """ """ function allocate_matrix(a::Assembler,cellidsrows,cellidscols) @@ -143,5 +187,7 @@ function test_assembler(a::Assembler,matvecdata,matdata,vecdata) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) @test num_free_dofs(test_fesp) == length(b) + strategy = get_assembly_strategy(a) + @test isa(strategy,AssemblyStrategy) end diff --git a/src/FESpaces/FEOperatorsFromTerms.jl b/src/FESpaces/FEOperatorsFromTerms.jl index 2747f1820..436909d92 100644 --- a/src/FESpaces/FEOperatorsFromTerms.jl +++ b/src/FESpaces/FEOperatorsFromTerms.jl @@ -18,12 +18,17 @@ end """ """ function FEOperator(trial::FESpace,test::FESpace,terms::FETerm...) - assem = SparseMatrixAssembler(test,trial) + assem = SparseMatrixAssembler(trial,test) FEOperator(trial,test,assem,terms...) end function FEOperator(mat::Type{<:AbstractSparseMatrix},trial::FESpace,test::FESpace,terms::FETerm...) - assem = SparseMatrixAssembler(mat,test,trial) + assem = SparseMatrixAssembler(mat,trial,test) + FEOperator(trial,test,assem,terms...) +end + +function FEOperator(mat::Type{<:AbstractSparseMatrix},vec::Type{<:AbstractVector},trial::FESpace,test::FESpace,terms::FETerm...) + assem = SparseMatrixAssembler(mat,vec,trial,test) FEOperator(trial,test,assem,terms...) end diff --git a/src/FESpaces/SparseMatrixAssemblers.jl b/src/FESpaces/SparseMatrixAssemblers.jl index d665b5eb3..e595ef03a 100644 --- a/src/FESpaces/SparseMatrixAssemblers.jl +++ b/src/FESpaces/SparseMatrixAssemblers.jl @@ -1,51 +1,81 @@ -struct SparseMatrixAssembler{E} <: Assembler - matrix_type::Type{E} - test::SingleFieldFESpace +struct SparseMatrixAssembler{M,V} <: Assembler + matrix_type::Type{M} + vector_type::Type{V} trial::SingleFieldFESpace + test::SingleFieldFESpace + strategy::AssemblyStrategy +end + +function SparseMatrixAssembler(mat::Type,vec::Type,trial::SingleFieldFESpace,test::SingleFieldFESpace) + strategy = DefaultAssemblyStrategy() + SparseMatrixAssembler(mat,vec,trial,test,strategy) +end + +function SparseMatrixAssembler(mat::Type,trial::SingleFieldFESpace,test::SingleFieldFESpace) + strategy = DefaultAssemblyStrategy() + SparseMatrixAssembler(mat,Vector{Float64},trial,test,strategy) end """ """ -function SparseMatrixAssembler(test::SingleFieldFESpace,trial::SingleFieldFESpace) - SparseMatrixAssembler(SparseMatrixCSC{Float64,Int},test,trial) +function SparseMatrixAssembler(trial::SingleFieldFESpace,test::SingleFieldFESpace) + matrix_type = SparseMatrixCSC{Float64,Int} + vector_type = Vector{Float64} + strategy = DefaultAssemblyStrategy() + SparseMatrixAssembler(matrix_type,vector_type,trial,test,strategy) end get_test(a::SparseMatrixAssembler) = a.test get_trial(a::SparseMatrixAssembler) = a.trial +get_assembly_strategy(a::SparseMatrixAssembler) = a.strategy + function allocate_vector(a::SparseMatrixAssembler,term_to_cellidsrows) - zero_free_values(a.test) + n = num_free_dofs(a.test) + allocate_vector(a.vector_type,n) end function assemble_vector!(b,a::SparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) + fill_entries!(b,zero(eltype(b))) + assemble_vector_add!(b,a,term_to_cellvec,term_to_cellidsrows) +end + +#TODO corresponding abstract method +function assemble_vector_add!(b,a::SparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) celldofs = get_cell_dofs(a.test) - fill!(b,zero(eltype(b))) for (cellvec, cellids) in zip(term_to_cellvec,term_to_cellidsrows) rows = reindex(celldofs,cellids) vals = apply_constraints_vector(a.test,cellvec,cellids) rows_cache = array_cache(rows) vals_cache = array_cache(vals) - _assemble_vector!(b,vals_cache,rows_cache,vals,rows) + _assemble_vector!(b,vals_cache,rows_cache,vals,rows,a.strategy) end b end -function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows) +function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows,strategy) @assert length(cell_vals) == length(cell_rows) for cell in 1:length(cell_rows) rows = getindex!(rows_cache,cell_rows,cell) vals = getindex!(vals_cache,cell_vals,cell) for (i,gid) in enumerate(rows) - if gid > 0 - vec[gid] += vals[i] + if gid > 0 && row_mask(strategy,gid) + _gid = row_map(strategy,gid) + add_entry!(vec,vals[i],_gid) end end end end -function allocate_matrix(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +#TODO +function count_matrix_nnz_coo(a::Assembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) + count_matrix_nnz_coo(a,term_to_cellidsrows, term_to_cellidscols) +end + +#TODO +function count_matrix_nnz_coo(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) n = 0 @@ -55,33 +85,24 @@ function allocate_matrix(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_c rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) @assert length(cell_cols) == length(cell_rows) - n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols) + n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) end - I, J, V = allocate_coo_vectors(a.matrix_type,n) - nini = 0 - for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols) - end - num_rows = num_free_dofs(a.test) - num_cols = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,num_rows,num_cols) - sparse_from_coo(a.matrix_type,I,J,V,num_rows,num_cols) + + n end -@noinline function _count_matrix_entries(::Type{M},rows_cache,cols_cache,cell_rows,cell_cols) where M +@noinline function _count_matrix_entries(::Type{M},rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M n = 0 for cell in 1:length(cell_cols) rows = getindex!(rows_cache,cell_rows,cell) cols = getindex!(cols_cache,cell_cols,cell) for gidcol in cols - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for gidrow in rows - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) n += 1 end end @@ -92,19 +113,40 @@ end n end -@noinline function _allocate_matrix!(a::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols) where M +#TODO +function fill_matrix_coo_symbolic!(I,J,a::Assembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) + fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows, term_to_cellidscols) +end + +#TODO +function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) + celldofs_rows = get_cell_dofs(a.test) + celldofs_cols = get_cell_dofs(a.trial) + nini = 0 + for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) + cell_rows = reindex(celldofs_rows,cellidsrows) + cell_cols = reindex(celldofs_cols,cellidscols) + rows_cache = array_cache(cell_rows) + cols_cache = array_cache(cell_cols) + nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) + end +end + +@noinline function _allocate_matrix!(a::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M n = nini for cell in 1:length(cell_cols) rows = getindex!(rows_cache,cell_rows,cell) cols = getindex!(cols_cache,cell_cols,cell) for gidcol in cols - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for gidrow in rows - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) n += 1 - @inbounds I[n] = gidrow - @inbounds J[n] = gidcol + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol end end end @@ -114,6 +156,17 @@ end n end +# TODO Perhapswe can move this to the abstract interface +function allocate_matrix(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) + n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + I,J,V = allocate_coo_vectors(a.matrix_type,n) + fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows,term_to_cellidscols) + m = num_free_dofs(a.test) + n = num_free_dofs(a.trial) + finalize_coo!(a.matrix_type,I,J,V,m,n) + sparse_from_coo(a.matrix_type,I,J,V,m,n) +end + function assemble_matrix!( mat,a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) z = zero(eltype(mat)) @@ -123,6 +176,7 @@ end function assemble_matrix_add!( mat,a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) @@ -133,12 +187,12 @@ function assemble_matrix_add!( rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) vals_cache = array_cache(cellmat) - _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) + _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols,a.strategy) end mat end -function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) +function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) @assert length(cell_cols) == length(cell_rows) @assert length(cell_vals) == length(cell_rows) for cell in 1:length(cell_cols) @@ -146,11 +200,13 @@ function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_r cols = getindex!(cols_cache,cell_cols,cell) vals = getindex!(vals_cache,cell_vals,cell) for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) v = vals[i,j] - add_entry!(mat,v,gidrow,gidcol) + add_entry!(mat,v,_gidrow,_gidcol) end end end @@ -160,53 +216,56 @@ end function assemble_matrix( a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + + n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + I,J,V = allocate_coo_vectors(a.matrix_type,n) + fill_matrix_coo_numeric!(I,J,V,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + m = num_free_dofs(a.test) + n = num_free_dofs(a.trial) + finalize_coo!(a.matrix_type,I,J,V,m,n) + sparse_from_coo(a.matrix_type,I,J,V,m,n) +end + +function fill_matrix_coo_numeric!( + I,J,V,a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) + + nini = n celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - n = 0 - for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - @assert length(cell_cols) == length(cell_rows) - n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols) - end - I, J, V = allocate_coo_vectors(a.matrix_type,n) - nini = 0 for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) - cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) + cell_vals = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) - vals_cache = array_cache(cellmat) - @assert length(cell_cols) == length(cell_rows) - @assert length(cellmat) == length(cell_rows) - nini = _assemble_matrix_fill!(a.matrix_type,nini,I,J,V,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) + vals_cache = array_cache(cell_vals) + nini = _fill_matrix!( + a.matrix_type,nini,I,J,V,rows_cache,cols_cache,vals_cache,cell_rows,cell_cols,cell_vals,a.strategy) end - num_rows = num_free_dofs(a.test) - num_cols = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,num_rows,num_cols) - sparse_from_coo(a.matrix_type,I,J,V,num_rows,num_cols) + + nini end -@noinline function _assemble_matrix_fill!(::Type{M},nini,I,J,V,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) where M +@noinline function _fill_matrix!( + a::Type{M},nini,I,J,V,rows_cache,cols_cache,vals_cache,cell_rows,cell_cols,cell_vals,strategy) where M + n = nini for cell in 1:length(cell_cols) rows = getindex!(rows_cache,cell_rows,cell) cols = getindex!(cols_cache,cell_cols,cell) vals = getindex!(vals_cache,cell_vals,cell) for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) n += 1 - @inbounds v = vals[i,j] - @inbounds I[n] = gidrow - @inbounds J[n] = gidcol - @inbounds V[n] = v + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol + @inbounds V[n] = vals[i,j] end end end @@ -217,9 +276,15 @@ end end function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, matvecdata, matdata, vecdata) - z = zero(eltype(A)) - fill_entries!(A,z) - fill!(b,zero(eltype(b))) + fill_entries!(A,zero(eltype(A))) + fill_entries!(b,zero(eltype(b))) + assemble_matrix_and_vector_add!(A,b,a,matvecdata, matdata, vecdata) + A, b +end + +function assemble_matrix_and_vector_add!( + A,b,a::SparseMatrixAssembler, matvecdata, matdata, vecdata) + celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) @@ -231,32 +296,14 @@ function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, matvecdata, m rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) vals_cache = array_cache(cellmatvec) - _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols) + _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - - for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) - cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - vals_cache = array_cache(cellmat) - _assemble_matrix!(A,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) - end - - for (cellvec, cellids) in zip(vecdata...) - rows = reindex(celldofs_rows,cellids) - vals = apply_constraints_vector(a.test,cellvec,cellids) - rows_cache = array_cache(rows) - vals_cache = array_cache(vals) - _assemble_vector!(b,vals_cache,rows_cache,vals,rows) - end - + assemble_matrix_add!(A,a,matdata...) + assemble_vector_add!(b,a,vecdata...) A, b end -function _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) +function _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) @assert length(cell_cols) == length(cell_rows) @assert length(cell_vals) == length(cell_rows) for cell in 1:length(cell_cols) @@ -265,43 +312,50 @@ function _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cell_ vals = getindex!(vals_cache,cell_vals,cell) matvals, vecvals = vals for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) v = matvals[i,j] - add_entry!(A,v,gidrow,gidcol) + add_entry!(A,v,_gidrow,_gidcol) end end end end for (i,gidrow) in enumerate(rows) - if gidrow > 0 + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) bi = vecvals[i] - b[gidrow] += bi + b[_gidrow] += bi end end end end function assemble_matrix_and_vector( a::SparseMatrixAssembler, matvecdata, matdata, vecdata) - celldofs_rows = get_cell_dofs(a.test) - celldofs_cols = get_cell_dofs(a.trial) term_to_cellidsrows, term_to_cellidscols, = _rearange_cell_ids(matvecdata,matdata,vecdata) + n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + I,J,V = allocate_coo_vectors(a.matrix_type,n) + b = allocate_vector(a,vecdata...) - n = 0 - for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - @assert length(cell_cols) == length(cell_rows) - n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols) - end + fill_matrix_and_vector_coo_numeric!(I,J,V,b,a, matvecdata, matdata, vecdata) - I, J, V = allocate_coo_vectors(a.matrix_type,n) - b = zero_free_values(a.test) - nini = 0 + m = num_free_dofs(a.test) + n = num_free_dofs(a.trial) + finalize_coo!(a.matrix_type,I,J,V,m,n) + A = sparse_from_coo(a.matrix_type,I,J,V,m,n) + + A, b +end + +function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) + + nini = n + + celldofs_rows = get_cell_dofs(a.test) + celldofs_cols = get_cell_dofs(a.trial) for (cellmatvec_rc,cellidsrows,cellidscols) in zip(matvecdata...) cell_rows = reindex(celldofs_rows,cellidsrows) @@ -313,39 +367,18 @@ function assemble_matrix_and_vector( a::SparseMatrixAssembler, matvecdata, matda vals_cache = array_cache(cellmatvec) @assert length(cell_cols) == length(cell_rows) @assert length(cellmatvec) == length(cell_rows) - nini = _assemble_matrix_and_vector_fill!(a.matrix_type,nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols) - end - - for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) - cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - vals_cache = array_cache(cellmat) - @assert length(cell_cols) == length(cell_rows) - @assert length(cellmat) == length(cell_rows) - nini = _assemble_matrix_fill!(a.matrix_type,nini,I,J,V,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) + nini = _assemble_matrix_and_vector_fill!( + a.matrix_type,nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - for (cellvec, cellids) in zip(vecdata...) - rows = reindex(celldofs_rows,cellids) - vals = apply_constraints_vector(a.test,cellvec,cellids) - rows_cache = array_cache(rows) - vals_cache = array_cache(vals) - _assemble_vector!(b,vals_cache,rows_cache,vals,rows) - end - - num_rows = num_free_dofs(a.test) - num_cols = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,num_rows,num_cols) - A = sparse_from_coo(a.matrix_type,I,J,V,num_rows,num_cols) + nini = fill_matrix_coo_numeric!(I,J,V,a,matdata...,nini) + assemble_vector_add!(b,a,vecdata...) - (A, b) + nini end -@noinline function _assemble_matrix_and_vector_fill!(::Type{M},nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) where M +@noinline function _assemble_matrix_and_vector_fill!( + ::Type{M},nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) where M n = nini for cell in 1:length(cell_cols) rows = getindex!(rows_cache,cell_rows,cell) @@ -353,14 +386,16 @@ end vals = getindex!(vals_cache,cell_vals,cell) matvals, vecvals = vals for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) if is_entry_stored(M,gidrow,gidcol) n += 1 @inbounds v = matvals[i,j] - @inbounds I[n] = gidrow - @inbounds J[n] = gidcol + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol @inbounds V[n] = v end end @@ -368,11 +403,13 @@ end end end for (i,gidrow) in enumerate(rows) - if gidrow > 0 + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) bi = vecvals[i] - b[gidrow] += bi + b[_gidrow] += bi end end end n end + diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl index 3ef825d28..fc12ce1d2 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -27,7 +27,7 @@ end function _setup_vols(V,trian,quad) U = TrialFESpace(V) - assem = SparseMatrixAssembler(V,U) + assem = SparseMatrixAssembler(U,V) bh = get_cell_basis(V) bh_trian = restrict(bh,trian) cellvec = integrate(bh_trian,trian,quad) diff --git a/test/FESpacesTests/AffineFEOperatorsTests.jl b/test/FESpacesTests/AffineFEOperatorsTests.jl index 623983810..a344271f7 100644 --- a/test/FESpacesTests/AffineFEOperatorsTests.jl +++ b/test/FESpacesTests/AffineFEOperatorsTests.jl @@ -37,7 +37,7 @@ cellmat = integrate(∇(v)*∇(u),trian,quad) cellvec = integrate(v*f,trian,quad) cellids = collect(1:num_cells(trian)) -assem = SparseMatrixAssembler(V,U) +assem = SparseMatrixAssembler(U,V) A = assemble_matrix(assem,[cellmat],[cellids],[cellids]) b = assemble_vector(assem,[cellvec],[cellids]) diff --git a/test/FESpacesTests/FESolversTests.jl b/test/FESpacesTests/FESolversTests.jl index 49526d560..907206802 100644 --- a/test/FESpacesTests/FESolversTests.jl +++ b/test/FESpacesTests/FESolversTests.jl @@ -36,7 +36,7 @@ cellmat = integrate(∇(v)*∇(u),trian,quad) cellvec = integrate(v*f,trian,quad) cellids = collect(1:num_cells(trian)) -assem = SparseMatrixAssembler(V,U) +assem = SparseMatrixAssembler(U,V) A = assemble_matrix(assem,[cellmat],[cellids],[cellids]) b = assemble_vector(assem,[cellvec],[cellids]) x = A \ b diff --git a/test/FESpacesTests/FETermsTests.jl b/test/FESpacesTests/FETermsTests.jl index 2381ccda0..3d3223a27 100644 --- a/test/FESpacesTests/FETermsTests.jl +++ b/test/FESpacesTests/FETermsTests.jl @@ -35,7 +35,7 @@ strian = SkeletonTriangulation(model) degree = 2 squad = CellQuadrature(strian,degree) -assem = SparseMatrixAssembler(V,U) +assem = SparseMatrixAssembler(U,V) uhd = zero(U) v = get_cell_basis(V) diff --git a/test/FESpacesTests/SparseMatrixAssemblersTests.jl b/test/FESpacesTests/SparseMatrixAssemblersTests.jl index e012a15fa..01988f836 100644 --- a/test/FESpacesTests/SparseMatrixAssemblersTests.jl +++ b/test/FESpacesTests/SparseMatrixAssemblersTests.jl @@ -74,7 +74,7 @@ mtypes = [ for T in mtypes - assem = SparseMatrixAssembler(T,V,U) + assem = SparseMatrixAssembler(T,Vector{Float64},U,V) test_assembler(assem,matvecdata,matdata,vecdata) mat = assemble_matrix(assem,[cellmat],[cellids],[cellids]) From cfb46ed790c637dde3edfaa34c02fb2e70ecf182 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 6 May 2020 11:07:39 +0200 Subject: [PATCH 09/17] Added SparseMatrixAssembler interface --- src/FESpaces/Assemblers.jl | 141 +++++++++++++++++++++++++ src/FESpaces/FESpaces.jl | 9 ++ src/FESpaces/SparseMatrixAssemblers.jl | 110 ++++--------------- src/MultiField/MultiField.jl | 8 ++ 4 files changed, 179 insertions(+), 89 deletions(-) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index cb5188862..1e4aa2536 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -100,18 +100,35 @@ function assemble_matrix!(A,a::Assembler,cellmat,cellidsrows,cellidscols) @abstractmethod end +""" +""" +function assemble_matrix_add!( + mat,a::Assembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + @abstractmethod +end + """ """ function assemble_vector!(b,a::Assembler,cellvec,cellids) @abstractmethod end +""" +""" +function assemble_vector_add!(b,a::Assembler,term_to_cellvec,term_to_cellidsrows) + @abstractmethod +end + """ """ function assemble_matrix_and_vector!(A,b,a::Assembler,matvecdata,matdata,vecdata) @abstractmethod end +function assemble_matrix_and_vector_add!(A,b,a::Assembler, matvecdata, matdata, vecdata) + @abstractmethod +end + function assemble_matrix_and_vector!(A,b,a::Assembler,matvecdata) matdata = ([],[],[]) vecdata = ([],[]) @@ -170,16 +187,19 @@ function test_assembler(a::Assembler,matvecdata,matdata,vecdata) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) assemble_matrix!(A,a,matdata...) + assemble_matrix_add!(A,a,matdata...) A = assemble_matrix(a,matdata...) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) b = allocate_vector(a,vecdata...) @test num_free_dofs(test_fesp) == length(b) assemble_vector!(b,a,vecdata...) + assemble_vector_add!(b,a,vecdata...) b = assemble_vector(a,vecdata...) @test num_free_dofs(test_fesp) == length(b) A, b = allocate_matrix_and_vector(a,matvecdata,matdata,vecdata) assemble_matrix_and_vector!(A,b,a,matvecdata,matdata,vecdata) + assemble_matrix_and_vector_add!(A,b,a,matvecdata,matdata,vecdata) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) @test num_free_dofs(test_fesp) == length(b) @@ -191,3 +211,124 @@ function test_assembler(a::Assembler,matvecdata,matdata,vecdata) @test isa(strategy,AssemblyStrategy) end +# This is an extended interface that only make sense for assemblers that build sparse matrices +# (e.g. not for matrix free assemblers) + +""" +""" +abstract type SparseMatrixAssembler <: Assembler end + +""" +""" +function get_matrix_type(a::SparseMatrixAssembler) + @abstractmethod +end + +""" +""" +function get_vector_type(a::SparseMatrixAssembler) + @abstractmethod +end + +function allocate_vector(a::SparseMatrixAssembler,term_to_cellidsrows) + n = num_free_dofs(a.test) + allocate_vector(get_vector_type(a),n) +end + +function assemble_vector!(b,a::SparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) + fill_entries!(b,zero(eltype(b))) + assemble_vector_add!(b,a,term_to_cellvec,term_to_cellidsrows) +end + +""" +""" +function count_matrix_nnz_coo(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) + @abstractmethod +end + +function count_matrix_nnz_coo(a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) + count_matrix_nnz_coo(a,term_to_cellidsrows, term_to_cellidscols) +end + +""" +""" +function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) + @abstractmethod +end + +function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) + fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows, term_to_cellidscols) +end + +function allocate_matrix(a::SparseMatrixAssembler, term_to_cellidsrows, term_to_cellidscols) + n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + I,J,V = allocate_coo_vectors(get_matrix_type(a),n) + fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows,term_to_cellidscols) + m = num_free_dofs(a.test) + n = num_free_dofs(a.trial) + finalize_coo!(get_matrix_type(a),I,J,V,m,n) + sparse_from_coo(get_matrix_type(a),I,J,V,m,n) +end + +function assemble_matrix!( + mat,a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + z = zero(eltype(mat)) + fill_entries!(mat,z) + assemble_matrix_add!(mat,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) +end + +""" +""" +function fill_matrix_coo_numeric!( + I,J,V,a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) + @abstractmethod +end + +function assemble_matrix( + a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + + n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + I,J,V = allocate_coo_vectors(get_matrix_type(a),n) + fill_matrix_coo_numeric!(I,J,V,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + m = num_free_dofs(a.test) + n = num_free_dofs(a.trial) + finalize_coo!(get_matrix_type(a),I,J,V,m,n) + sparse_from_coo(get_matrix_type(a),I,J,V,m,n) +end + +function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, matvecdata, matdata, vecdata) + fill_entries!(A,zero(eltype(A))) + fill_entries!(b,zero(eltype(b))) + assemble_matrix_and_vector_add!(A,b,a,matvecdata, matdata, vecdata) + A, b +end + +""" +""" +function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) + @abstractmethod +end + +function assemble_matrix_and_vector(a::SparseMatrixAssembler, matvecdata, matdata, vecdata) + + term_to_cellidsrows, term_to_cellidscols, = _rearange_cell_ids(matvecdata,matdata,vecdata) + n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + I,J,V = allocate_coo_vectors(get_matrix_type(a),n) + b = allocate_vector(a,vecdata...) + + fill_matrix_and_vector_coo_numeric!(I,J,V,b,a, matvecdata, matdata, vecdata) + + m = num_free_dofs(a.test) + n = num_free_dofs(a.trial) + finalize_coo!(get_matrix_type(a),I,J,V,m,n) + A = sparse_from_coo(get_matrix_type(a),I,J,V,m,n) + + A, b +end + +function test_sparse_matrix_assembler(a::SparseMatrixAssembler,matvecdata,matdata,vecdata) + test_assembler(a,matvecdata,matdata,vecdata) + _ = get_matrix_type(a) + _ = get_vector_type(a) +end + diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index 267e094ea..097257111 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -109,13 +109,22 @@ export assemble_matrix! export assemble_matrix_add! export assemble_matrix export assemble_vector! +export assemble_vector_add! export assemble_vector export assemble_matrix_and_vector! +export assemble_matrix_and_vector_add! export assemble_matrix_and_vector export allocate_vector export allocate_matrix export allocate_matrix_and_vector export test_assembler +export get_matrix_type +export get_vector_type +export count_matrix_nnz_coo +export fill_matrix_coo_symbolic! +export fill_matrix_coo_numeric! +export fill_matrix_and_vector_coo_numeric! +export test_sparse_matrix_assembler export SingleFieldFESpace export num_dirichlet_dofs diff --git a/src/FESpaces/SparseMatrixAssemblers.jl b/src/FESpaces/SparseMatrixAssemblers.jl index e595ef03a..777335c10 100644 --- a/src/FESpaces/SparseMatrixAssemblers.jl +++ b/src/FESpaces/SparseMatrixAssemblers.jl @@ -1,5 +1,5 @@ -struct SparseMatrixAssembler{M,V} <: Assembler +struct SingleFieldSparseMatrixAssembler{M,V} <: SparseMatrixAssembler matrix_type::Type{M} vector_type::Type{V} trial::SingleFieldFESpace @@ -7,14 +7,19 @@ struct SparseMatrixAssembler{M,V} <: Assembler strategy::AssemblyStrategy end +function SparseMatrixAssembler( + mat::Type,vec::Type,trial::SingleFieldFESpace,test::SingleFieldFESpace,strategy::AssemblyStrategy) + SingleFieldSparseMatrixAssembler(mat,vec,trial,test,strategy) +end + function SparseMatrixAssembler(mat::Type,vec::Type,trial::SingleFieldFESpace,test::SingleFieldFESpace) strategy = DefaultAssemblyStrategy() - SparseMatrixAssembler(mat,vec,trial,test,strategy) + SingleFieldSparseMatrixAssembler(mat,vec,trial,test,strategy) end function SparseMatrixAssembler(mat::Type,trial::SingleFieldFESpace,test::SingleFieldFESpace) strategy = DefaultAssemblyStrategy() - SparseMatrixAssembler(mat,Vector{Float64},trial,test,strategy) + SingleFieldSparseMatrixAssembler(mat,Vector{Float64},trial,test,strategy) end """ @@ -23,27 +28,20 @@ function SparseMatrixAssembler(trial::SingleFieldFESpace,test::SingleFieldFESpac matrix_type = SparseMatrixCSC{Float64,Int} vector_type = Vector{Float64} strategy = DefaultAssemblyStrategy() - SparseMatrixAssembler(matrix_type,vector_type,trial,test,strategy) + SingleFieldSparseMatrixAssembler(matrix_type,vector_type,trial,test,strategy) end -get_test(a::SparseMatrixAssembler) = a.test +get_test(a::SingleFieldSparseMatrixAssembler) = a.test -get_trial(a::SparseMatrixAssembler) = a.trial +get_trial(a::SingleFieldSparseMatrixAssembler) = a.trial -get_assembly_strategy(a::SparseMatrixAssembler) = a.strategy +get_matrix_type(a::SingleFieldSparseMatrixAssembler) = a.matrix_type -function allocate_vector(a::SparseMatrixAssembler,term_to_cellidsrows) - n = num_free_dofs(a.test) - allocate_vector(a.vector_type,n) -end +get_vector_type(a::SingleFieldSparseMatrixAssembler) = a.vector_type -function assemble_vector!(b,a::SparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) - fill_entries!(b,zero(eltype(b))) - assemble_vector_add!(b,a,term_to_cellvec,term_to_cellidsrows) -end +get_assembly_strategy(a::SingleFieldSparseMatrixAssembler) = a.strategy -#TODO corresponding abstract method -function assemble_vector_add!(b,a::SparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) +function assemble_vector_add!(b,a::SingleFieldSparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) celldofs = get_cell_dofs(a.test) for (cellvec, cellids) in zip(term_to_cellvec,term_to_cellidsrows) rows = reindex(celldofs,cellids) @@ -69,13 +67,7 @@ function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows,strateg end end -#TODO -function count_matrix_nnz_coo(a::Assembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) - count_matrix_nnz_coo(a,term_to_cellidsrows, term_to_cellidscols) -end - -#TODO -function count_matrix_nnz_coo(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function count_matrix_nnz_coo(a::SingleFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) n = 0 @@ -113,13 +105,7 @@ end n end -#TODO -function fill_matrix_coo_symbolic!(I,J,a::Assembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) - fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows, term_to_cellidscols) -end - -#TODO -function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function fill_matrix_coo_symbolic!(I,J,a::SingleFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) nini = 0 @@ -156,26 +142,8 @@ end n end -# TODO Perhapswe can move this to the abstract interface -function allocate_matrix(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) - n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) - I,J,V = allocate_coo_vectors(a.matrix_type,n) - fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows,term_to_cellidscols) - m = num_free_dofs(a.test) - n = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,m,n) - sparse_from_coo(a.matrix_type,I,J,V,m,n) -end - -function assemble_matrix!( - mat,a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) - z = zero(eltype(mat)) - fill_entries!(mat,z) - assemble_matrix_add!(mat,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) -end - function assemble_matrix_add!( - mat,a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + mat,a::SingleFieldSparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) @@ -214,20 +182,8 @@ function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_r end end -function assemble_matrix( - a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) - - n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) - I,J,V = allocate_coo_vectors(a.matrix_type,n) - fill_matrix_coo_numeric!(I,J,V,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) - m = num_free_dofs(a.test) - n = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,m,n) - sparse_from_coo(a.matrix_type,I,J,V,m,n) -end - function fill_matrix_coo_numeric!( - I,J,V,a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) + I,J,V,a::SingleFieldSparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) nini = n celldofs_rows = get_cell_dofs(a.test) @@ -275,15 +231,8 @@ end n end -function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, matvecdata, matdata, vecdata) - fill_entries!(A,zero(eltype(A))) - fill_entries!(b,zero(eltype(b))) - assemble_matrix_and_vector_add!(A,b,a,matvecdata, matdata, vecdata) - A, b -end - function assemble_matrix_and_vector_add!( - A,b,a::SparseMatrixAssembler, matvecdata, matdata, vecdata) + A,b,a::SingleFieldSparseMatrixAssembler, matvecdata, matdata, vecdata) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) @@ -333,24 +282,7 @@ function _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cell_ end end -function assemble_matrix_and_vector( a::SparseMatrixAssembler, matvecdata, matdata, vecdata) - - term_to_cellidsrows, term_to_cellidscols, = _rearange_cell_ids(matvecdata,matdata,vecdata) - n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) - I,J,V = allocate_coo_vectors(a.matrix_type,n) - b = allocate_vector(a,vecdata...) - - fill_matrix_and_vector_coo_numeric!(I,J,V,b,a, matvecdata, matdata, vecdata) - - m = num_free_dofs(a.test) - n = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,m,n) - A = sparse_from_coo(a.matrix_type,I,J,V,m,n) - - A, b -end - -function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) +function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SingleFieldSparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) nini = n diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 88a4432a5..28ec90d3b 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -61,17 +61,25 @@ import Gridap.FESpaces: get_test import Gridap.FESpaces: get_trial import Gridap.FESpaces: allocate_vector import Gridap.FESpaces: assemble_vector! +import Gridap.FESpaces: assemble_vector_add! import Gridap.FESpaces: allocate_matrix import Gridap.FESpaces: assemble_matrix! import Gridap.FESpaces: assemble_matrix_add! import Gridap.FESpaces: assemble_matrix import Gridap.FESpaces: allocate_matrix_and_vector import Gridap.FESpaces: assemble_matrix_and_vector! +import Gridap.FESpaces: assemble_matrix_and_vector_add! import Gridap.FESpaces: assemble_matrix_and_vector import Gridap.FESpaces: SparseMatrixAssembler import Gridap.FESpaces: AffineFEOperator import Gridap.FESpaces: FEOperator import Gridap.FESpaces: EvaluationFunction +export Gridap.FESpaces: get_matrix_type +export Gridap.FESpaces: get_vector_type +export Gridap.FESpaces: count_matrix_nnz_coo +export Gridap.FESpaces: fill_matrix_coo_symbolic! +export Gridap.FESpaces: fill_matrix_coo_numeric! +export Gridap.FESpaces: fill_matrix_and_vector_coo_numeric! import Base: +, - From 9fcaca44b249af3a394831e3d0cf6853953b08e0 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 6 May 2020 15:29:48 +0200 Subject: [PATCH 10/17] Refactored MultiFieldSparseMatrixAssembler --- src/FESpaces/Assemblers.jl | 12 +- src/FESpaces/FESpaces.jl | 6 + src/MultiField/MultiField.jl | 13 +- .../MultiFieldSparseMatrixAssemblers.jl | 327 ++++++++++-------- 4 files changed, 198 insertions(+), 160 deletions(-) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index 1e4aa2536..149337094 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -289,7 +289,12 @@ function assemble_matrix( n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) I,J,V = allocate_coo_vectors(get_matrix_type(a),n) - fill_matrix_coo_numeric!(I,J,V,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + + nmax = fill_matrix_coo_numeric!(I,J,V,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + resize!(I,nmax) + resize!(J,nmax) + resize!(V,nmax) + m = num_free_dofs(a.test) n = num_free_dofs(a.trial) finalize_coo!(get_matrix_type(a),I,J,V,m,n) @@ -316,7 +321,10 @@ function assemble_matrix_and_vector(a::SparseMatrixAssembler, matvecdata, matdat I,J,V = allocate_coo_vectors(get_matrix_type(a),n) b = allocate_vector(a,vecdata...) - fill_matrix_and_vector_coo_numeric!(I,J,V,b,a, matvecdata, matdata, vecdata) + nmax = fill_matrix_and_vector_coo_numeric!(I,J,V,b,a, matvecdata, matdata, vecdata) + resize!(I,nmax) + resize!(J,nmax) + resize!(V,nmax) m = num_free_dofs(a.test) n = num_free_dofs(a.trial) diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index 097257111..e59f6e467 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -103,6 +103,12 @@ export apply_constraints_matrix_and_vector_rows export test_fe_space export Assembler +export AssemblyStrategy +export row_map +export col_map +export row_mask +export col_mask +export DefaultAssemblyStrategy export get_test export get_trial export assemble_matrix! diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 28ec90d3b..821fbd6ce 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -74,12 +74,13 @@ import Gridap.FESpaces: SparseMatrixAssembler import Gridap.FESpaces: AffineFEOperator import Gridap.FESpaces: FEOperator import Gridap.FESpaces: EvaluationFunction -export Gridap.FESpaces: get_matrix_type -export Gridap.FESpaces: get_vector_type -export Gridap.FESpaces: count_matrix_nnz_coo -export Gridap.FESpaces: fill_matrix_coo_symbolic! -export Gridap.FESpaces: fill_matrix_coo_numeric! -export Gridap.FESpaces: fill_matrix_and_vector_coo_numeric! +import Gridap.FESpaces: get_matrix_type +import Gridap.FESpaces: get_vector_type +import Gridap.FESpaces: get_assembly_strategy +import Gridap.FESpaces: count_matrix_nnz_coo +import Gridap.FESpaces: fill_matrix_coo_symbolic! +import Gridap.FESpaces: fill_matrix_coo_numeric! +import Gridap.FESpaces: fill_matrix_and_vector_coo_numeric! import Base: +, - diff --git a/src/MultiField/MultiFieldSparseMatrixAssemblers.jl b/src/MultiField/MultiFieldSparseMatrixAssemblers.jl index a397efb6a..a8a745a38 100644 --- a/src/MultiField/MultiFieldSparseMatrixAssemblers.jl +++ b/src/MultiField/MultiFieldSparseMatrixAssemblers.jl @@ -1,84 +1,82 @@ """ - struct MultiFieldSparseMatrixAssembler{E} <: Assembler - matrix_type::Type{E} - test::MultiFieldFESpace - trial::MultiFieldFESpace + struct MultiFieldSparseMatrixAssembler{E} <: SparseMatrixAssembler + # private fields end """ -struct MultiFieldSparseMatrixAssembler{E} <: Assembler - matrix_type::Type{E} +struct MultiFieldSparseMatrixAssembler{M,V} <: SparseMatrixAssembler + matrix_type::Type{M} + vector_type::Type{V} test::MultiFieldFESpace trial::MultiFieldFESpace + strategy::AssemblyStrategy end -""" - MultiFieldSparseMatrixAssembler( - matrix_type::Type{<:AbstractSparseMatrix}, - test::Vector{<:SingleFieldFESpace}, - trial::Vector{<:SingleFieldFESpace}) -""" -function MultiFieldSparseMatrixAssembler( +function SparseMatrixAssembler( matrix_type::Type{<:AbstractSparseMatrix}, - test::Vector{<:SingleFieldFESpace}, - trial::Vector{<:SingleFieldFESpace}) - _test = MultiFieldFESpace(test) - _trial = MultiFieldFESpace(trial) - MultiFieldSparseMatrixAssembler(matrix_type,_test,_trial) + vector_type::Type{<:AbstractVector}, + test::MultiFieldFESpace, + trial::MultiFieldFESpace) + strategy = DefaultAssemblyStrategy() + MultiFieldSparseMatrixAssembler(matrix_type,vector_type,test,trial,strategy) end function SparseMatrixAssembler( matrix_type::Type{<:AbstractSparseMatrix}, test::MultiFieldFESpace, trial::MultiFieldFESpace) - MultiFieldSparseMatrixAssembler(matrix_type,test,trial) + SparseMatrixAssembler(matrix_type,Vector{Float64},test,trial) end function SparseMatrixAssembler( test::MultiFieldFESpace, trial::MultiFieldFESpace) matrix_type = SparseMatrixCSC{Float64,Int} - MultiFieldSparseMatrixAssembler(matrix_type,test,trial) + SparseMatrixAssembler(matrix_type,test,trial) end function SparseMatrixAssembler( matrix_type::Type{<:AbstractSparseMatrix}, test::Vector{<:SingleFieldFESpace}, trial::Vector{<:SingleFieldFESpace}) - MultiFieldSparseMatrixAssembler(matrix_type,test,trial) + _test = MultiFieldFESpace(test) + _trial = MultiFieldFESpace(trial) + SparseMatrixAssembler(matrix_type,_test,_trial) end function SparseMatrixAssembler( test::Vector{<:SingleFieldFESpace}, trial::Vector{<:SingleFieldFESpace}) + matrix_type = SparseMatrixCSC{Float64,Int} - SparseMatrixAssembler(matrix_type,test,trial) + _test = MultiFieldFESpace(test) + _trial = MultiFieldFESpace(trial) + SparseMatrixAssembler(matrix_type,_test,_trial) end get_test(a::MultiFieldSparseMatrixAssembler) = a.test get_trial(a::MultiFieldSparseMatrixAssembler) = a.trial -function allocate_vector(a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows) - zero_free_values(a.test) -end +get_matrix_type(a::MultiFieldSparseMatrixAssembler) = a.matrix_type + +get_vector_type(a::MultiFieldSparseMatrixAssembler) = a.vector_type -function assemble_vector!(b,a::MultiFieldSparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) +get_assembly_strategy(a::MultiFieldSparseMatrixAssembler) = a.strategy + +function assemble_vector_add!(b,a::MultiFieldSparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) celldofs = get_cell_dofs(a.test) - fill!(b,zero(eltype(b))) - @assert length(b) == num_free_dofs(a.test) - offsets = compute_field_offsets(a.test) for (cellvec, cellids) in zip(term_to_cellvec,term_to_cellidsrows) rows = reindex(celldofs,cellids) vals = apply_constraints_vector(a.test,cellvec,cellids) rows_cache = array_cache(rows) vals_cache = array_cache(vals) - _assemble_vector!(b,vals_cache,rows_cache,vals,rows) + _assemble_vector!(b,vals_cache,rows_cache,vals,rows,a.strategy) end b end -function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows) +function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows,strategy) @assert length(cell_vals) == length(cell_rows) for cell in 1:length(cell_rows) _rows = getindex!(rows_cache,cell_rows,cell) @@ -89,15 +87,16 @@ function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows) vals = _vals.blocks[block] rows = _rows.blocks[field] for (i,gid) in enumerate(rows) - if gid > 0 - vec[gid] += vals[i] + if gid > 0 && row_mask(strategy,gid) + _gid = row_map(strategy,gid) + vec[_gid] += vals[i] end end end end end -function allocate_matrix(a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function count_matrix_nnz_coo(a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) n = 0 @@ -107,24 +106,13 @@ function allocate_matrix(a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows, rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) @assert length(cell_cols) == length(cell_rows) - n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols) - end - I, J, V = allocate_coo_vectors(a.matrix_type,n) - nini = 0 - for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols) + n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) end - num_rows = num_free_dofs(a.test) - num_cols = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,num_rows,num_cols) - sparse_from_coo(a.matrix_type,I,J,V,num_rows,num_cols) + + n end -@noinline function _count_matrix_entries(::Type{M},rows_cache,cols_cache,cell_rows,cell_cols) where M +@noinline function _count_matrix_entries(::Type{M},rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M n = 0 for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) @@ -136,10 +124,12 @@ end for field_row in 1:nfields_rows rows = _rows.blocks[field_row] for gidcol in cols - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for gidrow in rows - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) n += 1 end end @@ -152,7 +142,20 @@ end n end -@noinline function _allocate_matrix!(::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols) where M +function fill_matrix_coo_symbolic!(I,J,a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) + celldofs_rows = get_cell_dofs(a.test) + celldofs_cols = get_cell_dofs(a.trial) + nini = 0 + for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) + cell_rows = reindex(celldofs_rows,cellidsrows) + cell_cols = reindex(celldofs_cols,cellidscols) + rows_cache = array_cache(cell_rows) + cols_cache = array_cache(cell_cols) + nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) + end +end + +@noinline function _allocate_matrix!(::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M n = nini for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) @@ -164,13 +167,15 @@ end for field_row in 1:nfields_rows @inbounds rows = _rows.blocks[field_row] for gidcol in cols - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for gidrow in rows - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) n += 1 - @inbounds I[n] = gidrow - @inbounds J[n] = gidcol + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol end end end @@ -182,15 +187,9 @@ end n end -function assemble_matrix!( - mat,a::MultiFieldSparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) - z = zero(eltype(mat)) - fill_entries!(mat,z) - assemble_matrix_add!(mat,a,term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) -end - function assemble_matrix_add!( mat,a::MultiFieldSparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) + celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) @@ -201,14 +200,12 @@ function assemble_matrix_add!( rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) vals_cache = array_cache(cellmat) - @assert length(cell_cols) == length(cell_rows) - @assert length(cellmat) == length(cell_rows) - _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) + _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols,a.strategy) end mat end -function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) +function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) _cols = getindex!(cols_cache,cell_cols,cell) @@ -220,11 +217,13 @@ function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_r rows = _rows.blocks[field_row] cols = _cols.blocks[field_col] for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) v = vals[i,j] - add_entry!(mat,v,gidrow,gidcol) + add_entry!(mat,v,_gidrow,_gidcol) end end end @@ -233,45 +232,29 @@ function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_r end end -function assemble_matrix( - a::MultiFieldSparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) +function fill_matrix_coo_numeric!( + I,J,V,a::MultiFieldSparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) + + nini = n celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - n = 0 - for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) - cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - vals_cache = array_cache(cellmat) - @assert length(cell_cols) == length(cell_rows) - @assert length(cellmat) == length(cell_rows) - n += _count_matrix_entries_ijv(a.matrix_type,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) - end - I, J, V = allocate_coo_vectors(a.matrix_type,n) - nini = 0 for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) - cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) + cell_vals = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) - vals_cache = array_cache(cellmat) - @assert length(cell_cols) == length(cell_rows) - @assert length(cellmat) == length(cell_rows) - nini = _assemble_matrix_ijv!(a.matrix_type,nini,I,J,V,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) + vals_cache = array_cache(cell_vals) + nini = _fill_matrix!( + a.matrix_type,nini,I,J,V,rows_cache,cols_cache,vals_cache,cell_rows,cell_cols,cell_vals,a.strategy) end - num_rows = num_free_dofs(a.test) - num_cols = num_free_dofs(a.trial) - finalize_coo!(a.matrix_type,I,J,V,num_rows,num_cols) - sparse_from_coo(a.matrix_type,I,J,V,num_rows,num_cols) + + nini end -@noinline function _count_matrix_entries_ijv(::Type{M},vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) where M - n = 0 +@noinline function _fill_matrix!(::Type{M},nini,I,J,V,rows_cache,cols_cache,vals_cache,cell_rows,cell_cols,cell_vals,strategy) where M + n = nini for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) _cols = getindex!(cols_cache,cell_cols,cell) @@ -279,14 +262,21 @@ end nblocks = length(_vals.blocks) for block in 1:nblocks field_row, field_col = _vals.coordinates[block] + vals = _vals.blocks[block] rows = _rows.blocks[field_row] cols = _cols.blocks[field_col] for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) n += 1 + @inbounds v = vals[i,j] + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol + @inbounds V[n] = v end end end @@ -297,43 +287,78 @@ end n end -@noinline function _assemble_matrix_ijv!(::Type{M},nini,I,J,V,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) where M - n = nini +function assemble_matrix_and_vector_add!( + A,b,a::MultiFieldSparseMatrixAssembler, matvecdata, matdata, vecdata) + + celldofs_rows = get_cell_dofs(a.test) + celldofs_cols = get_cell_dofs(a.trial) + + for (cellmatvec_rc,cellidsrows,cellidscols) in zip(matvecdata...) + cell_rows = reindex(celldofs_rows,cellidsrows) + cell_cols = reindex(celldofs_cols,cellidscols) + cellmatvec_r = apply_constraints_matrix_and_vector_cols(a.trial,cellmatvec_rc,cellidscols) + cellmatvec = apply_constraints_matrix_and_vector_rows(a.test,cellmatvec_r,cellidsrows) + rows_cache = array_cache(cell_rows) + cols_cache = array_cache(cell_cols) + vals_cache = array_cache(cellmatvec) + _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) + end + assemble_matrix_add!(A,a,matdata...) + assemble_vector_add!(b,a,vecdata...) + A, b +end + +@noinline function _assemble_matrix_and_vector!( + mat,vec,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) + for cell in 1:length(cell_cols) + _rows = getindex!(rows_cache,cell_rows,cell) _cols = getindex!(cols_cache,cell_cols,cell) _vals = getindex!(vals_cache,cell_vals,cell) - nblocks = length(_vals.blocks) + _valsmat, _valsvec = _vals + + nblocks = length(_valsmat.blocks) for block in 1:nblocks - field_row, field_col = _vals.coordinates[block] - vals = _vals.blocks[block] + field_row, field_col = _valsmat.coordinates[block] + valsmat = _valsmat.blocks[block] rows = _rows.blocks[field_row] cols = _cols.blocks[field_col] for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 - if is_entry_stored(M,gidrow,gidcol) - n += 1 - v = vals[i,j] - @inbounds I[n] = gidrow - @inbounds J[n] = gidcol - @inbounds V[n] = v - end + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + v = valsmat[i,j] + add_entry!(mat,v,_gidrow,_gidcol) end end end end end + + nblocks = length(_valsvec.blocks) + for block in 1:nblocks + field, = _valsvec.coordinates[block] + valsvec = _valsvec.blocks[block] + rows = _rows.blocks[field] + for (i,gid) in enumerate(rows) + if gid > 0 && row_mask(strategy,gid) + _gid = row_map(strategy,gid) + vec[_gid] += valsvec[i] + end + end + end + end - n end +function fill_matrix_and_vector_coo_numeric!( + I,J,V,b,a::MultiFieldSparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) + + nini = n -function assemble_matrix_and_vector!(A,b, a::MultiFieldSparseMatrixAssembler, matvecdata, matdata, vecdata) - z = zero(eltype(A)) - fill_entries!(A,z) - fill!(b,zero(eltype(b))) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) @@ -345,34 +370,22 @@ function assemble_matrix_and_vector!(A,b, a::MultiFieldSparseMatrixAssembler, ma rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) vals_cache = array_cache(cellmatvec) - _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols) - end - - for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) - cell_rows = reindex(celldofs_rows,cellidsrows) - cell_cols = reindex(celldofs_cols,cellidscols) - cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) - cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) - rows_cache = array_cache(cell_rows) - cols_cache = array_cache(cell_cols) - vals_cache = array_cache(cellmat) @assert length(cell_cols) == length(cell_rows) - @assert length(cellmat) == length(cell_rows) - _assemble_matrix!(A,vals_cache,rows_cache,cols_cache,cellmat,cell_rows,cell_cols) + @assert length(cellmatvec) == length(cell_rows) + nini = _assemble_matrix_and_vector_fill!( + a.matrix_type,nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - for (cellvec, cellids) in zip(vecdata...) - rows = reindex(celldofs_rows,cellids) - vals = apply_constraints_vector(a.test,cellvec,cellids) - rows_cache = array_cache(rows) - vals_cache = array_cache(vals) - _assemble_vector!(b,vals_cache,rows_cache,vals,rows) - end + nini = fill_matrix_coo_numeric!(I,J,V,a,matdata...,nini) + assemble_vector_add!(b,a,vecdata...) - A, b + nini end -function _assemble_matrix_and_vector!(mat,vec,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols) + +@noinline function _assemble_matrix_and_vector_fill!( + ::Type{M},nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) where M + n = nini for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) @@ -383,15 +396,22 @@ function _assemble_matrix_and_vector!(mat,vec,vals_cache,rows_cache,cols_cache,c nblocks = length(_valsmat.blocks) for block in 1:nblocks field_row, field_col = _valsmat.coordinates[block] - vals = _valsmat.blocks[block] + valsmat = _valsmat.blocks[block] rows = _rows.blocks[field_row] cols = _cols.blocks[field_col] for (j,gidcol) in enumerate(cols) - if gidcol > 0 + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) for (i,gidrow) in enumerate(rows) - if gidrow > 0 - v = vals[i,j] - add_entry!(mat,v,gidrow,gidcol) + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) + n += 1 + @inbounds v = valsmat[i,j] + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol + @inbounds V[n] = v + end end end end @@ -401,14 +421,17 @@ function _assemble_matrix_and_vector!(mat,vec,vals_cache,rows_cache,cols_cache,c nblocks = length(_valsvec.blocks) for block in 1:nblocks field, = _valsvec.coordinates[block] - vals = _valsvec.blocks[block] + valsvec = _valsvec.blocks[block] rows = _rows.blocks[field] for (i,gid) in enumerate(rows) - if gid > 0 - vec[gid] += vals[i] + if gid > 0 && row_mask(strategy,gid) + _gid = row_map(strategy,gid) + b[_gid] += valsvec[i] end end end end + n end + From 5f4f63b434aaeabce17c1f894c515f9fbaa57d92 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 7 May 2020 12:07:02 +0200 Subject: [PATCH 11/17] Simplified Assembler interface --- src/FESpaces/AffineFEOperators.jl | 8 +- src/FESpaces/Assemblers.jl | 209 ++++++++---------- src/FESpaces/FEOperatorsFromTerms.jl | 20 +- src/FESpaces/SparseMatrixAssemblers.jl | 54 +++-- src/FESpaces/ZeroMeanFESpaces.jl | 3 +- src/MultiField/MultiField.jl | 2 + .../MultiFieldSparseMatrixAssemblers.jl | 155 ++++++++----- test/FESpacesTests/AffineFEOperatorsTests.jl | 6 +- test/FESpacesTests/FESolversTests.jl | 6 +- test/FESpacesTests/FETermsTests.jl | 46 ++-- .../SparseMatrixAssemblersTests.jl | 30 ++- .../MultiFieldSparseMatrixAssemblersTests.jl | 13 +- 12 files changed, 296 insertions(+), 256 deletions(-) diff --git a/src/FESpaces/AffineFEOperators.jl b/src/FESpaces/AffineFEOperators.jl index e6ca9f6e9..1a2163484 100644 --- a/src/FESpaces/AffineFEOperators.jl +++ b/src/FESpaces/AffineFEOperators.jl @@ -28,13 +28,13 @@ function AffineFEOperator(trial::FESpace,test::FESpace,assem::Assembler,terms::A uhd = zero(trial) - matvecdata, matdata, vecdata = collect_cell_matrix_and_vector(uhd,u,v,terms) - A,b = assemble_matrix_and_vector(assem,matvecdata,matdata,vecdata) + data = collect_cell_matrix_and_vector(uhd,u,v,terms) + A,b = assemble_matrix_and_vector(assem,data) #matdata = collect_cell_matrix(u,v,terms) #vecdata = collect_cell_vector(uhd,v,terms) - #A = assemble_matrix(assem,matdata...) - #b = assemble_vector(assem,vecdata...) + #A = assemble_matrix(assem,matdata) + #b = assemble_vector(assem,vecdata) AffineFEOperator(trial,test,A,b) end diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index 149337094..854917d4f 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -61,149 +61,106 @@ end """ """ -function allocate_matrix(a::Assembler,cellidsrows,cellidscols) +function allocate_matrix(a::Assembler,matdata) @abstractmethod end -function allocate_matrix(a::Assembler,cellmat,cellidsrows,cellidscols) - allocate_matrix(a,cellidsrows,cellidscols) -end - """ """ -function allocate_vector(a::Assembler,cellidsrows) +function allocate_vector(a::Assembler,vecdata) @abstractmethod end -function allocate_vector(a::Assembler,cellvec,cellidsrows) - allocate_vector(a,cellidsrows) -end - """ """ -function allocate_matrix_and_vector(a::Assembler,matvecdata,matdata,vecdata) - matrows, matcols, vecrows = _rearange_cell_ids(matvecdata,matdata,vecdata) - A = allocate_matrix(a,matrows,matcols) - b = allocate_vector(a,vecrows) - (A,b) -end - -function allocate_matrix_and_vector(a::Assembler,matvecdata) - matdata = ([],[],[]) - vecdata = ([],[]) - allocate_matrix_and_vector(a,matvecdata,matdata,vecdata) +function allocate_matrix_and_vector(a::Assembler,data) + @abstractmethod end """ """ -function assemble_matrix!(A,a::Assembler,cellmat,cellidsrows,cellidscols) +function assemble_matrix!(A,a::Assembler,matdata) @abstractmethod end """ """ -function assemble_matrix_add!( - mat,a::Assembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) +function assemble_matrix_add!(mat,a::Assembler, matdata) @abstractmethod end """ """ -function assemble_vector!(b,a::Assembler,cellvec,cellids) +function assemble_vector!(b,a::Assembler,vecdata) @abstractmethod end """ """ -function assemble_vector_add!(b,a::Assembler,term_to_cellvec,term_to_cellidsrows) +function assemble_vector_add!(b,a::Assembler,vecdata) @abstractmethod end """ """ -function assemble_matrix_and_vector!(A,b,a::Assembler,matvecdata,matdata,vecdata) +function assemble_matrix_and_vector!(A,b,a::Assembler, data) @abstractmethod end -function assemble_matrix_and_vector_add!(A,b,a::Assembler, matvecdata, matdata, vecdata) +function assemble_matrix_and_vector_add!(A,b,a::Assembler, data) @abstractmethod end -function assemble_matrix_and_vector!(A,b,a::Assembler,matvecdata) - matdata = ([],[],[]) - vecdata = ([],[]) - assemble_matrix_and_vector!(A,b,a,matvecdata,matdata,vecdata) -end - """ """ -function assemble_matrix(a::Assembler,cellmat,cellidsrows,cellidscols) - A = allocate_matrix(a,cellidsrows,cellidscols) - assemble_matrix!(A,a,cellmat,cellidsrows,cellidscols) +function assemble_matrix(a::Assembler,matdata) + A = allocate_matrix(a,matdata) + assemble_matrix!(A,a,matdata) A end """ """ -function assemble_vector(a::Assembler,cellvec,cellids) - b = allocate_vector(a,cellids) - assemble_vector!(b,a,cellvec,cellids) +function assemble_vector(a::Assembler,vecdata) + b = allocate_vector(a,vecdata) + assemble_vector!(b,a,vecdata) b end """ """ -function assemble_matrix_and_vector(a::Assembler,matvecdata,matdata,vecdata) - A, b = allocate_matrix_and_vector(a,matvecdata,matdata,vecdata) - assemble_matrix_and_vector!(A,b,a,matvecdata,matdata,vecdata) +function assemble_matrix_and_vector(a::Assembler,data) + A, b = allocate_matrix_and_vector(a,data) + assemble_matrix_and_vector!(A,b,a,data) (A, b) end -function assemble_matrix_and_vector(a::Assembler,matvecdata) - matdata = ([],[],[]) - vecdata = ([],[]) - assemble_matrix_and_vector(a,matvecdata,matdata,vecdata) -end - -function _rearange_cell_ids(matvecdata,matdata,vecdata) - - matvec1, rows1, cols1 = matvecdata - mat2, rows2, cols2 = matdata - vec3, rows3 = vecdata - - matrows = vcat(rows1,rows2) - matcols = vcat(cols1,cols2) - vecrows = vcat(rows1,rows3) - - (matrows, matcols, vecrows) -end - """ """ -function test_assembler(a::Assembler,matvecdata,matdata,vecdata) +function test_assembler(a::Assembler,matdata,vecdata,data) trial_fesp = get_trial(a) test_fesp = get_test(a) - A = allocate_matrix(a,matdata...) + A = allocate_matrix(a,matdata) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) - assemble_matrix!(A,a,matdata...) - assemble_matrix_add!(A,a,matdata...) - A = assemble_matrix(a,matdata...) + assemble_matrix!(A,a,matdata) + assemble_matrix_add!(A,a,matdata) + A = assemble_matrix(a,matdata) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) - b = allocate_vector(a,vecdata...) + b = allocate_vector(a,vecdata) @test num_free_dofs(test_fesp) == length(b) - assemble_vector!(b,a,vecdata...) - assemble_vector_add!(b,a,vecdata...) - b = assemble_vector(a,vecdata...) + assemble_vector!(b,a,vecdata) + assemble_vector_add!(b,a,vecdata) + b = assemble_vector(a,vecdata) @test num_free_dofs(test_fesp) == length(b) - A, b = allocate_matrix_and_vector(a,matvecdata,matdata,vecdata) - assemble_matrix_and_vector!(A,b,a,matvecdata,matdata,vecdata) - assemble_matrix_and_vector_add!(A,b,a,matvecdata,matdata,vecdata) + A, b = allocate_matrix_and_vector(a,data) + assemble_matrix_and_vector!(A,b,a,data) + assemble_matrix_and_vector_add!(A,b,a,data) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) @test num_free_dofs(test_fesp) == length(b) - A, b = assemble_matrix_and_vector(a,matvecdata,matdata,vecdata) + A, b = assemble_matrix_and_vector(a,data) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) @test num_free_dofs(test_fesp) == length(b) @@ -211,8 +168,8 @@ function test_assembler(a::Assembler,matvecdata,matdata,vecdata) @test isa(strategy,AssemblyStrategy) end -# This is an extended interface that only make sense for assemblers that build sparse matrices -# (e.g. not for matrix free assemblers) +# This is an extended interface that only make sense for assemblers that build (sequential) sparse matrices +# (e.g. not for matrix free assemblers or for distributed assemblers) """ """ @@ -230,112 +187,122 @@ function get_vector_type(a::SparseMatrixAssembler) @abstractmethod end -function allocate_vector(a::SparseMatrixAssembler,term_to_cellidsrows) - n = num_free_dofs(a.test) +function allocate_vector(a::SparseMatrixAssembler,vecdata) + n = num_free_dofs(get_test(a)) allocate_vector(get_vector_type(a),n) end -function assemble_vector!(b,a::SparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) +function assemble_vector!(b,a::SparseMatrixAssembler,vecdata) fill_entries!(b,zero(eltype(b))) - assemble_vector_add!(b,a,term_to_cellvec,term_to_cellidsrows) + assemble_vector_add!(b,a,vecdata) end """ """ -function count_matrix_nnz_coo(a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function count_matrix_nnz_coo(a::SparseMatrixAssembler,vecdata) @abstractmethod end -function count_matrix_nnz_coo(a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) - count_matrix_nnz_coo(a,term_to_cellidsrows, term_to_cellidscols) +""" +""" +function count_matrix_and_vector_nnz_coo(a::SparseMatrixAssembler,data) + @abstractmethod end """ """ -function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,matdata,n=0) @abstractmethod end -function fill_matrix_coo_symbolic!(I,J,a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols) - fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows, term_to_cellidscols) +function fill_matrix_and_vector_coo_symbolic!(I,J,a::SparseMatrixAssembler,data,n=0) + @abstractmethod end -function allocate_matrix(a::SparseMatrixAssembler, term_to_cellidsrows, term_to_cellidscols) - n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) +function allocate_matrix(a::SparseMatrixAssembler,matdata) + n = count_matrix_nnz_coo(a,matdata) I,J,V = allocate_coo_vectors(get_matrix_type(a),n) - fill_matrix_coo_symbolic!(I,J,a,term_to_cellidsrows,term_to_cellidscols) - m = num_free_dofs(a.test) - n = num_free_dofs(a.trial) + fill_matrix_coo_symbolic!(I,J,a,matdata) + m = num_free_dofs(get_test(a)) + n = num_free_dofs(get_trial(a)) finalize_coo!(get_matrix_type(a),I,J,V,m,n) sparse_from_coo(get_matrix_type(a),I,J,V,m,n) end -function assemble_matrix!( - mat,a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) +function assemble_matrix!(mat,a::SparseMatrixAssembler,matdata) z = zero(eltype(mat)) fill_entries!(mat,z) - assemble_matrix_add!(mat,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + assemble_matrix_add!(mat,a,matdata) end """ """ -function fill_matrix_coo_numeric!( - I,J,V,a::SparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) +function fill_matrix_coo_numeric!(I,J,V,a::SparseMatrixAssembler,matdata,n=0) @abstractmethod end -function assemble_matrix( - a::SparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) +function assemble_matrix(a::SparseMatrixAssembler,matdata) - n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + n = count_matrix_nnz_coo(a,matdata) I,J,V = allocate_coo_vectors(get_matrix_type(a),n) - nmax = fill_matrix_coo_numeric!(I,J,V,a,term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) - resize!(I,nmax) - resize!(J,nmax) - resize!(V,nmax) + fill_matrix_coo_numeric!(I,J,V,a,matdata) - m = num_free_dofs(a.test) - n = num_free_dofs(a.trial) + m = num_free_dofs(get_test(a)) + n = num_free_dofs(get_trial(a)) finalize_coo!(get_matrix_type(a),I,J,V,m,n) sparse_from_coo(get_matrix_type(a),I,J,V,m,n) end -function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, matvecdata, matdata, vecdata) + +function allocate_matrix_and_vector(a::SparseMatrixAssembler,data) + + n = count_matrix_and_vector_nnz_coo(a,data) + + I,J,V = allocate_coo_vectors(get_matrix_type(a),n) + fill_matrix_and_vector_coo_symbolic!(I,J,a,data) + m = num_free_dofs(get_test(a)) + n = num_free_dofs(get_trial(a)) + finalize_coo!(get_matrix_type(a),I,J,V,m,n) + A = sparse_from_coo(get_matrix_type(a),I,J,V,m,n) + + b = allocate_vector(get_vector_type(a),m) + + A,b +end + +function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, data) fill_entries!(A,zero(eltype(A))) fill_entries!(b,zero(eltype(b))) - assemble_matrix_and_vector_add!(A,b,a,matvecdata, matdata, vecdata) + assemble_matrix_and_vector_add!(A,b,a,data) A, b end """ """ -function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) +function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SparseMatrixAssembler,data,n=0) @abstractmethod end -function assemble_matrix_and_vector(a::SparseMatrixAssembler, matvecdata, matdata, vecdata) +function assemble_matrix_and_vector(a::SparseMatrixAssembler, data) - term_to_cellidsrows, term_to_cellidscols, = _rearange_cell_ids(matvecdata,matdata,vecdata) - n = count_matrix_nnz_coo(a,term_to_cellidsrows,term_to_cellidscols) + n = count_matrix_and_vector_nnz_coo(a,data) I,J,V = allocate_coo_vectors(get_matrix_type(a),n) - b = allocate_vector(a,vecdata...) + n = num_free_dofs(get_test(a)) + b = allocate_vector(get_vector_type(a),n) - nmax = fill_matrix_and_vector_coo_numeric!(I,J,V,b,a, matvecdata, matdata, vecdata) - resize!(I,nmax) - resize!(J,nmax) - resize!(V,nmax) + fill_matrix_and_vector_coo_numeric!(I,J,V,b,a,data) - m = num_free_dofs(a.test) - n = num_free_dofs(a.trial) + m = num_free_dofs(get_test(a)) + n = num_free_dofs(get_trial(a)) finalize_coo!(get_matrix_type(a),I,J,V,m,n) A = sparse_from_coo(get_matrix_type(a),I,J,V,m,n) A, b end -function test_sparse_matrix_assembler(a::SparseMatrixAssembler,matvecdata,matdata,vecdata) - test_assembler(a,matvecdata,matdata,vecdata) +function test_sparse_matrix_assembler(a::SparseMatrixAssembler,matdata,vecdata,data) + test_assembler(a,matdata,vecdata,data) _ = get_matrix_type(a) _ = get_vector_type(a) end diff --git a/src/FESpaces/FEOperatorsFromTerms.jl b/src/FESpaces/FEOperatorsFromTerms.jl index 436909d92..b69e1e25d 100644 --- a/src/FESpaces/FEOperatorsFromTerms.jl +++ b/src/FESpaces/FEOperatorsFromTerms.jl @@ -43,15 +43,15 @@ end function allocate_residual(op::FEOperatorFromTerms,uh) @assert is_a_fe_function(uh) v = get_cell_basis(op.test) - _, cellids = collect_cell_residual(uh,v,op.terms) - allocate_vector(op.assem, cellids) + vecdata = collect_cell_residual(uh,v,op.terms) + allocate_vector(op.assem, vecdata) end function residual!(b::AbstractVector,op::FEOperatorFromTerms,uh) @assert is_a_fe_function(uh) v = get_cell_basis(op.test) - cellvecs, cellids = collect_cell_residual(uh,v,op.terms) - assemble_vector!(b,op.assem, cellvecs, cellids) + vecdata = collect_cell_residual(uh,v,op.terms) + assemble_vector!(b,op.assem, vecdata) b end @@ -59,16 +59,16 @@ function allocate_jacobian(op::FEOperatorFromTerms,uh) @assert is_a_fe_function(uh) du = get_cell_basis(op.trial) v = get_cell_basis(op.test) - _, cellidsrows, cellidscols = collect_cell_jacobian(uh,du,v,op.terms) - allocate_matrix(op.assem, cellidsrows, cellidscols) + matdata = collect_cell_jacobian(uh,du,v,op.terms) + allocate_matrix(op.assem, matdata) end function jacobian!(A::AbstractMatrix,op::FEOperatorFromTerms,uh) @assert is_a_fe_function(uh) du = get_cell_basis(op.trial) v = get_cell_basis(op.test) - cellmats, cellidsrows, cellidscols = collect_cell_jacobian(uh,du,v,op.terms) - assemble_matrix!(A,op.assem, cellmats, cellidsrows, cellidscols) + matdata = collect_cell_jacobian(uh,du,v,op.terms) + assemble_matrix!(A,op.assem,matdata) A end @@ -77,7 +77,7 @@ function residual_and_jacobian!(b::AbstractVector,A::AbstractMatrix,op::FEOperat du = get_cell_basis(op.trial) v = get_cell_basis(op.test) data = collect_cell_jacobian_and_residual(uh,du,v,op.terms) - assemble_matrix_and_vector!(A, b, op.assem,data...) + assemble_matrix_and_vector!(A, b, op.assem, data) (b,A) end @@ -86,6 +86,6 @@ function residual_and_jacobian(op::FEOperatorFromTerms,uh) du = get_cell_basis(op.trial) v = get_cell_basis(op.test) data = collect_cell_jacobian_and_residual(uh,du,v,op.terms) - A, b = assemble_matrix_and_vector(op.assem,data...) + A, b = assemble_matrix_and_vector(op.assem, data) (b, A) end diff --git a/src/FESpaces/SparseMatrixAssemblers.jl b/src/FESpaces/SparseMatrixAssemblers.jl index 777335c10..283aa50da 100644 --- a/src/FESpaces/SparseMatrixAssemblers.jl +++ b/src/FESpaces/SparseMatrixAssemblers.jl @@ -41,9 +41,9 @@ get_vector_type(a::SingleFieldSparseMatrixAssembler) = a.vector_type get_assembly_strategy(a::SingleFieldSparseMatrixAssembler) = a.strategy -function assemble_vector_add!(b,a::SingleFieldSparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) +function assemble_vector_add!(b,a::SingleFieldSparseMatrixAssembler,vecdata) celldofs = get_cell_dofs(a.test) - for (cellvec, cellids) in zip(term_to_cellvec,term_to_cellidsrows) + for (cellvec, cellids) in zip(vecdata...) rows = reindex(celldofs,cellids) vals = apply_constraints_vector(a.test,cellvec,cellids) rows_cache = array_cache(rows) @@ -67,7 +67,8 @@ function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows,strateg end end -function count_matrix_nnz_coo(a::SingleFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function count_matrix_nnz_coo(a::SingleFieldSparseMatrixAssembler,matdata) + _,term_to_cellidsrows, term_to_cellidscols = matdata celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) n = 0 @@ -105,10 +106,18 @@ end n end -function fill_matrix_coo_symbolic!(I,J,a::SingleFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function count_matrix_and_vector_nnz_coo(a::SingleFieldSparseMatrixAssembler,data) + matvecdata, matdata, vecdata = data + n = count_matrix_nnz_coo(a,matvecdata) + n += count_matrix_nnz_coo(a,matdata) + n +end + +function fill_matrix_coo_symbolic!(I,J,a::SingleFieldSparseMatrixAssembler,matdata,n=0) + _,term_to_cellidsrows, term_to_cellidscols = matdata celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - nini = 0 + nini = n for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) @@ -116,6 +125,7 @@ function fill_matrix_coo_symbolic!(I,J,a::SingleFieldSparseMatrixAssembler,term_ cols_cache = array_cache(cell_cols) nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) end + nini end @noinline function _allocate_matrix!(a::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M @@ -142,12 +152,18 @@ end n end -function assemble_matrix_add!( - mat,a::SingleFieldSparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) +function fill_matrix_and_vector_coo_symbolic!(I,J,a::SingleFieldSparseMatrixAssembler,data,n=0) + matvecdata, matdata, vecdata = data + nini = fill_matrix_coo_symbolic!(I,J,a,matvecdata,n) + nini = fill_matrix_coo_symbolic!(I,J,a,matdata,nini) + nini +end + +function assemble_matrix_add!(mat,a::SingleFieldSparseMatrixAssembler,matdata) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) @@ -182,13 +198,12 @@ function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_r end end -function fill_matrix_coo_numeric!( - I,J,V,a::SingleFieldSparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) +function fill_matrix_coo_numeric!(I,J,V,a::SingleFieldSparseMatrixAssembler,matdata,n=0) nini = n celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) @@ -231,9 +246,9 @@ end n end -function assemble_matrix_and_vector_add!( - A,b,a::SingleFieldSparseMatrixAssembler, matvecdata, matdata, vecdata) +function assemble_matrix_and_vector_add!(A,b,a::SingleFieldSparseMatrixAssembler, data) + matvecdata, matdata, vecdata = data celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) @@ -247,8 +262,8 @@ function assemble_matrix_and_vector_add!( vals_cache = array_cache(cellmatvec) _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - assemble_matrix_add!(A,a,matdata...) - assemble_vector_add!(b,a,vecdata...) + assemble_matrix_add!(A,a,matdata) + assemble_vector_add!(b,a,vecdata) A, b end @@ -282,8 +297,9 @@ function _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cell_ end end -function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SingleFieldSparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) - +function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SingleFieldSparseMatrixAssembler,data,n=0) + + matvecdata, matdata, vecdata = data nini = n celldofs_rows = get_cell_dofs(a.test) @@ -303,8 +319,8 @@ function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::SingleFieldSparseMatrixA a.matrix_type,nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - nini = fill_matrix_coo_numeric!(I,J,V,a,matdata...,nini) - assemble_vector_add!(b,a,vecdata...) + nini = fill_matrix_coo_numeric!(I,J,V,a,matdata,nini) + assemble_vector_add!(b,a,vecdata) nini end diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl index fc12ce1d2..dbd02abd8 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -32,7 +32,8 @@ function _setup_vols(V,trian,quad) bh_trian = restrict(bh,trian) cellvec = integrate(bh_trian,trian,quad) cellids = get_cell_id(trian) - vol_i = assemble_vector(assem,[cellvec],[cellids]) + vecdata = ([cellvec],[cellids]) + vol_i = assemble_vector(assem,vecdata) vol = sum(vol_i) (vol_i, vol) end diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 821fbd6ce..27577d11a 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -81,6 +81,8 @@ import Gridap.FESpaces: count_matrix_nnz_coo import Gridap.FESpaces: fill_matrix_coo_symbolic! import Gridap.FESpaces: fill_matrix_coo_numeric! import Gridap.FESpaces: fill_matrix_and_vector_coo_numeric! +import Gridap.FESpaces: count_matrix_and_vector_nnz_coo +import Gridap.FESpaces: fill_matrix_and_vector_coo_symbolic! import Base: +, - diff --git a/src/MultiField/MultiFieldSparseMatrixAssemblers.jl b/src/MultiField/MultiFieldSparseMatrixAssemblers.jl index a8a745a38..c26b24210 100644 --- a/src/MultiField/MultiFieldSparseMatrixAssemblers.jl +++ b/src/MultiField/MultiFieldSparseMatrixAssemblers.jl @@ -64,9 +64,9 @@ get_vector_type(a::MultiFieldSparseMatrixAssembler) = a.vector_type get_assembly_strategy(a::MultiFieldSparseMatrixAssembler) = a.strategy -function assemble_vector_add!(b,a::MultiFieldSparseMatrixAssembler,term_to_cellvec,term_to_cellidsrows) +function assemble_vector_add!(b,a::MultiFieldSparseMatrixAssembler,vecdata) celldofs = get_cell_dofs(a.test) - for (cellvec, cellids) in zip(term_to_cellvec,term_to_cellidsrows) + for (cellvec, cellids) in zip(vecdata...) rows = reindex(celldofs,cellids) vals = apply_constraints_vector(a.test,cellvec,cellids) rows_cache = array_cache(rows) @@ -96,42 +96,63 @@ function _assemble_vector!(vec,vals_cache,rows_cache,cell_vals,cell_rows,strateg end end -function count_matrix_nnz_coo(a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function count_matrix_nnz_coo(a::MultiFieldSparseMatrixAssembler,matdata) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) n = 0 - for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) + for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) + cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) + cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) @assert length(cell_cols) == length(cell_rows) - n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) + if length(cellmat) > 0 + coords = first(cellmat).coordinates + n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy,coords) + end end + n +end +function count_matrix_and_vector_nnz_coo(a::MultiFieldSparseMatrixAssembler,data) + matvecdata, matdata, vecdata = data + n = count_matrix_nnz_coo(a,matdata) + celldofs_rows = get_cell_dofs(a.test) + celldofs_cols = get_cell_dofs(a.trial) + for (cellmatvec_rc,cellidsrows,cellidscols) in zip(matvecdata...) + cell_rows = reindex(celldofs_rows,cellidsrows) + cell_cols = reindex(celldofs_cols,cellidscols) + cellmatvec_r = apply_constraints_matrix_and_vector_cols(a.trial,cellmatvec_rc,cellidscols) + cellmatvec = apply_constraints_matrix_and_vector_rows(a.test,cellmatvec_r,cellidsrows) + rows_cache = array_cache(cell_rows) + cols_cache = array_cache(cell_cols) + @assert length(cell_cols) == length(cell_rows) + if length(cellmatvec) > 0 + coords = first(cellmatvec)[1].coordinates + n += _count_matrix_entries(a.matrix_type,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy,coords) + end + end n end -@noinline function _count_matrix_entries(::Type{M},rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M +@noinline function _count_matrix_entries(::Type{M},rows_cache,cols_cache,cell_rows,cell_cols,strategy,coords) where M n = 0 for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) _cols = getindex!(cols_cache,cell_cols,cell) - nfields_rows = length(_rows.blocks) - nfields_cols = length(_cols.blocks) - for field_col in 1:nfields_cols + for (field_row,field_col) in coords cols = _cols.blocks[field_col] - for field_row in 1:nfields_rows - rows = _rows.blocks[field_row] - for gidcol in cols - if gidcol > 0 && col_mask(strategy,gidcol) - _gidcol = col_map(strategy,gidcol) - for gidrow in rows - if gidrow > 0 && row_mask(strategy,gidrow) - _gidrow = row_map(strategy,gidrow) - if is_entry_stored(M,_gidrow,_gidcol) - n += 1 - end + rows = _rows.blocks[field_row] + for gidcol in cols + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) + for gidrow in rows + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) + n += 1 end end end @@ -142,41 +163,66 @@ end n end -function fill_matrix_coo_symbolic!(I,J,a::MultiFieldSparseMatrixAssembler,term_to_cellidsrows, term_to_cellidscols) +function fill_matrix_coo_symbolic!(I,J,a::MultiFieldSparseMatrixAssembler,matdata,n=0) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - nini = 0 - for (cellidsrows,cellidscols) in zip(term_to_cellidsrows,term_to_cellidscols) + nini = n + for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) + cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) + cellmat = apply_constraints_matrix_rows(a.test,cellmat_r,cellidsrows) rows_cache = array_cache(cell_rows) cols_cache = array_cache(cell_cols) - nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy) + if length(cellmat) > 0 + coords = first(cellmat).coordinates + nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy,coords) + end end + nini end -@noinline function _allocate_matrix!(::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,strategy) where M +function fill_matrix_and_vector_coo_symbolic!(I,J,a::MultiFieldSparseMatrixAssembler,data,n=0) + matvecdata, matdata, vecdata = data + + nini = fill_matrix_coo_symbolic!(I,J,a,matdata,n) + + celldofs_rows = get_cell_dofs(a.test) + celldofs_cols = get_cell_dofs(a.trial) + for (cellmat_rc,cellidsrows,cellidscols) in zip(matvecdata...) + cell_rows = reindex(celldofs_rows,cellidsrows) + cell_cols = reindex(celldofs_cols,cellidscols) + cellmat_r = apply_constraints_matrix_and_vector_cols(a.trial,cellmat_rc,cellidscols) + cellmat = apply_constraints_matrix_and_vector_rows(a.test,cellmat_r,cellidsrows) + rows_cache = array_cache(cell_rows) + cols_cache = array_cache(cell_cols) + if length(cellmat) > 0 + coords = first(cellmat)[1].coordinates + nini = _allocate_matrix!(a.matrix_type,nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,a.strategy,coords) + end + end + nini +end + +@noinline function _allocate_matrix!( + ::Type{M},nini,I,J,rows_cache,cols_cache,cell_rows,cell_cols,strategy,coords) where M n = nini for cell in 1:length(cell_cols) _rows = getindex!(rows_cache,cell_rows,cell) _cols = getindex!(cols_cache,cell_cols,cell) - nfields_rows = length(_rows.blocks) - nfields_cols = length(_cols.blocks) - for field_col in 1:nfields_cols + for (field_row,field_col) in coords @inbounds cols = _cols.blocks[field_col] - for field_row in 1:nfields_rows - @inbounds rows = _rows.blocks[field_row] - for gidcol in cols - if gidcol > 0 && col_mask(strategy,gidcol) - _gidcol = col_map(strategy,gidcol) - for gidrow in rows - if gidrow > 0 && row_mask(strategy,gidrow) - _gidrow = row_map(strategy,gidrow) - if is_entry_stored(M,_gidrow,_gidcol) - n += 1 - @inbounds I[n] = _gidrow - @inbounds J[n] = _gidcol - end + @inbounds rows = _rows.blocks[field_row] + for gidcol in cols + if gidcol > 0 && col_mask(strategy,gidcol) + _gidcol = col_map(strategy,gidcol) + for gidrow in rows + if gidrow > 0 && row_mask(strategy,gidrow) + _gidrow = row_map(strategy,gidrow) + if is_entry_stored(M,_gidrow,_gidcol) + n += 1 + @inbounds I[n] = _gidrow + @inbounds J[n] = _gidcol end end end @@ -187,12 +233,11 @@ end n end -function assemble_matrix_add!( - mat,a::MultiFieldSparseMatrixAssembler, term_to_cellmat, term_to_cellidsrows, term_to_cellidscols) +function assemble_matrix_add!(mat,a::MultiFieldSparseMatrixAssembler, matdata) celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) @@ -232,13 +277,12 @@ function _assemble_matrix!(mat,vals_cache,rows_cache,cols_cache,cell_vals,cell_r end end -function fill_matrix_coo_numeric!( - I,J,V,a::MultiFieldSparseMatrixAssembler,term_to_cellmat,term_to_cellidsrows, term_to_cellidscols,n=0) +function fill_matrix_coo_numeric!(I,J,V,a::MultiFieldSparseMatrixAssembler,matdata,n=0) nini = n celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) - for (cellmat_rc,cellidsrows,cellidscols) in zip(term_to_cellmat,term_to_cellidsrows,term_to_cellidscols) + for (cellmat_rc,cellidsrows,cellidscols) in zip(matdata...) cell_rows = reindex(celldofs_rows,cellidsrows) cell_cols = reindex(celldofs_cols,cellidscols) cellmat_r = apply_constraints_matrix_cols(a.trial,cellmat_rc,cellidscols) @@ -287,9 +331,9 @@ end n end -function assemble_matrix_and_vector_add!( - A,b,a::MultiFieldSparseMatrixAssembler, matvecdata, matdata, vecdata) +function assemble_matrix_and_vector_add!(A,b,a::MultiFieldSparseMatrixAssembler, data) + matvecdata, matdata, vecdata = data celldofs_rows = get_cell_dofs(a.test) celldofs_cols = get_cell_dofs(a.trial) @@ -303,8 +347,8 @@ function assemble_matrix_and_vector_add!( vals_cache = array_cache(cellmatvec) _assemble_matrix_and_vector!(A,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - assemble_matrix_add!(A,a,matdata...) - assemble_vector_add!(b,a,vecdata...) + assemble_matrix_add!(A,a,matdata) + assemble_vector_add!(b,a,vecdata) A, b end @@ -354,9 +398,9 @@ end end end -function fill_matrix_and_vector_coo_numeric!( - I,J,V,b,a::MultiFieldSparseMatrixAssembler,matvecdata, matdata, vecdata,n=0) +function fill_matrix_and_vector_coo_numeric!(I,J,V,b,a::MultiFieldSparseMatrixAssembler,data,n=0) + matvecdata, matdata, vecdata = data nini = n celldofs_rows = get_cell_dofs(a.test) @@ -376,13 +420,12 @@ function fill_matrix_and_vector_coo_numeric!( a.matrix_type,nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cellmatvec,cell_rows,cell_cols,a.strategy) end - nini = fill_matrix_coo_numeric!(I,J,V,a,matdata...,nini) - assemble_vector_add!(b,a,vecdata...) + nini = fill_matrix_coo_numeric!(I,J,V,a,matdata,nini) + assemble_vector_add!(b,a,vecdata) nini end - @noinline function _assemble_matrix_and_vector_fill!( ::Type{M},nini,I,J,V,b,vals_cache,rows_cache,cols_cache,cell_vals,cell_rows,cell_cols,strategy) where M n = nini diff --git a/test/FESpacesTests/AffineFEOperatorsTests.jl b/test/FESpacesTests/AffineFEOperatorsTests.jl index a344271f7..2eca0c554 100644 --- a/test/FESpacesTests/AffineFEOperatorsTests.jl +++ b/test/FESpacesTests/AffineFEOperatorsTests.jl @@ -38,8 +38,10 @@ cellvec = integrate(v*f,trian,quad) cellids = collect(1:num_cells(trian)) assem = SparseMatrixAssembler(U,V) -A = assemble_matrix(assem,[cellmat],[cellids],[cellids]) -b = assemble_vector(assem,[cellvec],[cellids]) +matdata = ([cellmat],[cellids],[cellids]) +vecdata = ([cellvec],[cellids]) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) op = AffineFEOperator(U,V,A,b) @test A === get_matrix(op) diff --git a/test/FESpacesTests/FESolversTests.jl b/test/FESpacesTests/FESolversTests.jl index 907206802..47c70a3f1 100644 --- a/test/FESpacesTests/FESolversTests.jl +++ b/test/FESpacesTests/FESolversTests.jl @@ -37,8 +37,10 @@ cellvec = integrate(v*f,trian,quad) cellids = collect(1:num_cells(trian)) assem = SparseMatrixAssembler(U,V) -A = assemble_matrix(assem,[cellmat],[cellids],[cellids]) -b = assemble_vector(assem,[cellvec],[cellids]) +matdata = ([cellmat],[cellids],[cellids]) +vecdata = ([cellvec],[cellids]) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ b x0 = zeros(length(x)) diff --git a/test/FESpacesTests/FETermsTests.jl b/test/FESpacesTests/FETermsTests.jl index 3d3223a27..0de3bc6d4 100644 --- a/test/FESpacesTests/FETermsTests.jl +++ b/test/FESpacesTests/FETermsTests.jl @@ -61,24 +61,24 @@ t_affine = AffineFETerm(a,l,trian,quad) matdata = collect_cell_matrix(u,v,[t_affine, t_linear]) vecdata = collect_cell_vector(uhd,v,[t_affine, t_linear, t_source]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ b @test x ≈ get_free_values(uh) data = collect_cell_matrix_and_vector(uhd,u,v,[t_affine, t_linear, t_source]) -A, b = allocate_matrix_and_vector(assem,data...) -assemble_matrix_and_vector!(A,b,assem,data...) +A, b = allocate_matrix_and_vector(assem,data) +assemble_matrix_and_vector!(A,b,assem,data) x = A \ b @test x ≈ get_free_values(uh) -A, b = assemble_matrix_and_vector(assem,data...) +A, b = assemble_matrix_and_vector(assem,data) x = A \ b @test x ≈ get_free_values(uh) matdata = collect_cell_jacobian(uh,u,v,[t_affine]) vecdata = collect_cell_residual(uh,v,[t_affine]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ -b @test (x.+1) ≈ ones(length(x)) @@ -89,15 +89,15 @@ t_source = FESource(l,trian,quad) matdata = collect_cell_matrix(u,v,[t_linear, t_source]) vecdata = collect_cell_vector(uhd,v,[t_linear, t_source]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ b @test x ≈ get_free_values(uh) matdata = collect_cell_jacobian(uh,u,v,[t_linear, t_source]) vecdata = collect_cell_residual(uh,v,[t_linear, t_source]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ -b @test (x.+1) ≈ ones(length(x)) @@ -107,14 +107,14 @@ t_nonlinear = FETerm(r,j,trian,quad) matdata = collect_cell_jacobian(uh,u,v,[t_nonlinear]) vecdata = collect_cell_residual(uh,v,[t_nonlinear]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ -b @test (x.+1) ≈ ones(length(x)) data = collect_cell_jacobian_and_residual(uh,u,v,[t_nonlinear]) -A, b = allocate_matrix_and_vector(assem,data...) -assemble_matrix_and_vector!(A,b,assem,data...) +A, b = allocate_matrix_and_vector(assem,data) +assemble_matrix_and_vector!(A,b,assem,data) x = A \ -b @test (x.+1) ≈ ones(length(x)) @@ -154,14 +154,14 @@ t_matvec_Ω = AffineFETermFromCellMatVec(matvecfun,trian) matdata = collect_cell_matrix(u,v,[t_matvec_Ω,]) vecdata = collect_cell_vector(uhd,v,[t_matvec_Ω,]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ b @test x ≈ get_free_values(uh) data = collect_cell_matrix_and_vector(uhd,u,v,[t_matvec_Ω,]) -A, b = allocate_matrix_and_vector(assem,data...) -assemble_matrix_and_vector!(A,b,assem,data...) +A, b = allocate_matrix_and_vector(assem,data) +assemble_matrix_and_vector!(A,b,assem,data) x = A \ b @test x ≈ get_free_values(uh) @@ -195,14 +195,14 @@ t_jacres_Ω = FETermFromCellJacRes(jacresfun,trian) matdata = collect_cell_jacobian(uh,u,v,[t_jacres_Ω]) vecdata = collect_cell_residual(uh,v,[t_jacres_Ω]) -A = assemble_matrix(assem,matdata...) -b = assemble_vector(assem,vecdata...) +A = assemble_matrix(assem,matdata) +b = assemble_vector(assem,vecdata) x = A \ -b @test (x.+1) ≈ ones(length(x)) data = collect_cell_jacobian_and_residual(uh,u,v,[t_jacres_Ω]) -A, b = allocate_matrix_and_vector(assem,data...) -assemble_matrix_and_vector!(A,b,assem,data...) +A, b = allocate_matrix_and_vector(assem,data) +assemble_matrix_and_vector!(A,b,assem,data) x = A \ -b @test (x.+1) ≈ ones(length(x)) diff --git a/test/FESpacesTests/SparseMatrixAssemblersTests.jl b/test/FESpacesTests/SparseMatrixAssemblersTests.jl index 01988f836..8e38657a7 100644 --- a/test/FESpacesTests/SparseMatrixAssemblersTests.jl +++ b/test/FESpacesTests/SparseMatrixAssemblersTests.jl @@ -55,9 +55,6 @@ term_to_cellvec = [cellvec, bcellvec] term_to_cellids = [cellids, bcellids] term_to_cellmatvec = [ cellmatvec, bcellmatvec ] -matvecdata = ( term_to_cellmatvec , term_to_cellids, term_to_cellids) -matdata = (term_to_cellmat,term_to_cellids,term_to_cellids) -vecdata = (term_to_cellvec,term_to_cellids) mtypes = [ SparseMatrixCSC, @@ -74,16 +71,24 @@ mtypes = [ for T in mtypes + matvecdata = ( term_to_cellmatvec , term_to_cellids, term_to_cellids) + matdata = (term_to_cellmat,term_to_cellids,term_to_cellids) + vecdata = (term_to_cellvec,term_to_cellids) + data = (matvecdata,matdata,vecdata) + assem = SparseMatrixAssembler(T,Vector{Float64},U,V) - test_assembler(assem,matvecdata,matdata,vecdata) + test_assembler(assem,matdata,vecdata,data) - mat = assemble_matrix(assem,[cellmat],[cellids],[cellids]) - vec = assemble_vector(assem,[cellvec],[cellids]) + matdata = ([cellmat],[cellids],[cellids]) + vecdata = ([cellvec],[cellids]) + + mat = assemble_matrix(assem,matdata) + vec = assemble_vector(assem,vecdata) x = mat \ vec - assemble_matrix!(mat,assem,[cellmat],[cellids],[cellids]) - assemble_vector!(vec,assem,[cellvec],[cellids]) + assemble_matrix!(mat,assem,matdata) + assemble_vector!(vec,assem,vecdata) x2 = mat \ vec @test x ≈ x2 @@ -97,9 +102,10 @@ for T in mtypes @test mat[2, 3] ≈ -0.33333333333333 @test mat[3, 3] ≈ 1.333333333333333 - mat, vec = allocate_matrix_and_vector(assem,([cellmatvec],[cellids],[cellids])) - assemble_matrix_and_vector!(mat,vec,assem,([cellmatvec],[cellids],[cellids])) - assemble_matrix_and_vector!(mat,vec,assem,([cellmatvec],[cellids],[cellids])) + data = (([cellmatvec],[cellids],[cellids]),([],[],[]),([],[])) + mat, vec = allocate_matrix_and_vector(assem,data) + assemble_matrix_and_vector!(mat,vec,assem,data) + assemble_matrix_and_vector!(mat,vec,assem,data) @test vec ≈ [0.0625, 0.125, 0.0625] @test mat[1, 1] ≈ 1.333333333333333 @@ -108,7 +114,7 @@ for T in mtypes x3 = mat \ vec @test x ≈ x3 - mat, vec = assemble_matrix_and_vector(assem,([cellmatvec],[cellids],[cellids])) + mat, vec = assemble_matrix_and_vector(assem,data) x4 = mat \ vec @test x ≈ x4 diff --git a/test/MultiFieldTests/MultiFieldSparseMatrixAssemblersTests.jl b/test/MultiFieldTests/MultiFieldSparseMatrixAssemblersTests.jl index f02755ed4..3a764fd17 100644 --- a/test/MultiFieldTests/MultiFieldSparseMatrixAssemblersTests.jl +++ b/test/MultiFieldTests/MultiFieldSparseMatrixAssemblersTests.jl @@ -43,19 +43,20 @@ cellvec = integrate(dv*2,trian,quad) cellids = get_cell_id(trian) cellmatvec = pair_arrays(cellmat,cellvec) -assem = SparseMatrixAssembler(SparseMatrixCSR{0},Y,X) +assem = SparseMatrixAssembler(SparseMatrixCSR{0},X,Y) matvecdata = ([cellmatvec],[cellids],[cellids]) matdata = ([cellmat],[cellids],[cellids]) vecdata = ([cellvec],[cellids]) +data = (matvecdata,matdata,vecdata) -test_assembler(assem,matvecdata,matdata,vecdata) +test_assembler(assem,matdata,vecdata,data) -A = assemble_matrix(assem,matdata...) +A = assemble_matrix(assem,matdata) -A = allocate_matrix(assem,matdata...) +A = allocate_matrix(assem,matdata) -assem = SparseMatrixAssembler(Y,X) -test_assembler(assem,matvecdata,matdata,vecdata) +assem = SparseMatrixAssembler(X,Y) +test_assembler(assem,matdata,vecdata,data) end # module From 048f404e7f0a4abed2a6087fef673721617ae819 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 7 May 2020 15:55:28 +0200 Subject: [PATCH 12/17] Moved inner to outer constructor in FEOperatorFromTerms --- src/FESpaces/Assemblers.jl | 2 -- src/FESpaces/FEOperatorsFromTerms.jl | 5 +---- src/FESpaces/FESpaces.jl | 3 +++ 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index 854917d4f..da1b2d65f 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -164,8 +164,6 @@ function test_assembler(a::Assembler,matdata,vecdata,data) @test num_free_dofs(trial_fesp) == size(A,2) @test num_free_dofs(test_fesp) == size(A,1) @test num_free_dofs(test_fesp) == length(b) - strategy = get_assembly_strategy(a) - @test isa(strategy,AssemblyStrategy) end # This is an extended interface that only make sense for assemblers that build (sequential) sparse matrices diff --git a/src/FESpaces/FEOperatorsFromTerms.jl b/src/FESpaces/FEOperatorsFromTerms.jl index b69e1e25d..8c6a88cfc 100644 --- a/src/FESpaces/FEOperatorsFromTerms.jl +++ b/src/FESpaces/FEOperatorsFromTerms.jl @@ -4,15 +4,12 @@ struct FEOperatorFromTerms <: FEOperator test::FESpace assem::Assembler terms - function FEOperatorFromTerms(trial::FESpace,test::FESpace,assem::Assembler,terms::FETerm...) - new(trial,test,assem,terms) - end end """ """ function FEOperator(trial::FESpace,test::FESpace,assem::Assembler,terms::FETerm...) - FEOperatorFromTerms(trial,test,assem,terms...) + FEOperatorFromTerms(trial,test,assem,terms) end """ diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index e59f6e467..971584c25 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -127,7 +127,9 @@ export test_assembler export get_matrix_type export get_vector_type export count_matrix_nnz_coo +export count_matrix_and_vector_nnz_coo export fill_matrix_coo_symbolic! +export fill_matrix_and_vector_coo_symbolic! export fill_matrix_coo_numeric! export fill_matrix_and_vector_coo_numeric! export test_sparse_matrix_assembler @@ -189,6 +191,7 @@ export SparseMatrixAssembler export FEOperator export test_fe_operator export AffineFEOperator +export FEOperatorFromTerms export get_algebraic_operator export FESolver From 1a101e8e2938eadbf73e7d1ac99d0ec450cf72ce Mon Sep 17 00:00:00 2001 From: amartin Date: Tue, 12 May 2020 10:51:59 +1000 Subject: [PATCH 13/17] Update NEWS.md --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index 7cb380a98..009838fb1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,11 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added + - Added inner constructor to CartesianDiscreteModel. Since PR [#245](https://github.com/gridap/Gridap.jl/pull/245) + ## [0.9.2] - 2020-4-26 ### Added From 61597536239cab2815fce1a771f35c99a79ada76 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 14 May 2020 12:50:47 +0200 Subject: [PATCH 14/17] Updating FESpacesWithLinearConstraints to new interface --- src/FESpaces/FESpacesWithLinearConstraints.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/FESpaces/FESpacesWithLinearConstraints.jl b/src/FESpaces/FESpacesWithLinearConstraints.jl index fd68bde86..1c2430c90 100644 --- a/src/FESpaces/FESpacesWithLinearConstraints.jl +++ b/src/FESpaces/FESpacesWithLinearConstraints.jl @@ -350,7 +350,7 @@ function get_dirichlet_values(f::FESpaceWithLinearConstraints) end function scatter_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,fmdof_to_val,dmdof_to_val) - fdof_to_val = zero_free_values(eltype(fmdof_to_val),f.space) + fdof_to_val = zero_free_values(f.space) ddof_to_val = zero_dirichlet_values(f.space) _setup_dof_to_val!( fdof_to_val, @@ -406,7 +406,7 @@ end function gather_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,cell_to_ludof_to_val) fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ludof_to_val) - fmdof_to_val = zero_free_values(eltype(fdof_to_val),f) + fmdof_to_val = zero_free_values(f) dmdof_to_val = zero_dirichlet_values(f) _setup_mdof_to_val!( fmdof_to_val, @@ -458,8 +458,8 @@ function get_cell_basis(f::FESpaceWithLinearConstraints) get_cell_basis(f.space) end -function zero_free_values(::Type{T},f::FESpaceWithLinearConstraints) where T - zeros(T,num_free_dofs(f)) +function zero_free_values(f::FESpaceWithLinearConstraints) where T + zeros(num_free_dofs(f)) end constraint_style(::Type{<:FESpaceWithLinearConstraints}) = Val{true}() From 5d38d2e1e68ad1e533d4a92351c13464df8fea42 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 14 May 2020 13:10:01 +0200 Subject: [PATCH 15/17] Extending TriangulationPortion to boundary case --- src/Geometry/SkeletonTriangulations.jl | 3 ++- src/Geometry/TriangulationPortions.jl | 12 ++++++++++++ .../TriangulationPortionsTests.jl | 19 +++++++++++++++++-- 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/src/Geometry/SkeletonTriangulations.jl b/src/Geometry/SkeletonTriangulations.jl index d883aaf53..97500c669 100644 --- a/src/Geometry/SkeletonTriangulations.jl +++ b/src/Geometry/SkeletonTriangulations.jl @@ -10,9 +10,10 @@ The inner constructor enforces `B<:BoundaryTriangulation` struct SkeletonTriangulation{Dc,Dp,B} <: Triangulation{Dc,Dp} left::B right::B - function SkeletonTriangulation(left::B,right::B) where B<:BoundaryTriangulation + function SkeletonTriangulation(left::B,right::B) where B<:Triangulation Dc = num_cell_dims(left) Dp = num_point_dims(left) + @assert Dc + 1 == Dp new{Dc,Dp,B}(left,right) end end diff --git a/src/Geometry/TriangulationPortions.jl b/src/Geometry/TriangulationPortions.jl index ffdf766ca..79562c214 100644 --- a/src/Geometry/TriangulationPortions.jl +++ b/src/Geometry/TriangulationPortions.jl @@ -27,3 +27,15 @@ function get_cell_map(trian::TriangulationPortion) cell_map = get_cell_map(trian.oldtrian) reindex(cell_map,trian.cell_to_oldcell) end + +function get_cell_id(trian::TriangulationPortion) + reindex(get_cell_id(trian.oldtrian),trian.cell_to_oldcell) +end + +function restrict(f::AbstractArray,trian::TriangulationPortion) + reindex(f,trian) +end + +function get_normal_vector(trian::TriangulationPortion) + reindex(get_normal_vector(trian.oldtrian),trian.cell_to_oldcell) +end diff --git a/test/GeometryTests/TriangulationPortionsTests.jl b/test/GeometryTests/TriangulationPortionsTests.jl index f808e3fc0..e85deb1ce 100644 --- a/test/GeometryTests/TriangulationPortionsTests.jl +++ b/test/GeometryTests/TriangulationPortionsTests.jl @@ -2,14 +2,29 @@ module TriangulationPortionsTests using Gridap.ReferenceFEs using Gridap.Geometry +using Test domain = (0,1,0,1) partition = (10,10) -oldgrid = CartesianGrid(domain,partition) +oldmodel = CartesianDiscreteModel(domain,partition) +oldtrian = get_triangulation(oldmodel) cell_to_oldcell = collect(1:34) -trian = TriangulationPortion(oldgrid,cell_to_oldcell) +trian = TriangulationPortion(oldtrian,cell_to_oldcell) test_triangulation(trian) +oldbtrian = BoundaryTriangulation(oldmodel) + +bface_to_oldbface = collect(1:32) +btrian = TriangulationPortion(oldbtrian,bface_to_oldbface) +test_triangulation(btrian) + +nb = get_normal_vector(btrian) +@test length(nb) == num_cells(btrian) + +#using Gridap.Visualization +#writevtk(oldtrian,"oldtrian") +#writevtk(btrian,"btrian",cellfields=["normal"=>nb],celldata=["oldcell"=>get_cell_id(btrian)]) + end # module From 6c32cca87edfcb67287dde19b97472fa052d4c0a Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 14 May 2020 13:20:32 +0200 Subject: [PATCH 16/17] Extended TriangulationPortion to skeleton case --- src/Geometry/SkeletonTriangulations.jl | 6 ++++++ test/GeometryTests/TriangulationPortionsTests.jl | 10 ++++++++++ 2 files changed, 16 insertions(+) diff --git a/src/Geometry/SkeletonTriangulations.jl b/src/Geometry/SkeletonTriangulations.jl index 97500c669..c607b803a 100644 --- a/src/Geometry/SkeletonTriangulations.jl +++ b/src/Geometry/SkeletonTriangulations.jl @@ -165,6 +165,12 @@ function get_normal_vector(trian::SkeletonTriangulation) get_normal_vector(trian.left) end +function TriangulationPortion(oldtrian::SkeletonTriangulation,cell_to_oldcell::Vector{Int}) + left = TriangulationPortion(oldtrian.left,cell_to_oldcell) + right = TriangulationPortion(oldtrian.right,cell_to_oldcell) + SkeletonTriangulation(left,right) +end + # Specific API """ diff --git a/test/GeometryTests/TriangulationPortionsTests.jl b/test/GeometryTests/TriangulationPortionsTests.jl index e85deb1ce..a033b0756 100644 --- a/test/GeometryTests/TriangulationPortionsTests.jl +++ b/test/GeometryTests/TriangulationPortionsTests.jl @@ -22,9 +22,19 @@ test_triangulation(btrian) nb = get_normal_vector(btrian) @test length(nb) == num_cells(btrian) +sface_to_oldsface = collect(10:45) +oldstrian = SkeletonTriangulation(oldmodel) +strian = TriangulationPortion(oldstrian,sface_to_oldsface) +test_triangulation(strian) + +ns = get_normal_vector(strian) +@test length(ns) == num_cells(strian) + #using Gridap.Visualization #writevtk(oldtrian,"oldtrian") #writevtk(btrian,"btrian",cellfields=["normal"=>nb],celldata=["oldcell"=>get_cell_id(btrian)]) +#writevtk(strian,"strian",cellfields=["normal"=>ns], +# celldata=["oldcell_left"=>get_cell_id(strian).left,"oldcell_right"=>get_cell_id(strian).right]) end # module From 1b6161f0b425aedf4a022bc33ccb26290dd79c76 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 14 May 2020 13:38:44 +0200 Subject: [PATCH 17/17] Updated NEWS.md --- NEWS.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/NEWS.md b/NEWS.md index eb4f9f20d..767cc3998 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added + - Extended support of `TriangulationPortion` to boundary and skeleton triangulations. Since PR ... - Added `FESpaceWithLinearConstraints`. Since PR [#247](https://github.com/gridap/Gridap.jl/pull/247). - Added inner constructor to `CartesianDiscreteModel` allowing to build a model that represents a subgrid of a larger grid. Since PR [#245](https://github.com/gridap/Gridap.jl/pull/245). @@ -15,6 +16,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - The part associated with the imposition of constraints in the `FESpace` interface has changed slightly. Since PR [#247](https://github.com/gridap/Gridap.jl/pull/247). + - Simplified the signature of `zero_free_values(::FESpace)`. Since PR ... + - Major refactoring in the `Assembler` interface. + **Important change:** assembly-related functions take the data returned by functions like + `collect_cell_matrix` as it is. Example: the old user code `assemble_matrix(assembler,collect_cell_matrix(du,dv,terms)...)` + now is written simply as `assemble_matrix(assembler,collect_cell_matrix(du,dv,terms))`, i.e., the unpack of the last argument is not + used anymore. With new new assembler interface, it is also possible to customize the assembly process + via a so-called `AssemblerStrategy` object. Since PR ... - Change the types of the sizes and partition fields of CartesianDescriptor to tuples instead of points. Since PR [#246](https://github.com/gridap/Gridap.jl/pull/246).