diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs
new file mode 100644
index 00000000..d76fd1f8
--- /dev/null
+++ b/examples/rlst_interpolative_decomposition.rs
@@ -0,0 +1,125 @@
+//! Demo interpolative decomposition of a matrix.
+
+use rlst::dense::linalg::interpolative_decomposition::{Accuracy, MatrixIdDecomposition};
+pub use rlst::prelude::*;
+
+//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>) {
+    //Obtain n equally distributed angles 0<phi<pi and n equally distributed angles 0<theta<2pi
+    let pi: f64 = std::f64::consts::PI;
+
+    let mut angles1 = rlst_dynamic_array2!(f64, [n, 1]);
+    angles1.fill_from_seed_equally_distributed(0);
+    angles1.scale_inplace(pi);
+
+    let mut angles2 = rlst_dynamic_array2!(f64, [n, 1]);
+    angles2.fill_from_seed_equally_distributed(1);
+    angles2.scale_inplace(2.0 * pi);
+
+    //Calculate n points on a sphere given phi and theta
+    let mut points = rlst_dynamic_array2!(f64, [n, 3]);
+
+    for i in 0..n {
+        let phi: f64 = angles1.get_value([i, 0]).unwrap();
+        let theta: f64 = angles2.get_value([i, 0]).unwrap();
+        *points.get_mut([i, 0]).unwrap() = phi.sin() * theta.cos();
+        *points.get_mut([i, 1]).unwrap() = phi.sin() * theta.sin();
+        *points.get_mut([i, 2]).unwrap() = phi.cos();
+    }
+
+    for i in 0..arr.shape()[0] {
+        let point1 = points.r().into_subview([i, 0], [1, 3]);
+        for j in 0..n {
+            if i != j {
+                let point2 = points.r().into_subview([j, 0], [1, 3]);
+                //Calculate distance between all combinations of different points.
+                let res = point1.r().sub(point2.r());
+                //Calculate kernel e^{-|points1-points2|^2}
+                *arr.get_mut([i, j]).unwrap() = 1.0 / ((res.view_flat().norm_2().pow(2.0)).exp());
+            } else {
+                //If points are equal, set the value to 1
+                *arr.get_mut([i, j]).unwrap() = 1.0;
+            }
+        }
+    }
+}
+
+fn test_fat_matrix(slice: usize, n: usize, tol: f64) {
+    let mut arr: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, n]);
+    low_rank_matrix(n, &mut arr);
+
+    let res: IdDecomposition<f64> = arr
+        .r_mut()
+        .into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans)
+        .unwrap();
+
+    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);
+
+    //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.r_mut());
+    let perm_arr: DynamicArray<f64, 2> =
+        empty_array::<f64, 2>().simple_mult_into_resize(perm_mat.r_mut(), arr.r());
+    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.r(), res.skel);
+
+    let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2();
+    println!("Interpolative Decomposition L2 absolute error: {}", error);
+}
+
+fn test_skinny_matrix(slice: usize, n: usize, tol: f64) {
+    let mut arr: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, n]);
+    low_rank_matrix(n, &mut arr);
+
+    let mut arr_trans = empty_array();
+    arr_trans.fill_from_resize(arr.r().transpose());
+
+    let res: IdDecomposition<f64> = arr_trans
+        .r_mut()
+        .into_id_alloc(Accuracy::Tol(tol), TransMode::Trans)
+        .unwrap();
+
+    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);
+
+    //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.r_mut());
+    let perm_arr: DynamicArray<f64, 2> =
+        empty_array::<f64, 2>().simple_mult_into_resize(perm_mat.r_mut(), arr.r());
+    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.r(), res.skel);
+
+    let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2();
+    println!("Interpolative Decomposition L2 absolute error: {}", error);
+}
+
+pub fn main() {
+    let n: usize = 100;
+    let slice: usize = 50;
+    //Tolerance given for the
+    let tol: f64 = 1e-4;
+
+    //Create a low rank matrix
+    test_fat_matrix(slice, n, tol);
+    test_skinny_matrix(slice, n, tol);
+}
diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs
new file mode 100644
index 00000000..91e2eb9e
--- /dev/null
+++ b/examples/rlst_null_space.rs
@@ -0,0 +1,32 @@
+//! Demo the null space of a matrix.
+use rlst::dense::linalg::null_space::Method;
+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 tol = 1e-15;
+    let null_res = arr.r_mut().into_null_alloc(tol, Method::Svd).unwrap();
+    let res: Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2> =
+        empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r());
+    //Result
+    println!(
+        "Nullification via SVD. Value of |A*B|_2, where B=null(A): {}",
+        res.view_flat().norm_2()
+    );
+
+    let mut arr = rlst_dynamic_array2!(f64, [4, 3]);
+    arr.fill_from_seed_equally_distributed(0);
+    let shape = arr.shape();
+    let mut arr_qr = rlst_dynamic_array2!(f64, [shape[1], shape[0]]);
+    arr_qr.fill_from(arr.r().transpose());
+    let null_res = arr_qr.r_mut().into_null_alloc(tol, Method::Qr).unwrap();
+    let res: Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2> =
+        empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r());
+    //Result
+    println!(
+        "Nullification via QR. Value of |A*B|_2, where B=null(A): {}",
+        res.view_flat().norm_2()
+    );
+}
diff --git a/proc-macro/Cargo.toml b/proc-macro/Cargo.toml
index 73a4adab..39ec5d94 100644
--- a/proc-macro/Cargo.toml
+++ b/proc-macro/Cargo.toml
@@ -19,7 +19,7 @@ categories = ["mathematics", "science"]
 proc-macro = true
 
 [dependencies]
-syn = { version = "1.0", features = ["full"]}
+syn = { version = "2.0.100", features = ["full"]}
 quote = "1.0"
 darling = "0.20"
 
diff --git a/src/dense/linalg.rs b/src/dense/linalg.rs
index bbde1707..0dc24457 100644
--- a/src/dense/linalg.rs
+++ b/src/dense/linalg.rs
@@ -1,14 +1,19 @@
 //! Linear algebra routines
 
-use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixSvd, RlstScalar};
+use crate::{
+    MatrixId, MatrixInverse, MatrixNull, MatrixPseudoInverse, MatrixQr, MatrixSvd, RlstScalar,
+};
 
 use self::lu::MatrixLu;
 
+pub mod interpolative_decomposition;
 pub mod inverse;
 pub mod lu;
+pub mod null_space;
 pub mod pseudo_inverse;
 pub mod qr;
 pub mod svd;
