From fa22bf234574a1b367ec8a0bfc1407d99c5aee85 Mon Sep 17 00:00:00 2001
From: Tim Besard <tim.besard@gmail.com>
Date: Mon, 19 Oct 2026 19:35:06 +0200
Subject: [PATCH 1/3] Fall back from oneMKL for arrays it can't address as a
 vector

The BLAS-1 wrappers take any strided array and pass it to oneMKL as
`length(x)` elements with increment `stride(x, 1)`. For a view of a
matrix that skips rows that layout is wrong, so `norm`, `rmul!`,
`axpy!`, `axpby!`, `rotate!` and `reflect!` returned wrong results and
wrote outside the view. Limit the wrappers to strided vectors and
contiguous arrays, and have the LinearAlgebra methods fall back to
broadcasts over the flattened arrays for anything else, since BLAS
pairs elements linearly regardless of shape.

The ComplexF16 variants of `scal!` and `axpy!` compute in ComplexF32
and wrote the result back through `convert`, which only works for a
dense array; use a broadcast instead. `dotc` and `dotu` likewise widen
with a broadcast now, so they accept views.
---
 PR.md                    | 20 ++++++++++++
 lib/mkl/linalg.jl        | 42 +++++++++++++++++++-----
 lib/mkl/wrappers_blas.jl | 69 +++++++++++++++++++++-------------------
 test/onemkl.jl           | 43 +++++++++++++++++++++++++
 4 files changed, 134 insertions(+), 40 deletions(-)
 create mode 100644 PR.md

diff --git a/PR.md b/PR.md
new file mode 100644
index 000000000..f1b62a0c3
--- /dev/null
+++ b/PR.md
@@ -0,0 +1,20 @@
+The BLAS-1 methods (`norm`, `dot`, `BLAS.asum`, `axpy!`, `axpby!`, `rmul!`, `rotate!`, `reflect!`) and the storage-level `mul!` and triangular methods pass any strided view to oneMKL with `incx = stride(x, 1)` and `lda = max(1, stride(A, 2))`, without checking that oneMKL can address that layout. Reversed views get a negative increment, which BLAS walks backwards from the *lowest* address while the view's pointer is at its *highest* element, so the result is silently wrong. Views that skip rows get passed as `length(x)` elements with a single increment, which reads and writes outside the view. Reversed columns give a negative `lda`, and oneMKL's `invalid_argument` escapes the support library's C interface and aborts the process:
+
+```julia
+x = oneArray(Float32[3, 4]); r = view(x, 2:-1:1)
+norm(r)                    # 4.0, expected 5.0
+rmul!(r, 2f0); x           # [3, 8], expected [6, 8]
+mul!(C, 'N', 'N', view(A, :, 4:-1:1), B, true, false)
+# terminate called after throwing an instance of 'oneapi::mkl::invalid_argument'
+```
+
+The triangular case is reachable through public API, e.g. `lmul!(UpperTriangular(A), view(x, n:-1:1))`.
+
+This ports JuliaGPU/CUDA.jl#3411 and JuliaGPU/CUDA.jl#3409. The level-1 wrappers now take strided vectors or dense arrays only, and pass a pointer to the lowest address together with the signed increment, like LinearAlgebra's `vec_pointer_stride`; nrm2, asum, scal, iamax and iamin use the absolute increment since the order doesn't matter (oneMKL rejects a negative one for iamax/iamin, so their index is flipped back). The LinearAlgebra BLAS-1 methods fall back to broadcasts for strided arrays that aren't vectors or dense. Storage-level `mul!` only calls oneMKL when every matrix has a unit row stride and a leading dimension of at least its row count, and every vector a positive stride, falling back to `GPUArrays.generic_matmatmul!` otherwise as it already did for unsupported element types. The triangular methods run on contiguous copies and copy the result back.
+
+Not covered: the level-2 wrappers (`gemv!`, `symv!`, ...) still take raw strides when called directly, as in CUDA.jl, and the support library should turn `oneapi::mkl::exception` into an error code instead of aborting, which needs a JLL rebuild.
+
+Tested on an Iris Xe (Julia 1.12): the reproducers from both issues fail on main and pass here, and the new tests cover reversed and row-skipping views for every BLAS-1 method, the storage-level `mul!` paths (gemm, gemv, symm, hemm, symv, hemv) and the triangular products and solves.
+
+Fixes #712
+Fixes #700
diff --git a/lib/mkl/linalg.jl b/lib/mkl/linalg.jl
index c65147a79..03ec30d36 100644
--- a/lib/mkl/linalg.jl
+++ b/lib/mkl/linalg.jl
@@ -11,14 +11,15 @@ using LinearAlgebra: Transpose, Adjoint,
 # BLAS 1
 #
 
+# oneMKL can only scale arrays that it can address as a vector
 LinearAlgebra.rmul!(x::oneStridedVecOrMat{T}, k::Number) where T<:Union{onemklHalf,onemklFloat} =
-    scal!(length(x), T(k), x)
+    x isa oneStridedVecOrDenseMat ? scal!(length(x), T(k), x) : (x .*= k)
 
 # Work around ambiguity with GPUArrays wrapper
 LinearAlgebra.rmul!(x::oneStridedVecOrMat{<:onemklFloat}, k::Real) =
 	invoke(LinearAlgebra.rmul!, Tuple{typeof(x), Number}, x, k)
 
-LinearAlgebra.norm(x::oneStridedVecOrMat{<:Union{Float16,ComplexF16,onemklFloat}}) = nrm2(length(x), x)
+LinearAlgebra.norm(x::oneStridedVecOrDenseMat{<:Union{Float16,ComplexF16,onemklFloat}}) = nrm2(length(x), x)
 
 function LinearAlgebra.dot(x::oneStridedVector{T}, y::oneStridedVector{T}) where T<:Union{Float16, Float32, Float64}
     n = length(x)
