jsandham · jsandham · Aug 5, 2026 · Aug 5, 2026
diff --git a/clients/common/utility.cpp b/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

diff --git a/clients/testing/test_functions_dot_product.cpp b/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();

diff --git a/clients/testing/test_functions_krylov.cpp b/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();

diff --git a/clients/testing/tests/test_BICGSTAB.yaml b/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]

diff --git a/clients/testing/tests/test_CG.yaml b/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]

diff --git a/clients/testing/tests/test_dot_product.yaml b/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]
diff --git a/library/include/iterative_solvers/krylov/bicgstab.h b/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.

diff --git a/library/include/iterative_solvers/krylov/cg.h b/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

diff --git a/library/include/iterative_solvers/krylov/gmres.h b/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.

diff --git a/library/include/linalg_buffers.h b/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
diff --git a/library/include/linalg_math.h b/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.

diff --git a/library/src/CMakeLists.txt b/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

diff --git a/library/src/backend/device/cuda/cuda_axpy.cu b/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);
diff --git a/library/src/backend/device/cuda/cuda_axpy.h b/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
diff --git a/library/src/backend/device/device_axpy.cpp b/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;

diff --git a/library/src/backend/device/device_axpy.h b/library/src/backend/device/device_axpy.h
@@ -26,6 +26,7 @@
 #ifndef DEVICE_AXPY_H
 #define DEVICE_AXPY_H
 
+#include "linalg_buffers.h"
 #include "vector.h"
 
 namespace linalg
@@ -38,7 +39,9 @@ namespace linalg
                            const vector<double>& y,
                            double                gamma,
                            vector<double>&       z);
-    double device_dot_product(const vector<double>& x, const vector<double>& y);
+    double device_dot_product(const vector<double>&  x,
+                              const vector<double>&  y,
+                              dp_opt_buffer<double>& buffer);
 }
 
 #endif
diff --git a/library/src/backend/device/device_math.cpp b/library/src/backend/device/device_math.cpp
@@ -36,11 +36,11 @@
 #include "cuda/cuda_math.h"
 #endif
 
-double linalg::device_norm_euclid(const vector<double>& array)
+double linalg::device_norm_euclid(const vector<double>& array, dp_opt_buffer<double>& buffer)
 {
     ROUTINE_TRACE("linalg::device_norm_euclid");
 
-    return std::sqrt(device_dot_product(array, array));
+    return std::sqrt(device_dot_product(array, array, buffer));
 }
 
 double linalg::device_norm_inf(const vector<double>& array)

diff --git a/library/src/backend/device/device_math.h b/library/src/backend/device/device_math.h
@@ -43,6 +43,7 @@
 #include "device_ssor.h"
 #include "device_tridiagonal.h"
 
+#include "linalg_buffers.h"
 #include "linalg_export.h"
 
 /*! \file
@@ -51,7 +52,7 @@
 namespace linalg
 {
     // Euclidean norm
-    double device_norm_euclid(const vector<double>& array);
+    double device_norm_euclid(const vector<double>& array, dp_opt_buffer<double>& buffer);
 
     // Infinity norm
     double device_norm_inf(const vector<double>& array);

diff --git a/library/src/backend/host/host_axpy.cpp b/library/src/backend/host/host_axpy.cpp
@@ -120,7 +120,9 @@ void linalg::host_axpbypgz(double                alpha,
     host_axpbypgz_impl(x.get_size(), alpha, x.get_vec(), beta, y.get_vec(), gamma, z.get_vec());
 }
 
-double linalg::host_dot_product(const vector<double>& x, const vector<double>& y)
+double linalg::host_dot_product(const vector<double>&  x,
+                                const vector<double>&  y,
+                                dp_opt_buffer<double>& buffer)
 {
     ROUTINE_TRACE("linalg::host_dot_product");
     return host_dot_product_impl(x.get_vec(), y.get_vec(), x.get_size());

diff --git a/library/src/backend/host/host_axpy.h b/library/src/backend/host/host_axpy.h
@@ -27,6 +27,7 @@
 #ifndef HOST_AXPY_H
 #define HOST_AXPY_H
 
+#include "linalg_buffers.h"
 #include "vector.h"
 
 namespace linalg
@@ -39,7 +40,9 @@ namespace linalg
                          const vector<double>& y,
                          double                gamma,
                          vector<double>&       z);
-    double host_dot_product(const vector<double>& x, const vector<double>& y);
+    double host_dot_product(const vector<double>&  x,
+                            const vector<double>&  y,
+                            dp_opt_buffer<double>& buffer);
 }
 
 #endif