replace Eigen with OpenBLAS

i can finally see the light
samayala22 · Feb 4, 2024 · 6597804 · 6597804
1 parent c5932a0
commit 6597804
Show file tree

Hide file tree

Showing 5 changed files with 16 additions and 20 deletions.
diff --git a/vlm/backends/avx2/include/vlm_backend_avx2.hpp b/vlm/backends/avx2/include/vlm_backend_avx2.hpp
@@ -7,11 +7,9 @@ namespace vlm {
 
 class BackendAVX2 : public Backend {
     public:
-        struct linear_solver_t;
-        std::unique_ptr<linear_solver_t> solver;
-
         std::vector<f32> lhs;
         std::vector<f32> rhs;
+        std::vector<i32> ipiv;
         std::vector<f32> gamma;
         std::vector<f32> delta_gamma;
 

diff --git a/vlm/backends/avx2/src/vlm_backend_avx2.cpp b/vlm/backends/avx2/src/vlm_backend_avx2.cpp
@@ -16,23 +16,19 @@
 #include <oneapi/tbb/blocked_range.h>
 #include <oneapi/tbb/parallel_for.h>
 
-#include <Eigen/Dense>
-#include <Eigen/IterativeLinearSolvers>
+#include <lapacke.h>
+#include <cblas.h>
 
 using namespace vlm;
 
-struct BackendAVX2::linear_solver_t {
-    Eigen::PartialPivLU<Eigen::Ref<Eigen::MatrixXf>> lu;
-    linear_solver_t(Eigen::Map<Eigen::MatrixXf>& A) : lu(A) {};
-    ~linear_solver_t() = default;
-};
 
 BackendAVX2::~BackendAVX2() = default; // Destructor definition
 
 BackendAVX2::BackendAVX2(Mesh& mesh) : Backend(mesh) {
     //tbb::global_control global_limit(oneapi::tbb::global_control::max_allowed_parallelism, 1);
     lhs.resize((u64)mesh.nb_panels_wing() * (u64)mesh.nb_panels_wing());
     rhs.resize(mesh.nb_panels_wing());
+    ipiv.resize(mesh.nb_panels_wing());
     gamma.resize(mesh.nb_panels_wing());
     delta_gamma.resize(mesh.nb_panels_wing());
 }
@@ -332,17 +328,15 @@ void BackendAVX2::compute_rhs(const FlowData& flow, const std::vector<f32>& sect
 void BackendAVX2::lu_factor() {
     SimpleTimer timer("Factor");
     const u32 n = mesh.nb_panels_wing();
-    Eigen::Map<Eigen::MatrixXf> A(lhs.data(), n, n);
-    solver = std::make_unique<linear_solver_t>(A);
+    LAPACKE_sgetrf(LAPACK_COL_MAJOR, n, n, lhs.data(), n, ipiv.data());
 }
 
 void BackendAVX2::lu_solve() {
     SimpleTimer timer("Solve");
     const u32 n = mesh.nb_panels_wing();
-    Eigen::Map<Eigen::VectorXf> x(gamma.data(), n);
-    Eigen::Map<Eigen::VectorXf> b(rhs.data(), n);
-
-    x = solver->lu.solve(b);
+    std::copy(rhs.begin(), rhs.end(), gamma.begin());
+
+    LAPACKE_sgetrs(LAPACK_COL_MAJOR, 'N', n, 1, lhs.data(), n, ipiv.data(), gamma.data(), n);
 }
 
 f32 BackendAVX2::compute_coefficient_cl(const FlowData& flow, const f32 area,

diff --git a/vlm/backends/avx2/xmake.lua b/vlm/backends/avx2/xmake.lua
@@ -1,11 +1,13 @@
-add_requires("tbb", "openmp", "eigen")
+add_requires("tbb")
+add_requires("openblas")
 
 target("backend-avx2")
     set_kind("static")
     set_default(false)
     add_vectorexts("avx2", "fma")
-    add_packages("openmp", {public = true})
-    add_packages("tbb", "eigen")
+    add_packages("tbb")
+    add_defines("HAVE_LAPACK_CONFIG_H")
+    add_packages("openblas")
 
     add_includedirs("../../include")
     add_files("src/*.cpp")

diff --git a/vlm/backends/cuda/src/vlm_backend_cuda.cu b/vlm/backends/cuda/src/vlm_backend_cuda.cu
@@ -124,6 +124,8 @@ BackendCUDA::~BackendCUDA() {
 
     CHECK_CUDA(cudaFree(d_lhs));
     CHECK_CUDA(cudaFree(d_rhs));
+    CHECK_CUDA(cudaFree(d_gamma));
+    CHECK_CUDA(cudaFree(d_delta_gamma));
 }
 
 // For the moment, cuda backend just falls back to AVX2

diff --git a/vlm/xmake.lua b/vlm/xmake.lua
@@ -21,7 +21,7 @@ target("vlm")
     set_runargs({"-i"}, {"../../../../config/elliptic.vlm"}, {"-m"}, {"../../../../mesh/elliptic_64x64.x"}, {"-o"}, {"../../../../results/elliptic.vtu"})
     add_files("dev/main.cpp")
 
--- xmake run vlm -i ../../../../config/swept.vlm -m ../../../../mesh/rectangular_4x4.x -o ../../../../results/rec.vtu
+-- xmake run vlm -i ../../../../config/elliptic.vlm -m ../../../../mesh/elliptic_128x128.x -o ../../../../results/rec.vtu
 
 -- target("vlm-tests")
 --     set_kind("binary")