diff --git a/Cargo.lock b/Cargo.lock index 03f9dd9..9f0e999 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1148,6 +1148,7 @@ dependencies = [ "anyhow", "approx", "arcstr", + "bindgen 0.72.1", "cfgrammar", "const_format", "derive-where", @@ -1157,9 +1158,13 @@ dependencies = [ "indexmap 2.12.1", "itertools 0.14.0", "klayout-lyp", + "libc", "lrlex", "lrpar", "nalgebra", + "nalgebra-sparse", + "rand 0.8.5", + "rayon", "regex", "rgb", "serde", @@ -2403,9 +2408,9 @@ dependencies = [ [[package]] name = "generic-array" -version = "0.14.7" +version = "0.14.9" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "85649ca51fd72272d7821adaf274ad91c288277713d9c18820d8499a7ff69e9a" +checksum = "4bb6743198531e02858aeaea5398fcc883e71851fcbcb5a2f773e2fb6cb1edf2" dependencies = [ "typenum", "version_check", @@ -4011,6 +4016,7 @@ dependencies = [ "num-complex", "num-rational", "num-traits", + "rayon", "simba", "typenum", ] diff --git a/core/compiler/.gitignore b/core/compiler/.gitignore index d163863..dc84959 100644 --- a/core/compiler/.gitignore +++ b/core/compiler/.gitignore @@ -1 +1,2 @@ -build/ \ No newline at end of file +build/ + diff --git a/core/compiler/Cargo.toml b/core/compiler/Cargo.toml index 1b6574b..01d53e5 100644 --- a/core/compiler/Cargo.toml +++ b/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 } diff --git a/core/compiler/build.rs b/core/compiler/build.rs index 8a4bf93..acbd490 100644 --- a/core/compiler/build.rs +++ b/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"); } diff --git a/core/compiler/src/eigen_qr.cpp b/core/compiler/src/eigen_qr.cpp new file mode 100644 index 0000000..cf5d634 --- /dev/null +++ b/core/compiler/src/eigen_qr.cpp @@ -0,0 +1,140 @@ +#include +#include +#include +#include + +using namespace Eigen; +using namespace std; + +struct EigenQR { + SparseQR, COLAMDOrdering> qr_a; + SparseQR, COLAMDOrdering> 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> triplets(nnz); + + for (int i = 0; i < nnz; i++) { + triplets[i] = Triplet(rows[i], cols[i], vals[i]); + } + + SparseMatrix 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 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 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 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 in_vec(input, size); + Map 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; + } +} diff --git a/core/compiler/src/lib.rs b/core/compiler/src/lib.rs index de7c4c6..1fb715c 100644 --- a/core/compiler/src/lib.rs +++ b/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"); diff --git a/core/compiler/src/solver.rs b/core/compiler/src/solver.rs index 3ba5d1b..cee3201 100644 --- a/core/compiler/src/solver.rs +++ b/core/compiler/src/solver.rs @@ -1,7 +1,11 @@ -use approx::relative_eq; +use std::collections::HashMap; + use indexmap::{IndexMap, IndexSet}; use itertools::{Either, Itertools}; -use nalgebra::{DMatrix, DVector}; +use nalgebra::DVector; +use nalgebra_sparse::{CooMatrix, CsrMatrix}; +use rayon::iter::IntoParallelRefIterator; +use rayon::prelude::*; use serde::{Deserialize, Serialize}; const EPSILON: f64 = 1e-8; @@ -10,6 +14,7 @@ const INV_ROUND_STEP: f64 = 1. / ROUND_STEP; #[derive(Clone, Copy, Debug, Hash, PartialEq, Eq, Serialize, Deserialize, Ord, PartialOrd)] pub struct Var(u64); +use crate::spqr_eigen::SpqrFactorization; #[derive(Clone, Default)] pub struct Solver { @@ -91,36 +96,101 @@ impl Solver { if n_vars == 0 || self.constraints.is_empty() { return; } - let a = DMatrix::from_row_iterator( - self.constraints.len(), - n_vars, - self.constraints - .iter() - .flat_map(|c| c.expr.coeff_vec(n_vars)), - ); - let b = DVector::from_iterator( - self.constraints.len(), - self.constraints.iter().map(|c| -c.expr.constant), - ); - - let svd = a.clone().svd(true, true); - let vt = svd.v_t.as_ref().expect("No V^T matrix"); - let r = svd.rank(EPSILON); - if r == 0 { - return; + + let old_triplets: Vec<(usize, usize, f64)> = self + .constraints + .par_iter() + .enumerate() + .flat_map(|(c_index, c)| { + c.expr + .coeff_vec(n_vars) + .into_par_iter() + .enumerate() + .filter_map(move |(v_index, v)| { + if v != 0.0 { + Some((c_index, v_index, v)) + } else { + None + } + }) + }) + .collect(); + + let mut used = vec![false; n_vars]; + for (_, v_index, _) in &old_triplets { + used[*v_index] = true; } - let vt_recons = vt.rows(0, r); - let sol = svd.solve(&b, EPSILON).unwrap(); - for i in 0..self.next_var { - let recons = (vt_recons.transpose() * vt_recons.column(i as usize))[((i as usize), 0)]; - if !self.solved_vars.contains_key(&Var(i)) - && relative_eq!(recons, 1., epsilon = EPSILON) - { - let val = round(sol[(i as usize, 0)]); - self.solved_vars.insert(Var(i), val); + let mut var_map = vec![usize::MAX; n_vars]; + let mut rev_var_map = Vec::with_capacity(n_vars); + let mut new_index = 0; + + for (old_index, &is_used) in used.iter().enumerate() { + if is_used { + var_map[old_index] = new_index; + rev_var_map.push(old_index); + new_index += 1; + } + } + + let n = new_index; + let m = self.constraints.len(); + + let triplets: Vec<(usize, usize, f64)> = old_triplets + .into_par_iter() + .map(|(c_index, v_index, val)| { + let new_index = var_map[v_index]; + (c_index, new_index, val) + }) + .collect(); + + let temp_b: Vec = self + .constraints + .par_iter() + .map(|c| -c.expr.constant) + .collect(); + + let b = DVector::from_iterator(self.constraints.len(), temp_b); + + let temp_a_constraind_ids: Vec = self.constraints.par_iter().map(|c| c.id).collect(); + let a_constraint_ids = Vec::from_iter(temp_a_constraind_ids); + + let mut a_coo = CooMatrix::new(m, n); + for (i, j, v) in triplets.iter() { + a_coo.push(*i, *j, *v); + } + let a_sparse: CsrMatrix = CsrMatrix::from(&a_coo); + + let qr = SpqrFactorization::from_triplets(&triplets, m, n).unwrap(); + + let x = qr.solve(&b).unwrap(); + + let residual = &b - &a_sparse * &x; + + let tolerance = 1e-10; + + for i in 0..residual.nrows() { + let r = residual[(i, 0)]; + if r.abs() > tolerance { + self.inconsistent_constraints.insert(a_constraint_ids[i]); } } + + let null_space_mags = qr.get_free_indices().unwrap(); + + let par_solved_vars: HashMap = (0..n) + .into_par_iter() + .filter(|&i| null_space_mags[i] < tolerance) + .map(|i| { + let actual_val = x[(i, 0)]; + let actual_var = rev_var_map[i]; + (Var(actual_var as u64), round(actual_val)) + }) + .collect(); + + self.solved_vars.extend(par_solved_vars); + + for constraint in self.constraints.iter_mut() { substitute_expr(&self.solved_vars, &mut constraint.expr); if constraint.expr.coeffs.is_empty() @@ -300,4 +370,129 @@ mod tests { assert_relative_eq!(*solver.solved_vars.get(&y).unwrap(), 5., epsilon = EPSILON); assert!(!solver.solved_vars.contains_key(&z)); } + + #[test] + fn linear_constraints_solved_correctly_two() { + let mut solver = Solver::new(); + let x = solver.new_var(); + let y = solver.new_var(); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(1., x)], + constant: -3., + }); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(1., y)], + constant: -5., + }); + solver.solve(); + assert_relative_eq!(*solver.solved_vars.get(&x).unwrap(), 3., epsilon = EPSILON); + assert_relative_eq!(*solver.solved_vars.get(&y).unwrap(), 5., epsilon = EPSILON); + } + + #[test] + fn linear_constraints_solved_correctly_three() { + let mut solver = Solver::new(); + let x = solver.new_var(); + let y = solver.new_var(); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(1., x), (1., y)], + constant: -3., + }); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(1., y)], + constant: -5., + }); + solver.solve(); + assert_relative_eq!(*solver.solved_vars.get(&x).unwrap(), -2., epsilon = EPSILON); + assert_relative_eq!(*solver.solved_vars.get(&y).unwrap(), 5., epsilon = EPSILON); + } + + #[test] + fn linear_constraints_solved_correctly_four() { + let mut solver = Solver::new(); + let x = solver.new_var(); + let y = solver.new_var(); + let z = solver.new_var(); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(2., x), (1., y), (1., z)], + constant: -3., + }); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(1., x), (2., y), (1., z)], + constant: -5., + }); + solver.constrain_eq0(LinearExpr { + coeffs: vec![(1., x), (1., y), (2., z)], + constant: -8., + }); + solver.solve(); + assert_relative_eq!(*solver.solved_vars.get(&x).unwrap(), -1., epsilon = EPSILON); + assert_relative_eq!(*solver.solved_vars.get(&y).unwrap(), 1., epsilon = EPSILON); + assert_relative_eq!(*solver.solved_vars.get(&z).unwrap(), 4., epsilon = EPSILON); + } + + //big matrix + #[test] + fn big_num_stab_test() { + //make a graph laplacian for grid, want to see if QR = A + fn generate_bad_value() -> f64 { + use rand::Rng; + use rand::seq::SliceRandom; + use rand::thread_rng; + let mut rng = thread_rng(); + let exponents = [-8.0, 0.0, 8.0]; + let exp = exponents.choose(&mut rng).unwrap(); + let magnitude = 10.0_f64.powf(*exp); + let sign = if rng.gen_bool(0.5) { 1.0 } else { -1.0 }; + let noise = rng.gen_range(0.5..1.5); + magnitude * sign * noise + } + let size = 25; + + use nalgebra::DMatrix; + let mut a_dense = DMatrix::zeros(size * size, size * size); + let mut triplets: Vec<(usize, usize, f64)> = Vec::new(); + + for i in 0..(size) { + for j in 0..(size) { + let curr_node = i * size + j; + let mut bad_val = generate_bad_value(); + a_dense[(curr_node, curr_node)] = bad_val; + triplets.push((curr_node, curr_node, bad_val)); + if (j + 1) < size { + bad_val = generate_bad_value(); + a_dense[(curr_node, curr_node + 1)] = bad_val; + triplets.push((curr_node, curr_node + 1, bad_val)); + bad_val = generate_bad_value(); + a_dense[(curr_node + 1, curr_node)] = bad_val; + triplets.push((curr_node + 1, curr_node, bad_val)); + } + if (i + 1) < size { + let next_row_node = size * (i + 1) + j; + bad_val = generate_bad_value(); + a_dense[(curr_node, next_row_node)] = bad_val; + triplets.push((curr_node, next_row_node, bad_val)); + bad_val = generate_bad_value(); + a_dense[(next_row_node, curr_node)] = bad_val; + triplets.push((next_row_node, curr_node, bad_val)); + } + } + } + let qr = SpqrFactorization::from_triplets(&triplets, size * size, size * size).unwrap(); + let q = qr.qa_matrix().unwrap(); + let r = qr.ra_matrix().unwrap(); + + let p_indices = qr.permutation_a().unwrap(); + + let mut ap_dense = DMatrix::zeros(size * size, size * size); + for (new_col_idx, &old_col_idx) in p_indices.iter().enumerate() { + let col = a_dense.column(old_col_idx); + ap_dense.set_column(new_col_idx, &col); + } + let a_norm = ap_dense.norm(); + let resid = ap_dense - q * r; + let err = resid.norm(); + let relative_err = err / a_norm; + assert!(relative_err < 1e-12, "not num stab"); + } } diff --git a/core/compiler/src/spqr_eigen.rs b/core/compiler/src/spqr_eigen.rs new file mode 100644 index 0000000..8bf48bf --- /dev/null +++ b/core/compiler/src/spqr_eigen.rs @@ -0,0 +1,250 @@ +use nalgebra::{DMatrix, DVector}; +use std::ffi::c_void; + +unsafe extern "C" { + fn eigen_create(m: i32, n: i32) -> *mut c_void; + fn eigen_free(ptr: *mut c_void); + fn eigen_compute_qr( + ptr: *mut c_void, + rows: *const i32, + cols: *const i32, + vals: *const f64, + nnz: i32, + ) -> i32; + fn eigen_get_rank(ptr: *mut c_void) -> i32; + fn eigen_get_q_dense(ptr: *mut c_void, output: *mut f64, mode: i32); + fn eigen_get_r_dense(ptr: *mut c_void, output: *mut f64, mode: i32); + fn eigen_get_permutation(ptr: *mut c_void, output: *mut i32, mode: i32); + fn eigen_get_free_indices(ptr: *mut c_void, output: *mut f64); + fn eigen_apply_q(ptr: *mut c_void, input: *const f64, output: *mut f64, mode_mat_factor: i32, mode_transpose: i32); +} + +pub struct SpqrFactorization { + ptr: *mut c_void, + m: usize, + n: usize, + rank: usize, +} + +unsafe impl Send for SpqrFactorization {} +unsafe impl Sync for SpqrFactorization {} + +impl SpqrFactorization { + pub fn from_triplets( + triplet: &Vec<(usize, usize, f64)>, + m: usize, + n: usize, + ) -> Result { + unsafe { + if m == 0 || n == 0 { + + let ptr = eigen_create(m as i32, n as i32); + return Ok(SpqrFactorization { + ptr: ptr, + m: 0, + n: 0, + rank: 0 as usize, + }); + } + + let ptr = eigen_create(m as i32, n as i32); + + let nnz = triplet.len(); + let mut rows = Vec::with_capacity(nnz); + let mut cols = Vec::with_capacity(nnz); + let mut vals = Vec::with_capacity(nnz); + + for (r, c, v) in triplet { + rows.push(*r as i32); + cols.push(*c as i32); + vals.push(*v); + } + + let res = + eigen_compute_qr(ptr, rows.as_ptr(), cols.as_ptr(), vals.as_ptr(), nnz as i32); + + if res == 0 { + eigen_free(ptr); + } + + let rank = eigen_get_rank(ptr); + + Ok(SpqrFactorization { + ptr: ptr, + m: m, + n: n, + rank: rank as usize, + }) + } + } + + pub fn rank(&self) -> usize { + self.rank + } + + fn get_dense_matrix(&self, num_rows: usize, num_cols: usize, mode_matrix: i32, mode_factor: i32) -> DMatrix { + let mut output = vec![0.0; num_rows * num_cols]; + unsafe { + if mode_factor == 0 { + eigen_get_q_dense(self.ptr, output.as_mut_ptr(), mode_matrix); + } else { + eigen_get_r_dense(self.ptr, output.as_mut_ptr(), mode_matrix); + } + } + DMatrix::from_vec(num_rows, num_cols, output) + } + + ///returns the Dense R matrix from the QR of A + pub fn ra_matrix(&self) -> Result, String> { + Ok(self.get_dense_matrix(self.m, self.n, 0, 1)) + } + + ///Returns the Dense Q matrix from the QR of A + pub fn qa_matrix(&self) -> Result, String> { + Ok(self.get_dense_matrix(self.m, self.m, 0, 0)) + } + + ///Returns the Dense R (FULL) matrix from the QR of AT + pub fn rat_matrix(&self) -> Result, String> { + Ok(self.get_dense_matrix(self.n, self.m, 1, 1)) + } + + ///Returns the Dense Q (FULL) matrix from the QR of AT + pub fn qat_matrix(&self) -> Result, String> { + Ok(self.get_dense_matrix(self.n, self.n, 1, 0)) + } + + + pub fn permutation_a(&self) -> Result, String> { + let mut indices = vec![0_i32; self.n]; + unsafe { + eigen_get_permutation(self.ptr, indices.as_mut_ptr(), 0); + } + Ok(indices.into_iter().map(|x| x as usize).collect()) + } + + pub fn permutation_at(&self) -> Result, String> { + let mut indices = vec![0_i32; self.m]; + unsafe { + eigen_get_permutation(self.ptr, indices.as_mut_ptr(), 1); + } + Ok(indices.into_iter().map(|x| x as usize).collect()) + } + + fn apply_q_internal(&self, vec: &DVector, mode_mat_factor: i32, mode_transpose: i32) -> DVector { + let size = if mode_mat_factor == 0 { self.m } else { self.n }; + let mut output = vec![0.0_f64; size]; + + unsafe { + eigen_apply_q( + self.ptr, + vec.as_ptr(), + output.as_mut_ptr(), + mode_mat_factor, + mode_transpose, + ); + } + + DVector::from_vec(output) + } + + pub fn apply_qt(&self, b: &DVector) -> DVector { + self.apply_q_internal(b, 0, 0) + } + + + pub fn apply_q(&self, b: &DVector) -> DVector { + self.apply_q_internal(b, 0, 1) + } + + + ///Complete solve function; determines what path of action to take depending on dimensions of matrix A + pub fn solve(&self, b: &DVector) -> Result, String> { + if self.n == 0 { + return Ok(DVector::zeros(0)); + } + if self.m >= self.n { + return self.solve_regular(b); + } else { + return self.solve_underconstrained(b); + } + } + + ///Solves the system and returns the least squares solution in the case where the matrix has m >= n + pub fn solve_regular(&self, b: &DVector) -> Result, String> { + let r = self.ra_matrix().unwrap(); + + let r_r = r.view((0, 0), (self.rank, self.rank)); + let perm_vec = self.permutation_a().unwrap(); + + let c = self.apply_qt(&b); + let c_r = c.rows(0, self.rank); + + let y_r = r_r.solve_upper_triangular(&c_r).expect("asdf \n"); + + let mut y_prime = DVector::zeros(self.n); + + y_prime.rows_mut(0, self.rank).copy_from(&y_r); + + let mut x_prime = DVector::zeros(self.n); + for i in 0..self.n { + x_prime[perm_vec[i]] = y_prime[i]; + } + + return Ok(x_prime); + + } + + ///Solves the system for the underconstrained case when m < n; uses the precomputed QR of AT + pub fn solve_underconstrained(&self, b: &DVector) -> Result, String> { + let rank = self.rank(); + let r = self.rat_matrix().unwrap(); + let perm_vec = self.permutation_at().unwrap(); + + let mut c = DVector::zeros(self.m); + for i in 0..self.m { + c[i] = b[perm_vec[i]]; + } + + let r_acc = r.view((0, 0), (rank, rank)); + let c_main = c.rows(0, rank); + + let y_small = r_acc.transpose().solve_lower_triangular(&c_main).unwrap(); + + let mut y = DVector::zeros(self.n); + y.rows_mut(0, rank).copy_from(&y_small); + + let x = self.apply_q_internal(&y, 1, 1); + + Ok(x) + } + + ///Returns the nullspace vectors of A. Uses the last n - r rows of Q from AT + pub fn get_nullspace(&self) -> Result, String> { + let q = &self.qat_matrix().unwrap(); + let null_space_vectors = q.columns(self.rank, self.n - self.rank).clone_owned(); + return Ok(null_space_vectors); + } + + pub fn get_free_indices(&self) -> Result, String> { + if self.n == 0 { + return Ok(DVector::zeros(0)); + } + if self.m == 0 { + return Ok(DVector::zeros(self.n)); + } + let mut indices_mags = vec![0.0_f64; self.n]; + unsafe { + eigen_get_free_indices(self.ptr, indices_mags.as_mut_ptr()); + } + return Ok(DVector::from_vec(indices_mags)); + } +} + +impl Drop for SpqrFactorization { + fn drop(&mut self) { + unsafe { + eigen_free(self.ptr); + } + } +} diff --git a/examples/argon_workspace/lib.ar b/examples/argon_workspace/lib.ar index 09cd30b..b9699f2 100644 --- a/examples/argon_workspace/lib.ar +++ b/examples/argon_workspace/lib.ar @@ -2,5 +2,5 @@ mod utils; mod nested; cell test() { - rect("met1", x0=0., y0=0., x1=nested::nested::test(), y1=crate::utils::test()); + rect("met1", x0=0., y0=0., x1=#scope0 nested::nested::test(), y1=#scope1 crate::utils::test()); }