diff --git a/src/abstractgbarray.jl b/src/abstractgbarray.jl
index cf5a4fc2..2ffa7447 100644
--- a/src/abstractgbarray.jl
+++ b/src/abstractgbarray.jl
@@ -177,7 +177,7 @@ function Base.copy(A::M) where {M<:AbstractGBArray}
     M(_copyGrBMat(A.p), A.fill)
 end
 
-function Base.Matrix(A::GBArrayOrTranspose)
+function Base.Matrix(A::AbstractGBArray)
     format = sparsitystatus(A)
     if format === Dense()
         T = unsafeunpack!(A, Dense())
@@ -195,7 +195,7 @@ function Base.Matrix(A::GBArrayOrTranspose)
     return M
 end
 
-function Base.Vector(A::GBVectorOrTranspose)
+function Base.Vector(A::AbstractGBVector)
     format = sparsitystatus(A)
     if format === Dense()
         T = unsafeunpack!(A, Dense())
@@ -213,6 +213,10 @@ function Base.Vector(A::GBVectorOrTranspose)
     return M
 end
 
+# The unpacking above needs an AbstractGBArray, so materialize the transpose first.
+Base.Matrix(A::Transpose{<:Any, <:AbstractGBArray}) = Matrix(copy(A))
+Base.Vector(A::Transpose{<:Any, <:AbstractGBVector}) = Vector(parent(A))
+
 function SparseArrays.SparseMatrixCSC(A::GBArrayOrTranspose)
     T = copy(A) # avoid changing sparsity of A and destroying it.
     return unsafeunpack!(T, SparseMatrixCSC; attachfinalizer = true)
diff --git a/src/operations/mul.jl b/src/operations/mul.jl
index cb0101e4..cd8b97fa 100644
--- a/src/operations/mul.jl
+++ b/src/operations/mul.jl
@@ -165,6 +165,44 @@ function Base.:*(
     return @_densepack B (*(A, B, op; mask, accum, desc))
 end
 
+# Dense outputs: C = α * A * B + β * C, so that mul! callers (Arpack, KrylovKit, ...) run in
+# GraphBLAS rather than the generic LinearAlgebra fallback. See #125.
+# pack needs a DenseArray: alias contiguous views (e.g. Arpack's views into workd), copy others.
+_asdense(B::DenseVecOrMat) = B
+_asdense(B::Transpose) = transpose(_asdense(parent(B)))
+function _asdense(B::StridedVecOrMat)
+    contiguous = stride(B, 1) == 1 && (B isa AbstractVector || stride(B, 2) == size(B, 1))
+    return contiguous ? unsafe_wrap(Array, pointer(B), size(B)) : copy(B)
+end
+
+function _densemul!(C, A, B, α, β)
+    size(A, 2) == size(B, 1) || throw(DimensionMismatch("size(A, 2) != size(B, 1)"))
+    size(A, 1) == size(C, 1) || throw(DimensionMismatch("size(A, 1) != size(C, 1)"))
+    size(B, 2) == size(C, 2) || throw(DimensionMismatch("size(B, 2) != size(C, 2)"))
+    T = GC.@preserve B A * pack(_asdense(B))
+    R = unsafeunpack!(T, Dense())
+    try
+        R′ = reshape(R, size(C))
+        iszero(β) ? (C .= α .* R′) : (C .= α .* R′ .+ β .* C)
+    finally
+        _jlfree(R)
+    end
+    return C
+end
+
+const _GBVecMatOrTrans = Union{GBMatrixOrTranspose, GBVectorOrTranspose}
+
+for (TC, TB) in (
+    (:StridedVector, :(StridedVector{<:valid_union})),
+    (:StridedMatrix, :(Union{StridedVecOrMat{<:valid_union}, Transpose{<:Any, <:StridedVecOrMat{<:valid_union}}})),
+)
+    @eval begin
+        LinearAlgebra.mul!(C::$TC, A::_GBVecMatOrTrans, B::$TB, α::Number, β::Number) =
+            _densemul!(C, A, B, α, β)
+        LinearAlgebra.mul!(C::$TC, A::_GBVecMatOrTrans, B::$TB) = _densemul!(C, A, B, true, false)
+    end
+end
+
 function Base.:*((⊕)::Monoid, (⊗)::Function)
     return function(A::GBArrayOrTranspose, B::GBArrayOrTranspose; mask=nothing, accum=nothing, desc=nothing)
         *(A, B, (⊕, ⊗); mask, accum, desc)
diff --git a/test/gbarray.jl b/test/gbarray.jl
index f18b4341..5bfe9f9b 100644
--- a/test/gbarray.jl
+++ b/test/gbarray.jl
@@ -38,6 +38,9 @@
             # #104: non-default fill shows up in show
             @test !occursin("fill value", sprint(show, MIME"text/plain"(), A))
             @test occursin("fill value:  17", sprint(show, MIME"text/plain"(), setfill(A, 17)))
+            A = GBMatrix(sprand(ComplexF64, 5, 4, 0.5))
+            @test Matrix(A') == Matrix(A)'
+            @test Matrix(transpose(A)) == transpose(Matrix(A))
         end
         @testset "sparse" begin
             #Construction and indexing correct for sparse matrices
diff --git a/test/operations/mul.jl b/test/operations/mul.jl
index d59e4d98..3bef20d9 100644
--- a/test/operations/mul.jl
+++ b/test/operations/mul.jl
@@ -21,4 +21,15 @@
     @test m' * n == Matrix(m)' * Vector(n)
     @test n' * m == Vector(n)' * Matrix(m)
     @test n' * n == GBVector([Vector(n)' * Vector(n)])
+
+    # mul! into dense outputs runs in GraphBLAS (#125)
+    S = sprand(ComplexF64, 20, 15, 0.3)
+    A = GBMatrix(S)
+    x, X, y, Y = rand(ComplexF64, 15), rand(ComplexF64, 15, 3), rand(ComplexF64, 20), rand(ComplexF64, 20, 3)
+    @test mul!(copy(y), A, x) ≈ S * x
+    @test mul!(copy(y), A, x, 2, 3) ≈ 2 * S * x + 3 * y
+    @test mul!(copy(Y), A, X) ≈ S * X
+    @test mul!(copy(Y), A, X, 2, 3) ≈ 2 * S * X + 3 * Y
+    @test mul!(zeros(ComplexF64, 15), A', y) ≈ S' * y
+    @test mul!(view(zeros(ComplexF64, 25), 3:22), A, view([x; x], 16:30)) ≈ S * x
 end