@@ -39,31 +40,56 @@ function LinearAlgebra.:(*)(transx::Transpose{<:Any,<:oneStridedVector{T}}, y::o
     oneMKL.dotu(n, x, y)
 end
 
-LinearAlgebra.BLAS.asum(x::oneStridedVecOrMat{<:onemklFloat}) = asum(length(x), x)
+LinearAlgebra.BLAS.asum(x::oneStridedVecOrDenseMat{<:onemklFloat}) = asum(length(x), x)
 
 function LinearAlgebra.axpy!(alpha::Number, x::oneStridedVecOrMat{T}, y::oneStridedVecOrMat{T}) where T<:Union{onemklHalf,onemklFloat}
     length(x)==length(y) || throw(DimensionMismatch("axpy arguments have lengths $(length(x)) and $(length(y))"))
-    axpy!(length(x), alpha, x, y)
+    if x isa oneStridedVecOrDenseMat && y isa oneStridedVecOrDenseMat
+        axpy!(length(x), alpha, x, y)
+    else
+        vec(y) .+= vec(x) .* alpha
+    end
+    y
 end
 
 function LinearAlgebra.axpby!(alpha::Number, x::oneStridedVecOrMat{T}, beta::Number, y::oneStridedVecOrMat{T}) where T<:onemklFloat
     length(x)==length(y) || throw(DimensionMismatch("axpby arguments have lengths $(length(x)) and $(length(y))"))
-    axpby!(length(x), alpha, x, beta, y)
+    if x isa oneStridedVecOrDenseMat && y isa oneStridedVecOrDenseMat
+        axpby!(length(x), alpha, x, beta, y)
+    else
+        vec(y) .= vec(x) .* alpha .+ vec(y) .* beta
+    end
+    y
 end
 
 function LinearAlgebra.rotate!(x::oneStridedVecOrMat{T}, y::oneStridedVecOrMat{T}, c::Number, s::Number) where T<:onemklFloat
     nx = length(x)
     ny = length(y)
     nx==ny || throw(DimensionMismatch("rotate arguments have lengths $nx and $ny"))
-    rot!(nx, x, y, c, s)
+    if x isa oneStridedVecOrDenseMat && y isa oneStridedVecOrDenseMat
+        rot!(nx, x, y, c, s)
+    else
+        xv, yv = vec(x), vec(y)
+        t = c .* xv .+ s .* yv
+        yv .= c .* yv .- conj(s) .* xv
+        xv .= t
+    end
+    x, y
 end
 
 function LinearAlgebra.reflect!(x::oneStridedVecOrMat{T}, y::oneStridedVecOrMat{T}, c::Number, s::Number) where T<:onemklFloat
     nx = length(x)
     ny = length(y)
     nx==ny || throw(DimensionMismatch("reflect arguments have lengths $nx and $ny"))
-    rot!(nx, x, y, c, s)
-    scal!(ny, -one(real(T)), y)
+    if x isa oneStridedVecOrDenseMat && y isa oneStridedVecOrDenseMat
+        rot!(nx, x, y, c, s)
+        scal!(ny, -one(real(T)), y)
+    else
+        xv, yv = vec(x), vec(y)
+        t = c .* xv .+ s .* yv
+        yv .= conj(s) .* xv .- c .* yv
+        xv .= t
+    end
     x, y
 end
 
diff --git a/lib/mkl/wrappers_blas.jl b/lib/mkl/wrappers_blas.jl
index afecf7ed8..d936ca9e9 100644
--- a/lib/mkl/wrappers_blas.jl
+++ b/lib/mkl/wrappers_blas.jl
@@ -1,3 +1,9 @@
+# oneMKL addresses a vector through a single increment, which fits any strided vector and
+# any contiguous array, but not a view of a matrix that skips rows. A SubArray with fast
+# linear indexing (the last type parameter) is contiguous, or a strided vector.
+const oneContiguousSubArray{T,N} = SubArray{T,N,<:oneArray,<:Tuple,true}
+const oneStridedVecOrDenseMat{T} = Union{oneStridedVector{T}, oneDenseArray{T}, oneContiguousSubArray{T}}
+
 ## (GE) general matrix-matrix multiplication batched
 for (fname, elty) in
         ((:onemklHgemm_batch, :Float16),
@@ -484,8 +490,8 @@ for (fname, elty) in
     @eval begin
         function axpy!(n::Integer,
                        alpha::Number,
-                       x::oneStridedArray{$elty},
-                       y::oneStridedArray{$elty})
+                       x::oneStridedVecOrDenseMat{$elty},
+                       y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             alpha = $elty(alpha)
             $fname(sycl_queue(queue), n, alpha, x, stride(x,1), y, stride(y,1))
@@ -503,9 +509,9 @@ for (fname, elty) in
     @eval begin
         function axpby!(n::Integer,
                         alpha::Number,
-                        x::oneStridedArray{$elty},
+                        x::oneStridedVecOrDenseMat{$elty},
                         beta::Number,
-                        y::oneStridedArray{$elty})
+                        y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             alpha = $elty(alpha)
             beta = $elty(beta)
@@ -524,8 +530,8 @@ for (fname, elty, cty, sty, supty) in ((:onemklSrot,:Float32,:Float32,:Float32,:
                                        (:onemklZDrot,:ComplexF64,:Float64,:Float64,:Real))
     @eval begin
         function rot!(n::Integer,
-                      x::oneStridedArray{$elty},
-                      y::oneStridedArray{$elty},
+                      x::oneStridedVecOrDenseMat{$elty},
+                      y::oneStridedVecOrDenseMat{$elty},
                       c::Real,
                       s::$supty)
             queue = global_queue(context(x), device(x))
@@ -539,13 +545,12 @@ end
 
 function axpy!(n::Integer,
             alpha::Number,
-            x::oneStridedArray{ComplexF16},
-            y::oneStridedArray{ComplexF16})
+            x::oneStridedVecOrDenseMat{ComplexF16},
+            y::oneStridedVecOrDenseMat{ComplexF16})
     wide_x = widen.(x)
     wide_y = widen.(y)
     axpy!(n, alpha, wide_x, wide_y)
-    thin_y = convert(typeof(y), wide_y)
-    copyto!(y, thin_y)
+    y .= wide_y
     return y
 end
 
@@ -559,7 +564,7 @@ for (fname, elty) in
     @eval begin
         function scal!(n::Integer,
                        alpha::$elty,
-                       x::oneStridedArray{$elty})
+                       x::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             $fname(sycl_queue(queue), n, alpha, x, stride(x,1))
             x
@@ -569,11 +574,10 @@ end
 
 function scal!(n::Integer,
             alpha::Number,
-            x::oneStridedArray{ComplexF16})
+            x::oneStridedVecOrDenseMat{ComplexF16})
     wide_x = widen.(x)
     scal!(n, convert(ComplexF32, alpha), wide_x)
-    thin_x = convert(typeof(x), wide_x)
-    copyto!(x, thin_x)
+    x .= wide_x
     return x
 end
 
@@ -585,7 +589,7 @@ for (fname, elty, ret_type) in
      (:onemklCnrm2, :ComplexF32,:Float32),
      (:onemklZnrm2, :ComplexF64,:Float64))
     @eval begin
