Fixed issue in mul! for block arrays

Author: JordiManyer

NEWS.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Fixed `BlockPMatrix{V}(::UndefInitializer, rows, cols)` constructor dropping the `cols` argument, causing a `MethodError` at runtime. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194).
 - Fixed `local_views(::BlockPMatrix, rows, cols)` indexing 1D block-range vectors with a 2D `CartesianIndex`, causing `BoundsError` for any multi-field problem with ≥2 fields. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194).
+- Fixed `mul!(y::BlockPVector, A::BlockPMatrix, x::BlockPVector, α, β)` computing `α*β*(A*x)` instead of `α*(A*x) + β*y`; the 3-arg `mul!` was also updated to correctly zero `y` before accumulating. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194).
 
 ## [0.4.11] - 2026-02-20
 

src/BlockPartitionedArrays.jl

@@ -329,19 +329,23 @@ end
 
 function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector)
   o = one(eltype(A))
-  mul!(y,A,x,o,o)
+  z = zero(eltype(A))
+  mul!(y,A,x,o,z)
 end
 
 function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector,α::Number,β::Number)
   yb, Ab, xb = blocks(y), blocks(A), blocks(x)
-  z = zero(eltype(y))
   o = one(eltype(A))
+  z = zero(eltype(y))
   for i in 1:blocksize(A,1)
-    fill!(yb[i],z)
+    if iszero(β)
+      fill!(yb[i],z)
+    else
+      rmul!(yb[i],β)
+    end
     for j in 1:blocksize(A,2)
       mul!(yb[i],Ab[i,j],xb[j],α,o)
     end
-    rmul!(yb[i],β)
   end
   return y
 end

test/BlockPartitionedArraysTests.jl

@@ -66,6 +66,15 @@ function main(distribute, parts)
 
   maximum(abs, v)
   minimum(abs, v)
+
+  # BUG4 fix: mul!(y,A,x,α,β) must compute α*(A*x) + β*y, not α*β*(A*x)
+  # With A=0: y ← α*(0*x) + β*y_old = β*y_old
+  fill!(v, 1.0)
+  v4 = similar(v, block_range)
+  fill!(v4, 3.0)
+  LinearAlgebra.fillstored!(m, 0.0)
+  mul!(v4, m, v, 2.0, 2.0)   # y ← 2*(0*x) + 2*[3,...] = [6,...]
+  @test norm(v4) ≈ 6.0 * sqrt(2*n_global)
 end
 
 end # module

test/BlockSparseMatrixAssemblersTests.jl

@@ -21,7 +21,7 @@ function is_same_vector(x::BlockPVector,y::PVector,Ub,U)
   res = map(1:num_fields(Ub)) do i
     xi = restrict_to_field(Ub,x_fespace,i)
     yi = restrict_to_field(U,y_fespace,i)
-    xi ≈ yi
+    norm(xi - yi) < 1.0e-10
   end
   return all(res)
 end