From 1e3391659ea447a15679182ec7b8d53dcf758b87 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Sat, 16 Sep 2023 18:02:18 -0400 Subject: [PATCH] generic tensor product map now uses generic matrix --- Project.toml | 4 +- .../MatrixFreeOperators.jl | 13 +- src/MatrixFreeOperators/tensor_product_2d.jl | 132 ----------- src/MatrixFreeOperators/tensor_product_3d.jl | 223 ------------------ test/runtests.jl | 2 +- test/scaling_test_euler_2d.sh | 14 -- 6 files changed, 10 insertions(+), 378 deletions(-) delete mode 100644 src/MatrixFreeOperators/tensor_product_2d.jl delete mode 100644 src/MatrixFreeOperators/tensor_product_3d.jl delete mode 100644 test/scaling_test_euler_2d.sh diff --git a/Project.toml b/Project.toml index 53364762..47026efa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StableSpectralElements" uuid = "fb992021-99c7-4c2d-a14b-5e48ac4045b2" authors = ["Tristan Montoya "] -version = "0.2.5" +version = "0.2.6" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" @@ -69,4 +69,4 @@ TetGen = "1" TimerOutputs = "0.5" Triangulate = "2" WriteVTK = "1" -julia = "1.8" +julia = "1.8" \ No newline at end of file diff --git a/src/MatrixFreeOperators/MatrixFreeOperators.jl b/src/MatrixFreeOperators/MatrixFreeOperators.jl index 6eb861ca..fee83885 100644 --- a/src/MatrixFreeOperators/MatrixFreeOperators.jl +++ b/src/MatrixFreeOperators/MatrixFreeOperators.jl @@ -9,10 +9,6 @@ module MatrixFreeOperators struct GenericMatrixAlgorithm <: AbstractOperatorAlgorithm end struct GenericTensorProductAlgorithm <: AbstractOperatorAlgorithm end - export TensorProductMap2D, TensorProductMap3D - include("tensor_product_2d.jl") - include("tensor_product_3d.jl") - export WarpedTensorProductMap2D, WarpedTensorProductMap3D include("warped_product_2d.jl") include("warped_product_3d.jl") @@ -64,9 +60,14 @@ module MatrixFreeOperators alg::GenericTensorProductAlgorithm) = vcat( [make_operator(block, alg) for block in map.maps]...) make_operator(map::LinearMaps.KroneckerMap{Float64, <:NTuple{2,LinearMap}}, - ::GenericTensorProductAlgorithm) = TensorProductMap2D(map.maps...) + ::GenericTensorProductAlgorithm) = + make_operator(map.maps[1], GenericMatrixAlgorithm()) ⊗ + make_operator(map.maps[2], GenericMatrixAlgorithm()) make_operator(map::LinearMaps.KroneckerMap{Float64, <:NTuple{3,LinearMap}}, - ::GenericTensorProductAlgorithm) = TensorProductMap3D(map.maps...) + ::GenericTensorProductAlgorithm) = + make_operator(map.maps[1], GenericMatrixAlgorithm()) ⊗ + make_operator(map.maps[2], GenericMatrixAlgorithm()) ⊗ + make_operator(map.maps[3], GenericMatrixAlgorithm()) # count adds, muls, and muladds function count_ops(map::LinearMap) diff --git a/src/MatrixFreeOperators/tensor_product_2d.jl b/src/MatrixFreeOperators/tensor_product_2d.jl deleted file mode 100644 index 1884026c..00000000 --- a/src/MatrixFreeOperators/tensor_product_2d.jl +++ /dev/null @@ -1,132 +0,0 @@ -struct TensorProductMap2D{A_type,B_type,σᵢ_type, - σₒ_type} <: LinearMaps.LinearMap{Float64} - A::A_type - B::B_type - σᵢ::σᵢ_type - σₒ::σₒ_type -end - -@inline Base.size(L::TensorProductMap2D) = (size(L.σₒ,1)*size(L.σₒ,2), - size(L.σᵢ,1)*size(L.σᵢ,2)) - -function TensorProductMap2D(A, B) - (M1,N1) = size(A) - (M2,N2) = size(B) - σₒ = SMatrix{M1,M2,Int}([M2*(α1-1) + α2 for α1 in 1:M1, α2 in 1:M2]) - σᵢ = SMatrix{N1,N2,Int}([N2*(β1-1) + β2 for β1 in 1:N1, β2 in 1:N2]) - - if A isa LinearMaps.UniformScalingMap{Bool} - A = I - elseif A isa Union{LinearMaps.WrappedMap,OctavianMap} - A = SMatrix{M1,N1,Float64}(A.lmap) - end - if B isa LinearMaps.UniformScalingMap{Bool} - B = I - elseif B isa Union{LinearMaps.WrappedMap,OctavianMap} - B = SMatrix{M2,N2,Float64}(B.lmap) - end - return TensorProductMap2D(A,B, σᵢ, σₒ) -end - -""" -Return the transpose using - -(A ⊗ B)ᵀ = Aᵀ ⊗ Bᵀ -""" -function LinearAlgebra.transpose(L::TensorProductMap2D) - return TensorProductMap2D(transpose(L.A), transpose(L.B), L.σₒ, L.σᵢ) -end - -""" -Multiply a vector of length N1*N2 by the M1*M2 x N1*N2 matrix - -L[σₒ[α1,α2], σᵢ[β1,β2]] = A[α1,β1] B[α2,β2] - -or L = A ⊗ B. The action of this matrix on a vector x is - -(Lx)[σₒ[α1,α2]] = ∑_{β1,β2} A[α1,β1] B[α2,β2] x[σᵢ[β1,β2]] - = ∑_{β1} A[α1,β1] (∑_{β2} B[α2,β2] x[σᵢ[β1,β2]]) - = ∑_{β1} A[α1,β1] Z[α2,β1] -""" -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap2D{<:AbstractMatrix{Float64},<:AbstractMatrix{Float64}}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, σᵢ, σₒ) = L - - Z = MMatrix{size(σᵢ,1), size(σₒ,2),Float64}(undef) - - @inbounds for α2 in axes(σₒ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β2 in axes(σᵢ,2) - @muladd temp = temp + B[α2,β2] * x[σᵢ[β1,β2]] - end - Z[β1,α2] = temp - end - - @inbounds for α1 in axes(σₒ,1), α2 in axes(σₒ,2) - temp = 0.0 - @simd for β1 in axes(σᵢ,1) - @muladd temp = temp + A[α1,β1] * Z[β1,α2] - end - y[σₒ[α1,α2]] = temp - end - - return y -end - -""" -Multiply a vector of length N1*N2 by the M1*M2 x N1*M2 matrix - -L[σₒ[α1,α2], σᵢ[β1,β2]] = A[α1,β1] δ_{α2,β2} - -or L = A ⊗ I_{M2}. The action of this matrix on a vector x is - -(Lx)[σₒ[α1,α2]] = ∑_{β1,β2} A[α1,β1] δ_{α2,β2} x[σᵢ[β1,β2]] - = ∑_{β1} A[α1,β1] x[σᵢ[β1,α2]] -""" -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap2D{<:AbstractMatrix{Float64},<:UniformScaling}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, σᵢ, σₒ) = L - - @inbounds for α1 in axes(σₒ,1), α2 in axes(σₒ,2) - temp = 0.0 - @simd for β1 in axes(σᵢ,1) - @muladd temp = temp + A[α1,β1] * x[σᵢ[β1,α2]] - end - y[σₒ[α1,α2]] = temp - end - - return lmul!(B,y) -end - -""" -Multiply a vector of length N1*N2 by the M1*M2 x M1*N2 matrix -L[σₒ[α1,α2], σᵢ[β1,β2]] = δ_{α1,β1} B[α2,β2] - -or L = I_{M1} ⊗ B. The action of this matrix on a vector x is - -(Lx)[σₒ[α1,α2]] = ∑_{β1,β2} δ_{α1,β1} B[α2,β2] x[σᵢ[β1,β2]] - = ∑_{β2} B[α2,β2] x[σᵢ[α1,β2]]) -""" -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap2D{<:UniformScaling,<:AbstractMatrix{Float64}}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, σᵢ, σₒ) = L - - @inbounds for α1 in axes(σₒ,1), α2 in axes(σₒ,2) - temp = 0.0 - @simd for β2 in axes(σᵢ,2) - @muladd temp = temp + B[α2,β2] * x[σᵢ[α1,β2]] - end - y[σₒ[α1,α2]] = temp - end - - return lmul!(A,y) -end \ No newline at end of file diff --git a/src/MatrixFreeOperators/tensor_product_3d.jl b/src/MatrixFreeOperators/tensor_product_3d.jl deleted file mode 100644 index 1a9256fb..00000000 --- a/src/MatrixFreeOperators/tensor_product_3d.jl +++ /dev/null @@ -1,223 +0,0 @@ -struct TensorProductMap3D{A_type,B_type,C_type,σᵢ_type, - σₒ_type} <: LinearMaps.LinearMap{Float64} - A::A_type - B::B_type - C::C_type - σᵢ::σᵢ_type - σₒ::σₒ_type -end - -@inline Base.size(L::TensorProductMap3D) = ( - size(L.σₒ,1)*size(L.σₒ,2)*size(L.σₒ,3), - size(L.σᵢ,1)*size(L.σᵢ,2)*size(L.σᵢ,3)) - -function TensorProductMap3D(A, B, C) - (M1,N1) = size(A) - (M2,N2) = size(B) - (M3,N3) = size(C) - σₒ = SArray{Tuple{M1,M2,M3},Int}( - [M2*M3*(α1-1) + M3*(α2-1) + α3 for α1 in 1:M1, α2 in 1:M2, α3 in 1:M3]) - σᵢ = SArray{Tuple{N1,N2,N3},Int}( - [N2*N3*(β1-1) + N3*(β2-1) + β3 for β1 in 1:N1, β2 in 1:N2, β3 in 1:N3]) - - if A isa LinearMaps.UniformScalingMap{Bool} - A = I - elseif A isa Union{LinearMaps.WrappedMap,OctavianMap} - A = SMatrix{M1,N1,Float64}(A.lmap) - end - if B isa LinearMaps.UniformScalingMap{Bool} - B = I - elseif B isa Union{LinearMaps.WrappedMap,OctavianMap} - B = SMatrix{M2,N2,Float64}(B.lmap) - end - if C isa LinearMaps.UniformScalingMap{Bool} - C = I - elseif C isa Union{LinearMaps.WrappedMap,OctavianMap} - C = SMatrix{M3,N3,Float64}(C.lmap) - end - return TensorProductMap3D(A, B, C, σᵢ, σₒ) -end - -""" -Compute transpose using -(A ⊗ B ⊗ C)ᵀ = Aᵀ ⊗ Bᵀ ⊗ Cᵀ -""" -function LinearAlgebra.transpose(L::TensorProductMap3D) - return TensorProductMap3D(transpose(L.A), transpose(L.B), - transpose(L.C), L.σₒ, L.σᵢ) -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:AbstractMatrix{Float64},<:AbstractMatrix{Float64},<:AbstractMatrix{Float64}}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - Z_3 = MArray{Tuple{size(σᵢ,1), size(σᵢ,2), size(σₒ,3)},Float64}(undef) - Z_2 = MArray{Tuple{size(σᵢ,1), size(σₒ,2), size(σₒ,3)},Float64}(undef) - - @inbounds for α3 in axes(σₒ,3), β2 in axes(σᵢ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β3 in axes(σᵢ,3) - @muladd temp = temp + C[α3,β3] * x[σᵢ[β1,β2,β3]] - end - Z_3[β1,β2,α3] = temp - end - - @inbounds for α3 in axes(σₒ,3), α2 in axes(σₒ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β2 in axes(σᵢ,2) - @muladd temp = temp + B[α2,β2] * Z_3[β1,β2,α3] - end - Z_2[β1,α2,α3] = temp - end - - @inbounds for α3 in axes(σₒ,3), α2 in axes(σₒ,2), α1 in axes(σₒ,1) - temp = 0.0 - @simd for β1 in axes(σᵢ,1) - @muladd temp = temp + A[α1,β1] * Z_2[β1,α2,α3] - end - y[σₒ[α1,α2,α3]] = temp - end - - return y -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:AbstractMatrix{Float64},<:UniformScaling,<:UniformScaling}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - @inbounds for α1 in axes(σₒ,1), α2 in axes(σₒ,2), α3 in axes(σₒ,3) - temp = 0.0 - @simd for β1 in axes(σᵢ,1) - @muladd temp = temp + A[α1,β1] * x[σᵢ[β1,α2,α3]] - end - y[σₒ[α1,α2,α3]] = temp - end - - return lmul!(B*C,y) -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:UniformScaling,<:AbstractMatrix{Float64},<:UniformScaling}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - @inbounds for α1 in axes(σₒ,1), α2 in axes(σₒ,2), α3 in axes(σₒ,3) - temp = 0.0 - @simd for β2 in axes(σᵢ,2) - @muladd temp = temp + B[α2,β2] * x[σᵢ[α1,β2,α3]] - end - y[σₒ[α1,α2,α3]] = temp - end - - return lmul!(A*C,y) -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:UniformScaling,<:UniformScaling,<:AbstractMatrix{Float64}}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - @inbounds for α1 in axes(σₒ,1), α2 in axes(σₒ,2), α3 in axes(σₒ,3) - temp = 0.0 - @simd for β3 in axes(σᵢ,3) - @muladd temp = temp + C[α3,β3] * x[σᵢ[α1,α2,β3]] - end - y[σₒ[α1,α2,α3]] = temp - end - - return lmul!(A*B,y) -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:AbstractMatrix{Float64},<:AbstractMatrix{Float64},<:UniformScaling}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - Z_2 = MArray{Tuple{size(σᵢ,1), size(σₒ,2), size(σₒ,3)},Float64}(undef) - - @inbounds for α3 in axes(σₒ,3), α2 in axes(σₒ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β2 in axes(σᵢ,2) - @muladd temp = temp + B[α2,β2] * x[σᵢ[β1,β2,α3]] - end - Z_2[β1,α2,α3] = temp - end - - @inbounds for α3 in axes(σₒ,3), α2 in axes(σₒ,2), α1 in axes(σₒ,1) - temp = 0.0 - @simd for β1 in axes(σᵢ,1) - @muladd temp = temp + A[α1,β1] * Z_2[β1,α2,α3] - end - y[σₒ[α1,α2,α3]] = temp - end - - return lmul!(C,y) -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:AbstractMatrix{Float64},<:UniformScaling,<:AbstractMatrix{Float64}}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - Z_3 = MArray{Tuple{size(σᵢ,1), size(σᵢ,2), size(σₒ,3)},Float64}(undef) - - @inbounds for α3 in axes(σₒ,3), β2 in axes(σᵢ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β3 in axes(σᵢ,3) - @muladd temp = temp + C[α3,β3] * x[σᵢ[β1,β2,β3]] - end - Z_3[β1,β2,α3] = temp - end - - @inbounds for α3 in axes(σₒ,3), α2 in axes(σₒ,2), α1 in axes(σₒ,1) - temp = 0.0 - @simd for β1 in axes(σᵢ,1) - @muladd temp = temp + A[α1,β1] * Z_3[β1,α2,α3] - end - y[σₒ[α1,α2,α3]] = temp - end - - return lmul!(B, y) -end - -@inline function LinearAlgebra.mul!(y::AbstractVector{Float64}, - L::TensorProductMap3D{<:UniformScaling,<:AbstractMatrix{Float64},<:AbstractMatrix{Float64}}, - x::AbstractVector{Float64}) - - LinearMaps.check_dim_mul(y, L, x) - (; A, B, C, σᵢ, σₒ) = L - - Z_3 = MArray{Tuple{size(σᵢ,1), size(σᵢ,2), size(σₒ,3)},Float64}(undef) - - @inbounds for α3 in axes(σₒ,3), β2 in axes(σᵢ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β3 in axes(σᵢ,3) - @muladd temp = temp + C[α3,β3] * x[σᵢ[β1,β2,β3]] - end - Z_3[β1,β2,α3] = temp - end - - @inbounds for α3 in axes(σₒ,3), α2 in axes(σₒ,2), β1 in axes(σᵢ,1) - temp = 0.0 - @simd for β2 in axes(σᵢ,2) - @muladd temp = temp + B[α2,β2] * Z_3[β1,β2,α3] - end - y[σₒ[β1,α2,α3]] = temp - end - - return lmul!(A,y) -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 2a071c90..7d4096e6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,7 +47,7 @@ end LinearAdvectionEquation((1.0,1.0)), InitialDataSine(1.0,(2*π, 2*π)), FluxDifferencingForm(), - ReferenceOperator(), DefaultOperatorAlgorithm(), + ReferenceOperator(), GenericTensorProductAlgorithm(), 1.0, 2, 1.0, 1.0/100.0, 0.1, "test_advection_2d_dgsem_fluxdiff") diff --git a/test/scaling_test_euler_2d.sh b/test/scaling_test_euler_2d.sh deleted file mode 100644 index d172fb5c..00000000 --- a/test/scaling_test_euler_2d.sh +++ /dev/null @@ -1,14 +0,0 @@ -#!/bin/bash -#SBATCH --nodes=1 -#SBATCH --cpus-per-task=40 -#SBATCH --time=6:00:00 -#SBATCH --job-name scaling_test - -module load NiaEnv/2019b -cd /scratch/z/zingg/tmontoya/ScalingTests/scripts -export OPENBLAS_NUM_THREADS=1 - -for nthreads in 1 2 4 5 8 10 20 40 80 -do - julia --project=.. --threads $nthreads --check-bounds=no -e "using StableSpectralElements; scaling_test_euler_2d(4,16,"./results/p4M16/", Threaded())" > screen.txt -done \ No newline at end of file