-        function nrm2(n::Integer, x::oneStridedArray{$elty})
+        function nrm2(n::Integer, x::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             result = scalar_result($ret_type)
             $fname(sycl_queue(queue), n, x, stride(x,1), result)
@@ -594,9 +598,9 @@ for (fname, elty, ret_type) in
     end
 end
 
-nrm2(x::oneStridedArray) = nrm2(length(x), x)
+nrm2(x::oneStridedVecOrDenseMat) = nrm2(length(x), x)
 
-function nrm2(n::Integer, x::oneStridedArray{ComplexF16})
+function nrm2(n::Integer, x::oneStridedVecOrDenseMat{ComplexF16})
     wide_x = widen.(x)
     nrm = nrm2(n, wide_x)
     return convert(Float16, nrm)
@@ -613,8 +617,8 @@ for (jname, fname, elty) in
          (:dotu, :onemklZdotu, :ComplexF64))
     @eval begin
         function $jname(n::Integer,
-                         x::oneStridedArray{$elty},
-                         y::oneStridedArray{$elty})
+                         x::oneStridedVecOrDenseMat{$elty},
+                         y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             result = scalar_result($elty)
             $fname(sycl_queue(queue), n, x, stride(x,1), y, stride(y,1), result)
@@ -623,12 +627,12 @@ for (jname, fname, elty) in
     end
 end
 
-function dotc(n::Integer, x::oneStridedArray{ComplexF16}, y::oneStridedArray{ComplexF16})
-    convert(ComplexF16, dotc(n, convert(oneArray{ComplexF32}, x), convert(oneArray{ComplexF32}, y)))
+function dotc(n::Integer, x::oneStridedVecOrDenseMat{ComplexF16}, y::oneStridedVecOrDenseMat{ComplexF16})
+    convert(ComplexF16, dotc(n, widen.(x), widen.(y)))
 end
 
-function dotu(n::Integer, x::oneStridedArray{ComplexF16}, y::oneStridedArray{ComplexF16})
-    convert(ComplexF16, dotu(n, convert(oneArray{ComplexF32}, x), convert(oneArray{ComplexF32}, y)))
+function dotu(n::Integer, x::oneStridedVecOrDenseMat{ComplexF16}, y::oneStridedVecOrDenseMat{ComplexF16})
+    convert(ComplexF16, dotu(n, widen.(x), widen.(y)))
 end
 
 # level 2
@@ -673,9 +677,10 @@ for (fname, elty, celty) in ((:onemklCSscal, :Float32, :ComplexF32),
     @eval begin
         function scal!(n::Integer,
                        alpha::$elty,
-                       x::oneStridedArray{$celty})
+                       x::oneStridedVecOrDenseMat{$celty})
             queue = global_queue(context(x), device(x))
             $fname(sycl_queue(queue), n, alpha, x, stride(x,1))
+            x
         end
     end
 end
@@ -782,8 +787,8 @@ for (fname, elty) in
          (:onemklCcopy,:ComplexF32))
     @eval begin
         function copy!(n::Integer,
-                       x::oneStridedArray{$elty},
-                       y::oneStridedArray{$elty})
+                       x::oneStridedVecOrDenseMat{$elty},
+                       y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             $fname(sycl_queue(queue), n, x, stride(x, 1), y, stride(y, 1))
             y
@@ -791,7 +796,7 @@ for (fname, elty) in
     end
 end
 
-function copy!(n::Integer, x::oneStridedArray{T}, y::oneStridedArray{T}) where {T <: Union{Float16, ComplexF16}}
+function copy!(n::Integer, x::oneStridedVecOrDenseMat{T}, y::oneStridedVecOrDenseMat{T}) where {T <: Union{Float16, ComplexF16}}
     copyto!(y,x)
 end
 
@@ -803,7 +808,7 @@ for (fname, elty, ret_type) in
      (:onemklZasum, :ComplexF64, :Float64))
     @eval begin
         function asum(n::Integer,
-                      x::oneStridedArray{$elty})
+                      x::oneStridedVecOrDenseMat{$elty})
             result = scalar_result($ret_type)
             queue = global_queue(context(x), device(x))
             $fname(sycl_queue(queue), n, x, stride(x, 1), result)
@@ -819,7 +824,7 @@ for (fname, elty) in
      (:onemklZiamax_64,:ComplexF64),
      (:onemklCiamax_64,:ComplexF32))
     @eval begin
-        function iamax(x::oneStridedArray{$elty})
+        function iamax(x::oneStridedVecOrDenseMat{$elty})
             n = length(x)
             queue = global_queue(context(x), device(x))
             result = scalar_result(Int64)
@@ -836,7 +841,7 @@ for (fname, elty) in
      (:onemklZiamin_64,:ComplexF64),
      (:onemklCiamin_64,:ComplexF32))
     @eval begin
-        function iamin(x::StridedArray{$elty})
+        function iamin(x::oneStridedVecOrDenseMat{$elty})
             n = length(x)
             result = scalar_result(Int64)
             queue = global_queue(context(x), device(x))
@@ -853,8 +858,8 @@ for (fname, elty) in ((:onemklSswap,:Float32),
                       (:onemklZswap,:ComplexF64))
     @eval begin
         function swap!(n::Integer,
-            x::oneStridedArray{$elty},
-            y::oneStridedArray{$elty})
+            x::oneStridedVecOrDenseMat{$elty},
+            y::oneStridedVecOrDenseMat{$elty})
             # Assuming both memory allocated on same device & context
             queue = global_queue(context(x), device(x))
             $fname(sycl_queue(queue), n, x, stride(x, 1), y, stride(y, 1))
