From 7c5fb562e64bc4a0ec89415d5b4682ac897f93a6 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Tue, 22 Jul 2025 13:50:56 +0200 Subject: [PATCH 1/7] feat: make Lanczos part of the repo --- CMakeLists.txt | 8 - include/lanczos/LanczosAlgorithm.cu | 233 ++++++++++++++++++ include/lanczos/LanczosAlgorithm.h | 45 ++++ include/lanczos/utils/MatrixDot.h | 29 +++ include/lanczos/utils/cublasDebug.h | 40 +++ include/lanczos/utils/cuda_lib_defines.h | 40 +++ include/lanczos/utils/debugTools.h | 47 ++++ include/lanczos/utils/defines.h | 14 ++ include/lanczos/utils/device_blas.h | 59 +++++ include/lanczos/utils/device_container.h | 51 ++++ .../lanczos/utils/lapack_and_blas_defines.h | 24 ++ 11 files changed, 582 insertions(+), 8 deletions(-) create mode 100644 include/lanczos/LanczosAlgorithm.cu create mode 100644 include/lanczos/LanczosAlgorithm.h create mode 100644 include/lanczos/utils/MatrixDot.h create mode 100644 include/lanczos/utils/cublasDebug.h create mode 100644 include/lanczos/utils/cuda_lib_defines.h create mode 100644 include/lanczos/utils/debugTools.h create mode 100644 include/lanczos/utils/defines.h create mode 100644 include/lanczos/utils/device_blas.h create mode 100644 include/lanczos/utils/device_container.h create mode 100644 include/lanczos/utils/lapack_and_blas_defines.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 1bbcb019..71ed791f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -50,14 +50,6 @@ FetchContent_Declare( EXCLUDE_FROM_ALL ) FetchContent_MakeAvailable(uammd) -FetchContent_Declare( - lanczos - GIT_REPOSITORY https://github.com/RaulPPelaez/LanczosAlgorithm - GIT_TAG v1.0.1 - EXCLUDE_FROM_ALL -) -FetchContent_MakeAvailable(lanczos) -include_directories(${lanczos_SOURCE_DIR}/include) set(BLA_VENDOR Generic) find_package(BLAS REQUIRED) diff --git a/include/lanczos/LanczosAlgorithm.cu b/include/lanczos/LanczosAlgorithm.cu new file mode 100644 index 00000000..990a8a5b --- /dev/null +++ b/include/lanczos/LanczosAlgorithm.cu @@ -0,0 +1,233 @@ +/*Raul P. Pelaez 2017-2022. Lanczos algorithm + +References: + [1] Krylov subspace methods for computing hydrodynamic interactions in +Brownian dynamics simulations. -http://dx.doi.org/10.1063/1.4742347 + +*/ +#include "LanczosAlgorithm.h" +#include "utils/device_blas.h" +#include "utils/device_container.h" +#include "utils/lapack_and_blas_defines.h" +#include +#include +#include +#ifdef CUDA_ENABLED +#include "utils/debugTools.h" +#endif + +#include + +namespace lanczos { + +namespace detail { + +/*See algorithm I in [1]*/ +class KrylovSubspace { + Blas blas; + device_container w; // size N, v in each iteration + device_container + V; // size Nxmax_iter; Krylov subspace base transformation matrix + // Mobility Matrix in the Krylov subspace + std::vector + P; // Transformation Matrix to diagonalize H, max_iter x max_iter + /*upper diagonal and diagonal of H*/ + std::vector hdiag, hsup, htemp; + device_container htempGPU; + int N; + int subSpaceDimension; + + real normz; + + real computeNorm(const real *v, int numberElements) { + real norm2; + blas.nrm2(numberElements, v, 1, &norm2); + return norm2; + } + + void diagonalizeSubSpace() { + int size = getSubSpaceSize(); + /**************LAPACKE********************/ + /*The tridiagonal matrix is stored only with its diagonal and subdiagonal*/ + /*Store both in a temporal array*/ + for (int i = 0; i < size; i++) { + htemp[i] = hdiag[i]; + htemp[i + size] = hsup[i]; + } + /*P = eigenvectors must be filled with zeros, I do not know why*/ + real *h_P = P.data(); + memset(h_P, 0, size * size * sizeof(real)); + /*Compute eigenvalues and eigenvectors of a triangular symmetric matrix*/ + auto info = LAPACKE_steqr(LAPACK_COL_MAJOR, 'I', size, &htemp[0], + &htemp[0] + size, h_P, size); + if (info != 0) { + throw std::runtime_error("[Lanczos] Could not diagonalize tridiagonal " + "krylov matrix, steqr failed with code " + + std::to_string(info)); + } + } + + real *computeSquareRoot() { + int size = getSubSpaceSize(); + diagonalizeSubSpace(); + /***Hdiag_temp = Hdiag·P·e1****/ + for (int j = 0; j < size; j++) { + htemp[j] = sqrt(htemp[j]) * P[size * j]; + } + /***** Htemp = H^1/2·e1 = Pt· hdiag_temp ****/ + /*Compute with blas*/ + real *h_P = P.data(); + real alpha = 1.0; + real beta = 0.0; + cblas_gemv(CblasColMajor, CblasNoTrans, size, size, alpha, h_P, size, + &htemp[0], 1, beta, &htemp[0] + size, 1); + detail::device_copy(htemp.begin() + size, htemp.begin() + 2 * size, + htempGPU.begin()); + return detail::getRawPointer(htempGPU); + } + + real *getTransformationMatrix() { return detail::getRawPointer(V); } + + void resize(int subSpaceSize) { +#ifdef CUDA_ENABLED + CudaSafeCall(cudaDeviceSynchronize()); +#endif + try { + w.resize((N + 1), real()); + V.resize(N * subSpaceSize, 0); + P.resize(subSpaceSize * subSpaceSize, 0); + hdiag.resize(subSpaceSize + 1, 0); + hsup.resize(subSpaceSize + 1, 0); + htemp.resize(2 * subSpaceSize, 0); + htempGPU.resize(2 * subSpaceSize, 0); + } catch (...) { + throw std::runtime_error("[KrylovSubspace] Could not allocate memory"); + } + } + +public: + KrylovSubspace(int N) : subSpaceDimension(0), N(N) { this->resize(1); } + + /************v[0] = z/||z||_2*****/ + void setFirstBasisVector(const real *z) { + /*1/norm(z)*/ + real *Vm = getTransformationMatrix(); + this->normz = computeNorm(z, N); + /*v[0] = v[0]*1/norm(z)*/ + real invz2 = 1.0 / normz; + detail::device_copy(z, z + N, Vm); + blas.scal(N, &invz2, Vm, 1); + } + + void nextIteration(lanczos::Dot &dot) { + int i = subSpaceDimension; + this->subSpaceDimension++; + resize(subSpaceDimension + 1); + auto d_V = detail::getRawPointer(V); + auto d_w = detail::getRawPointer(w); + /*w = D·vi*/ + dot(d_V + N * i, d_w); + if (i > 0) { + /*w = w-h[i-1][i]·vi*/ + real alpha = -hsup[i - 1]; + blas.axpy(N, &alpha, d_V + N * (i - 1), 1, d_w, 1); + } + /*h[i][i] = dot(w, vi)*/ + blas.dot(N, d_w, 1, d_V + N * i, 1, &(hdiag[i])); + /*w = w-h[i][i]·vi*/ + real alpha = -hdiag[i]; + blas.axpy(N, &alpha, d_V + N * i, 1, d_w, 1); + /*h[i+1][i] = h[i][i+1] = norm(w)*/ + blas.nrm2(N, (real *)d_w, 1, &(hsup[i])); + /*v_(i+1) = w·1/ norm(w)*/ + real tol = 1e-3 * hdiag[i] / normz; + if (hsup[i] < tol) + hsup[i] = real(0.0); + if (hsup[i] > real(0.0)) { + real invw2 = 1.0 / hsup[i]; + blas.scal(N, &invw2, d_w, 1); + } else { /*If norm(w) = 0 that means all elements of w are zero, so set w = + e1*/ + detail::device_fill(w.begin(), w.end(), real()); + w[0] = 1; + } + detail::device_copy(w.begin(), w.begin() + N, d_V + N * (i + 1)); + } + + int getSubSpaceSize() { return subSpaceDimension; } + + // Computes the current result guess sqrt(M)·v, stores in BdW + void computeCurrentResultEstimation(real *BdW) { + int m = getSubSpaceSize(); + /**** y = ||z||_2 * Vm · H^1/2 · e_1 *****/ + /**** H^1/2·e1 = Pt· first_column_of(sqrt(Hdiag)·P) ******/ + real *HhalfDotE1 = computeSquareRoot(); + /*y = ||z||_2 * Vm · H^1/2 · e1 = Vm · (z2·hdiag_temp)*/ + real *Vm = getTransformationMatrix(); + real beta = 0.0; + blas.gemv(N, m, &this->normz, Vm, N, HhalfDotE1, 1, &beta, BdW, 1); + } +}; +} // namespace detail + +Solver::Solver() : check_convergence_steps(3) {} + +int Solver::run(lanczos::Dot &dot, real *Bz, const real *z, real tolerance, + int N) { + oldBz.resize((N + 1), real()); + /*Lanczos iterations for Krylov decomposition*/ + detail::KrylovSubspace solver(N); + solver.setFirstBasisVector(z); + const int checkConvergenceSteps = + std::min(check_convergence_steps, iterationHardLimit - 2); + for (int i = 0; i < iterationHardLimit; i++) { + solver.nextIteration(dot); + if (i >= checkConvergenceSteps) { + solver.computeCurrentResultEstimation(Bz); + if (i > 0) { + auto currentResidual = computeError(Bz, N); + if (currentResidual <= tolerance) { + registerRequiredStepsForConverge(i); + return i; + } + } + // Store current estimation + detail::device_copy(Bz, Bz + N, oldBz.begin()); + } + } + throw std::runtime_error("[Lanczos] Could not converge"); +} + +real Solver::computeError(real *Bz, int N) { + /*Compute error as in eq 27 in [1] + Error = ||Bz_i - Bz_{i-1}||_2 / ||Bz_{i-1}||_2 + */ + real normResult_prev; + real *d_oldBz = detail::getRawPointer(oldBz); + blas.nrm2(N, d_oldBz, 1, &normResult_prev); + /*oldBz = Bz-oldBz*/ + real a = -1.0; + blas.axpy(N, &a, Bz, 1, d_oldBz, 1); + /*yy = ||Bz_i - Bz_{i-1}||_2*/ + real yy; + blas.nrm2(N, d_oldBz, 1, &yy); + // eq. 27 in [1] + real Error = abs(yy / normResult_prev); + if (std::isnan(Error)) { + throw std::runtime_error( + "[Lanczos] Unknown error (found NaN in result guess)"); + } + return Error; +} + +void Solver::registerRequiredStepsForConverge(int steps_needed) { + if (steps_needed - 2 > check_convergence_steps) { + check_convergence_steps += 1; + } + // Or check more often if I performed too many iterations + else { + check_convergence_steps = std::max(1, check_convergence_steps - 2); + } +} + +} // namespace lanczos diff --git a/include/lanczos/LanczosAlgorithm.h b/include/lanczos/LanczosAlgorithm.h new file mode 100644 index 00000000..fc819b00 --- /dev/null +++ b/include/lanczos/LanczosAlgorithm.h @@ -0,0 +1,45 @@ + +/*Raul P. Pelaez 2022. Lanczos Algotihm, + Computes the matrix-vector product sqrt(M)·v using a recursive algorithm. + For that, it requires a functor in which the () operator takes an output real* array and an input real* (both device memory) as: + inline void operator()(real* in_v, real * out_Mv); + This function must fill "out" with the result of performing the M·v dot product- > out = M·a_v. + If M has size NxN and the cost of the dot product is O(M). The total cost of the algorithm is O(m·M). Where m << N. + If M·v performs a dense M-V product, the cost of the algorithm would be O(m·N^2). + References: + [1] Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations + J. Chem. Phys. 137, 064106 (2012); doi: 10.1063/1.4742347 +Some notes: + + From what I have seen, this algorithm converges to an error of ~1e-3 in a few steps (<5) and from that point a lot of iterations are needed to lower the error. + It usually achieves machine precision in under 50 iterations. + + If the matrix does not have a sqrt (not positive definite, not symmetric...) it will usually be reflected as a nan in the current error estimation. An exception will be thrown in this case. +*/ + +#pragma once + +#include"utils/defines.h" +#include +#include"utils/device_container.h" +#include "utils/device_blas.h" +namespace lanczos{ + using Dot = std::function; + struct Solver{ + Solver(); + + int run(Dot &dot, real *Bv, const real* v, real tolerance, int N); + + void setIterationHardLimit(int newLimit){this->iterationHardLimit = newLimit;} + + private: + real computeError(real* Bz, int N); + void registerRequiredStepsForConverge(int steps_needed); + + Blas blas; + device_container oldBz; + int check_convergence_steps; + int iterationHardLimit = 200; + }; +} +#include"LanczosAlgorithm.cu" diff --git a/include/lanczos/utils/MatrixDot.h b/include/lanczos/utils/MatrixDot.h new file mode 100644 index 00000000..a76da941 --- /dev/null +++ b/include/lanczos/utils/MatrixDot.h @@ -0,0 +1,29 @@ +#ifndef LANCZOS_MATRIX_DOT_H +#define LANCZOS_MATRIX_DOT_H +#include "defines.h" +#include +namespace lanczos{ + + struct MatrixDot{ + void setSize(int newsize){this->m_size = newsize;} + //virtual void dot(real* v, real*Mv) = 0; + virtual void operator()(real* v, real*Mv) = 0; + protected: + int m_size; + }; + + //Transforms any callable into a MatrixDot valid to use with Lanczos + template + struct MatrixDotAdaptor: public lanczos::MatrixDot{ + Foo& foo; + MatrixDotAdaptor(Foo &&foo):foo(foo){} + void operator()(real* v, real* Mv) override{foo(v,Mv);} + }; + + template + auto createMatrixDotAdaptor(Foo &&foo){ + return MatrixDotAdaptor(foo); + } + +} +#endif diff --git a/include/lanczos/utils/cublasDebug.h b/include/lanczos/utils/cublasDebug.h new file mode 100644 index 00000000..c28cac62 --- /dev/null +++ b/include/lanczos/utils/cublasDebug.h @@ -0,0 +1,40 @@ +#ifndef CUBLAS_DEBUG_H +#define CUBLAS_DEBUG_H +#ifdef CUDA_ERROR_CHECK +#define CUBLAS_ERROR_CHECK +#endif + +#include +#define CublasSafeCall(err) __cublasSafeCall(err, __FILE__, __LINE__) + + + +const char* cublasGetErrorString(cublasStatus_t status){ + switch(status){ + case CUBLAS_STATUS_SUCCESS: return "CUBLAS_STATUS_SUCCESS"; + case CUBLAS_STATUS_NOT_INITIALIZED: return "CUBLAS_STATUS_NOT_INITIALIZED"; + case CUBLAS_STATUS_ALLOC_FAILED: return "CUBLAS_STATUS_ALLOC_FAILED"; + case CUBLAS_STATUS_INVALID_VALUE: return "CUBLAS_STATUS_INVALID_VALUE"; + case CUBLAS_STATUS_ARCH_MISMATCH: return "CUBLAS_STATUS_ARCH_MISMATCH"; + case CUBLAS_STATUS_MAPPING_ERROR: return "CUBLAS_STATUS_MAPPING_ERROR"; + case CUBLAS_STATUS_EXECUTION_FAILED: return "CUBLAS_STATUS_EXECUTION_FAILED"; + case CUBLAS_STATUS_INTERNAL_ERROR: return "CUBLAS_STATUS_INTERNAL_ERROR"; + default: return "Cublas Unknown error"; + } +} + +inline void __cublasSafeCall( cublasStatus_t err, const char *file, const int line ) +{ + #ifdef CUBLAS_ERROR_CHECK + if ( CUBLAS_STATUS_SUCCESS != err ) + { + fprintf( stderr, "cublasSafeCall() failed at %s:%i : %s - code: %i\n", + file, line, cublasGetErrorString( err ), err); + exit( -1 ); + } + #endif + + return; +} + +#endif diff --git a/include/lanczos/utils/cuda_lib_defines.h b/include/lanczos/utils/cuda_lib_defines.h new file mode 100644 index 00000000..80d23263 --- /dev/null +++ b/include/lanczos/utils/cuda_lib_defines.h @@ -0,0 +1,40 @@ +/*Raul P. Pelaez 2016-2022. Precision agnostic cublas/cusolver function defines*/ +#ifndef CUDA_LIB_DEFINES_H +#define CUDA_LIB_DEFINES_H +#include + +#if defined SINGLE_PRECISION +#define cusolverDnpotrf cusolverDnSpotrf +#define cusolverDnpotrf_bufferSize cusolverDnSpotrf_bufferSize +#define cublastrmv cublasStrmv +#define cublassymv cublasSsymv +#define cublasgemv cublasSgemv +#define cublasnrm2 cublasSnrm2 +#define cublasscal cublasSscal +#define cublasaxpy cublasSaxpy +#define cublasdot cublasSdot +#define cusolverDnsyevd cusolverDnSsyevd +#define cusolverDnsyevd_bufferSize cusolverDnSsyevd_bufferSize +#define cusolverDngesvd_bufferSize cusolverDnSgesvd_bufferSize +#define cusolverDngesvd cusolverDnSgesvd +#define cublasgemm cublasSgemm +#else +#define cusolverDnpotrf cusolverDnDpotrf +#define cusolverDnpotrf_bufferSize cusolverDnDpotrf_bufferSize +#define cublastrmv cublasDtrmv +#define curandGenerateNormal curandGenerateNormalDouble +#define cublassymv cublasDsymv +#define cublasgemv cublasDgemv +#define cublasnrm2 cublasDnrm2 +#define cublasscal cublasDscal +#define cublasaxpy cublasDaxpy +#define cublasdot cublasDdot +#define cusolverDnsyevd cusolverDnDsyevd +#define cusolverDnsyevd_bufferSize cusolverDnDsyevd_bufferSize +#define cusolverDngesvd_bufferSize cusolverDnDgesvd_bufferSize +#define cusolverDngesvd cusolverDnDgesvd +#define cublasgemm cublasDgemm +#endif + + +#endif diff --git a/include/lanczos/utils/debugTools.h b/include/lanczos/utils/debugTools.h new file mode 100644 index 00000000..0adf9265 --- /dev/null +++ b/include/lanczos/utils/debugTools.h @@ -0,0 +1,47 @@ +/*Raul P. Pelaez 2019-2022. Some utilities for debugging GPU code + */ +#ifndef DEBUGTOOLS_H +#define DEBUGTOOLS_H + +#define CUDA_ERROR_CHECK + +#ifdef LANCZOS_DEBUG +#define CUDA_ERROR_CHECK_SYNC +#endif +#define CudaSafeCall(err) __cudaSafeCall(err, __FILE__, __LINE__) +#define CudaCheckError() __cudaCheckError(__FILE__, __LINE__) + +#include +#include +#include + +inline void __cudaSafeCall(cudaError err, const char *file, const int line){ + #ifdef CUDA_ERROR_CHECK + if (cudaSuccess != err){ + cudaGetLastError(); //Reset CUDA error status + throw std::runtime_error("CudaSafeCall() failed at "+ + std::string(file) + ":" + std::to_string(line)+ + " with error " + std::to_string(err)); + } + #endif +} + +inline void __cudaCheckError(const char *file, const int line){ + cudaError err; +#ifdef CUDA_ERROR_CHECK_SYNC + err = cudaDeviceSynchronize(); + if(cudaSuccess != err){ + throw std::runtime_error("CudaCheckError() with sync failed at "+ + std::string(file) + ":" + std::to_string(line)+ + " with error " + std::to_string(err)); + } +#endif + err = cudaGetLastError(); + if(cudaSuccess != err){ + throw std::runtime_error("CudaSafeCall() failed at "+ + std::string(file) + ":" + std::to_string(line)+ + " with error " + std::to_string(err)); + } +} + +#endif diff --git a/include/lanczos/utils/defines.h b/include/lanczos/utils/defines.h new file mode 100644 index 00000000..12921470 --- /dev/null +++ b/include/lanczos/utils/defines.h @@ -0,0 +1,14 @@ +#ifndef LANCZOS_DEFINES_H +#define LANCZOS_DEFINES_H +#ifndef DOUBLE_PRECISION +#define SINGLE_PRECISION +#endif + +namespace lanczos{ +#ifndef DOUBLE_PRECISION + using real = float; +#else + using real = double; +#endif +} +#endif diff --git a/include/lanczos/utils/device_blas.h b/include/lanczos/utils/device_blas.h new file mode 100644 index 00000000..b0966852 --- /dev/null +++ b/include/lanczos/utils/device_blas.h @@ -0,0 +1,59 @@ +#pragma once +#ifdef CUDA_ENABLED +#include "cublasDebug.h" +#include "cuda_lib_defines.h" +#include "defines.h" +namespace lanczos { +struct Blas { + cublasHandle_t cublas_handle; + Blas() { CublasSafeCall(cublasCreate(&cublas_handle)); } + ~Blas() { CublasSafeCall(cublasDestroy(cublas_handle)); } + template void gemv(Args &&...args) { + CublasSafeCall( + cublasSgemv(cublas_handle, CUBLAS_OP_N, std::forward(args)...)); + } + + template void nrm2(Args &&...args) { + CublasSafeCall(cublasSnrm2(cublas_handle, std::forward(args)...)); + } + + template void axpy(Args &&...args) { + CublasSafeCall(cublasSaxpy(cublas_handle, std::forward(args)...)); + } + + template void dot(Args &&...args) { + CublasSafeCall(cublasSdot(cublas_handle, std::forward(args)...)); + } + + template void scal(Args &&...args) { + CublasSafeCall(cublasSscal(cublas_handle, std::forward(args)...)); + } +}; +} // namespace lanczos +#else +#include "lapack_and_blas_defines.h" +namespace lanczos { +struct Blas { + void gemv(int n, int m, real *alpha, real *A, int inca, real *B, int incb, + real *beta, real *C, int incc) { + cblas_gemv(CblasColMajor, CblasNoTrans, n, m, *alpha, A, inca, B, incb, + *beta, C, incc); + } + void nrm2(int n, const real *A, int inca, real *res) { + *res = cblas_nrm2(n, A, inca); + } + + void axpy(int n, real *alpha, real *A, int inca, real *B, int incb) { + cblas_axpy(n, *alpha, A, inca, B, incb); + } + + void dot(int n, real *A, int inca, real *B, int incb, real *alpha) { + *alpha = cblas_dot(n, A, inca, B, incb); + } + + void scal(int n, real *alpha, real *A, int inca) { + cblas_scal(n, *alpha, A, inca); + } +}; +} // namespace lanczos +#endif diff --git a/include/lanczos/utils/device_container.h b/include/lanczos/utils/device_container.h new file mode 100644 index 00000000..4b654357 --- /dev/null +++ b/include/lanczos/utils/device_container.h @@ -0,0 +1,51 @@ +#ifndef LANCZOS_DEVICE_CONTAINER_H +#define LANCZOS_DEVICE_CONTAINER_H +#include +namespace lanczos{ + namespace detail{ + template + auto getRawPointer(Container &vec){ + return vec.data(); + } + } +} +#ifdef CUDA_ENABLED +#include +namespace lanczos{ + template using device_container = thrust::device_vector; + namespace detail{ + template + auto getRawPointer(thrust::device_vector &vec){ + return thrust::raw_pointer_cast(vec.data()); + } + + template + void device_copy(Iter begin, Iter end, Iter2 out){ + thrust::copy(thrust::cuda::par, begin, end, out); + } + template + void device_fill(Iter begin, Iter end, T value){ + thrust::fill(thrust::cuda::par, begin, end, value); + } + + } +} +#else +namespace lanczos{ + template using device_container = std::vector; + namespace detail{ + template + void device_copy(Iter begin, Iter end, Iter2 out){ + std::copy(begin, end, out); + } + template + void device_fill(Iter begin, Iter end, T value){ + std::fill(begin, end, value); + } + + + } +} +#endif + +#endif diff --git a/include/lanczos/utils/lapack_and_blas_defines.h b/include/lanczos/utils/lapack_and_blas_defines.h new file mode 100644 index 00000000..9dc80f0f --- /dev/null +++ b/include/lanczos/utils/lapack_and_blas_defines.h @@ -0,0 +1,24 @@ +#ifndef LAPACK_AND_BLAS_DEFINES_H +#define LAPACK_AND_BLAS_DEFINES_H +#ifdef USE_MKL +#include +#else +#include +#include +#endif +#ifdef SINGLE_PRECISION +#define LAPACKE_steqr LAPACKE_ssteqr +#define cblas_gemv cblas_sgemv +#define cblas_axpy cblas_saxpy +#define cblas_scal cblas_sscal +#define cblas_nrm2 cblas_snrm2 +#define cblas_dot cblas_sdot +#else +#define LAPACKE_steqr LAPACKE_dsteqr +#define cblas_gemv cblas_dgemv +#define cblas_axpy cblas_daxpy +#define cblas_scal cblas_dscal +#define cblas_nrm2 cblas_dnrm2 +#define cblas_dot cblas_ddot +#endif +#endif From 08a813b2ce90f1bb4e3a66476bdc26e05dbe463f Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Tue, 22 Jul 2025 13:51:12 +0200 Subject: [PATCH 2/7] feat: compile lanczos in CUDA mode --- include/MobilityInterface/MobilityInterface.h | 27 +++++++------- include/MobilityInterface/lanczos.h | 35 +++++++++++++++---- 2 files changed, 42 insertions(+), 20 deletions(-) diff --git a/include/MobilityInterface/MobilityInterface.h b/include/MobilityInterface/MobilityInterface.h index 5252ec34..724e6207 100644 --- a/include/MobilityInterface/MobilityInterface.h +++ b/include/MobilityInterface/MobilityInterface.h @@ -7,10 +7,10 @@ #include "defines.h" #include "lanczos.h" #include "memory/container.h" +#include #include #include #include - namespace libmobility { enum class periodicity_mode { @@ -28,6 +28,7 @@ struct Parameters { real tolerance = 1e-4; // Tolerance for Lanczos fluctuations std::uint64_t seed = 0; bool includeAngular = false; + std::function lanczosCallback; }; // A list of parameters that cannot be changed by reinitializing a solver and/or @@ -47,7 +48,7 @@ class Mobility { std::uint64_t lanczosSeed; real lanczosTolerance; std::shared_ptr lanczos; - std::vector lanczosOutput; + thrust::device_vector lanczosOutput; bool includeAngular = false; std::mt19937 rng; @@ -128,8 +129,8 @@ class Mobility { throw std::runtime_error( "[libMobility] The number of particles is not set. Did you " "forget to call setPositions?"); - device_adapter linear(ilinear, device::cpu); - device_adapter angular(iangular, device::cpu); + device_adapter linear(ilinear, device::cuda); + device_adapter angular(iangular, device::cuda); if (linear.empty()) throw std::runtime_error( "[libMobility] This solver requires linear velocities"); @@ -145,8 +146,8 @@ class Mobility { this->lanczosTolerance, this->lanczosSeed); } lanczosOutput.resize(3 * numberElements); - std::fill(lanczosOutput.begin(), lanczosOutput.end(), 0); - auto dev = linear.dev; + thrust::fill(lanczosOutput.begin(), lanczosOutput.end(), 0); + auto dev = device::cuda; lanczos->sqrtMdotW( [this, dev, numberParticles](const real *f, real *mv) { // Torques are stored at the end of the force array @@ -160,14 +161,14 @@ class Mobility { device_span s_mv({mv, mv + 3 * N}, dev); Mdot(s_f, s_t, s_mv, s_mt); }, - lanczosOutput.data(), numberElements, prefactor); - std::transform(lanczosOutput.begin(), - lanczosOutput.begin() + 3 * numberParticles, linear.begin(), - linear.begin(), thrust::plus()); + lanczosOutput.data().get(), numberElements, prefactor); + thrust::transform(thrust::cuda::par, lanczosOutput.begin(), + lanczosOutput.begin() + 3 * numberParticles, + linear.begin(), linear.begin(), thrust::plus()); if (this->includeAngular) - std::transform(lanczosOutput.begin() + 3 * numberParticles, - lanczosOutput.end(), angular.begin(), angular.begin(), - thrust::plus()); + thrust::transform(thrust::cuda::par, lanczosOutput.begin() + 3 * numberParticles, + lanczosOutput.end(), angular.begin(), angular.begin(), + thrust::plus()); } // computes velocities according to the Langevin equation. diff --git a/include/MobilityInterface/lanczos.h b/include/MobilityInterface/lanczos.h index aa3e99a4..b54f3436 100644 --- a/include/MobilityInterface/lanczos.h +++ b/include/MobilityInterface/lanczos.h @@ -3,17 +3,31 @@ */ #ifndef LIBMOBILITY_LANCZOS_ADAPTOR_H #define LIBMOBILITY_LANCZOS_ADAPTOR_H -#include "LanczosAlgorithm.h" +#define CUDA_ENABLED +#include "lanczos/LanczosAlgorithm.h" +#include "third_party/saruprng.cuh" #include #include #include #include +namespace detail { +using real = lanczos::real; +struct SaruFill { + uint seed1, seed2; + __device__ real operator()(uint id) { + Saru prng(seed1, seed2, id); + return prng.gf(real(0), real(1.0)).x; + } +}; + +} // namespace detail // This class uses the LanczosAlgorithm library to compute fluctuations. class LanczosStochasticVelocities { using real = lanczos::real; lanczos::Solver lanczos; - std::vector lanczosNoise; + // std::vector lanczosNoise; + thrust::device_vector lanczosNoise; real lanczosTolerance; std::mt19937 engine; @@ -27,13 +41,20 @@ class LanczosStochasticVelocities { // dW). Where B is an operator that applies the square root of the provided // mobility. template - void sqrtMdotW(MobilityDot dot, real *result, int numberParticles, + void sqrtMdotW(MobilityDot idot, real *result, int numberParticles, real prefactor = 1) { - std::normal_distribution dist{0, 1}; - auto gen = [&]() { return dist(engine); }; lanczosNoise.resize(3 * numberParticles); - std::generate(lanczosNoise.begin(), lanczosNoise.end(), gen); - lanczos.run(dot, result, lanczosNoise.data(), lanczosTolerance, + // std::generate(lanczosNoise.begin(), lanczosNoise.end(), gen); + uint seed1 = std::uniform_int_distribution(0, UINT32_MAX)(engine); + uint seed2 = std::uniform_int_distribution(0, UINT32_MAX)(engine); + auto cit = thrust::make_counting_iterator(0); + thrust::transform(cit, cit + 3 * numberParticles, lanczosNoise.begin(), + detail::SaruFill{seed1, seed2}); + + std::function dot = [&](real *f, real *mv) { + idot(f, mv); + }; + lanczos.run(dot, result, lanczosNoise.data().get(), lanczosTolerance, 3 * numberParticles); } }; From 6f744fba8d6169695da379b085c52e37981413fb Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Tue, 22 Jul 2025 16:47:08 +0200 Subject: [PATCH 3/7] feat: remove dead code --- include/lanczos/utils/MatrixDot.h | 29 --------------------- include/lanczos/utils/defines.h | 12 ++++----- include/lanczos/utils/device_blas.h | 28 --------------------- include/lanczos/utils/device_container.h | 32 +----------------------- 4 files changed, 6 insertions(+), 95 deletions(-) delete mode 100644 include/lanczos/utils/MatrixDot.h diff --git a/include/lanczos/utils/MatrixDot.h b/include/lanczos/utils/MatrixDot.h deleted file mode 100644 index a76da941..00000000 --- a/include/lanczos/utils/MatrixDot.h +++ /dev/null @@ -1,29 +0,0 @@ -#ifndef LANCZOS_MATRIX_DOT_H -#define LANCZOS_MATRIX_DOT_H -#include "defines.h" -#include -namespace lanczos{ - - struct MatrixDot{ - void setSize(int newsize){this->m_size = newsize;} - //virtual void dot(real* v, real*Mv) = 0; - virtual void operator()(real* v, real*Mv) = 0; - protected: - int m_size; - }; - - //Transforms any callable into a MatrixDot valid to use with Lanczos - template - struct MatrixDotAdaptor: public lanczos::MatrixDot{ - Foo& foo; - MatrixDotAdaptor(Foo &&foo):foo(foo){} - void operator()(real* v, real* Mv) override{foo(v,Mv);} - }; - - template - auto createMatrixDotAdaptor(Foo &&foo){ - return MatrixDotAdaptor(foo); - } - -} -#endif diff --git a/include/lanczos/utils/defines.h b/include/lanczos/utils/defines.h index 12921470..3d0ad8e4 100644 --- a/include/lanczos/utils/defines.h +++ b/include/lanczos/utils/defines.h @@ -1,14 +1,12 @@ -#ifndef LANCZOS_DEFINES_H -#define LANCZOS_DEFINES_H +#pragma once #ifndef DOUBLE_PRECISION #define SINGLE_PRECISION #endif -namespace lanczos{ +namespace lanczos { #ifndef DOUBLE_PRECISION - using real = float; +using real = float; #else - using real = double; -#endif -} +using real = double; #endif +} // namespace lanczos diff --git a/include/lanczos/utils/device_blas.h b/include/lanczos/utils/device_blas.h index b0966852..5fc3e93a 100644 --- a/include/lanczos/utils/device_blas.h +++ b/include/lanczos/utils/device_blas.h @@ -1,5 +1,4 @@ #pragma once -#ifdef CUDA_ENABLED #include "cublasDebug.h" #include "cuda_lib_defines.h" #include "defines.h" @@ -30,30 +29,3 @@ struct Blas { } }; } // namespace lanczos -#else -#include "lapack_and_blas_defines.h" -namespace lanczos { -struct Blas { - void gemv(int n, int m, real *alpha, real *A, int inca, real *B, int incb, - real *beta, real *C, int incc) { - cblas_gemv(CblasColMajor, CblasNoTrans, n, m, *alpha, A, inca, B, incb, - *beta, C, incc); - } - void nrm2(int n, const real *A, int inca, real *res) { - *res = cblas_nrm2(n, A, inca); - } - - void axpy(int n, real *alpha, real *A, int inca, real *B, int incb) { - cblas_axpy(n, *alpha, A, inca, B, incb); - } - - void dot(int n, real *A, int inca, real *B, int incb, real *alpha) { - *alpha = cblas_dot(n, A, inca, B, incb); - } - - void scal(int n, real *alpha, real *A, int inca) { - cblas_scal(n, *alpha, A, inca); - } -}; -} // namespace lanczos -#endif diff --git a/include/lanczos/utils/device_container.h b/include/lanczos/utils/device_container.h index 4b654357..7258024b 100644 --- a/include/lanczos/utils/device_container.h +++ b/include/lanczos/utils/device_container.h @@ -1,15 +1,4 @@ -#ifndef LANCZOS_DEVICE_CONTAINER_H -#define LANCZOS_DEVICE_CONTAINER_H -#include -namespace lanczos{ - namespace detail{ - template - auto getRawPointer(Container &vec){ - return vec.data(); - } - } -} -#ifdef CUDA_ENABLED +#pragma once #include namespace lanczos{ template using device_container = thrust::device_vector; @@ -30,22 +19,3 @@ namespace lanczos{ } } -#else -namespace lanczos{ - template using device_container = std::vector; - namespace detail{ - template - void device_copy(Iter begin, Iter end, Iter2 out){ - std::copy(begin, end, out); - } - template - void device_fill(Iter begin, Iter end, T value){ - std::fill(begin, end, value); - } - - - } -} -#endif - -#endif From a7c7ec93a1378b5325455961cec9d07078d25b06 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Tue, 22 Jul 2025 16:49:29 +0200 Subject: [PATCH 4/7] style: formatting --- include/lanczos/utils/cublasDebug.h | 54 ++++++++++--------- include/lanczos/utils/cuda_lib_defines.h | 6 +-- include/lanczos/utils/debugTools.h | 40 +++++++------- include/lanczos/utils/device_container.h | 34 ++++++------ .../lanczos/utils/lapack_and_blas_defines.h | 26 ++++----- 5 files changed, 82 insertions(+), 78 deletions(-) diff --git a/include/lanczos/utils/cublasDebug.h b/include/lanczos/utils/cublasDebug.h index c28cac62..24154421 100644 --- a/include/lanczos/utils/cublasDebug.h +++ b/include/lanczos/utils/cublasDebug.h @@ -4,35 +4,41 @@ #define CUBLAS_ERROR_CHECK #endif -#include +#include #define CublasSafeCall(err) __cublasSafeCall(err, __FILE__, __LINE__) - - -const char* cublasGetErrorString(cublasStatus_t status){ - switch(status){ - case CUBLAS_STATUS_SUCCESS: return "CUBLAS_STATUS_SUCCESS"; - case CUBLAS_STATUS_NOT_INITIALIZED: return "CUBLAS_STATUS_NOT_INITIALIZED"; - case CUBLAS_STATUS_ALLOC_FAILED: return "CUBLAS_STATUS_ALLOC_FAILED"; - case CUBLAS_STATUS_INVALID_VALUE: return "CUBLAS_STATUS_INVALID_VALUE"; - case CUBLAS_STATUS_ARCH_MISMATCH: return "CUBLAS_STATUS_ARCH_MISMATCH"; - case CUBLAS_STATUS_MAPPING_ERROR: return "CUBLAS_STATUS_MAPPING_ERROR"; - case CUBLAS_STATUS_EXECUTION_FAILED: return "CUBLAS_STATUS_EXECUTION_FAILED"; - case CUBLAS_STATUS_INTERNAL_ERROR: return "CUBLAS_STATUS_INTERNAL_ERROR"; - default: return "Cublas Unknown error"; +const char *cublasGetErrorString(cublasStatus_t status) { + switch (status) { + case CUBLAS_STATUS_SUCCESS: + return "CUBLAS_STATUS_SUCCESS"; + case CUBLAS_STATUS_NOT_INITIALIZED: + return "CUBLAS_STATUS_NOT_INITIALIZED"; + case CUBLAS_STATUS_ALLOC_FAILED: + return "CUBLAS_STATUS_ALLOC_FAILED"; + case CUBLAS_STATUS_INVALID_VALUE: + return "CUBLAS_STATUS_INVALID_VALUE"; + case CUBLAS_STATUS_ARCH_MISMATCH: + return "CUBLAS_STATUS_ARCH_MISMATCH"; + case CUBLAS_STATUS_MAPPING_ERROR: + return "CUBLAS_STATUS_MAPPING_ERROR"; + case CUBLAS_STATUS_EXECUTION_FAILED: + return "CUBLAS_STATUS_EXECUTION_FAILED"; + case CUBLAS_STATUS_INTERNAL_ERROR: + return "CUBLAS_STATUS_INTERNAL_ERROR"; + default: + return "Cublas Unknown error"; } } -inline void __cublasSafeCall( cublasStatus_t err, const char *file, const int line ) -{ - #ifdef CUBLAS_ERROR_CHECK - if ( CUBLAS_STATUS_SUCCESS != err ) - { - fprintf( stderr, "cublasSafeCall() failed at %s:%i : %s - code: %i\n", - file, line, cublasGetErrorString( err ), err); - exit( -1 ); - } - #endif +inline void __cublasSafeCall(cublasStatus_t err, const char *file, + const int line) { +#ifdef CUBLAS_ERROR_CHECK + if (CUBLAS_STATUS_SUCCESS != err) { + fprintf(stderr, "cublasSafeCall() failed at %s:%i : %s - code: %i\n", file, + line, cublasGetErrorString(err), err); + exit(-1); + } +#endif return; } diff --git a/include/lanczos/utils/cuda_lib_defines.h b/include/lanczos/utils/cuda_lib_defines.h index 80d23263..3c64722e 100644 --- a/include/lanczos/utils/cuda_lib_defines.h +++ b/include/lanczos/utils/cuda_lib_defines.h @@ -1,7 +1,8 @@ -/*Raul P. Pelaez 2016-2022. Precision agnostic cublas/cusolver function defines*/ +/*Raul P. Pelaez 2016-2022. Precision agnostic cublas/cusolver function + * defines*/ #ifndef CUDA_LIB_DEFINES_H #define CUDA_LIB_DEFINES_H -#include +#include #if defined SINGLE_PRECISION #define cusolverDnpotrf cusolverDnSpotrf @@ -36,5 +37,4 @@ #define cublasgemm cublasDgemm #endif - #endif diff --git a/include/lanczos/utils/debugTools.h b/include/lanczos/utils/debugTools.h index 0adf9265..dac3381b 100644 --- a/include/lanczos/utils/debugTools.h +++ b/include/lanczos/utils/debugTools.h @@ -11,36 +11,36 @@ #define CudaSafeCall(err) __cudaSafeCall(err, __FILE__, __LINE__) #define CudaCheckError() __cudaCheckError(__FILE__, __LINE__) -#include -#include -#include +#include +#include +#include -inline void __cudaSafeCall(cudaError err, const char *file, const int line){ - #ifdef CUDA_ERROR_CHECK - if (cudaSuccess != err){ - cudaGetLastError(); //Reset CUDA error status - throw std::runtime_error("CudaSafeCall() failed at "+ - std::string(file) + ":" + std::to_string(line)+ - " with error " + std::to_string(err)); +inline void __cudaSafeCall(cudaError err, const char *file, const int line) { +#ifdef CUDA_ERROR_CHECK + if (cudaSuccess != err) { + cudaGetLastError(); // Reset CUDA error status + throw std::runtime_error("CudaSafeCall() failed at " + std::string(file) + + ":" + std::to_string(line) + " with error " + + std::to_string(err)); } - #endif +#endif } -inline void __cudaCheckError(const char *file, const int line){ +inline void __cudaCheckError(const char *file, const int line) { cudaError err; #ifdef CUDA_ERROR_CHECK_SYNC err = cudaDeviceSynchronize(); - if(cudaSuccess != err){ - throw std::runtime_error("CudaCheckError() with sync failed at "+ - std::string(file) + ":" + std::to_string(line)+ - " with error " + std::to_string(err)); + if (cudaSuccess != err) { + throw std::runtime_error("CudaCheckError() with sync failed at " + + std::string(file) + ":" + std::to_string(line) + + " with error " + std::to_string(err)); } #endif err = cudaGetLastError(); - if(cudaSuccess != err){ - throw std::runtime_error("CudaSafeCall() failed at "+ - std::string(file) + ":" + std::to_string(line)+ - " with error " + std::to_string(err)); + if (cudaSuccess != err) { + throw std::runtime_error("CudaSafeCall() failed at " + std::string(file) + + ":" + std::to_string(line) + " with error " + + std::to_string(err)); } } diff --git a/include/lanczos/utils/device_container.h b/include/lanczos/utils/device_container.h index 7258024b..fc0bdf56 100644 --- a/include/lanczos/utils/device_container.h +++ b/include/lanczos/utils/device_container.h @@ -1,21 +1,19 @@ #pragma once -#include -namespace lanczos{ - template using device_container = thrust::device_vector; - namespace detail{ - template - auto getRawPointer(thrust::device_vector &vec){ - return thrust::raw_pointer_cast(vec.data()); - } - - template - void device_copy(Iter begin, Iter end, Iter2 out){ - thrust::copy(thrust::cuda::par, begin, end, out); - } - template - void device_fill(Iter begin, Iter end, T value){ - thrust::fill(thrust::cuda::par, begin, end, value); - } +#include +namespace lanczos { +template using device_container = thrust::device_vector; +namespace detail { +template auto getRawPointer(thrust::device_vector &vec) { + return thrust::raw_pointer_cast(vec.data()); +} - } +template +void device_copy(Iter begin, Iter end, Iter2 out) { + thrust::copy(thrust::cuda::par, begin, end, out); } +template void device_fill(Iter begin, Iter end, T value) { + thrust::fill(thrust::cuda::par, begin, end, value); +} + +} // namespace detail +} // namespace lanczos diff --git a/include/lanczos/utils/lapack_and_blas_defines.h b/include/lanczos/utils/lapack_and_blas_defines.h index 9dc80f0f..b345a7ef 100644 --- a/include/lanczos/utils/lapack_and_blas_defines.h +++ b/include/lanczos/utils/lapack_and_blas_defines.h @@ -1,24 +1,24 @@ #ifndef LAPACK_AND_BLAS_DEFINES_H #define LAPACK_AND_BLAS_DEFINES_H #ifdef USE_MKL -#include +#include #else -#include -#include +#include +#include #endif #ifdef SINGLE_PRECISION #define LAPACKE_steqr LAPACKE_ssteqr -#define cblas_gemv cblas_sgemv -#define cblas_axpy cblas_saxpy -#define cblas_scal cblas_sscal -#define cblas_nrm2 cblas_snrm2 -#define cblas_dot cblas_sdot +#define cblas_gemv cblas_sgemv +#define cblas_axpy cblas_saxpy +#define cblas_scal cblas_sscal +#define cblas_nrm2 cblas_snrm2 +#define cblas_dot cblas_sdot #else #define LAPACKE_steqr LAPACKE_dsteqr -#define cblas_gemv cblas_dgemv -#define cblas_axpy cblas_daxpy -#define cblas_scal cblas_dscal -#define cblas_nrm2 cblas_dnrm2 -#define cblas_dot cblas_ddot +#define cblas_gemv cblas_dgemv +#define cblas_axpy cblas_daxpy +#define cblas_scal cblas_dscal +#define cblas_nrm2 cblas_dnrm2 +#define cblas_dot cblas_ddot #endif #endif From e4bd24c55c13b56c8bcd8b2e0880fd5355b69ca9 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Tue, 22 Jul 2025 16:49:36 +0200 Subject: [PATCH 5/7] feat: add an optional callback function for lanczos --- include/MobilityInterface/MobilityInterface.h | 7 +- include/MobilityInterface/lanczos.h | 5 +- include/MobilityInterface/pythonify.h | 17 +++- include/lanczos/LanczosAlgorithm.cu | 17 ++-- include/lanczos/LanczosAlgorithm.h | 81 ++++++++++--------- 5 files changed, 75 insertions(+), 52 deletions(-) diff --git a/include/MobilityInterface/MobilityInterface.h b/include/MobilityInterface/MobilityInterface.h index 724e6207..00d57571 100644 --- a/include/MobilityInterface/MobilityInterface.h +++ b/include/MobilityInterface/MobilityInterface.h @@ -51,6 +51,7 @@ class Mobility { thrust::device_vector lanczosOutput; bool includeAngular = false; std::mt19937 rng; + std::function lanczosCallback; protected: Mobility() {}; @@ -97,6 +98,7 @@ class Mobility { this->initialized = true; this->lanczosSeed = this->rng(); this->lanczosTolerance = par.tolerance; + this->lanczosCallback = par.lanczosCallback; this->includeAngular = par.includeAngular; } @@ -161,12 +163,13 @@ class Mobility { device_span s_mv({mv, mv + 3 * N}, dev); Mdot(s_f, s_t, s_mv, s_mt); }, - lanczosOutput.data().get(), numberElements, prefactor); + lanczosOutput.data().get(), numberElements, lanczosCallback, prefactor); thrust::transform(thrust::cuda::par, lanczosOutput.begin(), lanczosOutput.begin() + 3 * numberParticles, linear.begin(), linear.begin(), thrust::plus()); if (this->includeAngular) - thrust::transform(thrust::cuda::par, lanczosOutput.begin() + 3 * numberParticles, + thrust::transform(thrust::cuda::par, + lanczosOutput.begin() + 3 * numberParticles, lanczosOutput.end(), angular.begin(), angular.begin(), thrust::plus()); } diff --git a/include/MobilityInterface/lanczos.h b/include/MobilityInterface/lanczos.h index b54f3436..d3ab929d 100644 --- a/include/MobilityInterface/lanczos.h +++ b/include/MobilityInterface/lanczos.h @@ -41,7 +41,7 @@ class LanczosStochasticVelocities { // dW). Where B is an operator that applies the square root of the provided // mobility. template - void sqrtMdotW(MobilityDot idot, real *result, int numberParticles, + void sqrtMdotW(MobilityDot idot, real *result, int numberParticles, std::function callback, real prefactor = 1) { lanczosNoise.resize(3 * numberParticles); // std::generate(lanczosNoise.begin(), lanczosNoise.end(), gen); @@ -50,12 +50,11 @@ class LanczosStochasticVelocities { auto cit = thrust::make_counting_iterator(0); thrust::transform(cit, cit + 3 * numberParticles, lanczosNoise.begin(), detail::SaruFill{seed1, seed2}); - std::function dot = [&](real *f, real *mv) { idot(f, mv); }; lanczos.run(dot, result, lanczosNoise.data().get(), lanczosTolerance, - 3 * numberParticles); + 3 * numberParticles, callback); } }; diff --git a/include/MobilityInterface/pythonify.h b/include/MobilityInterface/pythonify.h index d71b0506..dd0ed44b 100644 --- a/include/MobilityInterface/pythonify.h +++ b/include/MobilityInterface/pythonify.h @@ -7,13 +7,17 @@ python (accompanied by the default documentation of the mobility interface. #ifndef MOBILITY_PYTHONIFY_H #include "MobilityInterface/MobilityInterface.h" #include "memory/python_tensor.h" +#include #include #include +#include #include #include #include #include +#include #include + namespace nb = nanobind; using namespace nb::literals; namespace py = nb; @@ -159,6 +163,8 @@ includeAngular : bool, optional Whether the solver will produce angular velocities. Needed if torques are given. Default is false. tolerance : float, optional Tolerance, used for approximate methods and also for Lanczos (default fluctuation computation). Default is 1e-4. +lanczos_callback : callable, optional + Callback function to be called during the Lanczos process. It should take two arguments: the current iteration (int) and the current error (float). Default is None, which means no callback will be used. )pbdoc"; template auto call_sqrtMdotW(Solver &myself, real prefactor) { @@ -234,12 +240,19 @@ array_like template void call_initialize(Solver &myself, real eta, real a, bool includeAngular, - real tol) { + real tol, std::optional i_lanczosCallback) { libmobility::Parameters par; par.viscosity = eta; par.hydrodynamicRadius = {a}; par.tolerance = tol; par.includeAngular = includeAngular; + std::function lanczosCallback; + if (i_lanczosCallback.has_value() && bool(i_lanczosCallback.value())) { + lanczosCallback = [i_lanczosCallback](int i, real err) { + i_lanczosCallback.value()(i, err); + }; + } + par.lanczosCallback = lanczosCallback; myself.initialize(par); } @@ -356,7 +369,7 @@ auto define_module_content( "periodicityY"_a, "periodicityZ"_a) .def("initialize", call_initialize, initialize_docstring, "viscosity"_a, "hydrodynamicRadius"_a, "includeAngular"_a = false, - "tolerance"_a = 1e-4) + "tolerance"_a = 1e-4, "lanczos_callback"_a = nb::none()) .def("setPositions", call_setPositions, "The module will compute the mobility according to this set of " "positions.", diff --git a/include/lanczos/LanczosAlgorithm.cu b/include/lanczos/LanczosAlgorithm.cu index 990a8a5b..fac28e07 100644 --- a/include/lanczos/LanczosAlgorithm.cu +++ b/include/lanczos/LanczosAlgorithm.cu @@ -47,19 +47,17 @@ class KrylovSubspace { void diagonalizeSubSpace() { int size = getSubSpaceSize(); - /**************LAPACKE********************/ /*The tridiagonal matrix is stored only with its diagonal and subdiagonal*/ /*Store both in a temporal array*/ for (int i = 0; i < size; i++) { htemp[i] = hdiag[i]; htemp[i + size] = hsup[i]; } - /*P = eigenvectors must be filled with zeros, I do not know why*/ - real *h_P = P.data(); - memset(h_P, 0, size * size * sizeof(real)); + std::fill_n(P.begin(), size * size, + 0); // P must be zero filled for steqr to work correctly /*Compute eigenvalues and eigenvectors of a triangular symmetric matrix*/ auto info = LAPACKE_steqr(LAPACK_COL_MAJOR, 'I', size, &htemp[0], - &htemp[0] + size, h_P, size); + &htemp[0] + size, P.data(), size); if (info != 0) { throw std::runtime_error("[Lanczos] Could not diagonalize tridiagonal " "krylov matrix, steqr failed with code " + @@ -173,7 +171,7 @@ public: Solver::Solver() : check_convergence_steps(3) {} int Solver::run(lanczos::Dot &dot, real *Bz, const real *z, real tolerance, - int N) { + int N, lanczos::Callback callback) { oldBz.resize((N + 1), real()); /*Lanczos iterations for Krylov decomposition*/ detail::KrylovSubspace solver(N); @@ -182,10 +180,13 @@ int Solver::run(lanczos::Dot &dot, real *Bz, const real *z, real tolerance, std::min(check_convergence_steps, iterationHardLimit - 2); for (int i = 0; i < iterationHardLimit; i++) { solver.nextIteration(dot); - if (i >= checkConvergenceSteps) { + if (i >= checkConvergenceSteps || bool(callback)) { solver.computeCurrentResultEstimation(Bz); if (i > 0) { - auto currentResidual = computeError(Bz, N); + const real currentResidual = computeError(Bz, N); + if (callback) { + callback(i, currentResidual); + } if (currentResidual <= tolerance) { registerRequiredStepsForConverge(i); return i; diff --git a/include/lanczos/LanczosAlgorithm.h b/include/lanczos/LanczosAlgorithm.h index fc819b00..2ca5de6f 100644 --- a/include/lanczos/LanczosAlgorithm.h +++ b/include/lanczos/LanczosAlgorithm.h @@ -1,45 +1,52 @@ /*Raul P. Pelaez 2022. Lanczos Algotihm, Computes the matrix-vector product sqrt(M)·v using a recursive algorithm. - For that, it requires a functor in which the () operator takes an output real* array and an input real* (both device memory) as: - inline void operator()(real* in_v, real * out_Mv); - This function must fill "out" with the result of performing the M·v dot product- > out = M·a_v. - If M has size NxN and the cost of the dot product is O(M). The total cost of the algorithm is O(m·M). Where m << N. - If M·v performs a dense M-V product, the cost of the algorithm would be O(m·N^2). - References: - [1] Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations - J. Chem. Phys. 137, 064106 (2012); doi: 10.1063/1.4742347 -Some notes: - - From what I have seen, this algorithm converges to an error of ~1e-3 in a few steps (<5) and from that point a lot of iterations are needed to lower the error. - It usually achieves machine precision in under 50 iterations. - - If the matrix does not have a sqrt (not positive definite, not symmetric...) it will usually be reflected as a nan in the current error estimation. An exception will be thrown in this case. + For that, it requires a functor in which the () operator takes an output real* +array and an input real* (both device memory) as: inline void operator()(real* +in_v, real * out_Mv); This function must fill "out" with the result of +performing the M·v dot product- > out = M·a_v. If M has size NxN and the cost of +the dot product is O(M). The total cost of the algorithm is O(m·M). Where m << +N. If M·v performs a dense M-V product, the cost of the algorithm would be +O(m·N^2). References: [1] Krylov subspace methods for computing hydrodynamic +interactions in Brownian dynamics simulations J. Chem. Phys. 137, 064106 (2012); +doi: 10.1063/1.4742347 Some notes: + + From what I have seen, this algorithm converges to an error of ~1e-3 in a few +steps (<5) and from that point a lot of iterations are needed to lower the +error. It usually achieves machine precision in under 50 iterations. + + If the matrix does not have a sqrt (not positive definite, not symmetric...) +it will usually be reflected as a nan in the current error estimation. An +exception will be thrown in this case. */ #pragma once -#include"utils/defines.h" -#include -#include"utils/device_container.h" +#include "utils/defines.h" #include "utils/device_blas.h" -namespace lanczos{ - using Dot = std::function; - struct Solver{ - Solver(); - - int run(Dot &dot, real *Bv, const real* v, real tolerance, int N); - - void setIterationHardLimit(int newLimit){this->iterationHardLimit = newLimit;} - - private: - real computeError(real* Bz, int N); - void registerRequiredStepsForConverge(int steps_needed); - - Blas blas; - device_container oldBz; - int check_convergence_steps; - int iterationHardLimit = 200; - }; -} -#include"LanczosAlgorithm.cu" +#include "utils/device_container.h" +#include +namespace lanczos { +using Dot = std::function; +using Callback = std::function; +struct Solver { + Solver(); + + int run(Dot &dot, real *Bv, const real *v, real tolerance, int N, + Callback callback); + + void setIterationHardLimit(int newLimit) { + this->iterationHardLimit = newLimit; + } + +private: + real computeError(real *Bz, int N); + void registerRequiredStepsForConverge(int steps_needed); + + Blas blas; + device_container oldBz; + int check_convergence_steps; + int iterationHardLimit = 200; +}; +} // namespace lanczos +#include "LanczosAlgorithm.cu" From e5a13d099532d56636476aa1661890b940210fa5 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Tue, 22 Jul 2025 17:29:20 +0200 Subject: [PATCH 6/7] style: formatting --- include/MobilityInterface/lanczos.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/MobilityInterface/lanczos.h b/include/MobilityInterface/lanczos.h index d3ab929d..014e2868 100644 --- a/include/MobilityInterface/lanczos.h +++ b/include/MobilityInterface/lanczos.h @@ -41,8 +41,8 @@ class LanczosStochasticVelocities { // dW). Where B is an operator that applies the square root of the provided // mobility. template - void sqrtMdotW(MobilityDot idot, real *result, int numberParticles, std::function callback, - real prefactor = 1) { + void sqrtMdotW(MobilityDot idot, real *result, int numberParticles, + std::function callback, real prefactor = 1) { lanczosNoise.resize(3 * numberParticles); // std::generate(lanczosNoise.begin(), lanczosNoise.end(), gen); uint seed1 = std::uniform_int_distribution(0, UINT32_MAX)(engine); From 23e162ff0f4e538bef74d437a85c13d6643f16e7 Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Fri, 1 Aug 2025 17:28:29 -0600 Subject: [PATCH 7/7] bugfix: fix calling a host->GPU mem copy with a CUDA execution policy that failed --- include/lanczos/LanczosAlgorithm.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/lanczos/LanczosAlgorithm.cu b/include/lanczos/LanczosAlgorithm.cu index fac28e07..036111b9 100644 --- a/include/lanczos/LanczosAlgorithm.cu +++ b/include/lanczos/LanczosAlgorithm.cu @@ -79,8 +79,8 @@ class KrylovSubspace { real beta = 0.0; cblas_gemv(CblasColMajor, CblasNoTrans, size, size, alpha, h_P, size, &htemp[0], 1, beta, &htemp[0] + size, 1); - detail::device_copy(htemp.begin() + size, htemp.begin() + 2 * size, - htempGPU.begin()); + thrust::copy(htemp.begin() + size, htemp.begin() + 2 * size, + htempGPU.begin()); return detail::getRawPointer(htempGPU); }