From 6597804d08051de53f0439b5eb40348eaedc6c9d Mon Sep 17 00:00:00 2001
From: Samuel Ayala <samayala22@gmail.com>
Date: Sun, 4 Feb 2024 15:18:42 -0500
Subject: [PATCH] replace Eigen with OpenBLAS

i can finally see the light
---
 .../avx2/include/vlm_backend_avx2.hpp         |  4 +---
 vlm/backends/avx2/src/vlm_backend_avx2.cpp    | 20 +++++++------------
 vlm/backends/avx2/xmake.lua                   |  8 +++++---
 vlm/backends/cuda/src/vlm_backend_cuda.cu     |  2 ++
 vlm/xmake.lua                                 |  2 +-
 5 files changed, 16 insertions(+), 20 deletions(-)

diff --git a/vlm/backends/avx2/include/vlm_backend_avx2.hpp b/vlm/backends/avx2/include/vlm_backend_avx2.hpp
index 8928f38..732c598 100644
--- 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
index 94c1019..df05a0f 100644
--- a/vlm/backends/avx2/src/vlm_backend_avx2.cpp
+++ b/vlm/backends/avx2/src/vlm_backend_avx2.cpp
@@ -16,16 +16,11 @@
 #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
 
@@ -33,6 +28,7 @@ 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
index 5b2c5c9..3cd2efa 100644
--- 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
index 37ab15b..38b47f1 100644
--- 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
index ff0b05f..eeef64f 100644
--- 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")