diff --git a/test/onemkl.jl b/test/onemkl.jl
index 357d9eca1..b1ce36927 100644
--- a/test/onemkl.jl
+++ b/test/onemkl.jl
@@ -139,6 +139,49 @@ end
     end
 end
 
+# BLAS-1 methods on a view of a matrix that skips rows cannot pass it to oneMKL as a
+# vector, and fall back to generic code (#712)
+@testset "views oneMKL can't address as a vector" begin
+    A, B = rand(Float32, 3, 4), rand(Float32, 3, 4)
+    dA, dB = oneArray(A), oneArray(B)
+    rows(X) = view(X, 1:2, :)
+    @test norm(rows(dA)) ≈ norm(rows(A))
+    for f! in ((x, y) -> rmul!(x, 2f0),
+               (x, y) -> axpy!(2f0, y, x),
+               (x, y) -> axpby!(2f0, y, 3f0, x),
+               (x, y) -> rotate!(x, y, 0.5f0, 0.25f0),
+               (x, y) -> reflect!(x, y, 0.5f0, 0.25f0))
+        a, b, da, db = copy(A), copy(B), copy(dA), copy(dB)
+        f!(vec(rows(a)), vec(rows(b)))
+        f!(rows(da), rows(db))
+        @test Array(da) ≈ a
+        @test Array(db) ≈ b
+    end
+
+    # BLAS pairs elements linearly, so the shapes needn't match: a contiguous view of a
+    # matrix with a vector, and a row-skipping view with a vector
+    C = rand(Float32, 2, 4); v = rand(Float32, 8)
+    dC, dv = oneArray(C), oneArray(v)
+    @test Array(axpy!(2f0, view(dC, :, :), copy(dv))) ≈ axpy!(2f0, vec(C), copy(v))
+    @test Array(axpy!(2f0, view(dA, 1:2, :), copy(dv))) ≈ axpy!(2f0, vec(A[1:2, :]), copy(v))
+    @test Array(axpby!(2f0, view(dA, 1:2, :), 3f0, copy(dv))) ≈ axpby!(2f0, vec(A[1:2, :]), 3f0, copy(v))
+    dst = view(copy(dA), 1:2, :)
+    @test axpy!(2f0, view(dC, 1:2, 1:4), dst) === dst
+    @test axpby!(2f0, view(dC, 1:2, 1:4), 3f0, dst) === dst
+
+    # ComplexF16 is computed in ComplexF32 and written back
+    if ComplexF16 in eltypes
+        x, y = rand(ComplexF16, m), rand(ComplexF16, m)
+        dx, dy = oneArray(x), oneArray(y)
+        oneMKL.axpy!(m ÷ 2, ComplexF16(2), view(dx, 1:2:m), view(dy, 1:2:m))
+        axpy!(ComplexF16(2), view(x, 1:2:m), view(y, 1:2:m))
+        @test Array(dy) ≈ y
+        oneMKL.scal!(m ÷ 2, 3f0, view(dx, 1:2:m))
+        rmul!(view(x, 1:2:m), 3f0)
+        @test Array(dx) ≈ x
+    end
+end
+
 @testset "level 2" begin
     @testset for T in intersect(eltypes, [Float32, Float64, ComplexF32, ComplexF64])
         alpha = rand(T)

From 8b42e50ba9dda1dcc601351345700cf2b78bba22 Mon Sep 17 00:00:00 2001
From: Tim Besard <tim.besard@gmail.com>
Date: Mon, 19 Oct 2026 19:36:31 +0200
Subject: [PATCH 2/3] Fix BLAS-1 calls on reversed views

oneMKL follows BLAS and walks a vector with a negative increment
backwards from its lowest address, while the wrappers passed the
pointer of a reversed view's first element, the highest one. As a result
`norm`, `dot`, `BLAS.asum`, `axpy!`, `axpby!`, `rmul!`, `rotate!` and
`reflect!` on `view(x, n:-1:1)` returned wrong results or accessed the
wrong memory. Pass a pointer to the lowest address together with the
increment, like LinearAlgebra's `vec_pointer_stride`, and use the
absolute increment for nrm2, asum and scal, where the order doesn't
matter. iamax and iamin report the first of several extrema, so they
run on a contiguous copy of a reversed vector instead. This also fixes
direct calls like `oneMKL.nrm2(r)`.

Fixes #712
---
 lib/mkl/wrappers_blas.jl | 70 +++++++++++++++++++++++++++++++++-------
 test/onemkl.jl           | 56 ++++++++++++++++++++++++++++++++
 2 files changed, 115 insertions(+), 11 deletions(-)

diff --git a/lib/mkl/wrappers_blas.jl b/lib/mkl/wrappers_blas.jl
index d936ca9e9..86269a1aa 100644
--- a/lib/mkl/wrappers_blas.jl
+++ b/lib/mkl/wrappers_blas.jl
@@ -4,6 +4,14 @@
 const oneContiguousSubArray{T,N} = SubArray{T,N,<:oneArray,<:Tuple,true}
 const oneStridedVecOrDenseMat{T} = Union{oneStridedVector{T}, oneDenseArray{T}, oneContiguousSubArray{T}}
 
