From d3b5a6f0e6622fc244ad4ff6ae179e4076f603fd Mon Sep 17 00:00:00 2001 From: Tim Gymnich Date: Fri, 19 Apr 2024 13:06:53 +0200 Subject: [PATCH] Add generic_trimatmul and generic_mattrimul --- src/host/linalg.jl | 133 +++++++++++++++++++++++++++++++++++++++ test/testsuite/linalg.jl | 91 +++++++++++++++++++++++++++ 2 files changed, 224 insertions(+) diff --git a/src/host/linalg.jl b/src/host/linalg.jl index a2a99019..f0bb2436 100644 --- a/src/host/linalg.jl +++ b/src/host/linalg.jl @@ -372,6 +372,127 @@ function generic_matmatmul!(C::AbstractArray{R}, A::AbstractArray{T}, B::Abstrac C end +function generic_trimatmul!(C::AbstractGPUVecOrMat{R}, uploc, isunitc, tfun::Function, A::AbstractGPUMatrix{T}, B::AbstractGPUVecOrMat{S}) where {T,S,R} + if size(A,2) != size(B,1) + throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))")) + end + if size(C,1) != size(A,1) || size(C,2) != size(B,2) + throw(DimensionMismatch("result C has dimensions $(size(C)), needs $((size(A,1),size(B,2)))")) + end + if isempty(A) || isempty(B) + return fill!(C, zero(R)) + end + + upper = uploc == 'U' + unit = isunitc == 'U' + + function trimatmul(ctx, C, A, B) + idx = @linearidx C + assume.(size(C) .> 0) + i, j = @inbounds Tuple(CartesianIndices(C)[idx])..., 1 + l, m, n = size(A, 1), size(B, 1), size(B, 2) + + @inbounds if i <= l && j <= n + z2 = zero(A[i,1] * B[1,j] + A[i,1] * B[1,j]) + Cij = convert(promote_type(R, typeof(z2)), z2) + Cij += (unit ? one(Cij) : A[i,i]) * B[i,j] + for k in (upper ? (i + 1) : 1):(upper ? m : (i - 1)) + Cij += A[i,k] * B[k,j] + end + C[i,j] += Cij + end + + return + end + + function trimatmul_t(ctx, C, A, B) + idx = @linearidx C + assume.(size(C) .> 0) + i, j = @inbounds Tuple(CartesianIndices(C)[idx])..., 1 + l, m, n = size(A, 1), size(B, 1), size(B, 2) + + @inbounds if i <= l && j <= n + z2 = zero(A[i,1] * B[1,j] + A[i,1] * B[1,j]) + Cij = convert(promote_type(R, typeof(z2)), z2) + Cij += (unit ? one(Cij) : A[i,i]) * B[i,j] + for k in (upper ? (i + 1) : 1):(upper ? m : (i - 1)) + Ctmp += A[k,i] * B[k,j] + end + C[i,j] += Cij + end + + return + end + + function trimatmul_c(ctx, C, A, B) + idx = @linearidx C + assume.(size(C) .> 0) + i, j = @inbounds Tuple(CartesianIndices(C)[idx])..., 1 + l, m, n = size(A, 1), size(B, 1), size(B, 2) + + @inbounds if i <= l && j <= n + z2 = zero(A[i,1] * B[1,j] + A[i,1] * B[1,j]) + Cij = convert(promote_type(R, typeof(z2)), z2) + Cij += (unit ? one(Cij) : conj(A[i,i])) * B[i,j] + for k in (upper ? (i + 1) : 1):(upper ? m : (i - 1)) + Cij += conj(A[k,i]) * B[k,j] + end + C[i,j] += Cij + end + + return + end + + if tfun === identity + gpu_call(trimatmul, C, A, B; name="trimatmul") + elseif tfun == transpose + gpu_call(trimatmul_t, C, A, B; name="trimatmul") + elseif tfun === adjoint + gpu_call(trimatmul_c, C, A, B; name="trimatmul") + else + error("Not supported") + end + + C +end + +function generic_mattrimul!(C::AbstractGPUVecOrMat{R}, uploc, isunitc, tfun::Function, A::AbstractGPUMatrix{T}, B::AbstractGPUVecOrMat{S}) where {T,S,R} + if size(A,2) != size(B,1) + throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))")) + end + if size(C,1) != size(A,1) || size(C,2) != size(B,2) + throw(DimensionMismatch("result C has dimensions $(size(C)), needs $((size(A,1),size(B,2)))")) + end + if isempty(A) || isempty(B) + return fill!(C, zero(R)) + end + + upper = uploc == 'U' + unit = isunitc == 'U' + + # tA = tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C' + + gpu_call(C, A, B; name="mattrimul") do ctx, C, A, B + idx = @linearidx C + assume.(size(C) .> 0) + i, j = @inbounds Tuple(CartesianIndices(C)[idx])..., 1 + + @inbounds if i <= size(A,1) && j <= size(B,2) + z2 = zero(A[i, 1]*B[1, j] + A[i, 1]*B[1, j]) + Ctmp = convert(promote_type(R, typeof(z2)), z2) + Ctmp += A[i, j] * (unit ? one(Ctmp) : B[j, j]) + for k in (upper ? 1 : (j + 1)):(upper ? (j - 1) : (size(A,2))) + Ctmp += A[i, k]*B[k, j] + end + C[i,j] += Ctmp + end + + return + end + + C +end + function LinearAlgebra.generic_matvecmul!(C::AbstractGPUVector, tA::AbstractChar, A::AbstractGPUMatrix, B::AbstractGPUVector, _add::MulAddMul = MulAddMul()) generic_matmatmul!(C, wrap(A, tA), B, _add.alpha, _add.beta) end @@ -380,6 +501,18 @@ function LinearAlgebra.generic_matmatmul!(C::AbstractGPUVecOrMat, tA, tB, A::Abs generic_matmatmul!(C, wrap(A, tA), wrap(B, tB), _add.alpha, _add.beta) end +if VERSION >= v"1.10-" + +function LinearAlgebra.generic_trimatmul!(C::AbstractGPUVecOrMat, uploc, isunitc, tfun::Function, A::AbstractGPUMatrix, B::AbstractGPUVecOrMat) + generic_trimatmul!(C, uploc, isunitc, tfun, A, B) +end + +function LinearAlgebra.generic_mattrimul!(C::AbstractGPUMatrix, uploc, isunitc, tfun::Function, A::AbstractGPUMatrix, B::AbstractGPUMatrix) + generic_mattrimul!(C, uploc, isunitc, tfun, A, B) +end + +end + if VERSION < v"1.10.0-DEV.1365" # catch other functions that are called by LinearAlgebra's mul! function LinearAlgebra.gemv!(C::AbstractGPUVector, tA::AbstractChar, A::AbstractGPUMatrix, B::AbstractGPUVector, a::Number, b::Number) diff --git a/test/testsuite/linalg.jl b/test/testsuite/linalg.jl index da858608..24549a48 100644 --- a/test/testsuite/linalg.jl +++ b/test/testsuite/linalg.jl @@ -132,6 +132,97 @@ @test istriu(A) == istriu(B) end end + + if VERSION >= v"1.10-" + @testset "trimatmul" begin + n = 128 + b = AT(rand(Float32, n)) + B = AT(rand(Float32, n, n)) + + At = UpperTriangular(AT(rand(Float32, n,n))) + A = AT(At) + Ct = AT(zeros(Float32, n)) + C = zeros(Float32, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, b) + mul!(C, A, b) + @test Ct ≈ C + + Ct = AT(zeros(Float32, n, n)) + C = zeros(Float32, n, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, B) + mul!(C, A, B) + @test Ct ≈ C + + At = UnitUpperTriangular(AT(rand(Float32, n,n))) + A = AT(At) + Ct = AT(zeros(Float32, n)) + C = zeros(Float32, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, b) + mul!(C, A, b) + @test Ct ≈ C + + Ct = AT(zeros(Float32, n, n)) + C = zeros(Float32, n, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, B) + mul!(C, A, B) + @test Ct ≈ C + + + At = LowerTriangular(AT(rand(Float32, n,n))) + A = AT(At) + Ct = AT(zeros(Float32, n)) + C = zeros(Float32, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, b) + mul!(C, A, b) + @test Ct ≈ C + + Ct = AT(zeros(Float32, n, n)) + C = zeros(Float32, n, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, B) + mul!(C, A, B) + @test Ct ≈ C + + At = UnitLowerTriangular(AT(rand(Float32, n,n))) + A = AT(At) + Ct = AT(zeros(Float32, n)) + C = zeros(Float32, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, b) + mul!(C, A, b) + @test Ct ≈ C + + Ct = AT(zeros(Float32, n, n)) + C = zeros(Float32, n, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, B) + mul!(C, A, B) + @test Ct ≈ C + + At = UnitLowerTriangular(AT(rand(ComplexF32, n,n))) + A = AT(At) + b = AT(rand(ComplexF32, n)) + B = AT(rand(ComplexF32, n, n)) + Ct = AT(zeros(ComplexF32, n)) + C = zeros(ComplexF32, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, b) + mul!(C, A, b) + @test Ct ≈ C + + Ct = AT(zeros(ComplexF32, n, n)) + C = zeros(ComplexF32, n, n) + + LinearAlgebra.generic_trimatmul!(Ct, At, B) + mul!(C, A, B) + @test Ct ≈ C + end + end end @testset "diagonal" begin