fixes to interpolative decomposition and null space

examples/rlst_interpolative_decomposition.rs

@@ -1,6 +1,7 @@
 //! Demo the inverse of a matrix
 
 pub use rlst::prelude::*;
+use rlst::dense::linalg::interpolative_decomposition::{MatrixIdDecomposition, Accuracy};
 
 //Function that creates a low rank matrix by calculating a kernel given a random point distribution on an unit sphere.
 fn low_rank_matrix(n: usize, arr: &mut Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>){
@@ -54,22 +55,31 @@ pub fn main() {
     let tol: f64 = 1e-5;
 
     //Create a low rank matrix
-    let mut arr = rlst_dynamic_array2!(f64, [slice, n]);
-    low_rank_matrix(n,&mut arr);
+    let mut arr: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, n]);
+    low_rank_matrix(n, &mut arr);
 
-    let mut res = arr.view_mut().into_id_alloc(tol, None).unwrap();
+    let res: IdDecomposition<f64> = arr.view_mut().into_id_alloc(Accuracy::Tol(tol)).unwrap();
 
-    println!("The permuted matrix is:");
-    res.arr.pretty_print();
+    println!("The skeleton of the matrix is given by");
+    res.skel.pretty_print();
 
     println!("The interpolative decomposition matrix is:");
     res.id_mat.pretty_print();
 
     println!("The rank of this matrix is {}\n", res.rank);
 
-    let a_rs_app = empty_array().simple_mult_into_resize(res.id_mat.view(), res.arr.view_mut().into_subview([0, 0], [res.rank, n]));
-    let a_rs = res.arr.view_mut().into_subview([res.rank, 0], [slice-res.rank, n]);
+    //We extract the residuals of the matrix
+    let mut perm_mat: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, slice]);
+    res.get_p(perm_mat.view_mut());
+    let perm_arr: DynamicArray<f64, 2> = empty_array::<f64, 2>()
+        .simple_mult_into_resize(perm_mat.view_mut(), arr.view());
 
+    let mut a_rs: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice-res.rank, n]);
+    a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice-res.rank, n]));
+
+    //We compute an approximation of the residual columns of the matrix
+    let a_rs_app: DynamicArray<f64, 2> = empty_array().simple_mult_into_resize(res.id_mat.view(), res.skel);
+
     let error: f64 = a_rs.view().sub(a_rs_app.view()).view_flat().norm_2();
     println!("Interpolative Decomposition L2 absolute error: {}", error);
 }

examples/rlst_null_space.rs

@@ -1,14 +1,13 @@
 //! Demo the inverse of a matrix
-use rlst::dense::linalg::null_space::NullSpaceType;
 pub use rlst::prelude::*;
 
 //Here we compute the null space (B) of the rowspace of a short matrix (A).
 pub fn main() {
     let mut arr = rlst_dynamic_array2!(f64, [3, 4]);
     arr.fill_from_seed_equally_distributed(0);
-
-    let null_res = arr.view_mut().into_null_alloc(NullSpaceType::Row).unwrap();
-    let res = empty_array().simple_mult_into_resize(arr.view_mut(), null_res.null_space_arr.view());
+    let tol = 1e-15;
+    let null_res = arr.view_mut().into_null_alloc(tol).unwrap();
+    let res: Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2> = empty_array().simple_mult_into_resize(arr.view_mut(), null_res.null_space_arr.view());
 
     println!("Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2());
 

src/dense/linalg.rs

@@ -1,7 +1,7 @@
 //! Linear algebra routines
 
 
-use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar, ElementaryMatrixData};
+use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar};
 
 use self::lu::MatrixLu;
 
@@ -11,7 +11,6 @@ pub mod pseudo_inverse;
 pub mod qr;
 pub mod svd;
 pub mod interpolative_decomposition;
-pub mod elementary_matrix;
 pub mod null_space;
 
 /// Return true if stride is column major as required by Lapack.
@@ -24,9 +23,9 @@ pub fn assert_lapack_stride(stride: [usize; 2]) {
 }
 
 /// Marker trait for objects that support Matrix decompositions.
-pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull +ElementaryMatrixData{}
+pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull{}
 
-impl<T: RlstScalar + MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull + ElementaryMatrixData> LinAlg
+impl<T: RlstScalar + MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull> LinAlg
     for T
 {
 }