oneMKL functions taking multiple buffers and their aliasing/overlapping

[jakub-homola]

Hi,

can we pass the same buffer as input and output argument of oneMKL functions? Is "buffer aliasing" allowed?

E.g. the function oneapi::mkl::sparse::trsv. It takes multiple buffers as arguments, namely the input right-hand-side vector and the output solution vector. Can I pass the same buffer for both of these arguments? (mimicking the oneapi::mkl::blas::row_major::trsv behaviour, where there is just one argument for that) (experimentally tested with intel onemkl and it works, but don't know if it is a guarantee).

In general, can buffers passed to oneMKL functions overlap, alias, or be the same? I could not find any information about this in the spec (https://spec.oneapi.io/versions/1.2-rev-1/elements/oneMKL/source/domains/spblas/trsv.html and other pages on that site). Is it specified somewhere?

This might also apply to other one* libraries, but I am not familiar with them.

Thanks,
Jakub

[jakub-homola]

So, according to this reply in the intel onemkl forum, it is not allowed. But I got that information quite randomly, and only for the intel's onemkl. Still not sure how it should work, and where to find that information.

[mmeterel]

@jakub-homola Thanks for the question. It is hard to answer for all oneMKL since each domain might have different requirements/restrictions. I pinged @spencerpatty for your question on sparse::trsv.

For BLAS side, it is ok to pass the same buffer for input and output - like matrix C in GEMM.

Regarding overlap or alias questions, I am pinging @petercad.

[petercad]

Hi @jakub-homola, thanks for raising this issue! You're right, the spec should explain which arguments (if any) are allowed to alias.

For (non-sparse) BLAS routines, no overlap is allowed between an output matrix/vector and any other matrix/vector arguments. This is the standard convention for BLAS, but we should spell that out. I had meant to add that to the spec but forgot, so thanks for the reminder :) We'll take care of that.

[spencerpatty]

Hi @jakub-homola this is an interesting question. So the TRSV operation is solving (lets focus on lower solve) for y in the system "L * y = x". You are asking if x and y can be the same array. That is, you are asking if you can triangular solve for y in "L * y = y".

To me, mathematically this operation no longer makes sense if L is not an identity matrix or y is not an eigenvector corresponding to an eigenvalue lambda = 1 to begin with.

Practically, you could put the same array into the function as x and y and it might do something, but lets think a bit about what it might actually be doing...

The traditional CSR matrix sparse TRSV (L = rowptr, colind, vals CSR arrays and an extra invDiag array) algorithm is a variant (adding parallelism respecting the L sparsity pattern) of the following steps

y[0] = x[0] * invDiag[0];
...
y[i] = ( x[i] - sum_{k in (rowptr[i], lastLower[i] ) } ( val[k] * y[colind[k]] )  ) * invDiag[i];
...

or more simply thinking in 2D dense matrix terms

y[0] = x[0] / A[0][0]
...
y[i] = (x[i] - sum_{k in (0...i-1)} (A[i][k] * y[k])  ) / A[i][i]
...

so it looks like if x[i] is being implicitly updated when we change y[i], then maybe it actually does nothing to affect the end result (y output) for this algorithm... so maybe it turns out ok in this particular case.

In general, we don't envision the input/output matrices or vectors being the same array or overlapped in sparse BLAS either. We should spell it out as well in our documentation.

[spencerpatty]

I updated my previous response with a better answer, but updates don't get sent out with notice, so I am triggering a notice with this comment, in case it helps :)

[jakub-homola]

Hi,
thanks for the answer @spencerpatty.
I don't mean solving $Ly=y$. I still mean to solve $Ly=x$, but use just one function parameter for both the input right hand side and output solution, just as dense blas trsv does it. The code sample you showed suggests that it actually makes sense (at least for that algorithm), and we would not even need to allocate another array.

Anyway, thanks everyone for the info that the arrays cannot alias/overlap/be the same. Looking forward to the spec/documentation update.