Try #1399:

.buildkite/pipeline.yml

@@ -528,6 +528,16 @@ steps:
           slurm_gpus: 1
           slurm_mem: 20GB
 
+      - label: "Unit: operator matrices (CPU)"
+        key: unit_operator_matrices_cpu
+        command: "julia --color=yes --check-bounds=yes --project=test test/MatrixFields/operator_matrices.jl"
+
+      - label: "Unit: operator matrices (GPU)"
+        key: unit_operator_matrices_gpu
+        command: "julia --color=yes --project=test test/MatrixFields/operator_matrices.jl"
+        agents:
+          slurm_gpus: 1
+
   - group: "Unit: Hypsography"
     steps:
 

src/MatrixFields/MatrixFields.jl

@@ -52,6 +52,8 @@ include("band_matrix_row.jl")
 include("rmul_with_projection.jl")
 include("matrix_shape.jl")
 include("matrix_multiplication.jl")
+include("lazy_operators.jl")
+include("operator_matrices.jl")
 include("field2arrays.jl")
 
 const ColumnwiseBandMatrixField{V, S} = Fields.Field{

src/MatrixFields/band_matrix_row.jl

@@ -12,6 +12,8 @@ several aliases for commonly used subtypes of `BandMatrixRow`:
 - `TridiagonalMatrixRow(entry_1, entry_2, entry_3)`
 - `QuaddiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4)`
 - `PentadiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4, entry_5)`
+It is also possible to construct band matrix rows from other band matrix rows,
+e.g., `QuaddiagonalMatrixRow(BidiagonalMatrixRow(entry_1, entry_2))`.
 """
 struct BandMatrixRow{ld, bw, T} # bw is the bandwidth (the number of diagonals)
     entries::NTuple{bw, T}
@@ -25,6 +27,13 @@ BandMatrixRow{ld}(entries::Vararg{Any, bw}) where {ld, bw} =
 BandMatrixRow{ld, bw}(entries::Vararg{Any, bw}) where {ld, bw} =
     BandMatrixRow{ld, bw, rpromote_type(map(typeof, entries)...)}(entries)
 
+BandMatrixRow{ld, bw}(row::BandMatrixRow) where {ld, bw} =
+    BandMatrixRow{ld, bw, eltype(row)}(row)
+BandMatrixRow{ld, bw, T}(row::BandMatrixRow) where {ld, bw, T} =
+    convert(BandMatrixRow{ld, bw, T}, row)
+
+BandMatrixRow{ld, 0}() where {ld} = BandMatrixRow{ld, 0, Nothing}(())
+
 const DiagonalMatrixRow{T} = BandMatrixRow{0, 1, T}
 const BidiagonalMatrixRow{T} = BandMatrixRow{-1 + half, 2, T}
 const TridiagonalMatrixRow{T} = BandMatrixRow{-1, 3, T}

src/MatrixFields/lazy_operators.jl

@@ -0,0 +1,84 @@
+abstract type AbstractLazyOperator end
+
+abstract type AbstractLazyOperatorStyle <: Base.Broadcast.BroadcastStyle end
+struct LazyOperatorStyle <: AbstractLazyOperatorStyle end
+struct LazyOperatorStyleWithSpaceAvailable <: AbstractLazyOperatorStyle end
+
+Base.Broadcast.broadcastable(lazy_op::AbstractLazyOperator) = lazy_op
+
+Base.Broadcast.BroadcastStyle(::Type{<:AbstractLazyOperator}) =
+    LazyOperatorStyle()
+
+# Broadcasting over an AbstractLazyOperator and either a Ref or a Tuple
+Base.Broadcast.BroadcastStyle(
+    style::AbstractLazyOperatorStyle,
+    ::Union{Base.Broadcast.AbstractArrayStyle{0}, Base.Broadcast.Style{Tuple}},
+) = style
+
+# Broadcasting over an AbstractLazyOperator and at least one Field
+Base.Broadcast.BroadcastStyle(
+    ::AbstractLazyOperatorStyle,
+    ::Fields.AbstractFieldStyle,
+) = LazyOperatorStyleWithSpaceAvailable()
+
+# Broadcasting over more than one AbstractLazyOperator
+Base.Broadcast.BroadcastStyle(
+    ::LazyOperatorStyle,
+    ::LazyOperatorStyleWithSpaceAvailable,
+) = LazyOperatorStyleWithSpaceAvailable()
+
+Base.Broadcast.materialize(::Base.Broadcast.Broadcasted{LazyOperatorStyle}) =
+    error("Cannot materialize lazy operator broadcast because it does not \
+           contain any Fields with spatial information")
+
+Base.Broadcast.materialize(
+    bc::Base.Broadcast.Broadcasted{LazyOperatorStyleWithSpaceAvailable},
+) = Base.Broadcast.materialize(remove_lazy_operators(largest_space(bc), bc))
+
+Base.Broadcast.materialize!(
+    ::AbstractLazyOperatorStyle,
+    dest::Fields.Field,
+    bc::Base.Broadcast.Broadcasted,
+) = Base.Broadcast.materialize!(dest, remove_lazy_operators(axes(dest), bc))
+
+remove_lazy_operators(_, arg) = arg
+remove_lazy_operators(space, bc::Base.Broadcast.Broadcasted) =
+    Base.Broadcast.broadcasted(
+        bc.f,
+        map(arg -> remove_lazy_operators(space, arg), bc.args)...,
+        # remove_lazy_operators_args(space, bc.args...)...,
+    )
+remove_lazy_operators(_, lazy_op::AbstractLazyOperator) =
+    error("Missing method: remove_lazy_operators(space, ::$(typeof(lazy_op)))")
+
+# remove_lazy_operators_args(_) = ()
+# remove_lazy_operators_args(space, arg, args...) = (
+#     remove_lazy_operators(space, arg),
+#     remove_lazy_operators_args(space, args...)...,
+# )
+
+largest_space(_) = nothing
+largest_space(field::Fields.Field) = axes(field)
+largest_space(bc::Base.Broadcast.Broadcasted) = largest_space(bc.args...)
+largest_space(arg, args...) =
+    larger_space(largest_space(arg), largest_space(args...))
+
+larger_space(space1, ::Nothing) = space1
+larger_space(::Nothing, space2) = space2
+larger_space(space1::S, ::S) where {S} = space1 # Neither space is larger.
+larger_space(
+    space1::Spaces.ExtrudedFiniteDifferenceSpace,
+    ::Spaces.ExtrudedFiniteDifferenceSpace,
+) = space1 # The staggering does not matter here, so neither space is larger.
+larger_space(
+    space1::Spaces.FiniteDifferenceSpace,
+    ::Spaces.FiniteDifferenceSpace,
+) = space1 # The staggering does not matter here, so neither space is larger.
+larger_space(
+    space1::Spaces.ExtrudedFiniteDifferenceSpace,
+    ::Spaces.AbstractSpectralElementSpace,
+) = space1 # The types indicate that space2 is a subspace of space1.
+larger_space(
+    ::Spaces.AbstractSpectralElementSpace,
+    space2::Spaces.ExtrudedFiniteDifferenceSpace,
+) = space2 # The types indicate that space1 is a subspace of space2.

src/MatrixFields/matrix_multiplication.jl

@@ -266,7 +266,7 @@ function multiply_matrix_at_index(
             if is_interior || ld1_or_boundary_ld1 <= d <= ud1_or_boundary_ud1
                 @inbounds Operators.getidx(space, matrix2, loc, idx + d, hidx)
             else
-                rzero(eltype(matrix2)) # This value is never used.
+                zero(eltype(matrix2)) # This value is never used.
             end
         end # The rows are precomputed to avoid recomputing them multiple times.
         matrix2_rows_wrapper = BandMatrixRow{ld1}(matrix2_rows...)