+# BLAS walks a vector with a negative increment backwards from its lowest address, while a
+# reversed view points at its first element, the highest one. Like LinearAlgebra.BLAS, pass
+# a pointer to the lowest address of the `n` elements accessed, along with the increment.
+function vec_pointer_stride(n::Integer, x::oneStridedVecOrDenseMat)
+    st = stride(x, 1)
+    pointer(x) + min(st, 0) * sizeof(eltype(x)) * (n - 1), st
+end
+
 ## (GE) general matrix-matrix multiplication batched
 for (fname, elty) in
         ((:onemklHgemm_batch, :Float16),
@@ -494,7 +502,11 @@ for (fname, elty) in
                        y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             alpha = $elty(alpha)
-            $fname(sycl_queue(queue), n, alpha, x, stride(x,1), y, stride(y,1))
+            GC.@preserve x y begin
+                px, incx = vec_pointer_stride(n, x)
+                py, incy = vec_pointer_stride(n, y)
+                $fname(sycl_queue(queue), n, alpha, px, incx, py, incy)
+            end
             y
         end
     end
@@ -515,7 +527,11 @@ for (fname, elty) in
             queue = global_queue(context(x), device(x))
             alpha = $elty(alpha)
             beta = $elty(beta)
-            $fname(sycl_queue(queue), n, alpha, x, stride(x,1), beta, y, stride(y,1))
+            GC.@preserve x y begin
+                px, incx = vec_pointer_stride(n, x)
+                py, incy = vec_pointer_stride(n, y)
+                $fname(sycl_queue(queue), n, alpha, px, incx, beta, py, incy)
+            end
             y
         end
     end
@@ -537,7 +553,11 @@ for (fname, elty, cty, sty, supty) in ((:onemklSrot,:Float32,:Float32,:Float32,:
             queue = global_queue(context(x), device(x))
             c = $cty(c)
             s = $sty(s)
-            $fname(sycl_queue(queue), n, x, stride(x, 1), y, stride(y, 1), c, s)
+            GC.@preserve x y begin
+                px, incx = vec_pointer_stride(n, x)
+                py, incy = vec_pointer_stride(n, y)
+                $fname(sycl_queue(queue), n, px, incx, py, incy, c, s)
+            end
             x, y
         end
     end
@@ -566,7 +586,10 @@ for (fname, elty) in
                        alpha::$elty,
                        x::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
-            $fname(sycl_queue(queue), n, alpha, x, stride(x,1))
+            GC.@preserve x begin
+                px, incx = vec_pointer_stride(n, x)
+                $fname(sycl_queue(queue), n, alpha, px, abs(incx))
+            end
             x
         end
     end
@@ -592,7 +615,10 @@ for (fname, elty, ret_type) in
         function nrm2(n::Integer, x::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             result = scalar_result($ret_type)
-            $fname(sycl_queue(queue), n, x, stride(x,1), result)
+            GC.@preserve x begin
+                px, incx = vec_pointer_stride(n, x)
+                $fname(sycl_queue(queue), n, px, abs(incx), result)
+            end
             return result[1]
         end
     end
@@ -621,7 +647,11 @@ for (jname, fname, elty) in
                          y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
             result = scalar_result($elty)
-            $fname(sycl_queue(queue), n, x, stride(x,1), y, stride(y,1), result)
+            GC.@preserve x y begin
+                px, incx = vec_pointer_stride(n, x)
+                py, incy = vec_pointer_stride(n, y)
+                $fname(sycl_queue(queue), n, px, incx, py, incy, result)
+            end
             return result[1]
         end
     end
@@ -679,7 +709,10 @@ for (fname, elty, celty) in ((:onemklCSscal, :Float32, :ComplexF32),
                        alpha::$elty,
                        x::oneStridedVecOrDenseMat{$celty})
             queue = global_queue(context(x), device(x))
-            $fname(sycl_queue(queue), n, alpha, x, stride(x,1))
+            GC.@preserve x begin
+                px, incx = vec_pointer_stride(n, x)
+                $fname(sycl_queue(queue), n, alpha, px, abs(incx))
+            end
             x
         end
     end
@@ -790,7 +823,11 @@ for (fname, elty) in
                        x::oneStridedVecOrDenseMat{$elty},
                        y::oneStridedVecOrDenseMat{$elty})
             queue = global_queue(context(x), device(x))
-            $fname(sycl_queue(queue), n, x, stride(x, 1), y, stride(y, 1))
+            GC.@preserve x y begin
+                px, incx = vec_pointer_stride(n, x)
+                py, incy = vec_pointer_stride(n, y)
+                $fname(sycl_queue(queue), n, px, incx, py, incy)
+            end
             y
         end
     end
@@ -811,7 +848,10 @@ for (fname, elty, ret_type) in
                       x::oneStridedVecOrDenseMat{$elty})
             result = scalar_result($ret_type)
             queue = global_queue(context(x), device(x))
-            $fname(sycl_queue(queue), n, x, stride(x, 1), result)
+            GC.@preserve x begin
+                px, incx = vec_pointer_stride(n, x)
+                $fname(sycl_queue(queue), n, px, abs(incx), result)
+            end
             return result[1]
         end
     end
@@ -825,6 +865,9 @@ for (fname, elty) in
      (:onemklCiamax_64,:ComplexF32))
     @eval begin
         function iamax(x::oneStridedVecOrDenseMat{$elty})
+            # oneMKL rejects a negative increment, and walking the elements backwards
+            # would change which of several extrema is reported first
+            stride(x, 1) < 0 && return iamax(copy(x))
             n = length(x)
             queue = global_queue(context(x), device(x))
             result = scalar_result(Int64)
@@ -842,10 +885,11 @@ for (fname, elty) in
      (:onemklCiamin_64,:ComplexF32))
     @eval begin
         function iamin(x::oneStridedVecOrDenseMat{$elty})
+            stride(x, 1) < 0 && return iamin(copy(x))
             n = length(x)
             result = scalar_result(Int64)
             queue = global_queue(context(x), device(x))
-            $fname(sycl_queue(queue),n, x, stride(x, 1), result, 'O')
+            $fname(sycl_queue(queue), n, x, stride(x, 1), result, 'O')
             return result[1]
         end
     end
@@ -862,7 +906,11 @@ for (fname, elty) in ((:onemklSswap,:Float32),
             y::oneStridedVecOrDenseMat{$elty})
             # Assuming both memory allocated on same device & context
             queue = global_queue(context(x), device(x))
-            $fname(sycl_queue(queue), n, x, stride(x, 1), y, stride(y, 1))
+            GC.@preserve x y begin
+                px, incx = vec_pointer_stride(n, x)
+                py, incy = vec_pointer_stride(n, y)
+                $fname(sycl_queue(queue), n, px, incx, py, incy)
+            end
             x, y
         end
     end
diff --git a/test/onemkl.jl b/test/onemkl.jl
index b1ce36927..d882cc732 100644
--- a/test/onemkl.jl
+++ b/test/onemkl.jl
@@ -182,6 +182,62 @@ end
     end
 end
 
