Fix BLAS calls on views oneMKL can't address

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

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