Fix block setindex! on Mat: wrong array rank, and row-major values

src/mat.jl

@@ -301,14 +301,23 @@ function Base.setindex!(m::AbstractPetscMat{PetscLib}, vals::AbstractMatrix, row
     # Convert to 0-based indexing for PETSc
     petsc_rows = PetscInt[r - 1 for r in rows]
     petsc_cols = PetscInt[c - 1 for c in cols]
-    petsc_vals = PetscScalar.(vals)
-
-    # Use MatSetValues for block of entries
-    nrows = PetscInt.(length(petsc_rows))
-    ncols = PetscInt.(length(petsc_cols))
-
+
+    size(vals) == (length(rows), length(cols)) || throw(
+        DimensionMismatch(
+            "block is $(size(vals)) but the index ranges are " *
+            "$(length(rows))x$(length(cols))",
+        ),
+    )
+
+    # MatSetValues reads its value array row by row, so the block is transposed
+    # before it is flattened: Julia lays a matrix out column by column.
+    petsc_vals = vec(permutedims(PetscScalar.(vals)))
+
+    nrows = PetscInt(length(petsc_rows))
+    ncols = PetscInt(length(petsc_cols))
+
     LibPETSc.MatSetValues(PetscLib, m, nrows, petsc_rows, ncols, petsc_cols, petsc_vals, LibPETSc.INSERT_VALUES)
-    
+
     return m
 end
 

test/mat.jl

@@ -234,9 +234,50 @@ end
         @test B[1, 1] == PetscScalar(10.0)
         @test B[1, 2] == PetscScalar(11.0)
         @test B[2, 1] == PetscScalar(20.0)
-        
+
         PETSc.destroy!(B)
-
+
+        PETSc.finalize(petsclib)
+    end
+end
+
+# Writing a whole block at once (JuliaParallel/PETSc.jl#248). The block is
+# asymmetric so a transposed write shows up, which is the failure mode to watch:
+# MatSetValues reads its value array row by row and Julia stores a matrix column
+# by column.
+@testset "PetscMat block setindex!" begin
+    for petsclib in PETSc.petsclibs
+        PETSc.initialize(petsclib)
+        PetscScalar = petsclib.PetscScalar
+
+        A = PETSc.PetscMat(petsclib, 4, 4, 4)
+        block = PetscScalar[1 2 3; 4 5 6]
+        A[1:2, 1:3] = block
+        PETSc.assemble!(A)
+        @test [A[i, j] for i in 1:2, j in 1:3] == block
+
+        # Rows and columns need not be contiguous or ordered.
+        B = PETSc.PetscMat(petsclib, 5, 5, 5)
+        scattered = PetscScalar[10 20; 30 40; 50 60]
+        B[[1, 3, 5], [2, 4]] = scattered
+        PETSc.assemble!(B)
+        @test [B[i, j] for i in [1, 3, 5], j in [2, 4]] == scattered
+
+        # A block whose shape disagrees with the indices would otherwise read
+        # past the end of the array inside PETSc.
+        @test_throws DimensionMismatch B[1:2, 1:2] = PetscScalar[1 2 3; 4 5 6]
+
+        # Single row and single column go through their own methods.
+        C = PETSc.PetscMat(petsclib, 4, 4, 4)
+        C[2, [1, 3]] = PetscScalar[7, 8]
+        C[[1, 4], 3] = PetscScalar[9, 11]
+        PETSc.assemble!(C)
+        @test C[2, 1] == PetscScalar(7) && C[2, 3] == PetscScalar(8)
+        @test C[1, 3] == PetscScalar(9) && C[4, 3] == PetscScalar(11)
+
+        PETSc.destroy!(A)
+        PETSc.destroy!(B)
+        PETSc.destroy!(C)
         PETSc.finalize(petsclib)
     end
 end