+# BLAS walks a negative increment from the lowest address (#712)
+@testset "negative strides $T" for T in intersect(eltypes, [Float16, ComplexF16, Float32, Float64, ComplexF32, ComplexF64])
+    rev(a) = view(a, lastindex(a):-2:1)
+    fwd(a) = view(a, 1:2:lastindex(a))
+    x, y = rand(T, m), rand(T, m)
+    dx, dy = oneArray(x), oneArray(y)
+
+    @test norm(rev(dx)) ≈ norm(rev(x))
+    @test dot(rev(dx), fwd(dy)) ≈ dot(rev(x), fwd(y))
+    @test dot(rev(dx), rev(dy)) ≈ dot(rev(x), rev(y))
+    if T <: oneMKL.onemklFloat
+        @test BLAS.asum(rev(dx)) ≈ BLAS.asum(rev(x))
+        # BLAS compares |re| + |im|
+        @test oneMKL.iamax(rev(dx)) == argmax(abs.(real.(rev(x))) .+ abs.(imag.(rev(x))))
+        @test oneMKL.iamin(rev(dx)) == argmin(abs.(real.(rev(x))) .+ abs.(imag.(rev(x))))
+        # the first of several extrema, in the view's order, and 0 for an empty vector
+        ties = oneArray(T[9, 1, 9, 1])
+        @test oneMKL.iamax(view(ties, 4:-1:1)) == 2
+        @test oneMKL.iamin(view(ties, 4:-1:1)) == 1
+        @test oneMKL.iamax(view(ties, 1:-1:2)) == 0
+        @test oneMKL.iamin(view(ties, 1:-1:2)) == 0
+    end
+
+    # the LinearAlgebra methods for axpy! and rmul! don't cover ComplexF16
+    fs = Any[]
+    T <: Union{Float16, oneMKL.onemklFloat} && push!(fs, (x, y) -> axpy!(T(2), rev(x), fwd(y)),
+                                                        (x, y) -> axpy!(T(2), fwd(x), rev(y)),
+                                                        (x, y) -> rmul!(rev(y), T(3)))
+    if T <: oneMKL.onemklFloat
+        c, s = rand(real(T)), rand(T)
+        push!(fs, (x, y) -> axpby!(T(2), rev(x), T(3), rev(y)),
+                  (x, y) -> rotate!(rev(y), fwd(x), c, s),
+                  (x, y) -> reflect!(rev(y), fwd(x), c, s))
+    end
+    for f! in fs
+        a, b, da, db = copy(x), copy(y), copy(dx), copy(dy)
+        f!(a, b)
+        f!(da, db)
+        @test Array(da) ≈ a
+        @test Array(db) ≈ b
+    end
+
+    if T <: oneMKL.onemklFloat
+        a, da = copy(y), copy(dy)
+        copyto!(fwd(a), rev(x))
+        oneMKL.copy!(m ÷ 2, rev(dx), fwd(da))
+        @test Array(da) == a
+
+        a, b, da, db = copy(x), copy(y), copy(dx), copy(dy)
+        a[1:2:end], b[end:-2:1] = b[end:-2:1], a[1:2:end]
+        oneMKL.swap!(m ÷ 2, fwd(da), rev(db))
+        @test Array(da) == a
+        @test Array(db) == b
+    end
+end
+
 @testset "level 2" begin
     @testset for T in intersect(eltypes, [Float32, Float64, ComplexF32, ComplexF64])
         alpha = rand(T)

From 804dd6bb43cf8f178da06501a4cbf70c7fdaff8c Mon Sep 17 00:00:00 2001
From: Tim Besard <tim.besard@gmail.com>
Date: Mon, 19 Oct 2026 19:37:49 +0200
Subject: [PATCH 3/3] Check operand layouts before calling oneMKL in mul! and
 triangular methods

The storage-level `mul!` methods passed any strided `SubArray` to
oneMKL with `lda = max(1, stride(A, 2))` and `incx = stride(x, 1)`,
and the triangular `lmul!`/`ldiv!`/`rmul!`/`rdiv!` methods checked
nothing. Views that skip rows therefore gave silently wrong results,
reversed vectors were walked from the wrong end, and reversed columns
made oneMKL throw `invalid_argument` for `lda`, which escapes the
support library's C interface and aborts the process. The triangular
methods were reachable through public API, e.g.
`lmul!(UpperTriangular(A), view(x, n:-1:1))`.

Matrices, and vectors passed as matrices, now need a leading dimension
of at least the number of rows, and vectors a positive stride. `mul!`
falls back to `GPUArrays.generic_matmatmul!` otherwise, as it already
did for other unsupported layouts. The matrix-vector `mul!` now also
accepts a strided output vector, which the layout check keeps away
from oneMKL when needed. The triangular methods run on contiguous
copies and copy the result back.

Fixes #700
---
 lib/mkl/linalg.jl | 57 ++++++++++++++++++++--------
 test/onemkl.jl    | 95 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 137 insertions(+), 15 deletions(-)

diff --git a/lib/mkl/linalg.jl b/lib/mkl/linalg.jl
index 03ec30d36..621b8ccf4 100644
--- a/lib/mkl/linalg.jl
+++ b/lib/mkl/linalg.jl
@@ -7,6 +7,19 @@ using LinearAlgebra: Transpose, Adjoint,
                      UpperTriangular, UnitUpperTriangular,
                      MulAddMul, wrap
 
+# oneMKL addresses a matrix, or a vector passed as one, through its leading dimension
+# `stride(X, 2)`, which must be at least the number of rows. Reversed views fail this.
+column_major(X::oneStridedVecOrMat) = stride(X, 1) == 1 && stride(X, 2) >= max(1, size(X, 1))
+
+# BLAS walks a vector with a negative increment backwards from its end, while a reversed
+# view points at its first element, so vector operands need a positive stride.
+blas_layout(x::oneStridedVector) = stride(x, 1) > 0
+blas_layout(A::oneStridedMatrix) = column_major(A)
+
+# run a oneMKL routine on contiguous copies of operands it can't address
+blas_input(X) = blas_layout(X) ? X : copy(X)
+blas_output(f, C) = blas_layout(C) ? f(C) : copyto!(C, f(similar(C)))
+
 #
 # BLAS 1
 #
