Merge #1424

1424: Fix matrix multiplication at boundaries r=dennisYatunin a=dennisYatunin

Peeled off from #1399. Fixes a bug in matrix-matrix multiplication (i.e., the `MultiplyColumnwiseBandMatrixField()` operator with two matrix fields as inputs) that was causing entries from outside the second input matrix to be used in the computation of entries outside the product matrix. These entries are used for padding, allowing us to store the nonzero entries of band matrices in `Field`s with rectangular parent arrays, so using them during the computation of matrix-matrix products leads to an unnecessary performance loss. Moreover, it breaks the assumption that these outside entries have no effect on program behavior.

In addition to fixing the bug, this PR
- Updates the documentation of `MultiplyColumnwiseBandMatrixField` to reflect the change.
- Adds a unit test which ensures that outside entries are never used during matrix-matrix multiplication.
- Simplifies some of the code for `MultiplyColumnwiseBandMatrixField`, making sure that it reflects the updated documentation as closely as possible.
- Introduces the `column_axes(matrix_field)` utility function, which can be used to obtain the space that corresponds to the columns of `matrix_field` (as opposed to `axes(matrix_field)`, which corresponds to the rows). This function simplifies the computation of interior/boundary indices during matrix multiplication.
- Fixes a type instability in `field2arrays` and updates the corresponding unit test. This also decreases the memory allocated by the new unit test for outside entries.
- Modifies the `show` method for matrix `Field`s so that it also displays the matrix shape.


Co-authored-by: Dennis Yatunin <[email protected]>

.buildkite/pipeline.yml

@@ -513,6 +513,10 @@ steps:
         key: unit_field2arrays
         command: "julia --color=yes --check-bounds=yes --project=test test/MatrixFields/field2arrays.jl"
 
+      - label: "Unit: matrix multiplication at boundaries"
+        key: unit_matrix_multiplication_at_boundaries
+        command: "julia --color=yes --check-bounds=yes --project=test test/MatrixFields/matrix_multiplication_at_boundaries.jl"
+
       - label: "Unit: matrix field broadcasting (CPU)"
         key: unit_matrix_field_broadcasting_cpu
         command: "julia --color=yes --check-bounds=yes --project=test test/MatrixFields/matrix_field_broadcasting.jl"

docs/src/matrix_fields.md

@@ -30,6 +30,7 @@ rmul_with_projection
 mul_return_type
 rmul_return_type
 matrix_shape
+column_axes
 ```
 
 ## Utilities

src/MatrixFields/MatrixFields.jl

@@ -68,10 +68,11 @@ const ColumnwiseBandMatrixField{V, S} = Fields.Field{
 function Base.show(io::IO, field::ColumnwiseBandMatrixField)
     print(io, eltype(field), "-valued Field")
     if eltype(eltype(field)) <: Number
+        shape = typeof(matrix_shape(field)).name.name
         if field isa Fields.FiniteDifferenceField
-            println(io, " that corresponds to the matrix")
+            println(io, " that corresponds to the $shape matrix")
         else
-            println(io, " whose first column corresponds to the matrix")
+            println(io, " whose first column corresponds to the $shape matrix")
         end
         column_field = Fields.column(field, 1, 1, 1)
         io = IOContext(io, :compact => true, :limit => true)

src/MatrixFields/field2arrays.jl

@@ -94,7 +94,9 @@ all_columns(space::Spaces.ExtrudedFiniteDifferenceSpace) =
 #     Spaces.all_nodes(Spaces.horizontal_space(axes(field)))
 
 column_map(f::F, field) where {F} =
-    map((((i, j), h),) -> f(Spaces.column(field, i, j, h)), all_columns(field))
+    Iterators.map(all_columns(field)) do ((i, j), h)
+        f(Spaces.column(field, i, j, h))
+    end
 
 """
     field2arrays(field)
@@ -104,7 +106,7 @@ Converts a field defined on a `FiniteDifferenceSpace` or on an
 corresponds to a column of the field. This is done by calling
 `column_field2array` on each of the field's columns.
 """
-field2arrays(field) = column_map(column_field2array, field)
+field2arrays(field) = collect(column_map(column_field2array, field))
 
 """
     field2arrays_view(field)