+pub mod triangular_arrays;
 
 /// Return true if stride is column major as required by Lapack.
 pub fn assert_lapack_stride(stride: [usize; 2]) {
@@ -20,10 +25,21 @@ pub fn assert_lapack_stride(stride: [usize; 2]) {
 }
 
 /// Marker trait for objects that support Matrix decompositions.
-pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse {}
+pub trait LinAlg:
+    MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull
+{
+}
 
-impl<T: RlstScalar + MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse> LinAlg
-    for T
+impl<
+        T: RlstScalar
+            + MatrixInverse
+            + MatrixQr
+            + MatrixSvd
+            + MatrixLu
+            + MatrixPseudoInverse
+            + MatrixId
+            + MatrixNull,
+    > LinAlg for T
 {
 }
 
diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs
new file mode 100644
index 00000000..93693ff2
--- /dev/null
+++ b/src/dense/linalg/interpolative_decomposition.rs
@@ -0,0 +1,301 @@
+//! Interpolative decomposition of a matrix.
+use crate::dense::array::Array;
+use crate::dense::linalg::qr::MatrixQrDecomposition;
+use crate::dense::traits::ResizeInPlace;
+use crate::dense::traits::{
+    MultIntoResize, RandomAccessByRef, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef,
+    UnsafeRandomAccessByValue, UnsafeRandomAccessMut,
+};
+use crate::dense::types::{c32, c64, RlstResult, RlstScalar};
+use crate::dense::types::{Side, TransMode, TriangularType}; // Import TransMode from the appropriate module
+use crate::DynamicArray;
+use crate::{empty_array, rlst_dynamic_array2, BaseArray, VectorContainer};
+use crate::{TriangularMatrix, TriangularOperations}; // Import TriangularType from the appropriate module
+/// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix.
+///
+/// The matrix interpolative decomposition is defined for a two dimensional 'long' array `arr` of
+/// shape `[m, n]`, where `n>m`.
+///
+/// # Example
+///
+/// The following command computes the interpolative decomposition of an array `a` for a given tolerance, tol.
+/// ```
+/// # use rlst::rlst_dynamic_array2;
+/// # use rlst::dense::types::TransMode;
+/// # use rlst::dense::linalg::interpolative_decomposition::Accuracy;
+/// # let tol: f64 = 1e-5;
+/// # let mut a = rlst_dynamic_array2!(f64, [50, 100]);
+/// # a.fill_from_seed_equally_distributed(0);
+/// # let res = a.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans).unwrap();
+/// ```
+/// This method allocates memory for the interpolative decomposition.
+pub trait MatrixId: RlstScalar {
+    ///This method allocates space for ID
+    fn into_id_alloc<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+            + UnsafeRandomAccessMut<2, Item = Self>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self>,
+    >(
+        arr: Array<Self, ArrayImpl, 2>,
+        rank_param: Accuracy<<Self as RlstScalar>::Real>,
+        trans_mode: TransMode,
+    ) -> RlstResult<IdDecomposition<Self>>;
+}
+
+macro_rules! implement_into_id {
+    ($scalar:ty) => {
+        impl MatrixId for $scalar {
+            fn into_id_alloc<
+                ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+                    + UnsafeRandomAccessMut<2, Item = Self>
+                    + Stride<2>
+                    + Shape<2>
+                    + RawAccessMut<Item = Self>,
+            >(
+                arr: Array<Self, ArrayImpl, 2>,
+                rank_param: Accuracy<<Self as RlstScalar>::Real>,
+                trans_mode: TransMode,
+            ) -> RlstResult<IdDecomposition<Self>> {
+                IdDecomposition::<$scalar>::new(arr, rank_param, trans_mode)
+            }
+        }
+    };
+}
+
+implement_into_id!(f32);
+implement_into_id!(f64);
+implement_into_id!(c32);
+implement_into_id!(c64);
+
+impl<
+        Item: RlstScalar + MatrixId,
+        ArrayImplId: UnsafeRandomAccessByValue<2, Item = Item>
+            + UnsafeRandomAccessMut<2, Item = Item>
+            + Stride<2>
+            + RawAccessMut<Item = Item>
+            + Shape<2>,
+    > Array<Item, ArrayImplId, 2>
+{
+    /// Compute the interpolative decomposition of a given 2-dimensional array.
+    pub fn into_id_alloc(
+        self,
+        rank_param: Accuracy<<Item as RlstScalar>::Real>,
+        trans_mode: TransMode,
+    ) -> RlstResult<IdDecomposition<Item>> {
+        <Item as MatrixId>::into_id_alloc(self, rank_param, trans_mode)
+    }
+}
+
+/// Compute the matrix interpolative decomposition
+pub trait MatrixIdDecomposition: Sized {
+    /// Item type
+    type Item: RlstScalar;
+    /// Create a new Interpolative Decomposition from a given array.
+    fn new<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + UnsafeRandomAccessMut<2, Item = Self::Item>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self::Item>,
+    >(
+        arr: Array<Self::Item, ArrayImpl, 2>,
+        rank_param: Accuracy<<Self::Item as RlstScalar>::Real>,
+        trans_mode: TransMode,
+    ) -> RlstResult<Self>;
+
+    ///Compute the permutation matrix associated to the Interpolative Decomposition
+    fn get_p<
+        ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + Shape<2>
+            + UnsafeRandomAccessMut<2, Item = Self::Item>
+            + UnsafeRandomAccessByRef<2, Item = Self::Item>,
+    >(
+        &self,
+        arr: Array<Self::Item, ArrayImplMut, 2>,
+    );
+}
+
+///Stores the relevant features regarding interpolative decomposition.
+pub struct IdDecomposition<Item: RlstScalar> {
+    /// skel: skeleton of the interpolative decomposition
+    pub skel: DynamicArray<Item, 2>,
+    /// perm: permutation associated to the pivoting indiced interpolative decomposition
+    pub perm: Vec<usize>,
+    /// rank: rank of the matrix associated to the interpolative decomposition for a given tolerance
+    pub rank: usize,
+    ///id_mat: interpolative matrix calculated for a given tolerance
+    pub id_mat: Array<Item, BaseArray<Item, VectorContainer<Item>, 2>, 2>,
+}
+
+#[derive(Debug, Clone)]
+///Options to decide the matrix rank
+pub enum Accuracy<T> {
+    /// Indicates that the rank of the decomposition will be computed from a given tolerance
+    Tol(T),
+    /// Indicates that the rank of the decomposition is given beforehand by the user
+    FixedRank(usize),
+    /// Computes the rank from the tolerance, and if this one is smaller than a user set range, then we stick to the user set range
+    MaxRank(T, usize),
+}
+
+macro_rules! impl_id {
+    ($scalar:ty) => {
+        impl MatrixIdDecomposition for IdDecomposition<$scalar> {
+            type Item = $scalar;
+
+            fn new<
+                ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar>
+                    + UnsafeRandomAccessMut<2, Item = Self::Item>
+                    + Stride<2>
+                    + Shape<2>
+                    + RawAccessMut<Item = $scalar>,
+            >(
+                arr: Array<$scalar, ArrayImpl, 2>,
+                rank_param: Accuracy<<$scalar as RlstScalar>::Real>,
+                trans_mode: TransMode,
+            ) -> RlstResult<Self> {
+                //We compute the QR decomposition using rlst QR decomposition
+                let mut arr_work = empty_array();
+                let mut u_tri = empty_array();
+
+                match trans_mode {
+                    TransMode::Trans => arr_work.fill_from_resize(arr.r()),
+                    TransMode::NoTrans => arr_work.fill_from_resize(arr.r().transpose()),
+                    TransMode::ConjNoTrans => arr_work.fill_from_resize(arr.r().conj()),
+                    TransMode::ConjTrans => arr_work.fill_from_resize(arr.r().transpose().conj()),
+                };
+
+                let shape = arr_work.shape();
+                u_tri.resize_in_place([shape[1], shape[1]]);
+                let dim = shape[1];
+
+                let qr = arr_work.r_mut().into_qr_alloc().unwrap();
+                //We obtain relevant parameters of the decomposition: the permutation induced by the pivoting and the R matrix
+                let perm = qr.get_perm();
+                qr.get_r(u_tri.r_mut());
+
+                //The maximum rank is given by the number of columns of the transposed matrix
+                let rank: usize;
+
+                //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter.
+                match rank_param {
+                    Accuracy::Tol(tol) => {
+                        rank = rank_from_tolerance(u_tri.r_mut(), tol);
+                    }
+                    Accuracy::FixedRank(k) => rank = k,
+                    Accuracy::MaxRank(tol, k) => {
+                        rank = std::cmp::max(k, rank_from_tolerance(u_tri.r_mut(), tol));
+                    }
+                }
+
+                let mut permutation = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]);
+                permutation.set_zero();
+
+                let mut view = permutation.r_mut();
+                for (index, &elem) in perm.iter().enumerate() {
+                    view[[index, elem]] = <$scalar as num::One>::one();
+                }
+
+                let mut perm_arr = empty_array::<$scalar, 2>();
+                perm_arr.r_mut().mult_into_resize(
+                    TransMode::NoTrans,
+                    trans_mode,
+                    num::One::one(),
+                    permutation.r_mut(),
+                    arr.r(),
+                    num::Zero::zero(),
+                );
+
+                //We permute arr to extract the columns belonging to the skeleton
+                let mut skel = empty_array();
+                skel.fill_from_resize(perm_arr.into_subview([0, 0], [rank, shape[0]]));
+                //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix.
+                //If not, compute the Interpolative Decomposition matrix.
+                if rank == 0 || rank >= dim {
+                    let mut id_mat = rlst_dynamic_array2!($scalar, [dim, dim]);
+                    id_mat.set_identity();
+                    Ok(Self {
+                        skel,
+                        perm,
+                        rank,
+                        id_mat,
+                    })
+                } else {
+                    let shape: [usize; 2] = [shape[1], shape[1]];
+                    let mut id_mat: DynamicArray<$scalar, 2> =
+                        rlst_dynamic_array2!($scalar, [dim - rank, rank]);
+                    let r11 = TriangularMatrix::<$scalar>::new(
+                        &u_tri.r_mut().into_subview([0, 0], [rank, rank]),
+                        TriangularType::Upper,
+                    )
+                    .unwrap();
+
+                    let mut r12 = u_tri
+                        .r_mut()
+                        .into_subview([0, rank], [rank, shape[1] - rank]);
+                    r11.solve(&mut r12, Side::Left, TransMode::NoTrans);
+
+                    id_mat.fill_from(r12.r().conj().transpose().r());
+                    Ok(Self {
+                        skel,
+                        perm,
+                        rank,
+                        id_mat,
+                    })
+                }
+            }
+
+            fn get_p<
+                ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar>
+                    + Shape<2>
+                    + UnsafeRandomAccessMut<2, Item = $scalar>
+                    + UnsafeRandomAccessByRef<2, Item = $scalar>,
+            >(
+                &self,
+                mut arr: Array<$scalar, ArrayImplMut, 2>,
+            ) {
+                arr.set_zero();
+                let mut view = arr.r_mut();
+
+                for (index, &elem) in self.perm.iter().enumerate() {
+                    view[[index, elem]] = <$scalar as num::One>::one();
+                }
+            }
+        }
+    };
+}
+
+impl_id!(f64);
+impl_id!(f32);
+impl_id!(c32);
+impl_id!(c64);
+
+/// Compute the rank of the decomposition from a given tolerance
+fn rank_from_tolerance<
+    Item: RlstScalar,
+    ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Item>
+        + Shape<2>
+        + UnsafeRandomAccessMut<2, Item = Item>
+        + UnsafeRandomAccessByRef<2, Item = Item>,
+>(
+    ut_mat: Array<Item, ArrayImplMut, 2>,
+    tol: <Item as RlstScalar>::Real,
+) -> usize {
+    let dim = ut_mat.shape()[0];
+    let max = ut_mat.get([0, 0]).unwrap().abs();
+
+    //We compute the rank of the matrix
+    if max.re() > num::Zero::zero() {
+        let mut rank = 0;
+        for i in 0..dim {
+            if ut_mat.get([i, i]).unwrap().abs() > tol * max {
+                rank += 1;
+            }
+        }
+        rank
+    } else {
+        dim
+    }
+}
diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs
new file mode 100644
index 00000000..27263c95
--- /dev/null
+++ b/src/dense/linalg/null_space.rs
@@ -0,0 +1,261 @@
+//! Null space.
+use crate::dense::array::Array;
+use crate::dense::traits::{
+    DefaultIterator, RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue,
+};
+use crate::dense::types::{RlstResult, RlstScalar};
+use crate::{
+    empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixQr,
+    MatrixQrDecomposition, MatrixSvd, QrDecomposition, SvdMode, VectorContainer,
+};
+use itertools::min;
+use num::Zero;
+
+/// Compute the matrix nullspace.
+///
+/// The matrix nullspace is defined for a two dimensional array `arr` of
+/// shape `[m, n]`.
+///
+/// # Example
+///
+/// The following command computes the nullspace of an array `a`.
+/// The nullspace is found in
+/// ```
+/// # use rlst::prelude::*;
+/// # use rlst::dense::linalg::null_space::{Method};
+/// # let mut a = rlst_dynamic_array2!(f64, [4, 3]);
+/// # let tol = 1e-10;
+/// # a.fill_from_seed_equally_distributed(0);
+/// # let null_res = a.r_mut().into_null_alloc(tol, Method::Qr).unwrap();
+/// ```
+/// This method allocates memory for the nullspace computation.
+pub trait MatrixNull: RlstScalar {
+    /// Compute the matrix null space
+    fn into_null_alloc<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self>,
+    >(
+        arr: Array<Self, ArrayImpl, 2>,
+        tol: <Self as RlstScalar>::Real,
+        method: Method,
+    ) -> RlstResult<NullSpace<Self>>
+    where
+        QrDecomposition<Self, ArrayImpl>: MatrixQrDecomposition<Item = Self>;
+}
+
+impl<T: RlstScalar + MatrixSvd + MatrixQr> MatrixNull for T {
+    fn into_null_alloc<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self>,
+    >(
+        arr: Array<Self, ArrayImpl, 2>,
+        tol: <Self as RlstScalar>::Real,
+        method: Method,
+    ) -> RlstResult<NullSpace<Self>>
+    where
+        QrDecomposition<Self, ArrayImpl>: MatrixQrDecomposition<Item = Self>,
+    {
+        NullSpace::<Self>::new(arr, tol, method)
+    }
+}
+
+impl<
+        Item: RlstScalar + MatrixNull,
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item>
+            + Stride<2>
+            + RawAccessMut<Item = Item>
+            + Shape<2>,
+    > Array<Item, ArrayImpl, 2>
+{
+    /// Compute the Column or Row nullspace of a given 2-dimensional array.
+    pub fn into_null_alloc(
+        self,
+        tol: <Item as RlstScalar>::Real,
+        method: Method,
+    ) -> RlstResult<NullSpace<Item>>
+    where
+        QrDecomposition<Item, ArrayImpl>: MatrixQrDecomposition<Item = Item>,
+    {
+        <Item as MatrixNull>::into_null_alloc(self, tol, method)
+    }
+}
+
+type RealScalar<T> = <T as RlstScalar>::Real;
+/// Null Space
+pub struct NullSpace<Item: RlstScalar> {
+    ///Computed null space
+    pub null_space_arr: Array<Item, BaseArray<Item, VectorContainer<Item>, 2>, 2>,
+}
+
+///Method defines the way to compute the null space
+pub enum Method {
+    ///SVD method
+    Svd,
+    ///QR method
+    Qr,
+}
+
+///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct
+pub trait NullSpaceComputation {
+    ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64)
+    type Item: RlstScalar;
+
+    ///We create the null space decomposition
+    fn new<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + Stride<2>
+            + RawAccessMut<Item = Self::Item>
+            + Shape<2>,
+    >(
+        arr: Array<Self::Item, ArrayImpl, 2>,
+        tol: RealScalar<Self::Item>,
+        method: Method,
+    ) -> RlstResult<Self>
+    where
+        Self: Sized,
+        QrDecomposition<Self::Item, ArrayImpl>: MatrixQrDecomposition<Item = Self::Item>;
+}
+
+impl<T: RlstScalar + MatrixSvd + MatrixQr> NullSpaceComputation for NullSpace<T> {
+    type Item = T;
+
+    fn new<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + Stride<2>
+            + RawAccessMut<Item = Self::Item>
+            + Shape<2>,
+    >(
+        arr: Array<Self::Item, ArrayImpl, 2>,
+        tol: RealScalar<Self::Item>,
+        method: Method,
+    ) -> RlstResult<Self>
+    where
+        QrDecomposition<Self::Item, ArrayImpl>: MatrixQrDecomposition<Item = Self::Item>,
+    {
+        match method {
+            Method::Svd => {
+                let null_space_arr = svd_nullification(arr, tol);
+                Ok(Self { null_space_arr })
+            }
+            Method::Qr => {
+                let null_space_arr = qr_nullification(arr, tol);
+                Ok(Self { null_space_arr })
+            }
+        }
+    }
+}
+
+fn svd_nullification<
+    Item: RlstScalar + MatrixSvd,
+    ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + RawAccessMut<Item = Item> + Shape<2>,
+>(
+    arr: Array<Item, ArrayImpl, 2>,
+    tol: RealScalar<Item>,
+) -> DynamicArray<Item, 2> {
+    let shape: [usize; 2] = arr.shape();
+    let dim: usize = min(shape).unwrap();
+    let mut singular_values: DynamicArray<RealScalar<Item>, 1> =
+        rlst_dynamic_array1!(RealScalar<Item>, [dim]);
+    let mode: SvdMode = SvdMode::Full;
+    let mut u: DynamicArray<Item, 2> = rlst_dynamic_array2!(Item, [shape[0], shape[0]]);
+    let mut vt: DynamicArray<Item, 2> = rlst_dynamic_array2!(Item, [shape[1], shape[1]]);
+
+    arr.into_svd_alloc(u.r_mut(), vt.r_mut(), singular_values.data_mut(), mode)
+        .unwrap();
+
+    //For a full rank rectangular matrix, then rank = dim.
+    //find_matrix_rank checks if the matrix is full rank and recomputes the rank.
+    let rank: usize = find_svd_rank::<Item>(&mut singular_values, dim, tol);
+
+    //The null space is given by the last shape[1]-rank columns of V
+    let mut null_space_arr: DynamicArray<Item, 2> = empty_array();
+    null_space_arr.fill_from_resize(
+        vt.conj()
+            .transpose()
+            .into_subview([0, rank], [shape[1], shape[1] - rank]),
+    );
+
+    null_space_arr
+}
+
+fn qr_nullification<
+    Item: RlstScalar + MatrixQr,
+    ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut<Item = Item>,
+>(
+    arr: Array<Item, ArrayImpl, 2>,
+    tol: RealScalar<Item>,
+) -> DynamicArray<Item, 2>
+where
+    QrDecomposition<Item, ArrayImpl>: MatrixQrDecomposition<Item = Item>,
+{
+    let shape = arr.shape();
+    let dim: usize = min(shape).unwrap();
+    let qr = arr.into_qr_alloc().unwrap();
+    //We compute the QR decomposition to find a linearly independent basis of the space.
+    let mut q = rlst_dynamic_array2!(Item, [shape[0], shape[0]]);
+    let _ = qr.get_q_alloc(q.r_mut());
+    let mut r_mat: DynamicArray<Item, 2> = rlst_dynamic_array2!(Item, [dim, dim]);
+    qr.get_r(r_mat.r_mut());
+    //For a full rank rectangular matrix, then rank = dim.
+    //find_matrix_rank checks if the matrix is full rank and recomputes the rank.
+    let rank: usize = find_qr_rank::<Item>(&r_mat, dim, tol);
+    let mut null_space_arr = empty_array();
+    null_space_arr.fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0] - rank]));
+    null_space_arr
+}
+
+fn find_svd_rank<Item: RlstScalar>(
+    singular_values: &mut DynamicArray<RealScalar<Item>, 1>,
+    dim: usize,
+    tol: <Item as RlstScalar>::Real,
+) -> usize {
+    //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R.
+    let max: RealScalar<Item> = singular_values
+        .r()
+        .iter()
+        .max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap())
+        .unwrap()
+        .abs();
+
+    let rank: usize = if max.re() > <RealScalar<Item> as Zero>::zero() {
+        let aux_vec: Vec<RealScalar<Item>> = singular_values
+            .iter()
+            .filter(|el| el.abs() > tol * max)
+            .collect();
+        aux_vec.len()
+    } else {
+        dim
+    };
+
+    rank
+}
+
+fn find_qr_rank<Item: RlstScalar>(
+    r_mat: &DynamicArray<Item, 2>,
+    dim: usize,
+    tol: <Item as RlstScalar>::Real,
+) -> usize {
+    //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R.
+    let mut r_diag = rlst_dynamic_array1!(Item, [dim]);
+    r_mat.get_diag(r_diag.r_mut());
+
+    let max: RealScalar<Item> = r_diag
+        .r()
+        .iter()
+        .max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap())
+        .unwrap()
+        .abs();
+
+    let rank: usize = if max.re() > <RealScalar<Item> as Zero>::zero() {
+        let aux_vec: Vec<Item> = r_diag.iter().filter(|el| el.abs() > tol * max).collect();
+        aux_vec.len()
+    } else {
+        dim
+    };
+
+    rank
+}
diff --git a/src/dense/linalg/pseudo_inverse.rs b/src/dense/linalg/pseudo_inverse.rs
index eec186b6..5b0cc40e 100644
--- a/src/dense/linalg/pseudo_inverse.rs
+++ b/src/dense/linalg/pseudo_inverse.rs
@@ -45,9 +45,9 @@ pub trait MatrixPseudoInverse: RlstScalar + MatrixSvd {
     ///
     /// # Parameters
     /// - `pinv`: Array to store the pseudo-inverse in. If `self` has shape `[m, n]` then
-    ///           `pinv` must have shape `[n, m]`.
+    ///  `pinv` must have shape `[n, m]`.
     /// - `tol`: The relative tolerance. Singular values smaller or equal `tol * s\[0\]` will be discarded,
-    ///          where s\[0\] is the largest singular value.
+    ///  where s\[0\] is the largest singular value.
     fn into_pseudo_inverse_alloc<
         ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
             + Stride<2>
@@ -103,9 +103,9 @@ impl<
     ///
     /// # Parameters
     /// - `pinv`: Array to store the pseudo-inverse in. If `self` has shape `[m, n]` then
-    ///           `pinv` must have shape `[n, m]`.
+    ///  `pinv` must have shape `[n, m]`.
     /// - `tol`: The relative tolerance. Singular values smaller or equal `tol * s\[0\]` will be discarded,
-    ///          where s\[0\] is the largest singular value.
+    ///  where s\[0\] is the largest singular value.
     pub fn into_pseudo_inverse_alloc<
         ArrayImplPInv: UnsafeRandomAccessByValue<2, Item = Item>
             + Stride<2>
diff --git a/src/dense/linalg/qr.rs b/src/dense/linalg/qr.rs
index f932b864..3c6e200f 100644
--- a/src/dense/linalg/qr.rs
+++ b/src/dense/linalg/qr.rs
@@ -90,6 +90,94 @@ pub struct QrDecomposition<
     jpvt: Vec<i32>,
 }
 
