Optimizing dot product and expanding testing

clients/common/utility.cpp

@@ -430,6 +430,8 @@ bool check_solution(const linalg::csr_matrix&     A,
                     double                        tol,
                     int                           norm_type)
 {
+    linalg::dp_opt_buffer<double> buffer;
+
     for(size_t i = 0; i < x.get_size(); i++)
     {
         if(std::isnan(x[i]) || std::isinf(x[i]))
@@ -448,7 +450,7 @@ bool check_solution(const linalg::csr_matrix&     A,
     }
     else
     {
-        initial_residual_norm = linalg::norm_euclid(initial_residual);
+        initial_residual_norm = linalg::norm_euclid(initial_residual, buffer);
     }
 
     linalg::vector<double> residual(A.get_m());
@@ -461,7 +463,7 @@ bool check_solution(const linalg::csr_matrix&     A,
     }
     else
     {
-        residual_norm = linalg::norm_euclid(residual);
+        residual_norm = linalg::norm_euclid(residual, buffer);
     }
 
     std::cout << "absolute residual: " << residual_norm

clients/testing/test_functions_dot_product.cpp

@@ -42,6 +42,9 @@ bool testing::test_dot_product(Arguments arg)
     x.fill(2.0);
     y.fill(3.0);
 
+    linalg::dp_opt_buffer<double> buffer;
+    buffer.allocate_buffer(size);
+
     if(arg.backend == backend::GPU)
     {
         x.move_to_device();
@@ -53,15 +56,15 @@ bool testing::test_dot_product(Arguments arg)
     // Warmup
     for(int i = 0; i < 4; i++)
     {
-        result = linalg::dot_product(x, y);
+        result = linalg::dot_product(x, y, buffer);
     }
     linalg::synchronize();
 
     // Timed solve
     auto t1 = std::chrono::high_resolution_clock::now();
     for(int i = 0; i < 100; i++)
     {
-        result = linalg::dot_product(x, y);
+        result = linalg::dot_product(x, y, buffer);
     }
     linalg::synchronize();
     auto t2 = std::chrono::high_resolution_clock::now();

clients/testing/test_functions_krylov.cpp

@@ -38,10 +38,6 @@ bool testing::test_krylov(krylov_solver solver_type, Arguments arg)
     linalg::csr_matrix mat_A;
     mat_A.read_mtx(arg.filename);
 
-    linalg::vector<double> D1(mat_A.get_m());
-    linalg::vector<double> D2(mat_A.get_m());
-    mat_A.apply_ruiz_scaling(D1, D2, 30, 1e-03);
-
     // Solution vector
     linalg::vector<double> vec_x(mat_A.get_m());
     vec_x.zeros();

clients/testing/tests/CMakeLists.txt

@@ -52,6 +52,7 @@ target_sources(test_main
     test_axpy.cpp
     test_axpby.cpp
     test_axpbypgz.cpp
+    test_cg_regression.cpp
     test_dot_product.cpp
     test_ruiz_scaling.cpp
     test_symmetric_ruiz_scaling.cpp

clients/testing/tests/test_BICGSTAB.yaml

@@ -1,6 +1,6 @@
 Tests:
   quick_ci:
-    precond: [none, jacobi, SOR]
+    precond: [jacobi, SOR]
     matrix_file: ["matrices/SPD/nos7/nos7.mtx"]
     max_iters: [400]
     backend: [CPU]

clients/testing/tests/test_CG.yaml

@@ -1,6 +1,6 @@
 Tests:
   quick_ci:
-    precond: [none, jacobi, SOR]
+    precond: [jacobi]
     matrix_file: ["matrices/SPD/nos7/nos7.mtx"]
     max_iters: [400]
     backend: [CPU]

clients/testing/tests/test_cg_regression.cpp

@@ -0,0 +1,32 @@
+#include <gtest/gtest.h>
+
+#include "linalg.h"
+
+TEST(CGRegression, RespectsMaxIterations)
+{
+    std::vector<int>    row_ptr = {0, 1, 2};
+    std::vector<int>    col_ind = {0, 1};
+    std::vector<double> values  = {4.0, 4.0};
+
+    linalg::csr_matrix A(row_ptr, col_ind, values, 2, 2, 2);
+
+    linalg::vector<double> x(2);
+    x.zeros();
+
+    std::vector<double>    b_data = {1.0, 2.0};
+    linalg::vector<double> b(b_data);
+
+    linalg::iter_control control;
+    control.max_iter = 1;
+    control.abs_tol  = 1e-14;
+    control.rel_tol  = 1e-14;
+
+    linalg::cg_solver solver;
+    solver.build(A);
+
+    int iterations = solver.solve(A, x, b, nullptr, control);
+
+    EXPECT_EQ(iterations, 1);
+    EXPECT_NEAR(x[0], 0.25, 1e-12);
+    EXPECT_NEAR(x[1], 0.5, 1e-12);
+}

clients/testing/tests/test_dot_product.yaml

@@ -8,9 +8,10 @@ Tests:
     backend: [CPU, GPU]
 
   medium:
-    m: [555, 678, 801, 978, 1024, 1436, 1867, 2048, 2345, 2567, 3001, 3456]
+    m: [555, 678, 801, 978, 1024, 1436, 1867, 2048, 2345, 2567, 3001, 3456, 4096, 8192]
     backend: [CPU, GPU]
 
   large:
-    m: [2048, 4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576]
+    m: [16384, 32768, 65536, 131072, 262144, 524288, 1048576, 2345918, 3456789, 4567890,
+        5678901, 6789012, 7890123, 8901234, 9012345]
     backend: [CPU, GPU]

library/include/iterative_solvers/krylov/bicgstab.h

@@ -27,6 +27,7 @@
 #ifndef BICGSTAB_H
 #define BICGSTAB_H
 
+#include "../../linalg_buffers.h"
 #include "../../linalg_export.h"
 
 #include "../iter_control.h"
@@ -216,6 +217,9 @@ namespace linalg
         /*! \brief Intermediate vector for preconditioning: \f$M^{-1} \mathbf{s}\f$. */
         vector<double> q;
 
+        dp_opt_buffer<double>
+            buffer; /*!< \brief Buffer for optimization data used in dot products. */
+
         /*! \brief Number of iterations after which the solver should restart.
      * A value of 0 or a very large number typically means no restart.
      * Restarts can help to avoid potential breakdowns or loss of orthogonality.

library/include/iterative_solvers/krylov/cg.h

@@ -27,6 +27,7 @@
 #ifndef CG_H
 #define CG_H
 
+#include "../../linalg_buffers.h"
 #include "../../linalg_export.h"
 
 #include "../iter_control.h"
@@ -200,6 +201,9 @@ namespace linalg
         /*! \brief Residual vector in the CG algorithm. */
         vector<double> res;
 
+        dp_opt_buffer<double>
+            buffer; /*!< \brief Buffer for optimization data used in dot products. */
+
         /*! \brief Number of iterations after which the solver should restart.
      * A value of 0 or a very large number typically means no restart.
      * For CG, restarts are usually not needed for exact arithmetic but can
@@ -252,10 +256,7 @@ namespace linalg
      * \param b The right-hand side vector.
      * \param control An `iter_control` object that manages the iteration process,
      * including convergence tolerance and maximum iterations.
-     * \return An integer status code:
-     * - `0` if the solver converged successfully within the specified tolerance.
-     * - `1` if the maximum number of iterations was reached without convergence.
-     * - Negative values might indicate issues like a non-positive definite matrix or division by zero.
+     * \return The number of iterations performed by the solver.
      */
         int solve_nonprecond(const csr_matrix&     A,
                              vector<double>&       x,
@@ -278,10 +279,7 @@ namespace linalg
      * of the preconditioner should have been called previously for matrix `A`.
      * \param control An `iter_control` object that manages the iteration process,
      * including convergence tolerance and maximum iterations.
-     * \return An integer status code:
-     * - `0` if the solver converged successfully within the specified tolerance.
-     * - `1` if the maximum number of iterations was reached without convergence.
-     * - Negative values might indicate issues with the matrix, preconditioner, or numerical stability.
+     * \return The number of iterations performed by the solver.
      */
         int solve_precond(const csr_matrix&     A,
                           vector<double>&       x,
@@ -302,7 +300,7 @@ namespace linalg
      * \param precond A pointer to a `preconditioner` object to be used. If `nullptr`, no preconditioning is applied.
      * \param control An `iter_control` object that manages the iteration process,
      * including convergence tolerance and maximum iterations.
-     * \return An integer status code, consistent with `solve_nonprecond` or `solve_precond`.
+     * \return The number of iterations performed by the solver.
      */
         int solve(const csr_matrix&     A,
                   vector<double>&       x,

library/include/iterative_solvers/krylov/gmres.h

@@ -27,6 +27,7 @@
 #ifndef GMRES_H
 #define GMRES_H
 
+#include "../../linalg_buffers.h"
 #include "../../linalg_export.h"
 
 #include "../iter_control.h"
@@ -244,6 +245,9 @@ namespace linalg
         /*! \brief Intermediate vector for preconditioning or other operations. */
         vector<double> z;
 
+        dp_opt_buffer<double>
+            buffer; /*!< \brief Buffer for device operations, used to optimize memory usage and performance on GPU backends. */
+
         /*! \brief The restart parameter `m` for GMRES(m).
      * \details This defines the maximum dimension of the krylov subspace before restarting.
      * A smaller `restart` value means less memory usage but potentially more restarts.

library/include/linalg_buffers.h

@@ -0,0 +1,52 @@
+//********************************************************************************
+//
+// MIT License
+//
+// Copyright(c) 2026 James Sandham
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this softwareand associated documentation files(the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions :
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//
+//********************************************************************************
+
+#ifndef OPTIMIZATION_BUFFERS_H
+#define OPTIMIZATION_BUFFERS_H
+
+#include <cstddef>
+
+namespace linalg
+{
+    template <typename T>
+    class dp_opt_buffer
+    {
+    private:
+        T*   data;
+        bool on_host;
+
+    public:
+        dp_opt_buffer();
+        ~dp_opt_buffer();
+
+        void allocate_buffer(size_t size);
+        void free_buffer();
+
+        T* get_buffer();
+    };
+}
+
+#endif // OPTIMIZATION_BUFFERS_H

library/include/linalg_enums.h

@@ -63,6 +63,7 @@ namespace linalg
     enum class csrmv_algorithm
     {
         default_algorithm, /*!< Default algorithm */
+        lrb,
         merge_path,
         rowsplit,
         nnzsplit

library/include/linalg_math.h

@@ -30,6 +30,7 @@
 #include <string>
 
 #include "csr_matrix.h"
+#include "linalg_buffers.h"
 #include "linalg_enums.h"
 #include "linalg_export.h"
 #include "linalg_types.h"
@@ -126,9 +127,12 @@ namespace linalg
      *
      * @param x The first input vector.
      * @param y The second input vector.
+     * @param buffer Buffer storing optimization data.
      * @return The double-precision floating-point result of the dot product.
      */
-    LINALGLIB_API double dot_product(const vector<double>& x, const vector<double>& y);
+    LINALGLIB_API double dot_product(const vector<double>&  x,
+                                     const vector<double>&  y,
+                                     dp_opt_buffer<double>& buffer);
 
     /**
      * @brief Computes the residual vector for a linear system: \f$res = b - A \cdot x\f$.
@@ -159,9 +163,10 @@ namespace linalg
      * The Euclidean norm is calculated as \f$\sqrt{\sum_{i=0}^{n-1} |array_i|^2}\f$.
      *
      * @param array The input vector.
+     * @param buffer Buffer storing optimization data.
      * @return The double-precision floating-point value of the Euclidean norm.
      */
-    LINALGLIB_API double norm_euclid(const vector<double>& array);
+    LINALGLIB_API double norm_euclid(const vector<double>& array, dp_opt_buffer<double>& buffer);
 
     /**
      * @brief Computes the infinity (maximum absolute value) norm of a vector.

library/src/CMakeLists.txt

@@ -29,6 +29,7 @@ target_sources(linalglib
     linalg_primitives.cpp
     linalg_math.cpp
     linalg_memory.cpp
+    linalg_buffers.cpp
     vector.cpp
     csr_matrix.cpp
     perfetto_trace.cpp

library/src/backend/device/cuda/cuda_axpy.cu

@@ -68,21 +68,17 @@ void linalg::cuda_axpbypgz(int size, T alpha, const T* x, T beta, const T* y, T
 // dot product z = x*y
 //-------------------------------------------------------------------------------
 template <typename T>
-T linalg::cuda_dot_product(const T* x, const T* y, int size)
+T linalg::cuda_dot_product(const T* x, const T* y, T* buffer, int size)
 {
     ROUTINE_TRACE("linalg::cuda_dot_product_impl");
-    T* workspace = nullptr;
-    CHECK_CUDA(cudaMalloc((void**)&workspace, sizeof(T) * 256));
-
-    dot_product_kernel_part1<256><<<256, 256>>>(size, x, y, workspace);
+    dot_product_kernel_part1<256><<<256, 256>>>(size, x, y, buffer);
     CHECK_CUDA_LAUNCH_ERROR();
 
-    dot_product_kernel_part2<256><<<1, 256>>>(workspace);
+    dot_product_kernel_part2<256><<<1, 256>>>(buffer);
     CHECK_CUDA_LAUNCH_ERROR();
 
     T result;
-    CHECK_CUDA(cudaMemcpy(&result, workspace, sizeof(T), cudaMemcpyDeviceToHost));
-    CHECK_CUDA(cudaFree(workspace));
+    CHECK_CUDA(cudaMemcpy(&result, buffer, sizeof(T), cudaMemcpyDeviceToHost));
 
     return result;
 }
@@ -95,5 +91,5 @@ template void linalg::cuda_axpbypgz<double>(
     int, double, const double*, double, const double*, double, double*);
 template void
     linalg::cuda_axpbypgz<float>(int, float, const float*, float, const float*, float, float*);
-template double linalg::cuda_dot_product<double>(const double*, const double*, int);
-template float  linalg::cuda_dot_product<float>(const float*, const float*, int);
+template double linalg::cuda_dot_product<double>(const double*, const double*, double*, int);
+template float  linalg::cuda_dot_product<float>(const float*, const float*, float*, int);

library/src/backend/device/cuda/cuda_axpy.h

@@ -35,7 +35,7 @@ namespace linalg
     template <typename T>
     void cuda_axpbypgz(int size, T alpha, const T* x, T beta, const T* y, T gamma, T* z);
     template <typename T>
-    T cuda_dot_product(const T* x, const T* y, int size);
+    T cuda_dot_product(const T* x, const T* y, T* buffer, int size);
 }
 
 #endif

library/src/backend/device/device_axpy.cpp

@@ -24,6 +24,8 @@
 //
 //********************************************************************************
 
+// #include "../../../include/linalg_buffers.h"
+
 #include "device_axpy.h"
 
 #include <iostream>
@@ -83,12 +85,15 @@ void linalg::device_axpbypgz(double                alpha,
     }
 }
 
-double linalg::device_dot_product(const vector<double>& x, const vector<double>& y)
+double linalg::device_dot_product(const vector<double>&  x,
+                                  const vector<double>&  y,
+                                  dp_opt_buffer<double>& buffer)
 {
     ROUTINE_TRACE("linalg::device_dot_product");
     if constexpr(is_cuda_available())
     {
-        return RETURN_CALL_CUDA(cuda_dot_product(x.get_vec(), y.get_vec(), x.get_size()));
+        return RETURN_CALL_CUDA(
+            cuda_dot_product(x.get_vec(), y.get_vec(), buffer.get_buffer(), x.get_size()));
     }
     std::cout << "Error: Not device backend available for the function " << __func__ << std::endl;
     return 0.0;