diff --git a/PR.md b/PR.md
new file mode 100644
index 00000000..f1b62a0c
--- /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 c65147a7..621b8ccf 100644
--- a/lib/mkl/linalg.jl
+++ b/lib/mkl/linalg.jl
@@ -7,18 +7,32 @@ 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
 #
 
+# 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,38 +53,63 @@ 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
 
 #
 # 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)
@@ -90,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
@@ -111,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
 
 
 #
@@ -150,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)
@@ -184,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)
@@ -200,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/lib/mkl/wrappers_blas.jl b/lib/mkl/wrappers_blas.jl
index afecf7ed..86269a1a 100644
--- a/lib/mkl/wrappers_blas.jl
+++ b/lib/mkl/wrappers_blas.jl
@@ -1,3 +1,17 @@
+# 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}}
+
+# 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),
@@ -484,11 +498,15 @@ 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))
+            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
@@ -503,13 +521,17 @@ 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)
-            $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
@@ -524,14 +546,18 @@ 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))
             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
@@ -539,13 +565,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,9 +584,12 @@ 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))
+            GC.@preserve x begin
+                px, incx = vec_pointer_stride(n, x)
+                $fname(sycl_queue(queue), n, alpha, px, abs(incx))
+            end
             x
         end
     end
@@ -569,11 +597,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,18 +612,21 @@ 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)
+            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
 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,22 +643,26 @@ 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)
+            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
 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 +707,13 @@ 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))
+            GC.@preserve x begin
+                px, incx = vec_pointer_stride(n, x)
+                $fname(sycl_queue(queue), n, alpha, px, abs(incx))
+            end
+            x
         end
     end
 end
@@ -782,16 +820,20 @@ 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))
+            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
 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,10 +845,13 @@ 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)
+            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
@@ -819,7 +864,10 @@ for (fname, elty) in
      (:onemklZiamax_64,:ComplexF64),
      (:onemklCiamax_64,:ComplexF32))
     @eval begin
-        function iamax(x::oneStridedArray{$elty})
+        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)
@@ -836,11 +884,12 @@ for (fname, elty) in
      (:onemklZiamin_64,:ComplexF64),
      (:onemklCiamin_64,:ComplexF32))
     @eval begin
-        function iamin(x::StridedArray{$elty})
+        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
@@ -853,11 +902,15 @@ 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))
+            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 357d9eca..fa8b4a83 100644
--- a/test/onemkl.jl
+++ b/test/onemkl.jl
@@ -139,6 +139,105 @@ 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
+
+# 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)
@@ -966,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)