[LAPACK][rocSOLVER] Fix incorrect QR results on repeated runs (#626)

qr_diag_repro.cpp

@@ -0,0 +1,98 @@
+// Standalone reproducer mirroring the geqrf+orgqr diagonal regression test
+// added for https://github.com/uxlfoundation/oneMath/issues/626.
+// Runs geqrf + orgqr repeatedly on a diagonal matrix (whose Q is identity)
+// and verifies Q's diagonal is 1 on every run.
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+#if __has_include(<sycl/sycl.hpp>)
+#include <sycl/sycl.hpp>
+#else
+#include <CL/sycl.hpp>
+#endif
+
+#include "oneapi/math.hpp"
+
+template <typename fp>
+bool run_case(sycl::queue& queue, std::int64_t n, std::int64_t lda, std::int64_t num_runs) {
+    const std::int64_t m = n;
+    const std::int64_t k = n;
+
+    const fp diag_value = fp(2.0);
+    std::vector<fp> A_initial(static_cast<size_t>(lda) * n, fp(0.0));
+    for (std::int64_t j = 0; j < n; ++j)
+        A_initial[j * lda + j] = diag_value;
+
+    const fp tol = fp(10.0) * n * std::numeric_limits<fp>::epsilon();
+
+    auto* A_dev = sycl::malloc_device<fp>(A_initial.size(), queue);
+    auto* tau_dev = sycl::malloc_device<fp>(k, queue);
+
+    const std::int64_t geqrf_ss = oneapi::math::lapack::geqrf_scratchpad_size<fp>(queue, m, n, lda);
+    const std::int64_t orgqr_ss =
+        oneapi::math::lapack::orgqr_scratchpad_size<fp>(queue, m, n, k, lda);
+    const std::int64_t ss = std::max(geqrf_ss, orgqr_ss);
+    auto* scratch_dev = sycl::malloc_device<fp>(ss, queue);
+
+    std::vector<fp> Q(A_initial.size());
+    bool ok = true;
+    for (std::int64_t run = 0; run < num_runs && ok; ++run) {
+        queue.memcpy(A_dev, A_initial.data(), A_initial.size() * sizeof(fp)).wait();
+
+        oneapi::math::lapack::geqrf(queue, m, n, A_dev, lda, tau_dev, scratch_dev, geqrf_ss);
+        oneapi::math::lapack::orgqr(queue, m, n, k, A_dev, lda, tau_dev, scratch_dev, orgqr_ss);
+        queue.wait_and_throw();
+
+        queue.memcpy(Q.data(), A_dev, Q.size() * sizeof(fp)).wait();
+
+        int bad = 0;
+        for (std::int64_t j = 0; j < n && bad < 5; ++j) {
+            for (std::int64_t i = 0; i < m && bad < 5; ++i) {
+                const fp expected = (i == j) ? fp(1.0) : fp(0.0);
+                const fp got = std::abs(Q[j * lda + i]);
+                if (std::abs(got - expected) > tol) {
+                    std::cout << "    run " << run << ": Q(" << i << "," << j << ")=" << got
+                              << " expected " << expected << "\n";
+                    ++bad;
+                    ok = false;
+                }
+            }
+        }
+    }
+
+    sycl::free(A_dev, queue);
+    sycl::free(tau_dev, queue);
+    sycl::free(scratch_dev, queue);
+    return ok;
+}
+
+int main() {
+    sycl::queue queue{ sycl::gpu_selector_v };
+    std::cout << "Device: " << queue.get_device().get_info<sycl::info::device::name>() << "\n";
+
+    struct Case {
+        std::int64_t n, lda, runs;
+    };
+    const std::vector<Case> cases = { { 256, 256, 6 }, { 512, 512, 6 }, { 1024, 1024, 6 } };
+
+    bool all_ok = true;
+    for (auto c : cases) {
+        std::cout << "[float ] n=" << c.n << " runs=" << c.runs << ": " << std::flush;
+        bool ok = run_case<float>(queue, c.n, c.lda, c.runs);
+        std::cout << (ok ? "PASS" : "FAIL") << "\n";
+        all_ok &= ok;
+
+        std::cout << "[double] n=" << c.n << " runs=" << c.runs << ": " << std::flush;
+        ok = run_case<double>(queue, c.n, c.lda, c.runs);
+        std::cout << (ok ? "PASS" : "FAIL") << "\n";
+        all_ok &= ok;
+    }
+
+    std::cout << (all_ok ? "ALL PASS" : "SOME FAILED") << "\n";
+    return all_ok ? 0 : 1;
+}

src/lapack/backends/cusolver/cusolver_helper.hpp

@@ -290,6 +290,16 @@ struct CudaEquivalentType<std::complex<double>> {
 
 /* devinfo */
 
+/* Allocates the device memory cuSOLVER writes its `devInfo` output to. A failed
+ * allocation must not be passed on to cuSOLVER as a null pointer, since that
+ * silently disables the routine's error reporting. */
+inline int* create_devinfo(sycl::queue& queue, const char* func_name, int dev_info_size = 1) {
+    int* devInfo = sycl::malloc_device<int>(dev_info_size, queue);
+    if (!devInfo)
+        throw oneapi::math::device_bad_alloc("lapack", func_name, queue.get_device());
+    return devInfo;
+}
+
 inline void get_cusolver_devinfo(sycl::queue& queue, sycl::buffer<int>& devInfo,
                                  std::vector<int>& dev_info_) {
     sycl::host_accessor<int, 1, sycl::access::mode::read> dev_info_acc{ devInfo };

src/lapack/backends/cusolver/cusolver_scope_handle.cpp

@@ -141,7 +141,13 @@ cusolverDnHandle_t CusolverScopedContextHandler::get_handle(const sycl::queue& q
 }
 
 CUstream CusolverScopedContextHandler::get_stream(const sycl::queue& queue) {
-    return sycl::get_native<sycl::backend::ext_oneapi_cuda>(queue);
+    // Return the CUDA stream backing the current host task / native command
+    // submission. When `ext_codeplay_enqueue_native_command` is used, the SYCL
+    // runtime tracks completion of the native command via this stream, so
+    // cuSOLVER work must be enqueued onto it (not the queue's default stream)
+    // to avoid the returned SYCL event signalling before the computation
+    // finishes. See uxlfoundation/oneMath#626.
+    return ih.get_native_queue<sycl::backend::ext_oneapi_cuda>();
 }
 sycl::context CusolverScopedContextHandler::get_context(const sycl::queue& queue) {
     return queue.get_context();

src/lapack/backends/rocsolver/rocsolver_helper.hpp

@@ -150,12 +150,6 @@ class hip_error : virtual public std::runtime_error {
         throw rocsolver_error(std::string(#name) + std::string(" : "), err); \
     }
 
-#define ROCSOLVER_ERROR_FUNC_T(name, func, err, ...)                        \
-    err = func(__VA_ARGS__);                                                \
-    if (err != rocblas_status_success) {                                    \
-        throw rocsolver_error(std::string(name) + std::string(" : "), err); \
-    }
-
 #define ROCSOLVER_ERROR_FUNC_T_SYNC(name, func, err, handle, ...)            \
     err = func(handle, __VA_ARGS__);                                         \
     if (err != rocblas_status_success) {                                     \
@@ -166,14 +160,18 @@ class hip_error : virtual public std::runtime_error {
     hipError_t hip_err;                                                      \
     HIP_ERROR_FUNC(hipStreamSynchronize, hip_err, currentStreamId);
 
+/* rocSOLVER calls only enqueue work on the HIP stream of the surrounding host
+ * task or native command. A submission depending on that work may be scheduled
+ * on a different stream of the queue's stream pool without being ordered against
+ * this one, and then reads a computation that has not finished yet. Unlike
+ * getrf/potrf/..., the Householder routines have no `info` output either, so they
+ * do not get an implicit host-side wait from `lapack_info_check`. Waiting for the
+ * stream keeps the enqueued work complete once the call returns.
+ * See uxlfoundation/oneMath#626. */
 template <class Func, class... Types>
 inline void rocsolver_native_named_func(const char* func_name, Func func, rocsolver_status err,
                                         rocsolver_handle handle, Types... args) {
-#ifdef SYCL_EXT_ONEAPI_ENQUEUE_NATIVE_COMMAND
-    ROCSOLVER_ERROR_FUNC_T(func_name, func, err, handle, args...)
-#else
     ROCSOLVER_ERROR_FUNC_T_SYNC(func_name, func, err, handle, args...)
-#endif
 };
 
 inline rocblas_eform get_rocsolver_itype(std::int64_t itype) {
@@ -257,6 +255,16 @@ struct RocmEquivalentType<std::complex<double>> {
 
 /* devinfo */
 
+/* Allocates the device memory rocSOLVER writes its `devInfo` output to. A failed
+ * allocation must not be passed on to rocSOLVER as a null pointer, since that
+ * silently disables the routine's error reporting. */
+inline int* create_devinfo(sycl::queue& queue, const char* func_name, int dev_info_size = 1) {
+    int* devInfo = sycl::malloc_device<int>(dev_info_size, queue);
+    if (!devInfo)
+        throw oneapi::math::device_bad_alloc("lapack", func_name, queue.get_device());
+    return devInfo;
+}
+
 inline int get_rocsolver_devinfo(sycl::queue& queue, sycl::buffer<int>& devInfo) {
     sycl::host_accessor<int, 1, sycl::access::mode::read> dev_info_{ devInfo };
     return dev_info_[0];

src/lapack/backends/rocsolver/rocsolver_lapack.cpp

@@ -1273,7 +1273,7 @@ inline sycl::event getrf(const char* func_name, Func func, sycl::queue& queue, s
     std::uint64_t ipiv_size = std::min(n, m);
     int* ipiv32 = (int*)malloc_device(sizeof(int) * ipiv_size, queue);
 
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1411,7 +1411,7 @@ inline sycl::event gesvd(const char* func_name, Func func, sycl::queue& queue,
     using rocmDataType_A = typename RocmEquivalentType<T_A>::Type;
     using rocmDataType_B = typename RocmEquivalentType<T_B>::Type;
     overflow_check(m, n, lda, ldu, ldvt, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1462,7 +1462,7 @@ inline sycl::event heevd(const char* func_name, Func func, sycl::queue& queue,
     using rocmDataType_A = typename RocmEquivalentType<T_A>::Type;
     using rocmDataType_B = typename RocmEquivalentType<T_B>::Type;
     overflow_check(n, lda, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1508,7 +1508,7 @@ inline sycl::event hegvd(const char* func_name, Func func, sycl::queue& queue, s
     using rocmDataType_A = typename RocmEquivalentType<T_A>::Type;
     using rocmDataType_B = typename RocmEquivalentType<T_B>::Type;
     overflow_check(n, lda, ldb, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1818,7 +1818,7 @@ inline sycl::event potrf(const char* func_name, Func func, sycl::queue& queue,
                          const std::vector<sycl::event>& dependencies) {
     using rocmDataType = typename RocmEquivalentType<T>::Type;
     overflow_check(n, lda, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1860,7 +1860,7 @@ inline sycl::event potri(const char* func_name, Func func, sycl::queue& queue,
                          const std::vector<sycl::event>& dependencies) {
     using rocmDataType = typename RocmEquivalentType<T>::Type;
     overflow_check(n, lda, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1944,7 +1944,7 @@ inline sycl::event syevd(const char* func_name, Func func, sycl::queue& queue,
                          const std::vector<sycl::event>& dependencies) {
     using rocmDataType = typename RocmEquivalentType<T>::Type;
     overflow_check(n, lda, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -1989,7 +1989,7 @@ inline sycl::event sygvd(const char* func_name, Func func, sycl::queue& queue, s
                          const std::vector<sycl::event>& dependencies) {
     using rocmDataType = typename RocmEquivalentType<T>::Type;
     overflow_check(n, lda, ldb, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
     auto done = queue.submit([&](sycl::handler& cgh) {
         int64_t num_events = dependencies.size();
         for (int64_t i = 0; i < num_events; i++) {
@@ -2074,7 +2074,7 @@ inline sycl::event sytrf(const char* func_name, Func func, sycl::queue& queue,
                          const std::vector<sycl::event>& dependencies) {
     using rocmDataType = typename RocmEquivalentType<T>::Type;
     overflow_check(n, lda, scratchpad_size);
-    int* devInfo = (int*)malloc_device(sizeof(int), queue);
+    int* devInfo = create_devinfo(queue, __func__);
 
     // rocsolver legacy api does not accept 64-bit ints.
     // To get around the limitation.

src/lapack/backends/rocsolver/rocsolver_scope_handle.cpp

@@ -143,7 +143,10 @@ rocblas_handle RocsolverScopedContextHandler::get_handle(const sycl::queue& queu
 }
 
 hipStream_t RocsolverScopedContextHandler::get_stream(const sycl::queue& queue) {
-    return sycl::get_native<sycl::backend::ext_oneapi_hip>(queue);
+    // Return the HIP stream backing the current host task / native command
+    // submission, which is the stream the SYCL runtime records the submission's
+    // completion event on. Matches the cuBLAS and cuSOLVER backends.
+    return ih.get_native_queue<sycl::backend::ext_oneapi_hip>();
 }
 sycl::context RocsolverScopedContextHandler::get_context(const sycl::queue& queue) {
     return queue.get_context();

tests/unit_tests/CMakeLists.txt

@@ -31,9 +31,7 @@ endif()
 if("lapack" IN_LIST TEST_TARGET_DOMAINS)
   find_package(LAPACKE)
   if(NOT LAPACKE_FOUND)
-    # TODO: add list of tests without Netlib dependency
-    message(WARNING "Netlib LAPACKE headers or libraries are not found, LAPACK unit tests will be skipped")
-    list(REMOVE_ITEM TEST_TARGET_DOMAINS "lapack")
+    message(WARNING "Netlib LAPACKE headers or libraries are not found, LAPACK unit tests comparing against a reference implementation will be skipped")
   endif()
 endif()
 
@@ -66,9 +64,7 @@ set(blas_TEST_LINK "")
 set(lapack_TEST_LIST
         lapack_source)
 
-if(LAPACKE_FOUND)
-  set(lapack_TEST_LINK "")
-endif()
+set(lapack_TEST_LINK "")
 
 # RNG config
 set(rng_TEST_LIST

tests/unit_tests/lapack/source/CMakeLists.txt

@@ -19,7 +19,9 @@
 
 #Build object from all test sources
 # TODO: add list of tests without Netlib dependency
-set(LAPACK_SOURCES)
+set(LAPACK_SOURCES
+    "geqrf_orgqr_diagonal.cpp"
+)
 
 set(LAPACK_SOURCES_W_LAPACKE
     "gebrd.cpp"