diff --git a/lib/mkl/wrappers_blas.jl b/lib/mkl/wrappers_blas.jl
index 86269a1a..067117ef 100644
--- a/lib/mkl/wrappers_blas.jl
+++ b/lib/mkl/wrappers_blas.jl
@@ -1214,6 +1214,9 @@ for (mmname_variant, smname_variant, elty) in
             lda = max(1,stride(A,2))
             ldb = max(1,stride(B,2))
             ldc = max(1,stride(C,2))
+            # the variant reads C regardless of beta (unlike gemm), and trsm leaks NaNs
+            # from an uninitialized C into the result
+            iszero(beta) && fill!(C, zero($elty))
             queue = global_queue(context(A), device(A))
             $mmname_variant(sycl_queue(queue), side, uplo, transa, diag, m, n, alpha, A, lda, B, ldb, beta, C, ldc)
             C
@@ -1236,6 +1239,9 @@ for (mmname_variant, smname_variant, elty) in
             lda = max(1,stride(A,2))
             ldb = max(1,stride(B,2))
             ldc = max(1,stride(C,2))
+            # the variant reads C regardless of beta (unlike gemm), and trsm leaks NaNs
+            # from an uninitialized C into the result
+            iszero(beta) && fill!(C, zero($elty))
             queue = global_queue(context(A), device(A))
             $smname_variant(sycl_queue(queue), side, uplo, transa, diag, m, n, alpha, A, lda, B, ldb, beta, C, ldc)
             C
diff --git a/test/onemkl.jl b/test/onemkl.jl
index 9d52cf22..c9e91022 100644
--- a/test/onemkl.jl
+++ b/test/onemkl.jl
@@ -782,6 +782,11 @@ end
             D = alpha*A*B + beta*C
             @test D ≈ h_C
             @test B == Array(dB)
+
+            # with beta = 0, the contents of C must not matter
+            dC = oneArray(fill(T(NaN), m, n))
+            oneMKL.trmm!('L','U','N','N',alpha,dA,dB,dC)
+            @test alpha*A*B ≈ Array(dC)
         end
 
         @testset "right trmm!" begin
@@ -833,6 +838,11 @@ end
             D = alpha*(A\B) + beta*C
             @test D ≈ h_C
             @test B == Array(dB)
+
+            # with beta = 0, the contents of C must not matter
+            dC = oneArray(fill(T(NaN), m, n))
+            oneMKL.trsm!('L','U','N','N',alpha,dA,dB,dC)
+            @test alpha*(A\B) ≈ Array(dC)
         end
 
         @testset "left trsm" begin
@@ -1118,6 +1128,16 @@ end
     @test all(iszero, Array(C))
 end
 
+@testset "triangular products and solves into uninitialized outputs" begin
+    S = rand(Float32, 4, 4) + 4I; B = rand(Float32, 4, 5)
+    dS, dB = oneArray(S), oneArray(B)
+    for (f, g) in ((mul!, *), (ldiv!, \))
+        C = oneArray(fill(NaN32, 4, 5))
+        f(C, UpperTriangular(dS), dB)
+        @test Array(C) ≈ g(UpperTriangular(S), B)
+    end
+end
+
 @testset "triangular products and solves with reversed views" begin
     S = rand(Float32, 4, 4) + 4I; B = rand(Float32, 4, 5); x = rand(Float32, 4)
     dS = oneArray(S)