Optimize LU factorization without pivoting

clients/gtest/lapack/getf2_getrf_gtest.cpp

@@ -73,9 +73,7 @@ const vector<vector<int>> large_matrix_size_range = {
     {1000, 1024, 0},
 };
 
-const vector<int> large_n_size_range = {
-    45, 64, 520, 1024, 2000,
-};
+const vector<int> large_n_size_range = {45, 64, 520, 1024, 2000};
 
 Arguments getrf_setup_arguments(getrf_tuple tup)
 {

library/src/CMakeLists.txt

@@ -260,6 +260,11 @@ set(rocsolver_specialized_source
   specialized/roclapack_potf2_specialized_kernels_d.cpp
   specialized/roclapack_potf2_specialized_kernels_c.cpp
   specialized/roclapack_potf2_specialized_kernels_z.cpp
+  # getf2_nopiv
+  specialized/roclapack_getf2_nopiv_specialized_kernels_s.cpp
+  specialized/roclapack_getf2_nopiv_specialized_kernels_d.cpp
+  specialized/roclapack_getf2_nopiv_specialized_kernels_c.cpp
+  specialized/roclapack_getf2_nopiv_specialized_kernels_z.cpp
 )
 
 if(OPTIMAL)

library/src/include/ideal_sizes.hpp

@@ -387,6 +387,21 @@
 #define GETRF_NPVT_BATCH_BLKSIZES_COMPLEX 0, -16, -32, -48, 64, 128
 #endif
 
+// ---------------------------------------------------------------
+// size of submatrix that can fit in 64 KBytes of LDS shared memory
+// ---------------------------------------------------------------
+#ifndef GETRF_NOPIV_BLOCKSIZE
+#define GETRF_NOPIV_BLOCKSIZE(T) ((sizeof(T) == 4) ? 128 : (sizeof(T) == 8) ? 90 : 64)
+#endif
+
+// --------------------------------------------
+// assume there is at least 8 MBytes of last level cache
+// terminate recursion if matrix can fit in cache
+// --------------------------------------------
+#ifndef GETRF_NOPIV_STOPPING_NB
+#define GETRF_NOPIV_STOPPING_NB(T) ((sizeof(T) == 4) ? 1408 : (sizeof(T) == 8) ? 1024 : 704)
+#endif
+
 /****************************** getri *****************************************
 *******************************************************************************/
 #ifndef GETRI_MAX_COLS

library/src/include/rocsolver_run_specialized_kernels.hpp

@@ -258,6 +258,16 @@ rocblas_status potf2_run_small(rocblas_handle handle,
                                rocblas_int* info,
                                const rocblas_int batch_count);
 
+template <typename T, typename U>
+rocblas_status getf2_nopiv_run_small(rocblas_handle handle,
+                                     const rocblas_int n,
+                                     U AA,
+                                     const rocblas_int shiftA,
+                                     const rocblas_int lda,
+                                     const rocblas_stride strideA,
+                                     rocblas_int* info,
+                                     const rocblas_int batch_count);
+
 #ifdef OPTIMAL
 
 template <typename T, typename I, typename INFO, typename U>

library/src/lapack/roclapack_getf2_nopiv.hpp

@@ -0,0 +1,120 @@
+/************************************************************************
+ * Derived from the BSD3-licensed
+ * LAPACK routine (version 3.7.0) --
+ *     Univ. of Tennessee, Univ. of California Berkeley,
+ *     Univ. of Colorado Denver and NAG Ltd..
+ *     December 2016
+ * Copyright (C) 2019-2024 Advanced Micro Devices, Inc. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ * *************************************************************************/
+
+#pragma once
+
+#include "rocblas.hpp"
+#include "rocsolver/rocsolver.h"
+#include "rocsolver_run_specialized_kernels.hpp"
+
+ROCSOLVER_BEGIN_NAMESPACE
+
+template <typename T>
+rocblas_status rocsolver_getf2_nopiv_potrf_argCheck(rocblas_handle handle,
+                                                    const rocblas_int n,
+                                                    const rocblas_int lda,
+                                                    T A,
+                                                    rocblas_int* info,
+                                                    const rocblas_int batch_count = 1)
+{
+    // order is important for unit tests:
+
+    // 1. invalid size
+    if(n < 0 || lda < n || batch_count < 0)
+        return rocblas_status_invalid_size;
+
+    // 2. skip pointer check if querying memory size
+    if(rocblas_is_device_memory_size_query(handle))
+        return rocblas_status_continue;
+
+    // 3. invalid pointers
+    if((n && !A) || (batch_count && !info))
+        return rocblas_status_invalid_pointer;
+
+    return rocblas_status_continue;
+}
+
+template <typename T, typename U, bool COMPLEX = rocblas_is_complex<T>>
+rocblas_status rocsolver_getf2_nopiv_template(rocblas_handle handle,
+                                              const rocblas_int n,
+                                              U A,
+                                              const rocblas_int shiftA,
+                                              const rocblas_int lda,
+                                              const rocblas_stride strideA,
+                                              rocblas_int* info,
+                                              const rocblas_int batch_count)
+{
+    ROCSOLVER_ENTER("getf2_nopiv", "n:", n, "shiftA:", shiftA, "lda:", lda, "bc:", batch_count);
+
+    // quick return if zero instances in batch
+    if(batch_count == 0)
+        return rocblas_status_success;
+
+    hipStream_t stream;
+    rocblas_get_stream(handle, &stream);
+
+    rocblas_int blocksReset = (batch_count - 1) / BS1 + 1;
+    dim3 gridReset(blocksReset, 1, 1);
+    dim3 threads(BS1, 1, 1);
+
+    // info=0 (starting with a positive definite matrix)
+    ROCSOLVER_LAUNCH_KERNEL(reset_info, gridReset, threads, 0, stream, info, batch_count, 0);
+
+    // quick return if no dimensions
+    if(n == 0)
+        return rocblas_status_success;
+
+    // everything must be executed with scalars on the device
+    rocblas_pointer_mode old_mode;
+    rocblas_get_pointer_mode(handle, &old_mode);
+    rocblas_set_pointer_mode(handle, rocblas_pointer_mode_device);
+
+    assert(n <= GETRF_NOPIV_BLOCKSIZE(T));
+
+    rocblas_status istat = rocblas_status_success;
+    if(n <= GETRF_NOPIV_BLOCKSIZE(T))
+    {
+        // ----------------------
+        // use specialized kernel
+        // ----------------------
+        getf2_nopiv_run_small<T>(handle, n, A, shiftA, lda, strideA, info, batch_count);
+    }
+    else
+    {
+        istat = rocblas_status_internal_error;
+    }
+
+    rocblas_set_pointer_mode(handle, old_mode);
+    return istat;
+}
+
+ROCSOLVER_END_NAMESPACE

library/src/lapack/roclapack_getrf.cpp

@@ -26,6 +26,7 @@
  * *************************************************************************/
 
 #include "roclapack_getrf.hpp"
+#include "roclapack_getrf_nopiv.hpp"
 
 ROCSOLVER_BEGIN_NAMESPACE
 
@@ -109,6 +110,84 @@ rocblas_status rocsolver_getrf_impl(rocblas_handle handle,
         optim_mem, pivot);
 }
 
