ENH: Implement GPU version of spectral decomposition as linear regression fallback

[david-cortes-intel]

Description

As a follow up from https://github.com/oneapi-src/oneDAL/pull/2930

This PR adds a GPU version of the fallback for linear regression based on spectral decomposition, which mirrors the CPU version. It reuses the same tests, but this time executing them with fp32 too.

A few things I wasn't sure about:

• Here I'm taking the minimum of the diagonal of the Cholesky factorization using a sycl reduction class, but I'm not sure if that's the usual way in oneDAL.
• I am also not 100% sure if it is guaranteed that this operation (minimum of the diagonal) with an 'at_device' at the end would force execution to wait until the operation is completed before returning the value.
• There's an operation here for which oneMKL doesn't provide the equivalent LAPACK functionality, which is: "column-major matrix divided by a broadcasted row vector" (i.e. multiply each row of the matrix by a scalar from the vector). I implemented it as a 2D for-loop with a multiplication by the inverse (which is faster but less precise), but I don't know if this is better than a 1D loop that would take a whole row, as rows are contiguous in memory.

---

PR should start as a draft, then move to ready for review state after CI is passed and all applicable checkboxes are closed. This approach ensures that reviewers don't spend extra time asking for regular requirements.

You can remove a checkbox as not applicable only if it doesn't relate to this PR in any way. For example, PR with docs update doesn't require checkboxes for performance while PR with any change in actual code should have checkboxes and justify how this code change is expected to affect performance (or justification should be self-evident).

Checklist to comply with before moving PR from draft:

PR completeness and readability

• [x] I have reviewed my changes thoroughly before submitting this pull request.
• [x] I have commented my code, particularly in hard-to-understand areas.
• [x] Git commit message contains an appropriate signed-off-by string (see CONTRIBUTING.md for details).
• [x] I have added a respective label(s) to PR if I have a permission for that.
• [x] I have resolved any merge conflicts that might occur with the base branch.

Testing

• [x] I have run it locally and tested the changes extensively.
• [x] All CI jobs are green or I have provided justification why they aren't.
• [x] I have extended testing suite if new functionality was introduced in this PR.

Performance

Not applicable (addition is a fallback for cases that would have failed previously).

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

Oddly, the first version of this function (minimum of the diagonal) that I wrote seems to succeed in some tests, but fail in others (particularly when it it called repeatedly, as in incremental linear regression):

template <typename Float>
Float diagonal_minimum(sycl::queue& queue,
                       const Float* Matrix,
                       const std::int64_t dim_matrix,
                       sycl::event& event_Matrix) {
    constexpr auto alloc = sycl::usm::alloc::device;
    auto diag_min_holder = array<Float>::empty(queue, 1, alloc);
    sycl::event diag_min_event = queue.submit([&](sycl::handler& h) {
        auto min_reduction = sycl::reduction(diag_min_holder.get_mutable_data(), sycl::minimum<>());
        h.depends_on(event_Matrix);
        h.parallel_for(dim_matrix, min_reduction, [=](const auto& i, auto& min_obj) {
            min_obj.combine(Matrix[i * (dim_matrix + 1)]);
        });
    });
    return ndview<Float, 1>::wrap(diag_min_holder).at_device(queue, 0, { diag_min_event });
}

(tried also adding event_Matrix to the dependencies of the return line, but it made no difference)

But this slower version doesn't seem to have any issue:

template <typename Float>
Float diagonal_minimum(sycl::queue& queue,
                       const Float* Matrix,
                       const std::int64_t dim_matrix,
                       sycl::event& event_Matrix) {
    constexpr auto alloc = sycl::usm::alloc::device;
    auto idx_min_holder = array<std::int64_t>::empty(queue, 1, alloc);
    sycl::event diag_min_event = mkl::blas::column_major::iamin(
        queue,
        dim_matrix,
        Matrix,
        dim_matrix + 1,
        idx_min_holder.get_mutable_data(),
        mkl::index_base::zero,
        { event_Matrix }
    );
    const std::int64_t idx_min = ndview<std::int64_t, 1>::wrap(idx_min_holder).at_device(queue, 0, { diag_min_event });
    return ndview<Float, 1>::wrap(Matrix, dim_matrix * dim_matrix).at_device(queue, idx_min * (dim_matrix + 1), { event_Matrix });
}

Any hints and what might be happening?

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

Another thing that I'm now not very sure about is: should calls to oneMKL lapack manually call queue.wait_and_throw() in order to catch runtime exceptions from them?

This is what the oneMKL docs mention:

oneMKL error handling relies on the mechanism of C++ exceptions. Should error occur, it will be propagated at the point of a function call where it is caught using standard C++ error handling mechanism.

... which sounds like it would let it wait until the exception can be thrown, but then it leaves me wondering how it handles queueing.

[david-cortes-intel]

After some investigation regarding the diagonal_minimum function, turns out the issue with the function for diagonal minimum was that it wasn't initializing the reducer object, despite what the docs mention: https://github.khronos.org/SYCL_Reference/iface/reduction-variables.html#known-identities

Adding a manual initialization did the trick.

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

/intelci: run

[Vika-F]

Another thing that I'm now not very sure about is: should calls to oneMKL lapack manually call queue.wait_and_throw() in order to catch runtime exceptions from them?

@david-cortes-intel Good point. I think yes, .wait_and_throw() is needed after the last oneMKL function call in try to catch the exception. Otherwise it can be missed as oneMKL calls run asynchronously on device.

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

Another thing that I'm now not very sure about is: should calls to oneMKL lapack manually call queue.wait_and_throw() in order to catch runtime exceptions from them?

@david-cortes-intel Good point. I think yes, .wait_and_throw() is needed after the last oneMKL function call in try to catch the exception. Otherwise it can be missed as oneMKL calls run asynchronously on device.

Thanks for confirming. Added a manual call to wait_and_throw.

[Vika-F]

/intelci: run

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

/intelci: run

[david-cortes-intel]

CI is passing now, will merge.