+/// Compute the QR decomposition of a matrix.
+///
+/// The QR Decomposition of an `(m,n)` matrix `A` is defined
+/// by `AP = QR`, where `P` is an `(n, n)` permutation matrix,
+/// `Q` is a `(m, m)` orthogonal matrix, and `R` is
+/// an `(m, n)` upper triangular matrix.
+pub trait MatrixQrDecomposition: Sized {
+    /// Item type
+    type Item: RlstScalar;
+    /// Array implementaion
+    type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+        + Stride<2>
+        + RawAccessMut<Item = Self::Item>
+        + Shape<2>;
+
+    /// Create a new QR Decomposition.
+    fn new(arr: Array<Self::Item, Self::ArrayImpl, 2>) -> RlstResult<Self>;
+
+    /// Return the permuation vector for the QR decomposition.
+    ///
+    /// If `perm[i] = j` then the ith column of QR corresponds
+    /// to the jth column of the original array.
+    fn get_perm(&self) -> Vec<usize>;
+
+    /// Return the R matrix of the QR decomposition.
+    ///
+    /// If `A` has dimension `(m, n)` then R has
+    /// dimension `(k, n)` with `k=min(m, n)`.
+    fn get_r<
+        ArrayImplR: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + UnsafeRandomAccessMut<2, Item = Self::Item>
+            + RawAccessMut<Item = Self::Item>
+            + Shape<2>
+            + Stride<2>,
+    >(
+        &self,
+        arr: Array<Self::Item, ArrayImplR, 2>,
+    );
+
+    /// Return the permuation matrix `P`.
+    ///
+    /// For `A` an `(m,n)` matrix `P` has dimension `(n, n)`.
+    fn get_p<
+        ArrayImplQ: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + UnsafeRandomAccessMut<2, Item = Self::Item>
+            + RawAccessMut<Item = Self::Item>
+            + Shape<2>
+            + Stride<2>,
+    >(
+        &self,
+        arr: Array<Self::Item, ArrayImplQ, 2>,
+    );
+
+    /// Return the Q matrix of the QR decomposition.
+    ///
+    /// If `A` has dimension `(m, n)` then `arr` needs
+    /// to be of dimension `(m, r)`, where `r<= m``
+    /// is the desired number of columns of `Q`.
+    ///
+    /// This method allocates temporary memory during execution.
+    fn get_q_alloc<
+        ArrayImplQ: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + UnsafeRandomAccessMut<2, Item = Self::Item>
+            + RawAccessMut<Item = Self::Item>
+            + Shape<2>
+            + Stride<2>,
+    >(
+        &self,
+        arr: Array<Self::Item, ArrayImplQ, 2>,
+    ) -> RlstResult<()>;
+
+    /// Apply Q to a given matrix.
+    ///
+    /// This method allocates temporary memory during execution.
+    fn apply_q_alloc<
+        ArrayImplQ: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + UnsafeRandomAccessMut<2, Item = Self::Item>
+            + RawAccessMut<Item = Self::Item>
+            + Shape<2>
+            + Stride<2>,
+    >(
+        &self,
+        arr: Array<Self::Item, ArrayImplQ, 2>,
+        side: ApplyQSide,
+        trans: ApplyQTrans,
+    ) -> RlstResult<()>;
+}
+
 macro_rules! implement_qr_real {
     ($scalar:ty, $geqp3:expr, $ormqr:expr) => {
         impl<
@@ -97,10 +185,12 @@ macro_rules! implement_qr_real {
                     + Stride<2>
                     + Shape<2>
                     + RawAccessMut<Item = $scalar>,
-            > QrDecomposition<$scalar, ArrayImpl>
+            > MatrixQrDecomposition for QrDecomposition<$scalar, ArrayImpl>
         {
-            /// Create a new QR Decomposition.
-            pub fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult<Self> {
+            type Item = $scalar;
+            type ArrayImpl = ArrayImpl;
+
+            fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult<Self> {
                 let stride = arr.stride();
                 let shape = arr.shape();
 
@@ -165,22 +255,14 @@ macro_rules! implement_qr_real {
                 }
             }
 
-            /// Return the permuation vector for the QR decomposition.
-            ///
-            /// If `perm[i] = j` then the ith column of QR corresponds
-            /// to the jth column of the original array.
-            pub fn get_perm(&self) -> Vec<usize> {
+            fn get_perm(&self) -> Vec<usize> {
                 self.jpvt
                     .iter()
                     .map(|&elem| elem as usize - 1)
                     .collect_vec()
             }
 
-            /// Return the R matrix of the QR decomposition.
-            ///
-            /// If `A` has dimension `(m, n)` then R has
-            /// dimension `(k, n)` with `k=min(m, n)`.
-            pub fn get_r<
+            fn get_r<
                 ArrayImplR: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -205,10 +287,7 @@ macro_rules! implement_qr_real {
                 }
             }
 
-            /// Return the permuation matrix `P`.
-            ///
-            /// For `A` an `(m,n)` matrix `P` has dimension `(n, n)`.
-            pub fn get_p<
+            fn get_p<
                 ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -226,14 +305,7 @@ macro_rules! implement_qr_real {
                 }
             }
 
-            /// Return the Q matrix of the QR decomposition.
-            ///
-            /// If `A` has dimension `(m, n)` then `arr` needs
-            /// to be of dimension `(m, r)`, where `r<= m``
-            /// is the desired number of columns of `Q`.
-            ///
-            /// This method allocates temporary memory during execution.
-            pub fn get_q_alloc<
+            fn get_q_alloc<
                 ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -249,10 +321,7 @@ macro_rules! implement_qr_real {
                 self.apply_q_alloc(arr, ApplyQSide::Left, ApplyQTrans::NoTrans)
             }
 
-            /// Apply Q to a given matrix.
-            ///
-            /// This method allocates temporary memory during execution.
-            pub fn apply_q_alloc<
+            fn apply_q_alloc<
                 ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -360,10 +429,12 @@ macro_rules! implement_qr_complex {
                     + Stride<2>
                     + Shape<2>
                     + RawAccessMut<Item = $scalar>,
-            > QrDecomposition<$scalar, ArrayImpl>
+            > MatrixQrDecomposition for QrDecomposition<$scalar, ArrayImpl>
         {
+            type Item = $scalar;
+            type ArrayImpl = ArrayImpl;
             /// Create a new QR Decomposition.
-            pub fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult<Self> {
+            fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult<Self> {
                 let stride = arr.stride();
                 let shape = arr.shape();
 
@@ -437,7 +508,7 @@ macro_rules! implement_qr_complex {
             ///
             /// If `perm[i] = j` then the ith column of QR corresponds
             /// to the jth column of the original array.
-            pub fn get_perm(&self) -> Vec<usize> {
+            fn get_perm(&self) -> Vec<usize> {
                 self.jpvt
                     .iter()
                     .map(|&elem| elem as usize - 1)
@@ -448,7 +519,7 @@ macro_rules! implement_qr_complex {
             ///
             /// If `A` has dimension `(m, n)` then R has
             /// dimension `(k, n)` with `k=min(m, n)`.
-            pub fn get_r<
+            fn get_r<
                 ArrayImplR: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -476,7 +547,7 @@ macro_rules! implement_qr_complex {
             /// Return the permuation matrix `P`.
             ///
             /// For `A` an `(m,n)` matrix `P` has dimension `(n, n)`.
-            pub fn get_p<
+            fn get_p<
                 ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -501,7 +572,7 @@ macro_rules! implement_qr_complex {
             /// is the desired number of columns of `Q`.
             ///
             /// This method allocates temporary memory during execution.
-            pub fn get_q_alloc<
+            fn get_q_alloc<
                 ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
@@ -520,7 +591,7 @@ macro_rules! implement_qr_complex {
             /// Apply Q to a given matrix.
             ///
             /// This method allocates temporary memory during execution.
-            pub fn apply_q_alloc<
+            fn apply_q_alloc<
                 ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar>
                     + UnsafeRandomAccessMut<2, Item = $scalar>
                     + RawAccessMut<Item = $scalar>
diff --git a/src/dense/linalg/svd.rs b/src/dense/linalg/svd.rs
index 7d4b9ea8..a29ab1a5 100644
--- a/src/dense/linalg/svd.rs
+++ b/src/dense/linalg/svd.rs
@@ -34,11 +34,11 @@ pub trait MatrixSvd: RlstScalar {
     /// # Parameters
     ///
     /// - `u` - Stores the matrix `U`. For the full SVD the shape
-    ///         needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`.
+    ///  needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`.
     /// - `vt` - Stores the matrix `Vt`. For the full SVD the shape needs to be `(n, n)`.
-    ///          For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores
-    ///          the complex conjugate transpose of the matrix of right singular vectors.
-    ///          Hence, the columns of `vt.transpose().conj()` will be the right singular vectors.
+    ///  For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores
+    ///  the complex conjugate transpose of the matrix of right singular vectors.
+    ///  Hence, the columns of `vt.transpose().conj()` will be the right singular vectors.
     /// - `singular_values` - Stores the `k` singular values of `A`.
     /// - `mode` - Choose between full SVD [SvdMode::Full] or reduced SVD [SvdMode::Reduced].
     ///
@@ -504,11 +504,11 @@ impl<
     /// # Parameters
     ///
     /// - `u` - Stores the matrix `U`. For the full SVD the shape
-    ///         needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`.
+    ///  needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`.
     /// - `vt` - Stores the matrix `Vt`. For the full SVD the shape needs to be `(n, n)`.
-    ///          For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores
-    ///          the complex conjugate transpose of the matrix of right singular vectors.
-    ///          Hence, the columns of `vt.transpose().conj()` will be the right singular vectors.
+    ///  For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores
+    ///  the complex conjugate transpose of the matrix of right singular vectors.
+    ///  Hence, the columns of `vt.transpose().conj()` will be the right singular vectors.
     /// - `singular_values` - Stores the `k` singular values of `A`.
     /// - `mode` - Choose between full SVD [SvdMode::Full] or reduced SVD [SvdMode::Reduced].
     ///
diff --git a/src/dense/linalg/triangular_arrays.rs b/src/dense/linalg/triangular_arrays.rs
new file mode 100644
index 00000000..766360e3
--- /dev/null
+++ b/src/dense/linalg/triangular_arrays.rs
@@ -0,0 +1,253 @@
+//! Interpolative decomposition of a matrix.
+use crate::dense::array::Array;
+use crate::dense::traits::{
+    RandomAccessByRef, RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue,
+};
+use crate::dense::types::{c32, c64, RlstResult, RlstScalar};
+use crate::{rlst_dynamic_array2, Side, TransMode, TriangularType};
+use crate::{DynamicArray, RawAccess};
+use blas::{ctrmm, ctrsm, dtrmm, dtrsm, strmm, strsm, ztrmm, ztrsm};
+
+///Interface to obtain the upper-triangular part of the matrix
+pub trait Triangular: RlstScalar {
+    ///Compute the upper triangular
+    fn into_triangular_alloc<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self>
+            + RandomAccessByRef<2, Item = Self>,
+    >(
+        arr: &Array<Self, ArrayImpl, 2>,
+        triangular_type: TriangularType,
+    ) -> RlstResult<TriangularMatrix<Self>>;
+}
+
+macro_rules! implement_into_triangular {
+    ($scalar:ty) => {
+        impl Triangular for $scalar {
+            fn into_triangular_alloc<
+                ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+                    + Stride<2>
+                    + Shape<2>
+                    + RawAccessMut<Item = Self>
+                    + RandomAccessByRef<2, Item = Self>,
+            >(
+                arr: &Array<Self, ArrayImpl, 2>,
+                triangular_type: TriangularType,
+            ) -> RlstResult<TriangularMatrix<Self>> {
+                TriangularMatrix::<$scalar>::new(arr, triangular_type)
+            }
+        }
+    };
+}
+
+implement_into_triangular!(f32);
+implement_into_triangular!(f64);
+implement_into_triangular!(c32);
+implement_into_triangular!(c64);
+
+/// A struct representing an upper triangular matrix.
+pub struct TriangularMatrix<Item: RlstScalar> {
+    /// The upper triangular part of the matrix.
+    pub tri: DynamicArray<Item, 2>,
+    /// Defines if the matrix is lower or upper triangular
+    pub triangular_type: TriangularType,
+}
+
+///Solver for upper-triangular matrix.
+pub trait TriangularOperations: Sized {
+    ///Defines an abstract type for the matrix to solve
+    type Item: RlstScalar;
+    ///Extract the upper triangular and save it
+    fn new<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + Stride<2>
+            + Shape<2>
+            + RandomAccessByRef<2, Item = Self::Item>,
+    >(
+        arr: &Array<Self::Item, ArrayImpl, 2>,
+        triangular_type: TriangularType,
+    ) -> RlstResult<Self>;
+    ///Solves the upper-triangular system
+    fn solve<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self::Item>,
+    >(
+        &self,
+        arr_b: &mut Array<Self::Item, ArrayImpl, 2>,
+        side: Side,
+        trans: TransMode,
+    );
+
+    /// Multiplies the triangular matrix by a regular matrix
+    fn mul<
+        ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+            + Stride<2>
+            + Shape<2>
+            + RawAccessMut<Item = Self::Item>,
+    >(
+        &self,
+        arr_b: &mut Array<Self::Item, ArrayImpl, 2>,
+        side: Side,
+        trans: TransMode,
+    );
+}
+
+macro_rules! implement_solve_upper_triangular {
+    ($scalar:ty, $trsm:expr, $trmm:expr) => {
+        impl TriangularOperations for TriangularMatrix<$scalar> {
+            type Item = $scalar;
+
+            fn new<
+                ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+                    + Stride<2>
+                    + Shape<2>
+                    + RandomAccessByRef<2, Item = Self::Item>,
+            >(
+                arr: &Array<Self::Item, ArrayImpl, 2>,
+                triangular_type: TriangularType,
+            ) -> RlstResult<Self> {
+                let shape = arr.shape();
+                let mut tri = rlst_dynamic_array2!(Self::Item, shape);
+                let mut view = tri.r_mut();
+
+                match triangular_type {
+                    TriangularType::Upper => {
+                        for i in 0..shape[0] {
+                            for j in i..shape[1] {
+                                view[[i, j]] = arr[[i, j]];
+                            }
+                        }
+                    }
+                    TriangularType::Lower => {
+                        for i in 0..shape[0] {
+                            for j in 0..=i {
+                                view[[i, j]] = arr[[i, j]];
+                            }
+                        }
+                    }
+                }
+                Ok(Self {
+                    tri,
+                    triangular_type,
+                })
+            }
+
+            fn solve<
+                ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+                    + Stride<2>
+                    + Shape<2>
+                    + RawAccessMut<Item = Self::Item>,
+            >(
+                &self,
+                arr_b: &mut Array<$scalar, ArrayImpl, 2>,
+                side: Side,
+                trans: TransMode,
+            ) {
+                let lda = self.tri.stride()[1] as i32;
+                let m = arr_b.shape()[0] as i32;
+                let n = arr_b.shape()[1] as i32;
+                let ldb = arr_b.stride()[1] as i32;
+                let alpha = 1.0;
+
+                let t_type_char = match self.triangular_type {
+                    TriangularType::Upper => b'U', // A is on the left
+                    TriangularType::Lower => b'L', // A is on the right
+                };
+
+                let side_char = match side {
+                    Side::Left => b'L',  // A is on the left
+                    Side::Right => b'R', // A is on the right
+                };
+
+                let trans_char = match trans {
+                    TransMode::NoTrans => b'N',
+                    TransMode::ConjNoTrans => {
+                        panic!("TransMode::ConjNoTrans not supported for trsm implementation.")
+                    }
+                    TransMode::Trans => b'T',
+                    TransMode::ConjTrans => b'C',
+                };
+
+                unsafe {
+                    // dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
+                    $trsm(
+                        side_char,   // Side: A is on the left
+                        t_type_char, // Uplo: A is upper triangular
+                        trans_char,  // TransA: no transpose
+                        b'N',        // Diag: A is not unit triangular
+                        m,
+                        n,
+                        alpha.into(),
+                        self.tri.data(),
+                        lda,
+                        arr_b.data_mut(),
+                        ldb,
+                    );
+                }
+            }
+
+            fn mul<
+                ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+                    + Stride<2>
+                    + Shape<2>
+                    + RawAccessMut<Item = Self::Item>,
+            >(
+                &self,
+                arr_b: &mut Array<$scalar, ArrayImpl, 2>,
+                side: Side,
+                trans: TransMode,
+            ) {
+                let lda = self.tri.stride()[1] as i32;
+                let m = arr_b.shape()[0] as i32;
+                let n = arr_b.shape()[1] as i32;
+                let ldb = arr_b.stride()[1] as i32;
+                let alpha = 1.0;
+
+                let t_type_char = match self.triangular_type {
+                    TriangularType::Upper => b'U', // A is on the left
+                    TriangularType::Lower => b'L', // A is on the right
+                };
+
+                let side_char = match side {
+                    Side::Left => b'L',  // A is on the left
+                    Side::Right => b'R', // A is on the right
+                };
+
+                let trans_char = match trans {
+                    TransMode::NoTrans => b'N',
+                    TransMode::ConjNoTrans => {
+                        panic!("TransMode::ConjNoTrans not supported for trsm implementation.")
+                    }
+                    TransMode::Trans => b'T',
+                    TransMode::ConjTrans => b'C',
+                };
+
+                unsafe {
+                    // dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
+                    $trmm(
+                        side_char,   // Side: A is on the left
+                        t_type_char, // Uplo: A is upper triangular
+                        trans_char,  // TransA: no transpose
+                        b'N',        // Diag: A is not unit triangular
+                        m,
+                        n,
+                        alpha.into(),
+                        self.tri.data(),
+                        lda,
+                        arr_b.data_mut(),
+                        ldb,
+                    );
+                }
+            }
+        }
+    };
+}
+
+implement_solve_upper_triangular!(f32, strsm, strmm);
+implement_solve_upper_triangular!(f64, dtrsm, dtrmm);
+implement_solve_upper_triangular!(c32, ctrsm, ctrmm);
+implement_solve_upper_triangular!(c64, ztrsm, ztrmm);
diff --git a/src/dense/types.rs b/src/dense/types.rs
index 815ff0b4..9cafa02d 100644
--- a/src/dense/types.rs
+++ b/src/dense/types.rs
@@ -431,3 +431,21 @@ pub enum TransMode {
     /// Conjugate transpose of matrix.
     ConjTrans,
 }
+
+/// Defines if the matrix is lower or upper triangular
+#[derive(Clone, Copy, PartialEq)]
+pub enum TriangularType {
+    /// The matrix is upper-triangular
+    Upper,
+    /// The matrix is lower-triangular
+    Lower,
+}
+
+/// Defines the side through which the upper or lower triangular matrix acts
+#[derive(Clone, Copy, PartialEq)]
+pub enum Side {
+    /// Left side
+    Left,
+    /// Right side
+    Right,
+}
diff --git a/src/doc/dense_linear_algebra.rs b/src/doc/dense_linear_algebra.rs
index 4b1fe2b0..366b4684 100644
--- a/src/doc/dense_linear_algebra.rs
+++ b/src/doc/dense_linear_algebra.rs
@@ -487,7 +487,7 @@
 //! qr.get_q_alloc(q_mat.r_mut());
 //! qr.get_p(p_mat.r_mut());
 //! ````
-//! The content of `arr` is overwritten with the QR decomposition. The method [get_q_alloc](crate::dense::linalg::qr::QrDecomposition::get_q_alloc)
+//! The content of `arr` is overwritten with the QR decomposition. The method [get_q_alloc](crate::dense::linalg::qr::MatrixQrDecomposition::get_q_alloc)
 //! needs to allocate additional temporary memory on the heap. This is why it is annoted with `_alloc`.
 //!
 //! # Singular value decomposition
@@ -514,6 +514,30 @@
 //! ```
 //! To compute the full SVD use the parameter [SvdMode::Full](crate::SvdMode::Full).
 //!
+//! # Interpolative decomposition
+//!
+//! To compute the Interpolative Decomposition of a two-dimensional array (long or square) `arr`, for a given tolerance, use
+//!
+//! ```
+//! # use rlst::prelude::*;
+//! # use rlst::dense::linalg::interpolative_decomposition::Accuracy;
+//! let mut rand = rand::thread_rng();
+//! let mut arr = rlst_dynamic_array2!(f64, [5, 8]);
+//! let tol: f64 = 1e-5;
+//! let res = arr.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans).unwrap();
+//! ```
+//!
+//! You can also obtain the Interpolative Decomposition by giving the algorithm a fixed rank, and the tolerance parameter is ignored, like so
+//!
+//! ```
+//! # use rlst::prelude::*;
+//! # use rlst::dense::linalg::interpolative_decomposition::Accuracy;
+//! let mut rand = rand::thread_rng();
+//! let mut arr = rlst_dynamic_array2!(f64, [5, 8]);
+//! let k = 2;
+//! let res = arr.r_mut().into_id_alloc(Accuracy::FixedRank(k), TransMode::NoTrans).unwrap();
+//! ```
+//!
 //! # Other vector functions
 //!
 //! The following other functions for arrays with a single dimension (i.e. vectors) are provided.
diff --git a/src/prelude.rs b/src/prelude.rs
index 30243460..9238c7dd 100644
--- a/src/prelude.rs
+++ b/src/prelude.rs
@@ -51,7 +51,8 @@ pub use crate::dense::traits::{
 };
 
 pub use crate::dense::types::{
-    c32, c64, DataChunk, RlstBase, RlstError, RlstNum, RlstResult, RlstScalar, TransMode,
+    c32, c64, DataChunk, RlstBase, RlstError, RlstNum, RlstResult, RlstScalar, Side, TransMode,
+    TriangularType,
 };
 
 pub use crate::dense::base_array::BaseArray;
@@ -61,11 +62,16 @@ pub use crate::dense::data_container::{
 
 pub use crate::dense::array::empty_axis::AxisPosition;
 
+pub use crate::dense::linalg::interpolative_decomposition::{IdDecomposition, MatrixId};
 pub use crate::dense::linalg::inverse::MatrixInverse;
 pub use crate::dense::linalg::lu::{LuDecomposition, MatrixLuDecomposition};
+pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace};
 pub use crate::dense::linalg::pseudo_inverse::MatrixPseudoInverse;
-pub use crate::dense::linalg::qr::{MatrixQr, QrDecomposition};
+pub use crate::dense::linalg::qr::{MatrixQr, MatrixQrDecomposition, QrDecomposition};
 pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode};
+pub use crate::dense::linalg::triangular_arrays::{
+    Triangular, TriangularMatrix, TriangularOperations,
+};
 
 pub use crate::dense::array::rank1_array::Rank1Array;