diff --git a/clients/common/utility.cpp b/clients/common/utility.cpp index c4152eb..f84d795 100644 --- 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 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 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 index f464bd5..2c52943 100644 --- 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 buffer; + buffer.allocate_buffer(size); + if(arg.backend == backend::GPU) { x.move_to_device(); @@ -53,7 +56,7 @@ 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(); @@ -61,7 +64,7 @@ bool testing::test_dot_product(Arguments arg) 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 index 8f5702d..27c625c 100644 --- 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 D1(mat_A.get_m()); - linalg::vector D2(mat_A.get_m()); - mat_A.apply_ruiz_scaling(D1, D2, 30, 1e-03); - // Solution vector linalg::vector 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 index 5b005bf..000c06f 100644 --- 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 index c436a79..ffe8fcf 100644 --- 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 index 0487b9c..ab1edaf 100644 --- 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 index cc0867d..158703b 100644 --- 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 q; + dp_opt_buffer + 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 index a7db738..1814f7c 100644 --- 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 res; + dp_opt_buffer + 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 index efbaabd..cb008c6 100644 --- 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 z; + dp_opt_buffer + 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 new file mode 100644 index 0000000..ac7bcd5 --- /dev/null +++ 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 + +namespace linalg +{ + template + 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 index 6c361be..f3d6312 100644 --- a/library/include/linalg_math.h +++ b/library/include/linalg_math.h @@ -30,6 +30,7 @@ #include #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& x, const vector& y); + LINALGLIB_API double dot_product(const vector& x, + const vector& y, + dp_opt_buffer& 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& array); + LINALGLIB_API double norm_euclid(const vector& array, dp_opt_buffer& buffer); /** * @brief Computes the infinity (maximum absolute value) norm of a vector. diff --git a/library/src/CMakeLists.txt b/library/src/CMakeLists.txt index 4de2d65..4a47f4d 100644 --- 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 index 9c2db40..7a51f1d 100644 --- 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 -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( int, double, const double*, double, const double*, double, double*); template void linalg::cuda_axpbypgz(int, float, const float*, float, const float*, float, float*); -template double linalg::cuda_dot_product(const double*, const double*, int); -template float linalg::cuda_dot_product(const float*, const float*, int); +template double linalg::cuda_dot_product(const double*, const double*, double*, int); +template float linalg::cuda_dot_product(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 index bcb9f3d..237d1cd 100644 --- 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 void cuda_axpbypgz(int size, T alpha, const T* x, T beta, const T* y, T gamma, T* z); template - 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 index f641a9f..268c135 100644 --- 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 @@ -83,12 +85,15 @@ void linalg::device_axpbypgz(double alpha, } } -double linalg::device_dot_product(const vector& x, const vector& y) +double linalg::device_dot_product(const vector& x, + const vector& y, + dp_opt_buffer& 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 index 32c1965..f202578 100644 --- 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& y, double gamma, vector& z); - double device_dot_product(const vector& x, const vector& y); + double device_dot_product(const vector& x, + const vector& y, + dp_opt_buffer& buffer); } #endif diff --git a/library/src/backend/device/device_math.cpp b/library/src/backend/device/device_math.cpp index 99c0405..d11508e 100644 --- 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& array) +double linalg::device_norm_euclid(const vector& array, dp_opt_buffer& 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& array) diff --git a/library/src/backend/device/device_math.h b/library/src/backend/device/device_math.h index cc2a03c..080815e 100644 --- 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& array); + double device_norm_euclid(const vector& array, dp_opt_buffer& buffer); // Infinity norm double device_norm_inf(const vector& array); diff --git a/library/src/backend/host/host_axpy.cpp b/library/src/backend/host/host_axpy.cpp index a41a8df..7209399 100644 --- 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& x, const vector& y) +double linalg::host_dot_product(const vector& x, + const vector& y, + dp_opt_buffer& 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 index 5fd40d2..ad9d98f 100644 --- 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& y, double gamma, vector& z); - double host_dot_product(const vector& x, const vector& y); + double host_dot_product(const vector& x, + const vector& y, + dp_opt_buffer& buffer); } #endif diff --git a/library/src/backend/host/host_math.cpp b/library/src/backend/host/host_math.cpp index 0d602f9..02a031b 100644 --- a/library/src/backend/host/host_math.cpp +++ b/library/src/backend/host/host_math.cpp @@ -64,11 +64,11 @@ namespace linalg } } -double linalg::host_norm_euclid(const vector& array) +double linalg::host_norm_euclid(const vector& array, dp_opt_buffer& buffer) { ROUTINE_TRACE("linalg::host_norm_euclid"); - return std::sqrt(host_dot_product(array, array)); + return std::sqrt(host_dot_product(array, array, buffer)); } double linalg::host_norm_inf(const vector& array) diff --git a/library/src/backend/host/host_math.h b/library/src/backend/host/host_math.h index 302984f..02ad273 100644 --- a/library/src/backend/host/host_math.h +++ b/library/src/backend/host/host_math.h @@ -43,11 +43,12 @@ #include "host_ssor.h" #include "host_tridiagonal.h" +#include "linalg_buffers.h" #include "linalg_export.h" namespace linalg { - double host_norm_euclid(const vector& array); + double host_norm_euclid(const vector& array, dp_opt_buffer& buffer); double host_norm_inf(const vector& array); void host_jacobi_solve(const vector& rhs, const vector& diag, vector& x); diff --git a/library/src/iterative_solvers/krylov/bicgstab.cpp b/library/src/iterative_solvers/krylov/bicgstab.cpp index 7df3166..616a14d 100644 --- a/library/src/iterative_solvers/krylov/bicgstab.cpp +++ b/library/src/iterative_solvers/krylov/bicgstab.cpp @@ -59,6 +59,8 @@ void bicgstab_solver::build(const csr_matrix& A) t.resize(A.get_m()); z.resize(A.get_m()); q.resize(A.get_m()); + + buffer.allocate_buffer(A.get_m()); } int bicgstab_solver::solve_nonprecond(const csr_matrix& A, @@ -76,7 +78,7 @@ int bicgstab_solver::solve_nonprecond(const csr_matrix& A, // r0 = r r0.copy_from(r); - double rho = dot_product(r0, r); + double rho = dot_product(r0, r, buffer); // p = r p.copy_from(r); @@ -89,7 +91,7 @@ int bicgstab_solver::solve_nonprecond(const csr_matrix& A, // v = Ap A.multiply_by_vector(v, p); - double alpha = rho / dot_product(r0, v); + double alpha = rho / dot_product(r0, v, buffer); // r = r - alpha * v axpy(-1.0 * alpha, v, r); @@ -97,8 +99,8 @@ int bicgstab_solver::solve_nonprecond(const csr_matrix& A, // t = A * r A.multiply_by_vector(t, r); - double omega1 = dot_product(t, r); - double omega2 = dot_product(t, t); + double omega1 = dot_product(t, r, buffer); + double omega2 = dot_product(t, t, buffer); if(omega1 == 0.0 || omega2 == 0.0) { @@ -122,7 +124,7 @@ int bicgstab_solver::solve_nonprecond(const csr_matrix& A, } double rho_prev = rho; - rho = dot_product(r0, r); + rho = dot_product(r0, r, buffer); double beta = (rho / rho_prev) * (alpha / omega); // p = r + beta * (p - omega * v) @@ -158,7 +160,7 @@ int bicgstab_solver::solve_precond(const csr_matrix& A, // r0 = r r0.copy_from(r); - double rho = dot_product(r0, r); + double rho = dot_product(r0, r, buffer); // p = r p.copy_from(r); @@ -174,7 +176,7 @@ int bicgstab_solver::solve_precond(const csr_matrix& A, // q = A*z A.multiply_by_vector(q, z); - double alpha = rho / dot_product(r0, q); + double alpha = rho / dot_product(r0, q, buffer); // r = r - alpha * q axpy(-1.0 * alpha, q, r); @@ -185,8 +187,8 @@ int bicgstab_solver::solve_precond(const csr_matrix& A, // t = A * v A.multiply_by_vector(t, v); - double omega1 = dot_product(t, r); - double omega2 = dot_product(t, t); + double omega1 = dot_product(t, r, buffer); + double omega2 = dot_product(t, t, buffer); if(omega1 == 0.0 || omega2 == 0.0) { @@ -210,7 +212,7 @@ int bicgstab_solver::solve_precond(const csr_matrix& A, } double rho_prev = rho; - rho = dot_product(r0, r); + rho = dot_product(r0, r, buffer); double beta = (rho / rho_prev) * (alpha / omega); // p = r + beta * (p - omega * q) diff --git a/library/src/iterative_solvers/krylov/cg.cpp b/library/src/iterative_solvers/krylov/cg.cpp index 4a3418f..757e029 100644 --- a/library/src/iterative_solvers/krylov/cg.cpp +++ b/library/src/iterative_solvers/krylov/cg.cpp @@ -55,6 +55,8 @@ void cg_solver::build(const csr_matrix& A) p.resize(A.get_m()); z.resize(A.get_m()); res.resize(A.get_m()); + + buffer.allocate_buffer(A.get_m()); } int cg_solver::solve_nonprecond(const csr_matrix& A, @@ -77,7 +79,7 @@ int cg_solver::solve_nonprecond(const csr_matrix& A, // p = res p.copy_from(res); - gamma = dot_product(res, res); + gamma = dot_product(res, res, buffer); } auto t1 = std::chrono::high_resolution_clock::now(); @@ -94,12 +96,12 @@ int cg_solver::solve_nonprecond(const csr_matrix& A, // p = res p.copy_from(res); - gamma = dot_product(res, res); + gamma = dot_product(res, res, buffer); } // z = A * p and alpha = (r, r) / (A * p, p) A.multiply_by_vector(z, p); - double alpha = gamma / dot_product(z, p); + double alpha = gamma / dot_product(z, p, buffer); // update x = x + alpha * p axpy(alpha, p, x); @@ -116,7 +118,7 @@ int cg_solver::solve_nonprecond(const csr_matrix& A, // find beta double old_gamma = gamma; - gamma = dot_product(res, res); + gamma = dot_product(res, res, buffer); double beta = gamma / old_gamma; // update p = res + beta * p @@ -159,7 +161,7 @@ int cg_solver::solve_precond(const csr_matrix& A, // p = z p.copy_from(z); - gamma = dot_product(z, res); + gamma = dot_product(z, res, buffer); } auto t1 = std::chrono::high_resolution_clock::now(); @@ -179,12 +181,12 @@ int cg_solver::solve_precond(const csr_matrix& A, // p = z p.copy_from(z); - gamma = dot_product(z, res); + gamma = dot_product(z, res, buffer); } // z = A * p and alpha = (z, r) / (Ap, p) A.multiply_by_vector(z, p); - double alpha = gamma / dot_product(z, p); + double alpha = gamma / dot_product(z, p, buffer); // update x = x + alpha * p axpy(alpha, p, x); @@ -204,7 +206,7 @@ int cg_solver::solve_precond(const csr_matrix& A, // find beta double old_gamma = gamma; - gamma = dot_product(z, res); + gamma = dot_product(z, res, buffer); double beta = gamma / old_gamma; // update p = z + beta * p @@ -265,4 +267,4 @@ void cg_solver::move_to_device() bool cg_solver::is_on_host() const { return on_host; -} \ No newline at end of file +} diff --git a/library/src/iterative_solvers/krylov/gmres.cpp b/library/src/iterative_solvers/krylov/gmres.cpp index 504c555..ab105b7 100644 --- a/library/src/iterative_solvers/krylov/gmres.cpp +++ b/library/src/iterative_solvers/krylov/gmres.cpp @@ -33,7 +33,6 @@ #include #include - #include "../../trace.h" using namespace linalg; @@ -220,6 +219,8 @@ void gmres_solver::build(const csr_matrix& A, int restart) Q.resize(A.get_m() * (restart + 1)); c.resize(restart); s.resize(restart); + + buffer.allocate_buffer(A.get_m()); } int gmres_solver::solve_nonprecond(const csr_matrix& A, @@ -232,7 +233,7 @@ int gmres_solver::solve_nonprecond(const csr_matrix& A, // res = b - A * x compute_residual(A, x, b, res); - double res_norm = norm_euclid(res); + double res_norm = norm_euclid(res, buffer); double initial_res_norm = res_norm; // Check norm of residual against tolerance @@ -340,7 +341,7 @@ int gmres_solver::solve_nonprecond(const csr_matrix& A, // res = b - A * x compute_residual(A, x, b, res); - res_norm = norm_euclid(res); + res_norm = norm_euclid(res, buffer); // Check norm of residual against tolerance if(control.residual_converges(res_norm, initial_res_norm)) @@ -384,7 +385,7 @@ int gmres_solver::solve_precond(const csr_matrix& A, // z = (M^-1) * res precond->solve(res, z); - double res_norm = norm_euclid(z); + double res_norm = norm_euclid(z, buffer); double initial_res_norm = res_norm; // Check norm of residual against tolerance @@ -495,7 +496,7 @@ int gmres_solver::solve_precond(const csr_matrix& A, // z = (M^-1) * res precond->solve(res, z); - res_norm = norm_euclid(z); + res_norm = norm_euclid(z, buffer); // Check norm of residual against tolerance if(control.residual_converges(res_norm, initial_res_norm)) @@ -539,4 +540,4 @@ int gmres_solver::solve(const csr_matrix& A, { return solve_precond(A, x, b, precond, control); } -} \ No newline at end of file +} diff --git a/library/src/linalg_buffers.cpp b/library/src/linalg_buffers.cpp new file mode 100644 index 0000000..18c9a1d --- /dev/null +++ b/library/src/linalg_buffers.cpp @@ -0,0 +1,69 @@ +//******************************************************************************** +// +// 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. +// +//******************************************************************************** + +#include "../include/linalg_buffers.h" + +#include "backend/device/device_memory.h" + +using namespace linalg; + +template +dp_opt_buffer::dp_opt_buffer() + : data(nullptr) + , on_host(true) +{ +} + +template +dp_opt_buffer::~dp_opt_buffer() +{ + free_buffer(); +} + +template +void dp_opt_buffer::allocate_buffer(size_t size) +{ + device_allocate(&data, 256); +} + +template +void dp_opt_buffer::free_buffer() +{ + if(data) + { + device_free(data); + data = nullptr; + } +} + +template +T* dp_opt_buffer::get_buffer() +{ + return data; +} + +template class dp_opt_buffer; +template class dp_opt_buffer; diff --git a/library/src/linalg_math.cpp b/library/src/linalg_math.cpp index 0ff4101..3b05384 100644 --- a/library/src/linalg_math.cpp +++ b/library/src/linalg_math.cpp @@ -99,11 +99,14 @@ void linalg::transpose_matrix(const csr_matrix& A, csr_matrix& transposeA) } // Dot product -double linalg::dot_product(const vector& x, const vector& y) +double linalg::dot_product(const vector& x, + const vector& y, + dp_opt_buffer& buffer) { ROUTINE_TRACE("linalg::dot_product"); - return backend_dispatch("linalg::dot_product", host_dot_product, device_dot_product, x, y); + return backend_dispatch( + "linalg::dot_product", host_dot_product, device_dot_product, x, y, buffer); } // Compute residual @@ -127,11 +130,12 @@ void linalg::diagonal(const csr_matrix& A, vector& d) } // Euclidean norm -double linalg::norm_euclid(const vector& array) +double linalg::norm_euclid(const vector& array, dp_opt_buffer& buffer) { ROUTINE_TRACE("linalg::norm_euclid"); - return backend_dispatch("linalg::norm_euclid", host_norm_euclid, device_norm_euclid, array); + return backend_dispatch( + "linalg::norm_euclid", host_norm_euclid, device_norm_euclid, array, buffer); } // Infinity norm