Skip to content

Commit 2b1bb15

Browse files
michel2323michel2323
committed
Support CSC
1 parent 919a2b7 commit 2b1bb15

3 files changed

Lines changed: 424 additions & 80 deletions

File tree

‎lib/mkl/wrappers_sparse.jl‎

Lines changed: 314 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -121,6 +121,78 @@ for SparseMatrix in (:oneSparseMatrixCSR, :oneSparseMatrixCOO)
121121
end
122122
end
123123

124+
# Special handling for CSC matrices since they are stored as transposed CSR
125+
for (fname, elty) in ((:onemklSsparse_gemv, :Float32),
126+
(:onemklDsparse_gemv, :Float64),
127+
(:onemklCsparse_gemv, :ComplexF32),
128+
(:onemklZsparse_gemv, :ComplexF64))
129+
@eval begin
130+
function sparse_gemv!(trans::Char,
131+
alpha::Number,
132+
A::oneSparseMatrixCSC{$elty},
133+
x::oneStridedVector{$elty},
134+
beta::Number,
135+
y::oneStridedVector{$elty})
136+
137+
# CSC(A) is represented by storing CSR(A^T). Map operations accordingly:
138+
# - trans = 'N': want A*x -> use op(S)='T' with S=A^T.
139+
# - trans = 'T': want A^T*x -> use op(S)='N' with S=A^T.
140+
# - trans = 'C': want A^H*x.
141+
# * For real eltypes, A^H == A^T -> use op(S)='N'.
142+
# * For complex eltypes, we cannot express A^H using a single op(S).
143+
# Use identity: conj(y_new) = conj(alpha) * A * conj(x) + conj(beta) * conj(y)
144+
# and compute with op(S)='T' (since S^T = A), conjugating x and y around the call.
145+
146+
if trans == 'N'
147+
queue = global_queue(context(x), device())
148+
$fname(sycl_queue(queue), 'T', alpha, A.handle, x, beta, y)
149+
return y
150+
elseif trans == 'T'
151+
queue = global_queue(context(x), device())
152+
$fname(sycl_queue(queue), 'N', alpha, A.handle, x, beta, y)
153+
return y
154+
else
155+
# trans == 'C'
156+
if $elty <: Complex
157+
# Compute A^H*x via identity:
158+
# conj(y_new) = conj(alpha) * (A^T) * conj(x) + conj(beta) * conj(y)
159+
# Since S=A^T and op='N' computes S*x = A^T*x, we can realize this with one call.
160+
y .= conj.(y)
161+
x_conj = similar(x)
162+
x_conj .= conj.(x)
163+
164+
queue = global_queue(context(x), device())
165+
$fname(sycl_queue(queue), 'N', conj(alpha), A.handle, x_conj, conj(beta), y)
166+
167+
y .= conj.(y)
168+
return y
169+
else
170+
# real eltype: A^H == A^T
171+
queue = global_queue(context(x), device())
172+
$fname(sycl_queue(queue), 'N', alpha, A.handle, x, beta, y)
173+
return y
174+
end
175+
end
176+
end
177+
end
178+
end
179+
180+
function sparse_optimize_gemv!(trans::Char, A::oneSparseMatrixCSC)
181+
# For CSC matrices stored as transposed CSR, we need to optimize with the transposed operation
182+
if trans == 'N'
183+
actual_trans = 'T'
184+
elseif trans == 'T'
185+
actual_trans = 'N'
186+
else # trans == 'C'
187+
# complex 'C' case is implemented using op='N' on S=A^T with conjugation trick
188+
actual_trans = 'N'
189+
end
190+
191+
queue = global_queue(context(A.nzVal), device(A.nzVal))
192+
onemklXsparse_optimize_gemv(sycl_queue(queue), actual_trans, A.handle)
193+
return A
194+
end
195+
124196
for (fname, elty) in ((:onemklSsparse_gemm, :Float32),
125197
(:onemklDsparse_gemm, :Float64),
126198
(:onemklCsparse_gemm, :ComplexF32),
@@ -160,6 +232,115 @@ function sparse_optimize_gemm!(trans::Char, transB::Char, nrhs::Int, A::oneSpars
160232
return A
161233
end
162234

235+
# Special handling for CSC matrices since they are stored as transposed CSR (S = A^T)
236+
for (fname, elty) in ((:onemklSsparse_gemm, :Float32),
237+
(:onemklDsparse_gemm, :Float64),
238+
(:onemklCsparse_gemm, :ComplexF32),
239+
(:onemklZsparse_gemm, :ComplexF64))
240+
@eval begin
241+
function sparse_gemm!(transa::Char,
242+
transb::Char,
243+
alpha::Number,
244+
A::oneSparseMatrixCSC{$elty},
245+
B::oneStridedMatrix{$elty},
246+
beta::Number,
247+
C::oneStridedMatrix{$elty})
248+
249+
# Map op(A) to op(S) where S = A^T stored as CSR in the handle
250+
# transa: 'N' -> op(S)='T'; 'T' -> op(S)='N'; 'C' ->
251+
# real: op(S)='N' (since A^H == A^T)
252+
# complex: use conjugation identity on B and C with op(S)='N'
253+
254+
mB, nB = size(B)
255+
mC, nC = size(C)
256+
(nB != nC) && (transb == 'N') && throw(ArgumentError("B and C must have the same number of columns."))
257+
(mB != nC) && (transb != 'N') && throw(ArgumentError("Bᵀ and C must have the same number of columns."))
258+
nrhs = size(B, 2)
259+
ldb = max(1,stride(B,2))
260+
ldc = max(1,stride(C,2))
261+
262+
queue = global_queue(context(C), device())
263+
264+
if transa == 'N'
265+
# Want A * opB -> use S^T * opB
266+
$fname(sycl_queue(queue), 'C', 'T', transb, alpha, A.handle, B, nrhs, ldb, beta, C, ldc)
267+
return C
268+
elseif transa == 'T'
269+
# Want A^T * opB -> use S * opB
270+
$fname(sycl_queue(queue), 'C', 'N', transb, alpha, A.handle, B, nrhs, ldb, beta, C, ldc)
271+
return C
272+
else
273+
# transa == 'C'
274+
if $elty <: Complex
275+
# Use identity: conj(C_new) = conj(alpha) * S * conj(opB(B)) + conj(beta) * conj(C)
276+
# Prepare conj(C) in-place and conj(B) into a temporary if needed
277+
C .= conj.(C)
278+
279+
# Determine how to supply opB under conjugation
280+
# - transb == 'N': pass transb='N' and use conj(B)
281+
# - transb == 'T': pass transb='T' and use conj(B)
282+
# - transb == 'C': since conj(B^H) = B^T, pass transb='T' and use B as-is
283+
local transb_eff::Char
284+
local Beff
285+
if transb == 'N'
286+
transb_eff = 'N'
287+
Beff = similar(B)
288+
Beff .= conj.(B)
289+
elseif transb == 'T'
290+
transb_eff = 'T'
291+
Beff = similar(B)
292+
Beff .= conj.(B)
293+
else
294+
# transb == 'C'
295+
transb_eff = 'T'
296+
Beff = B
297+
end
298+
299+
$fname(sycl_queue(queue), 'C', 'N', transb_eff, conj(alpha), A.handle, Beff, nrhs, ldb, conj(beta), C, ldc)
300+
301+
# Undo conjugation to obtain C_new
302+
C .= conj.(C)
303+
return C
304+
else
305+
# real eltype: A^H == A^T -> use S * opB
306+
$fname(sycl_queue(queue), 'C', 'N', transb, alpha, A.handle, B, nrhs, ldb, beta, C, ldc)
307+
return C
308+
end
309+
end
310+
end
311+
end
312+
end
313+
314+
function sparse_optimize_gemm!(trans::Char, A::oneSparseMatrixCSC)
315+
# Map requested op(A) to op(S) for S = A^T
316+
actual_trans = if trans == 'N'
317+
'T'
318+
elseif trans == 'T'
319+
'N'
320+
else
321+
# 'C' case: complex handled via conjugation with op(S)='N'; real reduces to 'T' which maps to 'N'
322+
'N'
323+
end
324+
queue = global_queue(context(A.nzVal), device(A.nzVal))
325+
onemklXsparse_optimize_gemm(sycl_queue(queue), actual_trans, A.handle)
326+
return A
327+
end
328+
329+
function sparse_optimize_gemm!(trans::Char, transB::Char, nrhs::Int, A::oneSparseMatrixCSC)
330+
# Map as in the basic optimize, and adjust transB for the complex 'C' case if needed at runtime.
331+
# We don't know eltype here; choose conservative mapping for A like above.
332+
actual_trans = if trans == 'N'
333+
'T'
334+
elseif trans == 'T'
335+
'N'
336+
else
337+
'N'
338+
end
339+
queue = global_queue(context(A.nzVal), device(A.nzVal))
340+
onemklXsparse_optimize_gemm_advanced(sycl_queue(queue), 'C', actual_trans, transB, A.handle, nrhs)
341+
return A
342+
end
343+
163344
for (fname, elty) in ((:onemklSsparse_symv, :Float32),
164345
(:onemklDsparse_symv, :Float64),
165346
(:onemklCsparse_symv, :ComplexF32),
@@ -179,6 +360,32 @@ for (fname, elty) in ((:onemklSsparse_symv, :Float32),
179360
end
180361
end
181362

363+
# Special handling for CSC matrices since they are stored as transposed CSR
364+
for (fname, elty) in ((:onemklSsparse_symv, :Float32),
365+
(:onemklDsparse_symv, :Float64),
366+
(:onemklCsparse_symv, :ComplexF32),
367+
(:onemklZsparse_symv, :ComplexF64))
368+
@eval begin
369+
function sparse_symv!(uplo::Char,
370+
alpha::Number,
371+
A::oneSparseMatrixCSC{$elty},
372+
x::oneStridedVector{$elty},
373+
beta::Number,
374+
y::oneStridedVector{$elty})
375+
376+
# CSC(A) is represented by storing CSR(A^T). For symmetric matrices:
377+
# If A is symmetric, then A^T is also symmetric.
378+
# For uplo mapping: 'U' (upper) becomes 'L' (lower) and vice versa
379+
# because the upper triangle of A becomes the lower triangle of A^T
380+
actual_uplo = uplo == 'U' ? 'L' : 'U'
381+
382+
queue = global_queue(context(y), device())
383+
$fname(sycl_queue(queue), actual_uplo, alpha, A.handle, x, beta, y)
384+
y
385+
end
386+
end
387+
end
388+
182389
for (fname, elty) in ((:onemklSsparse_trmv, :Float32),
183390
(:onemklDsparse_trmv, :Float64),
184391
(:onemklCsparse_trmv, :ComplexF32),
@@ -200,6 +407,31 @@ for (fname, elty) in ((:onemklSsparse_trmv, :Float32),
200407
end
201408
end
202409

410+
# Special handling for CSC matrices since they are stored as transposed CSR
411+
for (fname, elty) in ((:onemklSsparse_trmv, :Float32),
412+
(:onemklDsparse_trmv, :Float64),
413+
(:onemklCsparse_trmv, :ComplexF32),
414+
(:onemklZsparse_trmv, :ComplexF64))
415+
@eval begin
416+
function sparse_trmv!(uplo::Char,
417+
trans::Char,
418+
diag::Char,
419+
alpha::Number,
420+
A::oneSparseMatrixCSC{$elty},
421+
x::oneStridedVector{$elty},
422+
beta::Number,
423+
y::oneStridedVector{$elty})
424+
425+
# Intel oneAPI sparse trmv only supports nontrans operations.
426+
# Since CSC(A) is stored as CSR(A^T), we cannot map CSC operations
427+
# to CSR operations for triangular operations without transpose support.
428+
throw(ArgumentError("sparse_trmv! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
429+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
430+
"Convert to oneSparseMatrixCSR format instead."))
431+
end
432+
end
433+
end
434+
203435
function sparse_optimize_trmv!(uplo::Char, trans::Char, diag::Char, A::oneSparseMatrixCSR)
204436
queue = global_queue(context(A.nzVal), device(A.nzVal))
205437
onemklXsparse_optimize_trmv(sycl_queue(queue), uplo, trans, diag, A.handle)
@@ -226,6 +458,30 @@ for (fname, elty) in ((:onemklSsparse_trsv, :Float32),
226458
end
227459
end
228460

461+
# Special handling for CSC matrices since they are stored as transposed CSR
462+
for (fname, elty) in ((:onemklSsparse_trsv, :Float32),
463+
(:onemklDsparse_trsv, :Float64),
464+
(:onemklCsparse_trsv, :ComplexF32),
465+
(:onemklZsparse_trsv, :ComplexF64))
466+
@eval begin
467+
function sparse_trsv!(uplo::Char,
468+
trans::Char,
469+
diag::Char,
470+
alpha::Number,
471+
A::oneSparseMatrixCSC{$elty},
472+
x::oneStridedVector{$elty},
473+
y::oneStridedVector{$elty})
474+
475+
# Intel oneAPI sparse trsv only supports nontrans operations.
476+
# Since CSC(A) is stored as CSR(A^T), we cannot map CSC operations
477+
# to CSR operations for triangular solve operations without transpose support.
478+
throw(ArgumentError("sparse_trsv! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
479+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
480+
"Convert to oneSparseMatrixCSR format instead."))
481+
end
482+
end
483+
end
484+
229485
function sparse_optimize_trsv!(uplo::Char, trans::Char, diag::Char, A::oneSparseMatrixCSR)
230486
queue = global_queue(context(A.nzVal), device(A.nzVal))
231487
onemklXsparse_optimize_trsv(sycl_queue(queue), uplo, trans, diag, A.handle)
@@ -262,6 +518,31 @@ for (fname, elty) in ((:onemklSsparse_trsm, :Float32),
262518
end
263519
end
264520

521+
# Special handling for CSC matrices since they are stored as transposed CSR
522+
for (fname, elty) in ((:onemklSsparse_trsm, :Float32),
523+
(:onemklDsparse_trsm, :Float64),
524+
(:onemklCsparse_trsm, :ComplexF32),
525+
(:onemklZsparse_trsm, :ComplexF64))
526+
@eval begin
527+
function sparse_trsm!(uplo::Char,
528+
transA::Char,
529+
transX::Char,
530+
diag::Char,
531+
alpha::Number,
532+
A::oneSparseMatrixCSC{$elty},
533+
X::oneStridedMatrix{$elty},
534+
Y::oneStridedMatrix{$elty})
535+
536+
# Intel oneAPI sparse trsm only supports nontrans operations for the matrix A.
537+
# Since CSC(A) is stored as CSR(A^T), we cannot map CSC operations
538+
# to CSR operations for triangular solve operations without transpose support.
539+
throw(ArgumentError("sparse_trsm! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
540+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
541+
"Convert to oneSparseMatrixCSR format instead."))
542+
end
543+
end
544+
end
545+
265546
function sparse_optimize_trsm!(uplo::Char, trans::Char, diag::Char, A::oneSparseMatrixCSR)
266547
queue = global_queue(context(A.nzVal), device(A.nzVal))
267548
onemklXsparse_optimize_trsm(sycl_queue(queue), uplo, trans, diag, A.handle)
@@ -273,3 +554,36 @@ function sparse_optimize_trsm!(uplo::Char, trans::Char, diag::Char, nrhs::Int, A
273554
onemklXsparse_optimize_trsm_advanced(sycl_queue(queue), 'C', uplo, trans, diag, A.handle, nrhs)
274555
return A
275556
end
557+
558+
# CSC optimization functions
559+
function sparse_optimize_trmv!(uplo::Char, trans::Char, diag::Char, A::oneSparseMatrixCSC)
560+
# Intel oneAPI sparse trmv only supports nontrans operations.
561+
# Since CSC(A) is stored as CSR(A^T), triangular operations are not supported for CSC format.
562+
throw(ArgumentError("sparse_optimize_trmv! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
563+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
564+
"Convert to oneSparseMatrixCSR format instead."))
565+
end
566+
567+
function sparse_optimize_trsv!(uplo::Char, trans::Char, diag::Char, A::oneSparseMatrixCSC)
568+
# Intel oneAPI sparse trsv only supports nontrans operations.
569+
# Since CSC(A) is stored as CSR(A^T), triangular operations are not supported for CSC format.
570+
throw(ArgumentError("sparse_optimize_trsv! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
571+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
572+
"Convert to oneSparseMatrixCSR format instead."))
573+
end
574+
575+
function sparse_optimize_trsm!(uplo::Char, trans::Char, diag::Char, A::oneSparseMatrixCSC)
576+
# Intel oneAPI sparse trsm only supports nontrans operations.
577+
# Since CSC(A) is stored as CSR(A^T), triangular operations are not supported for CSC format.
578+
throw(ArgumentError("sparse_optimize_trsm! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
579+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
580+
"Convert to oneSparseMatrixCSR format instead."))
581+
end
582+
583+
function sparse_optimize_trsm!(uplo::Char, trans::Char, diag::Char, nrhs::Int, A::oneSparseMatrixCSC)
584+
# Intel oneAPI sparse trsm only supports nontrans operations.
585+
# Since CSC(A) is stored as CSR(A^T), triangular operations are not supported for CSC format.
586+
throw(ArgumentError("sparse_optimize_trsm! is not supported for oneSparseMatrixCSC due to Intel oneAPI limitations. " *
587+
"Intel sparse library only supports nontrans operations for triangular matrix operations. " *
588+
"Convert to oneSparseMatrixCSR format instead."))
589+
end

‎test/Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,4 +19,5 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1919
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
2020
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2121
libigc_jll = "94295238-5935-5bd7-bb0f-b00942e9bdd5"
22+
oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"
2223
oneAPI_Support_jll = "b049733a-a71d-5ed3-8eba-7d323ac00b36"

0 commit comments

Comments
 (0)