feat(solver): improve solver performance using Eigen

core/compiler/.gitignore

@@ -1 +1,2 @@
-build/
+build/
+

core/compiler/Cargo.toml

@@ -5,8 +5,12 @@ edition = "2024"
 
 [dependencies]
 derive-where = { version = "1", features = ["serde"] }
-nalgebra = "0.34"
+nalgebra = { version = "0.34", features = ["rayon"] }
+nalgebra-sparse = "0.11"
 klayout-lyp = "0.1.1"
+rayon = "1.10"
+libc = "0.2"
+rand = "0.8.5"
 gds21 = "0.2"
 
 anyhow = { workspace = true }
@@ -32,6 +36,7 @@ const_format = "0.2"
 
 [build-dependencies]
 lrpar = { version = "0.13", features = ["serde"] }
+cc = "1.2.51"
 
 cfgrammar = { workspace = true }
 lrlex = { workspace = true }

core/compiler/build.rs

@@ -30,4 +30,24 @@ fn main() {
         .mod_name("cell_l")
         .build()
         .unwrap();
+
+    let mut build = cc::Build::new();
+
+        build
+            .cpp(true)               // Use C++ compiler
+            .std("c++14")            // Eigen requires C++14 or newer
+            .file("src/eigen_qr.cpp") // Path to your C++ file
+            .flag("-O3");
+
+
+    if std::path::Path::new("/opt/homebrew/include/eigen3").exists() {
+        build.include("/opt/homebrew/include/eigen3");
+    } else {
+        build.include("/usr/local/include/eigen3");
+    }
+
+    build.compile("eigen_qr");
+
+    println!("cargo:rerun-if-changed=src/eigen_qr.cpp");
+    println!("cargo:rerun-if-changed=build.rs");
 }

core/compiler/src/eigen_qr.cpp