@@ -96,7 +109,7 @@ end
 #
 # BLAS 2
 
-function LinearAlgebra.mul!(Y::oneVector, tA::AbstractChar, A::oneStridedMatrix, B::oneStridedVector, a::Number, b::Number)
+function LinearAlgebra.mul!(Y::oneStridedVector, tA::AbstractChar, A::oneStridedMatrix, B::oneStridedVector, a::Number, b::Number)
     mA, nA = tA == 'N' ? size(A) : reverse(size(A))
 
     if nA != length(B)
@@ -116,7 +129,8 @@ function LinearAlgebra.mul!(Y::oneVector, tA::AbstractChar, A::oneStridedMatrix,
 
     alpha, beta = promote(a, b, zero(T))
     if alpha isa Union{Bool,T} && beta isa Union{Bool,T}
-        if T <: onemklFloat && eltype(A) == eltype(B) == T
+        if T <: onemklFloat && eltype(A) == eltype(B) == T &&
+           blas_layout(A) && blas_layout(B) && blas_layout(Y)
             if tA in ('N', 'T', 'C')
                 return gemv!(tA, alpha, A, B, beta, Y)
             elseif tA in ('S', 's') && T <: Real
@@ -137,10 +151,14 @@ end
 # triangular
 ## multiplication
 LinearAlgebra.generic_trimatmul!(c::oneStridedVector{T}, uploc, isunitc, tfun::Function, A::oneStridedMatrix{T}, b::oneStridedVector{T}) where {T<:onemklFloat} =
-    trmv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, c === b ? c : copyto!(c, b))
+    blas_output(c) do c
+        trmv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, blas_input(A), c === b ? c : copyto!(c, b))
+    end
 ## division
 LinearAlgebra.generic_trimatdiv!(C::oneStridedVector{T}, uploc, isunitc, tfun::Function, A::oneStridedMatrix{T}, B::oneStridedVector{T}) where {T<:onemklFloat} =
-    trsv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, C === B ? C : copyto!(C, B))
+    blas_output(C) do C
+        trsv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, blas_input(A), C === B ? C : copyto!(C, B))
+    end
 
 
 #
@@ -176,7 +194,8 @@ function LinearAlgebra.mul!(
         return iszero(beta) ? fill!(C, zero(T)) : rmul!(C, beta)
     end
 
-    if T <: Union{onemklFloat, onemklComplex, onemklHalf} &&
+    mkl_layout = column_major(A) && column_major(B) && column_major(C)
+    if mkl_layout && T <: Union{onemklFloat, onemklComplex, onemklHalf} &&
             alpha isa Union{Bool,T} && beta isa Union{Bool,T}
         # TODO: should the gemm part above be included in this branch?
         α, β = T(alpha), T(beta)
@@ -210,9 +229,9 @@ end
 # the old names, both the alpha/beta variants (1.12) and the ones taking a final MulAddMul
 # (1.10 and 1.11).
 @static if VERSION < v"1.13.0-rc4"
-    LinearAlgebra.generic_matvecmul!(Y::oneVector, tA::AbstractChar, A::oneStridedMatrix, B::oneStridedVector, alpha::Number, beta::Number) =
+    LinearAlgebra.generic_matvecmul!(Y::oneStridedVector, tA::AbstractChar, A::oneStridedMatrix, B::oneStridedVector, alpha::Number, beta::Number) =
         LinearAlgebra.mul!(Y, tA, A, B, alpha, beta)
-    LinearAlgebra.generic_matvecmul!(Y::oneVector, tA::AbstractChar, A::oneStridedMatrix, B::oneStridedVector, _add::MulAddMul) =
+    LinearAlgebra.generic_matvecmul!(Y::oneStridedVector, tA::AbstractChar, A::oneStridedMatrix, B::oneStridedVector, _add::MulAddMul) =
         LinearAlgebra.mul!(Y, tA, A, B, _add.alpha, _add.beta)
     LinearAlgebra.generic_matmatmul!(C::oneStridedVecOrMat, tA, tB, A::oneStridedVecOrMat, B::oneStridedVecOrMat, alpha::Number, beta::Number) =
         LinearAlgebra.mul!(C, tA, tB, A, B, alpha, beta)
@@ -226,17 +245,25 @@ end
 # operand; that case stays on the in-place oneMKL routines. Otherwise the out-of-place
 # "variants" write straight into C without first copying the operand over.
 LinearAlgebra.generic_trimatmul!(C::oneStridedMatrix{T}, uploc, isunitc, tfun::Function, A::oneStridedMatrix{T}, B::oneStridedMatrix{T}) where {T<:onemklFloat} =
-    C === B ? trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C) :
-              trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, B, C)
+    blas_output(C) do C
+        C === B ? trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(A), C) :
+                  trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(A), blas_input(B), C)
+    end
 LinearAlgebra.generic_mattrimul!(C::oneStridedMatrix{T}, uploc, isunitc, tfun::Function, A::oneStridedMatrix{T}, B::oneStridedMatrix{T}) where {T<:onemklFloat} =
-    C === A ? trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C) :
-              trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, A, C)
+    blas_output(C) do C
+        C === A ? trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(B), C) :
+                  trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(B), blas_input(A), C)
+    end
 LinearAlgebra.generic_trimatdiv!(C::oneStridedMatrix{T}, uploc, isunitc, tfun::Function, A::oneStridedMatrix{T}, B::oneStridedMatrix{T}) where {T<:onemklFloat} =
-    C === B ? trsm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C) :
-              trsm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, B, C)
+    blas_output(C) do C
+        C === B ? trsm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(A), C) :
+                  trsm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(A), blas_input(B), C)
+    end
 LinearAlgebra.generic_mattridiv!(C::oneStridedMatrix{T}, uploc, isunitc, tfun::Function, A::oneStridedMatrix{T}, B::oneStridedMatrix{T}) where {T<:onemklFloat} =
-    C === A ? trsm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C) :
-              trsm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, A, C)
+    blas_output(C) do C
+        C === A ? trsm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(B), C) :
+                  trsm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), blas_input(B), blas_input(A), C)
+    end
 
 #
 # BLAS extensions