+template <typename T, typename U>
+rocblas_status rocsolver_getrf_nopiv_impl(rocblas_handle handle,
+                                          const rocblas_int m,
+                                          const rocblas_int n,
+                                          U A,
+                                          const rocblas_int lda,
+                                          rocblas_int* info)
+{
+    const char* name = "getrf_nopiv";
+    ROCSOLVER_ENTER_TOP(name, "-m", m, "-n", n, "--lda", lda);
+
+    using S = decltype(std::real(T{}));
+
+    if(!handle)
+        return rocblas_status_invalid_handle;
+
+    // argument checking
+    {
+        bool const pivot = false;
+        rocblas_int* const ipiv = nullptr;
+
+        rocblas_status st = rocsolver_getf2_getrf_argCheck(handle, m, n, lda, A, ipiv, info, pivot);
+        if(st != rocblas_status_continue)
+            return st;
+    }
+
+    // working with unshifted arrays
+    rocblas_int shiftA = 0;
+
+    // normal (non-batched non-strided) execution
+    rocblas_stride strideA = 0;
+    rocblas_int batch_count = 1;
+
+    // memory workspace sizes:
+    //
+    // size of reusable workspace (and for calling TRSM)
+    bool optim_mem = false;
+    size_t size_work1 = 0;
+    size_t size_work2 = 0;
+    size_t size_work3 = 0;
+    size_t size_work4 = 0;
+
+    // size to store info about singularity of each subblock
+    size_t size_iinfo = 0;
+
+    constexpr bool is_batched = false;
+    constexpr bool is_strided = false;
+    rocsolver_getrf_nopiv_getMemorySize<is_batched, is_strided, T>(
+        m, n, batch_count, &size_work1, &size_work2, &size_work3, &size_work4, &size_iinfo,
+        &optim_mem);
+
+    if(rocblas_is_device_memory_size_query(handle))
+        return rocblas_set_optimal_device_memory_size(handle, size_work1, size_work2, size_work3,
+                                                      size_work4, size_iinfo);
+
+    // memory workspace allocation
+
+    void* work1 = nullptr;
+    void* work2 = nullptr;
+    void* work3 = nullptr;
+    void* work4 = nullptr;
+    void* iinfo = nullptr;
+    rocblas_device_malloc mem(handle, size_work1, size_work2, size_work3, size_work4, size_iinfo);
+
+    if(!mem)
+        return rocblas_status_memory_error;
+
+    work1 = mem[0];
+    work2 = mem[1];
+    work3 = mem[2];
+    work4 = mem[3];
+    iinfo = mem[4];
+
+    // execution
+    return rocsolver_getrf_nopiv_template<false, false, T>(handle, m, n, A, shiftA, lda, strideA,
+                                                           info, batch_count, work1, work2, work3,
+                                                           work4, (rocblas_int*)iinfo, optim_mem);
+}
 ROCSOLVER_END_NAMESPACE
 
 /*
@@ -165,6 +244,7 @@ rocblas_status rocsolver_zgetrf(rocblas_handle handle,
                                                                    true);
 }
 
+static constexpr bool use_getrf_nopiv = true;
 rocblas_status rocsolver_sgetrf_64(rocblas_handle handle,
                                    const int64_t m,
                                    const int64_t n,
@@ -218,8 +298,15 @@ rocblas_status rocsolver_sgetrf_npvt(rocblas_handle handle,
                                      const rocblas_int lda,
                                      rocblas_int* info)
 {
-    rocblas_int* ipiv = nullptr;
-    return rocsolver::rocsolver_getrf_impl<float>(handle, m, n, A, lda, ipiv, info, false);
+    if(use_getrf_nopiv)
+    {
+        return rocsolver::rocsolver_getrf_nopiv_impl<float>(handle, m, n, A, lda, info);
+    }
+    else
+    {
+        rocblas_int* ipiv = nullptr;
+        return rocsolver::rocsolver_getrf_impl<float>(handle, m, n, A, lda, ipiv, info, false);
+    }
 }
 
 rocblas_status rocsolver_dgetrf_npvt(rocblas_handle handle,
@@ -229,8 +316,15 @@ rocblas_status rocsolver_dgetrf_npvt(rocblas_handle handle,
                                      const rocblas_int lda,
                                      rocblas_int* info)
 {
-    rocblas_int* ipiv = nullptr;
-    return rocsolver::rocsolver_getrf_impl<double>(handle, m, n, A, lda, ipiv, info, false);
+    if(use_getrf_nopiv)
+    {
+        return rocsolver::rocsolver_getrf_nopiv_impl<double>(handle, m, n, A, lda, info);
+    }
+    else
+    {
+        rocblas_int* ipiv = nullptr;
+        return rocsolver::rocsolver_getrf_impl<double>(handle, m, n, A, lda, ipiv, info, false);
+    }
 }
 
 rocblas_status rocsolver_cgetrf_npvt(rocblas_handle handle,
@@ -240,9 +334,17 @@ rocblas_status rocsolver_cgetrf_npvt(rocblas_handle handle,
                                      const rocblas_int lda,
                                      rocblas_int* info)
 {
-    rocblas_int* ipiv = nullptr;
-    return rocsolver::rocsolver_getrf_impl<rocblas_float_complex>(handle, m, n, A, lda, ipiv, info,
-                                                                  false);
+    if(use_getrf_nopiv)
+    {
+        return rocsolver::rocsolver_getrf_nopiv_impl<rocblas_float_complex>(handle, m, n, A, lda,
+                                                                            info);
+    }
+    else
+    {
+        rocblas_int* ipiv = nullptr;
+        return rocsolver::rocsolver_getrf_impl<rocblas_float_complex>(handle, m, n, A, lda, ipiv,
+                                                                      info, false);
+    }
 }
 
 rocblas_status rocsolver_zgetrf_npvt(rocblas_handle handle,
@@ -252,9 +354,17 @@ rocblas_status rocsolver_zgetrf_npvt(rocblas_handle handle,
                                      const rocblas_int lda,
                                      rocblas_int* info)
 {
-    rocblas_int* ipiv = nullptr;
-    return rocsolver::rocsolver_getrf_impl<rocblas_double_complex>(handle, m, n, A, lda, ipiv, info,
-                                                                   false);
+    if(use_getrf_nopiv)
+    {
+        return rocsolver::rocsolver_getrf_nopiv_impl<rocblas_double_complex>(handle, m, n, A, lda,
+                                                                             info);
+    }
+    else
+    {
+        rocblas_int* ipiv = nullptr;
+        return rocsolver::rocsolver_getrf_impl<rocblas_double_complex>(handle, m, n, A, lda, ipiv,
+                                                                       info, false);
+    }
 }
 
 rocblas_status rocsolver_sgetrf_npvt_64(rocblas_handle handle,