@@ -0,0 +1,140 @@
+#include <Eigen/Sparse>
+#include <Eigen/SparseQR>
+#include <vector>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+struct EigenQR {
+    SparseQR<SparseMatrix<double>, COLAMDOrdering<int>> qr_a;
+    SparseQR<SparseMatrix<double>, COLAMDOrdering<int>> qr_at;
+
+    int m;
+    int n;
+
+    EigenQR(int num_rows, int num_cols) : m(num_rows), n(num_cols) {}
+};
+
+extern "C" {
+    void* eigen_create(int m, int n) {
+        return new EigenQR(m, n);
+    }
+
+    void* eigen_free(void* ptr) {
+        if (ptr) {
+            delete (EigenQR*) ptr;
+        }
+    }
+
+    int eigen_compute_qr(void* ptr, int* rows, int* cols, double* vals, int nnz) {
+        EigenQR* qr = (EigenQR*) ptr;
+
+        std::vector<Triplet<double>> triplets(nnz);
+
+        for (int i = 0; i < nnz; i++) {
+            triplets[i] = Triplet<double>(rows[i], cols[i], vals[i]);
+        }
+
+        SparseMatrix<double> A(qr->m, qr->n);
+        A.setFromTriplets(triplets.begin(), triplets.end());
+        A.makeCompressed();
+
+        qr->qr_a.compute(A);
+        if (qr->qr_a.info() != Eigen::Success) {
+            return 0;
+        }
+        qr->qr_at.compute(A.transpose());
+        if (qr->qr_at.info() != Eigen::Success) {
+            return 0;
+        }
+        return 1;
+    }
+
+    int eigen_get_rank(void* ptr) {
+        if (ptr) {
+            EigenQR* qr = (EigenQR*) ptr;
+            return qr->qr_a.rank();
+        }
+        return 0;
+    }
+
+    void eigen_get_q_dense(void* ptr, double* output, int mode) {
+        EigenQR* qr = (EigenQR*) ptr;
+        auto& qr_mat = (mode == 0) ? qr->qr_a : qr->qr_at;
+        int size = (mode == 0) ? qr->m : qr-> n;
+
+        MatrixXd I = MatrixXd::Identity(size, size);
+        MatrixXd Q_dense = qr_mat.matrixQ() * I;
+
+        Map<MatrixXd> ret(output, size, size);
+        ret = Q_dense;
+    }
+
+    void eigen_get_r_dense(void* ptr, double* output, int mode) {
+        EigenQR* qr = (EigenQR*) ptr;
+        auto& qr_mat = (mode == 0) ? qr->qr_a : qr->qr_at;
+        int r_rows = (mode == 0) ? qr->m : qr->n;
+        int r_cols = (mode == 0) ? qr->n : qr->m;
+
+        Map<MatrixXd> ret(output, r_rows, r_cols);
+        ret.setZero();
+
+        MatrixXd R_dense = MatrixXd(qr_mat.matrixR());
+        int h = std::min((int) R_dense.rows(), r_rows);
+        int w = std::min((int) R_dense.cols(), r_cols);
+        ret.block(0, 0, h, w) = R_dense.block(0, 0, h, w);
+    }
+
+    void eigen_get_permutation(void* ptr, int* output, int mode) {
+        EigenQR* qr = (EigenQR*) ptr;
+        auto& qr_mat = (mode == 0) ? qr->qr_a : qr->qr_at;
+        int size = (mode == 0) ? qr->n : qr->m;
+
+        auto& p = qr_mat.colsPermutation();
+        for (int i = 0; i < size; i++) {
+            output[i] = p.indices()(i);
+        }
+    }
+
+    //gives the summed absolute value of each null space basis vector along an index
+    void eigen_get_free_indices(void* ptr, double* output) {
+        EigenQR* qr = (EigenQR*) ptr;
+
+        Map<VectorXd> ret(output, qr->n);
+        ret.setZero();
+
+        if (qr->n - qr->qr_a.rank() <= 0) {
+            return;
+        }
+
+        VectorXd temp = VectorXd::Zero(qr->n);
+
+        for (int i = qr->qr_a.rank(); i < qr->n; i++) {
+            temp.setZero();
+            temp(i) = 1.0;
+
+            VectorXd col =  qr->qr_at.matrixQ() * temp;
+            ret = ret + col.cwiseAbs();
+        }
+    }
+
+    void eigen_apply_q(void* ptr, double* input, double* output, int mode_mat_factor, int mode_transpose) {
+        EigenQR* qr = (EigenQR*) ptr;
+        auto& qr_mat = (mode_mat_factor == 0) ? qr->qr_a : qr->qr_at;
+        int size = (mode_mat_factor == 0) ? qr->m : qr->n;
+
+        Map<VectorXd> in_vec(input, size);
+        Map<VectorXd> out_vec(output, size);
+
+        VectorXd temp;
+
+        if (mode_transpose == 0) {
+            temp = qr_mat.matrixQ().transpose() * in_vec;
+        } else {
+            temp = qr_mat.matrixQ() * in_vec;
+        }
+
+        out_vec = temp;
+    }
+}

core/compiler/src/lib.rs

@@ -5,10 +5,13 @@ pub mod gds;
 pub mod layer;
 pub mod parse;
 pub mod solver;
+//pub mod spqr;
+pub mod spqr_eigen;
 
 #[cfg(test)]
 mod tests {
 
+    use gds21::GdsUnits;
     use std::path::PathBuf;
 
     use crate::{
@@ -19,7 +22,7 @@ mod tests {
     use approx::assert_relative_eq;
     use const_format::concatcp;
 
-    use crate::compile::{compile, CellArg, CompileInput};
+    use crate::compile::{CellArg, CompileInput, compile};
     const EPSILON: f64 = 1e-10;
 
     const EXAMPLES_DIR: &str = concat!(env!("CARGO_MANIFEST_DIR"), "/../../examples");
@@ -439,6 +442,7 @@ mod tests {
         cells
             .to_gds(
                 GdsMap::from_lyp(SKY130_LYP).expect("failed to create GDS map"),
+                GdsUnits::new(1e-3, 1e-9),
                 work_dir.join("layout.gds"),
             )
             .expect("Failed to write to GDS");
@@ -681,6 +685,7 @@ mod tests {
         cells
             .to_gds(
                 GdsMap::from_lyp(SKY130_LYP).expect("failed to create GDS map"),
+                GdsUnits::new(1e-3, 1e-9),
                 work_dir.join("layout.gds"),
             )
             .expect("Failed to write to GDS");