diff --git a/test/onemkl.jl b/test/onemkl.jl
index d882cc732..fa8b4a83c 100644
--- a/test/onemkl.jl
+++ b/test/onemkl.jl
@@ -1065,6 +1065,101 @@ end
     end
 end
 
+# storage-level mul! only calls oneMKL for operands it can address (#700)
+@testset "storage-level mul! falls back where oneMKL can't go" begin
+    A = rand(Float32, 8, 6); B = rand(Float32, 6, 5); x = rand(Float32, 6)
+    dA, dB, dx = oneArray(A), oneArray(B), oneArray(x)
+
+    # columns that aren't contiguous (row step 2)
+    C = oneArray(zeros(Float32, 4, 5))
+    mul!(C, 'N', 'N', view(dA, 1:2:8, :), dB, true, false)
+    @test Array(C) ≈ A[1:2:8, :] * B
+    y = oneArray(zeros(Float32, 4))
+    mul!(y, 'N', view(dA, 1:2:8, :), dx, true, false)
+    @test Array(y) ≈ A[1:2:8, :] * x
+
+    # reversed views, which oneMKL would get as a negative leading dimension or increment
+    C = oneArray(zeros(Float32, 8, 5))
+    mul!(C, 'N', 'N', view(dA, :, 6:-1:1), dB, true, false)
+    @test Array(C) ≈ A[:, 6:-1:1] * B
+    y = oneArray(zeros(Float32, 8))
+    mul!(y, 'N', view(dA, :, 6:-1:1), dx, true, false)
+    @test Array(y) ≈ A[:, 6:-1:1] * x
+    mul!(y, 'N', dA, view(dx, 6:-1:1), true, false)
+    @test Array(y) ≈ A * reverse(x)
+    mul!(view(y, 8:-1:1), 'N', dA, dx, true, false)
+    @test Array(y) ≈ reverse(A * x)
+    S = rand(ComplexF32, 6, 6); dS = oneArray(S)
+    dBc, dxc = oneArray(ComplexF32.(B)), oneArray(ComplexF32.(x))
+    for (tA, wrapper) in (('S', Symmetric), ('H', Hermitian))
+        Cc = oneArray(zeros(ComplexF32, 6, 5))
+        mul!(Cc, tA, 'N', view(dS, :, 6:-1:1), dBc, true, false)
+        @test Array(Cc) ≈ wrapper(S[:, 6:-1:1]) * B
+        yc = oneArray(zeros(ComplexF32, 6))
+        mul!(yc, tA, dS, view(dxc, 6:-1:1), true, false)
+        @test Array(yc) ≈ wrapper(S) * reverse(x)
+    end
+
+    # element types oneMKL has no symm/hemm for
+    S16 = rand(Float16, 4, 4); B16 = rand(Float16, 2, 3)
+    C16 = oneArray(zeros(Float16, 2, 3))
+    mul!(C16, 'S', 'N', view(oneArray(S16), 1:2, 1:2), oneArray(B16), true, false)
+    @test Array(C16) ≈ Symmetric(S16[1:2, 1:2]) * B16
+    if Float64 in eltypes
+        S32 = rand(Float32, 4, 4); B64 = rand(Float64, 4, 3)
+        C64 = oneArray(zeros(Float64, 4, 3))
+        mul!(C64, 'S', 'N', oneArray(S32), oneArray(B64), true, false)
+        @test Array(C64) ≈ Symmetric(S32) * B64
+    end
+
+    # an empty inner dimension only scales C, which needn't be addressable by oneMKL
+    C = oneArray(fill(NaN32, 3, 2))
+    mul!(view(C, :, 2:-1:1), 'N', 'N', oneArray(zeros(Float32, 3, 0)), oneArray(zeros(Float32, 0, 2)), true, false)
+    @test all(iszero, Array(C))
+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)
+    U, dU = UpperTriangular(S[:, 4:-1:1]), UpperTriangular(view(dS, :, 4:-1:1))
+
+    # reversed vectors
+    z = oneArray(x)
+    lmul!(UpperTriangular(dS), view(z, 4:-1:1))
+    @test Array(z) ≈ reverse(UpperTriangular(S) * reverse(x))
+    z = oneArray(x)
+    ldiv!(UpperTriangular(dS), view(z, 4:-1:1))
+    @test Array(z) ≈ reverse(UpperTriangular(S) \ reverse(x))
+
+    # reversed triangular matrix
+    z = oneArray(x)
+    @test Array(lmul!(dU, z)) ≈ U * x
+    z = oneArray(x)
+    @test Array(ldiv!(dU, z)) ≈ U \ x
+    Z = oneArray(B)
+    @test Array(lmul!(dU, Z)) ≈ U * B
+    Z = oneArray(B)
+    @test Array(ldiv!(dU, Z)) ≈ U \ B
+    Z = oneArray(copy(B'))
+    @test Array(rmul!(Z, dU)) ≈ B' * U
+    Z = oneArray(copy(B'))
+    @test Array(rdiv!(Z, dU)) ≈ B' / U
+
+    # reversed right-hand side
+    Z = oneArray(B)
+    lmul!(UpperTriangular(dS), view(Z, :, 5:-1:1))
+    @test Array(Z) ≈ (UpperTriangular(S) * B[:, 5:-1:1])[:, 5:-1:1]
+    Z = oneArray(B)
+    ldiv!(UpperTriangular(dS), view(Z, :, 5:-1:1))
+    @test Array(Z) ≈ (UpperTriangular(S) \ B[:, 5:-1:1])[:, 5:-1:1]
+    Z = oneArray(copy(B'))
+    rmul!(view(Z, 5:-1:1, :), UpperTriangular(dS))
+    @test Array(Z) ≈ (B'[5:-1:1, :] * UpperTriangular(S))[5:-1:1, :]
+    Z = oneArray(copy(B'))
+    rdiv!(view(Z, 5:-1:1, :), UpperTriangular(dS))
+    @test Array(Z) ≈ (B'[5:-1:1, :] / UpperTriangular(S))[5:-1:1, :]
+end
+
 @testset "Batch Primitives" begin
     @testset for T in intersect(eltypes, [Float16, Float32, Float64, ComplexF32, ComplexF64])
         alpha = rand(T)