diff --git a/ndarray-linalg/src/assert.rs b/ndarray-linalg/src/assert.rs
index 74c1618b..7acdd46d 100644
--- a/ndarray-linalg/src/assert.rs
+++ b/ndarray-linalg/src/assert.rs
@@ -29,11 +29,9 @@ pub fn aclose<A: Scalar>(test: A, truth: A, atol: A::Real) {
 }
 
 /// check two arrays are close in maximum norm
-pub fn close_max<A, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, atol: A::Real)
+pub fn close_max<A, D>(test: &ArrayRef<A, D>, truth: &ArrayRef<A, D>, atol: A::Real)
 where
     A: Scalar + Lapack,
-    S1: Data<Elem = A>,
-    S2: Data<Elem = A>,
     D: Dimension,
     D::Pattern: PartialEq + Debug,
 {
@@ -48,11 +46,9 @@ where
 }
 
 /// check two arrays are close in L1 norm
-pub fn close_l1<A, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: A::Real)
+pub fn close_l1<A, D>(test: &ArrayRef<A, D>, truth: &ArrayRef<A, D>, rtol: A::Real)
 where
     A: Scalar + Lapack,
-    S1: Data<Elem = A>,
-    S2: Data<Elem = A>,
     D: Dimension,
     D::Pattern: PartialEq + Debug,
 {
@@ -67,11 +63,9 @@ where
 }
 
 /// check two arrays are close in L2 norm
-pub fn close_l2<A, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: A::Real)
+pub fn close_l2<A, D>(test: &ArrayRef<A, D>, truth: &ArrayRef<A, D>, rtol: A::Real)
 where
     A: Scalar + Lapack,
-    S1: Data<Elem = A>,
-    S2: Data<Elem = A>,
     D: Dimension,
     D::Pattern: PartialEq + Debug,
 {
diff --git a/ndarray-linalg/src/cholesky.rs b/ndarray-linalg/src/cholesky.rs
index 58cc5cee..73ecf672 100644
--- a/ndarray-linalg/src/cholesky.rs
+++ b/ndarray-linalg/src/cholesky.rs
@@ -166,13 +166,7 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A>,
 {
-    fn solvec_inplace<'a, Sb>(
-        &self,
-        b: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+    fn solvec_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>> {
         A::solve_cholesky(
             self.factor.square_layout()?,
             self.uplo,
@@ -225,10 +219,9 @@ pub trait CholeskyInplace {
     fn cholesky_inplace(&mut self, uplo: UPLO) -> Result<&mut Self>;
 }
 
-impl<A, S> Cholesky for ArrayBase<S, Ix2>
+impl<A> Cholesky for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type Output = Array2<A>;
 
@@ -251,10 +244,9 @@ where
     }
 }
 
-impl<A, S> CholeskyInplace for ArrayBase<S, Ix2>
+impl<A> CholeskyInplace for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: DataMut<Elem = A>,
 {
     fn cholesky_inplace(&mut self, uplo: UPLO) -> Result<&mut Self> {
         A::cholesky(self.square_layout()?, uplo, self.as_allocated_mut()?)?;
@@ -301,10 +293,9 @@ where
     }
 }
 
-impl<A, Si> FactorizeC<OwnedRepr<A>> for ArrayBase<Si, Ix2>
+impl<A> FactorizeC<OwnedRepr<A>> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    Si: Data<Elem = A>,
 {
     fn factorizec(&self, uplo: UPLO) -> Result<CholeskyFactorized<OwnedRepr<A>>> {
         Ok(CholeskyFactorized {
@@ -320,7 +311,7 @@ pub trait SolveC<A: Scalar> {
     /// Solves a system of linear equations `A * x = b` with Hermitian (or real
     /// symmetric) positive definite matrix `A`, where `A` is `self`, `b` is
     /// the argument, and `x` is the successful result.
-    fn solvec<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
+    fn solvec(&self, b: &ArrayRef<A, Ix1>) -> Result<Array1<A>> {
         let mut b = replicate(b);
         self.solvec_inplace(&mut b)?;
         Ok(b)
@@ -339,24 +330,14 @@ pub trait SolveC<A: Scalar> {
     /// symmetric) positive definite matrix `A`, where `A` is `self`, `b` is
     /// the argument, and `x` is the successful result. The value of `x` is
     /// also assigned to the argument.
-    fn solvec_inplace<'a, S: DataMut<Elem = A>>(
-        &self,
-        b: &'a mut ArrayBase<S, Ix1>,
-    ) -> Result<&'a mut ArrayBase<S, Ix1>>;
+    fn solvec_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>>;
 }
 
-impl<A, S> SolveC<A> for ArrayBase<S, Ix2>
+impl<A> SolveC<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
-    fn solvec_inplace<'a, Sb>(
-        &self,
-        b: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+    fn solvec_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>> {
         self.factorizec(UPLO::Upper)?.solvec_inplace(b)
     }
 }
@@ -377,10 +358,9 @@ pub trait InverseCInto {
     fn invc_into(self) -> Result<Self::Output>;
 }
 
-impl<A, S> InverseC for ArrayBase<S, Ix2>
+impl<A> InverseC for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type Output = Array2<A>;
 
@@ -435,10 +415,9 @@ pub trait DeterminantCInto {
     fn ln_detc_into(self) -> Self::Output;
 }
 
-impl<A, S> DeterminantC for ArrayBase<S, Ix2>
+impl<A> DeterminantC for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type Output = Result<<A as Scalar>::Real>;
 
diff --git a/ndarray-linalg/src/convert.rs b/ndarray-linalg/src/convert.rs
index c808211e..d588b461 100644
--- a/ndarray-linalg/src/convert.rs
+++ b/ndarray-linalg/src/convert.rs
@@ -46,33 +46,31 @@ where
     }
 }
 
-pub fn replicate<A, Sv, So, D>(a: &ArrayBase<Sv, D>) -> ArrayBase<So, D>
+pub fn replicate<A, S, D>(a: &ArrayRef<A, D>) -> ArrayBase<S, D>
 where
     A: Copy,
-    Sv: Data<Elem = A>,
-    So: DataOwned<Elem = A> + DataMut,
+    S: DataOwned<Elem = A> + DataMut,
     D: Dimension,
 {
     unsafe {
-        let ret = ArrayBase::<So, D>::build_uninit(a.dim(), |view| {
+        let ret = ArrayBase::<S, D>::build_uninit(a.dim(), |view| {
             a.assign_to(view);
         });
         ret.assume_init()
     }
 }
 
-fn clone_with_layout<A, Si, So>(l: MatrixLayout, a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
+fn clone_with_layout<A, S>(l: MatrixLayout, a: &ArrayRef<A, Ix2>) -> ArrayBase<S, Ix2>
 where
     A: Copy,
-    Si: Data<Elem = A>,
-    So: DataOwned<Elem = A> + DataMut,
+    S: DataOwned<Elem = A> + DataMut,
 {
     let shape_builder = match l {
         MatrixLayout::C { row, lda } => (row as usize, lda as usize).set_f(false),
         MatrixLayout::F { col, lda } => (lda as usize, col as usize).set_f(true),
     };
     unsafe {
-        let ret = ArrayBase::<So, _>::build_uninit(shape_builder, |view| {
+        let ret = ArrayBase::<S, _>::build_uninit(shape_builder, |view| {
             a.assign_to(view);
         });
         ret.assume_init()
@@ -119,10 +117,9 @@ where
 /// data in the triangular portion corresponding to `uplo`.
 ///
 /// ***Panics*** if `a` is not square.
-pub(crate) fn triangular_fill_hermitian<A, S>(a: &mut ArrayBase<S, Ix2>, uplo: UPLO)
+pub(crate) fn triangular_fill_hermitian<A>(a: &mut ArrayRef<A, Ix2>, uplo: UPLO)
 where
     A: Scalar + Lapack,
-    S: DataMut<Elem = A>,
 {
     assert!(a.is_square());
     match uplo {
diff --git a/ndarray-linalg/src/diagonal.rs b/ndarray-linalg/src/diagonal.rs
index a7adad4a..4006a1ef 100644
--- a/ndarray-linalg/src/diagonal.rs
+++ b/ndarray-linalg/src/diagonal.rs
@@ -24,7 +24,7 @@ impl<S: Data> IntoDiagonal<S> for ArrayBase<S, Ix1> {
     }
 }
 
-impl<A, S: Data<Elem = A>> AsDiagonal<A> for ArrayBase<S, Ix1> {
+impl<A> AsDiagonal<A> for ArrayRef<A, Ix1> {
     fn as_diagonal(&self) -> Diagonal<ViewRepr<&A>> {
         Diagonal { diag: self.view() }
     }
@@ -37,10 +37,7 @@ where
 {
     type Elem = A;
 
-    fn apply_mut<S>(&self, a: &mut ArrayBase<S, Ix1>)
-    where
-        S: DataMut<Elem = A>,
-    {
+    fn apply_mut(&self, a: &mut ArrayRef<A, Ix1>) {
         for (val, d) in a.iter_mut().zip(self.diag.iter()) {
             *val *= *d;
         }
diff --git a/ndarray-linalg/src/eig.rs b/ndarray-linalg/src/eig.rs
index 03e3ee03..ea61a232 100644
--- a/ndarray-linalg/src/eig.rs
+++ b/ndarray-linalg/src/eig.rs
@@ -39,10 +39,9 @@ pub trait Eig {
     fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)>;
 }
 
-impl<A, S> Eig for ArrayBase<S, Ix2>
+impl<A> Eig for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type EigVal = Array1<A::Complex>;
     type EigVec = Array2<A::Complex>;
@@ -65,10 +64,9 @@ pub trait EigVals {
     fn eigvals(&self) -> Result<Self::EigVal>;
 }
 
-impl<A, S> EigVals for ArrayBase<S, Ix2>
+impl<A> EigVals for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type EigVal = Array1<A::Complex>;
 
@@ -127,27 +125,71 @@ pub trait EigGeneralized {
     /// computing the eigenvalues as α/β. If `None`, no approximate comparisons to zero will be
     /// made.
     fn eig_generalized(
-        &self,
+        self,
         thresh_opt: Option<Self::Real>,
     ) -> Result<(Self::EigVal, Self::EigVec)>;
 }
 
-impl<A, S> EigGeneralized for (ArrayBase<S, Ix2>, ArrayBase<S, Ix2>)
+/// Turn arrays, references to arrays, and [`ArrayRef`]s into owned arrays
+pub trait MaybeOwnedMatrix {
+    type Elem;
+
+    /// Convert into an owned array, cloning only when necessary.
+    fn into_owned(self) -> Array2<Self::Elem>;
+}
+
+impl<S> MaybeOwnedMatrix for ArrayBase<S, Ix2>
 where
-    A: Scalar + Lapack,
-    S: Data<Elem = A>,
+    S: Data,
+    S::Elem: Clone,
 {
-    type EigVal = Array1<GeneralizedEigenvalue<A::Complex>>;
-    type EigVec = Array2<A::Complex>;
-    type Real = A::Real;
+    type Elem = S::Elem;
+
+    fn into_owned(self) -> Array2<S::Elem> {
+        ArrayBase::into_owned(self)
+    }
+}
+
+impl<S> MaybeOwnedMatrix for &ArrayBase<S, Ix2>
+where
+    S: Data,
+    S::Elem: Clone,
+{
+    type Elem = S::Elem;
+
+    fn into_owned(self) -> Array2<S::Elem> {
+        self.to_owned()
+    }
+}
+
+impl<A> MaybeOwnedMatrix for &ArrayRef2<A>
+where
+    A: Clone,
+{
+    type Elem = A;
+
+    fn into_owned(self) -> Array2<A> {
+        self.to_owned()
+    }
+}
+
+impl<T1, T2> EigGeneralized for (T1, T2)
+where
+    T1: MaybeOwnedMatrix,
+    T1::Elem: Lapack + Scalar,
+    T2: MaybeOwnedMatrix<Elem = T1::Elem>,
+{
+    type EigVal = Array1<GeneralizedEigenvalue<<T1::Elem as Scalar>::Complex>>;
+    type EigVec = Array2<<T1::Elem as Scalar>::Complex>;
+    type Real = <T1::Elem as Scalar>::Real;
 
     fn eig_generalized(
-        &self,
+        self,
         thresh_opt: Option<Self::Real>,
     ) -> Result<(Self::EigVal, Self::EigVec)> {
-        let (mut a, mut b) = (self.0.to_owned(), self.1.to_owned());
+        let (mut a, mut b) = (self.0.into_owned(), self.1.into_owned());
         let layout = a.square_layout()?;
-        let (s, t) = A::eig_generalized(
+        let (s, t) = T1::Elem::eig_generalized(
             true,
             layout,
             a.as_allocated_mut()?,
@@ -161,3 +203,25 @@ where
         ))
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use crate::MaybeOwnedMatrix;
+
+    #[test]
+    fn test_maybe_owned_matrix() {
+        let a = array![[1.0, 2.0], [3.0, 4.0]];
+        let a_ptr = a.as_ptr();
+        let a1 = MaybeOwnedMatrix::into_owned(a);
+        assert_eq!(a_ptr, a1.as_ptr());
+
+        let b = a1.clone();
+        let b1 = MaybeOwnedMatrix::into_owned(&b);
+        assert_eq!(b, b1);
+        assert_ne!(b.as_ptr(), b1.as_ptr());
+
+        let b2 = MaybeOwnedMatrix::into_owned(&*b);
+        assert_eq!(b, b2);
+        assert_ne!(b.as_ptr(), b2.as_ptr());
+    }
+}
diff --git a/ndarray-linalg/src/eigh.rs b/ndarray-linalg/src/eigh.rs
index 837f51f8..a171e9ad 100644
--- a/ndarray-linalg/src/eigh.rs
+++ b/ndarray-linalg/src/eigh.rs
@@ -117,10 +117,9 @@ where
     }
 }
 
-impl<A, S> EighInplace for ArrayBase<S, Ix2>
+impl<A> EighInplace for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: DataMut<Elem = A>,
 {
     type EigVal = Array1<A::Real>;
 
diff --git a/ndarray-linalg/src/generate.rs b/ndarray-linalg/src/generate.rs
index 5646c808..3397ea12 100644
--- a/ndarray-linalg/src/generate.rs
+++ b/ndarray-linalg/src/generate.rs
@@ -9,13 +9,12 @@ use super::qr::*;
 use super::types::*;
 
 /// Hermite conjugate matrix
-pub fn conjugate<A, Si, So>(a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
+pub fn conjugate<A, S>(a: &ArrayRef<A, Ix2>) -> ArrayBase<S, Ix2>
 where
     A: Scalar,
-    Si: Data<Elem = A>,
-    So: DataOwned<Elem = A> + DataMut,
+    S: DataOwned<Elem = A> + DataMut,
 {
-    let mut a: ArrayBase<So, Ix2> = replicate(&a.t());
+    let mut a: ArrayBase<S, Ix2> = replicate(&a.t());
     for val in a.iter_mut() {
         *val = val.conj();
     }
diff --git a/ndarray-linalg/src/inner.rs b/ndarray-linalg/src/inner.rs
index 87ca96cc..4ab122b9 100644
--- a/ndarray-linalg/src/inner.rs
+++ b/ndarray-linalg/src/inner.rs
@@ -9,18 +9,16 @@ pub trait InnerProduct {
     type Elem: Scalar;
 
     /// Inner product `(self.conjugate, rhs)
-    fn inner<S>(&self, rhs: &ArrayBase<S, Ix1>) -> Self::Elem
-    where
-        S: Data<Elem = Self::Elem>;
+    fn inner(&self, rhs: &ArrayRef<Self::Elem, Ix1>) -> Self::Elem;
 }
 
-impl<A, S> InnerProduct for ArrayBase<S, Ix1>
+impl<A> InnerProduct for ArrayRef<A, Ix1>
 where
     A: Scalar,
-    S: Data<Elem = A>,
 {
     type Elem = A;
-    fn inner<St: Data<Elem = A>>(&self, rhs: &ArrayBase<St, Ix1>) -> A {
+
+    fn inner(&self, rhs: &ArrayRef<A, Ix1>) -> A {
         assert_eq!(self.len(), rhs.len());
         Zip::from(self)
             .and(rhs)
diff --git a/ndarray-linalg/src/krylov/householder.rs b/ndarray-linalg/src/krylov/householder.rs
index c04d9049..88b0d5ed 100644
--- a/ndarray-linalg/src/krylov/householder.rs
+++ b/ndarray-linalg/src/krylov/householder.rs
@@ -8,10 +8,9 @@ use crate::{inner::*, norm::*};
 use num_traits::One;
 
 /// Calc a reflactor `w` from a vector `x`
-pub fn calc_reflector<A, S>(x: &mut ArrayBase<S, Ix1>)
+pub fn calc_reflector<A>(x: &mut ArrayRef<A, Ix1>)
 where
     A: Scalar + Lapack,
-    S: DataMut<Elem = A>,
 {
     assert!(!x.is_empty());
     let norm = x.norm_l2();
@@ -26,11 +25,9 @@ where
 /// Panic
 /// ------
 /// - if the size of `w` and `a` mismaches
-pub fn reflect<A, S1, S2>(w: &ArrayBase<S1, Ix1>, a: &mut ArrayBase<S2, Ix1>)
+pub fn reflect<A>(w: &ArrayRef<A, Ix1>, a: &mut ArrayRef<A, Ix1>)
 where
     A: Scalar + Lapack,
-    S1: Data<Elem = A>,
-    S2: DataMut<Elem = A>,
 {
     assert_eq!(w.len(), a.len());
     let n = a.len();
@@ -66,10 +63,7 @@ impl<A: Scalar + Lapack> Householder<A> {
     }
 
     /// Take a Reflection `P = I - 2ww^T`
-    fn fundamental_reflection<S>(&self, k: usize, a: &mut ArrayBase<S, Ix1>)
-    where
-        S: DataMut<Elem = A>,
-    {
+    fn fundamental_reflection(&self, k: usize, a: &mut ArrayRef<A, Ix1>) {
         assert!(k < self.v.len());
         assert_eq!(
             a.len(),
@@ -80,10 +74,7 @@ impl<A: Scalar + Lapack> Householder<A> {
     }
 
     /// Take forward reflection `P = P_l ... P_1`
-    pub fn forward_reflection<S>(&self, a: &mut ArrayBase<S, Ix1>)
-    where
-        S: DataMut<Elem = A>,
-    {
+    pub fn forward_reflection(&self, a: &mut ArrayRef<A, Ix1>) {
         assert!(a.len() == self.dim);
         let l = self.v.len();
         for k in 0..l {
@@ -92,10 +83,7 @@ impl<A: Scalar + Lapack> Householder<A> {
     }
 
     /// Take backward reflection `P = P_1 ... P_l`
-    pub fn backward_reflection<S>(&self, a: &mut ArrayBase<S, Ix1>)
-    where
-        S: DataMut<Elem = A>,
-    {
+    pub fn backward_reflection(&self, a: &mut ArrayRef<A, Ix1>) {
         assert!(a.len() == self.dim);
         let l = self.v.len();
         for k in (0..l).rev() {
@@ -104,10 +92,7 @@ impl<A: Scalar + Lapack> Householder<A> {
     }
 
     /// Compose coefficients array using reflected vector
-    fn compose_coefficients<S>(&self, a: &ArrayBase<S, Ix1>) -> Coefficients<A>
-    where
-        S: Data<Elem = A>,
-    {
+    fn compose_coefficients(&self, a: &ArrayRef<A, Ix1>) -> Coefficients<A> {
         let k = self.len();
         let res = a.slice(s![k..]).norm_l2();
         let mut c = Array1::zeros(k + 1);
@@ -122,10 +107,7 @@ impl<A: Scalar + Lapack> Householder<A> {
     }
 
     /// Construct the residual vector from reflected vector
-    fn construct_residual<S>(&self, a: &mut ArrayBase<S, Ix1>)
-    where
-        S: DataMut<Elem = A>,
-    {
+    fn construct_residual(&self, a: &mut ArrayRef<A, Ix1>) {
         let k = self.len();
         azip!((a in a.slice_mut(s![..k])) *a = A::zero());
         self.backward_reflection(a);
@@ -147,10 +129,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.tol
     }
 
-    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
-    where
-        S: DataMut<Elem = A>,
-    {
+    fn decompose(&self, a: &mut ArrayRef<A, Ix1>) -> Array1<A> {
         self.forward_reflection(a);
         let coef = self.compose_coefficients(a);
         self.construct_residual(a);
@@ -166,10 +145,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.compose_coefficients(&a)
     }
 
-    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
-    where
-        S: DataMut<Elem = A>,
-    {
+    fn div_append(&mut self, a: &mut ArrayRef<A, Ix1>) -> AppendResult<A> {
         assert_eq!(a.len(), self.dim);
         let k = self.len();
         self.forward_reflection(a);
diff --git a/ndarray-linalg/src/krylov/mgs.rs b/ndarray-linalg/src/krylov/mgs.rs
index 67806a2c..25fa45c2 100644
--- a/ndarray-linalg/src/krylov/mgs.rs
+++ b/ndarray-linalg/src/krylov/mgs.rs
@@ -42,10 +42,7 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
         self.tol
     }
 
-    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
-    where
-        S: DataMut<Elem = A>,
-    {
+    fn decompose(&self, a: &mut ArrayRef<A, Ix1>) -> Array1<A> {
         assert_eq!(a.len(), self.dim());
         let mut coef = Array1::zeros(self.len() + 1);
         for i in 0..self.len() {
@@ -77,10 +74,9 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
         self.div_append(&mut a)
     }
 
-    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
+    fn div_append(&mut self, a: &mut ArrayRef<A, Ix1>) -> AppendResult<A>
     where
         A: Lapack,
-        S: DataMut<Elem = A>,
     {
         let coef = self.decompose(a);
         let nrm = coef[coef.len() - 1].re();
diff --git a/ndarray-linalg/src/krylov/mod.rs b/ndarray-linalg/src/krylov/mod.rs
index 14c9f7ef..bc8e2a26 100644
--- a/ndarray-linalg/src/krylov/mod.rs
+++ b/ndarray-linalg/src/krylov/mod.rs
@@ -92,9 +92,7 @@ pub trait Orthogonalizer {
     /// - `a` becomes the tangent vector
     /// - The Coefficients to the current basis is returned.
     ///
-    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
-    where
-        S: DataMut<Elem = Self::Elem>;
+    fn decompose(&self, a: &mut ArrayRef<Self::Elem, Ix1>) -> Coefficients<Self::Elem>;
 
     /// Calculate the coefficient to the current basis basis
     ///
@@ -112,9 +110,7 @@ pub trait Orthogonalizer {
 
     /// Add new vector if the residual is larger than relative tolerance,
     /// and return the residual vector
-    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
-    where
-        S: DataMut<Elem = Self::Elem>;
+    fn div_append(&mut self, a: &mut ArrayRef<Self::Elem, Ix1>) -> AppendResult<Self::Elem>;
 
     /// Get Q-matrix of generated basis
     fn get_q(&self) -> Q<Self::Elem>;
diff --git a/ndarray-linalg/src/layout.rs b/ndarray-linalg/src/layout.rs
index d0d13585..5bc9c9c9 100644
--- a/ndarray-linalg/src/layout.rs
+++ b/ndarray-linalg/src/layout.rs
@@ -18,10 +18,7 @@ pub trait AllocatedArrayMut: AllocatedArray {
     fn as_allocated_mut(&mut self) -> Result<&mut [Self::Elem]>;
 }
 
-impl<A, S> AllocatedArray for ArrayBase<S, Ix2>
-where
-    S: Data<Elem = A>,
-{
+impl<A> AllocatedArray for ArrayRef<A, Ix2> {
     type Elem = A;
 
     fn layout(&self) -> Result<MatrixLayout> {
@@ -72,10 +69,7 @@ where
     }
 }
 
-impl<A, S> AllocatedArrayMut for ArrayBase<S, Ix2>
-where
-    S: DataMut<Elem = A>,
-{
+impl<A> AllocatedArrayMut for ArrayRef<A, Ix2> {
     fn as_allocated_mut(&mut self) -> Result<&mut [A]> {
         self.as_slice_memory_order_mut()
             .ok_or(LinalgError::MemoryNotCont)
diff --git a/ndarray-linalg/src/least_squares.rs b/ndarray-linalg/src/least_squares.rs
index f376f569..62eb8503 100644
--- a/ndarray-linalg/src/least_squares.rs
+++ b/ndarray-linalg/src/least_squares.rs
@@ -91,9 +91,8 @@ pub struct LeastSquaresResult<E: Scalar, I: Dimension> {
     pub residual_sum_of_squares: Option<Array<E::Real, I::Smaller>>,
 }
 /// Solve least squares for immutable references
-pub trait LeastSquaresSvd<D, E, I>
+pub trait LeastSquaresSvd<E, I>
 where
-    D: Data<Elem = E>,
     E: Scalar + Lapack,
     I: Dimension,
 {
@@ -104,7 +103,7 @@ where
     /// `A` and `rhs` must have the same layout, i.e. they must
     /// be both either row- or column-major format, otherwise a
     /// `IncompatibleShape` error is raised.
-    fn least_squares(&self, rhs: &ArrayBase<D, I>) -> Result<LeastSquaresResult<E, I>>;
+    fn least_squares(&self, rhs: &ArrayRef<E, I>) -> Result<LeastSquaresResult<E, I>>;
 }
 
 /// Solve least squares for owned matrices
@@ -127,9 +126,8 @@ where
 
 /// Solve least squares for mutable references, overwriting
 /// the input fields in the process
-pub trait LeastSquaresSvdInPlace<D, E, I>
+pub trait LeastSquaresSvdInPlace<E, I>
 where
-    D: Data<Elem = E>,
     E: Scalar + Lapack,
     I: Dimension,
 {
@@ -143,7 +141,7 @@ where
     /// `IncompatibleShape` error is raised.
     fn least_squares_in_place(
         &mut self,
-        rhs: &mut ArrayBase<D, I>,
+        rhs: &mut ArrayRef<E, I>,
     ) -> Result<LeastSquaresResult<E, I>>;
 }
 
@@ -151,11 +149,9 @@ where
 /// column vector as a right-hand side.
 /// `E` is one of `f32`, `f64`, `c32`, `c64`. `D1`, `D2` can be any
 /// valid representation for `ArrayBase` (over `E`).
-impl<E, D1, D2> LeastSquaresSvd<D2, E, Ix1> for ArrayBase<D1, Ix2>
+impl<E> LeastSquaresSvd<E, Ix1> for ArrayRef<E, Ix2>
 where
     E: Scalar + Lapack,
-    D1: Data<Elem = E>,
-    D2: Data<Elem = E>,
 {
     /// Solve a least squares problem of the form `Ax = rhs`
     /// by calling `A.least_squares(&rhs)`, where `rhs` is a
@@ -164,7 +160,7 @@ where
     /// `A` and `rhs` must have the same layout, i.e. they must
     /// be both either row- or column-major format, otherwise a
     /// `IncompatibleShape` error is raised.
-    fn least_squares(&self, rhs: &ArrayBase<D2, Ix1>) -> Result<LeastSquaresResult<E, Ix1>> {
+    fn least_squares(&self, rhs: &ArrayRef<E, Ix1>) -> Result<LeastSquaresResult<E, Ix1>> {
         let a = self.to_owned();
         let b = rhs.to_owned();
         a.least_squares_into(b)
@@ -175,11 +171,9 @@ where
 /// (=mulitipe vectors) as a right-hand side.
 /// `E` is one of `f32`, `f64`, `c32`, `c64`. `D1`, `D2` can be any
 /// valid representation for `ArrayBase` (over `E`).
-impl<E, D1, D2> LeastSquaresSvd<D2, E, Ix2> for ArrayBase<D1, Ix2>
+impl<E> LeastSquaresSvd<E, Ix2> for ArrayRef<E, Ix2>
 where
     E: Scalar + Lapack,
-    D1: Data<Elem = E>,
-    D2: Data<Elem = E>,
 {
     /// Solve a least squares problem of the form `Ax = rhs`
     /// by calling `A.least_squares(&rhs)`, where `rhs` is
@@ -188,7 +182,7 @@ where
     /// `A` and `rhs` must have the same layout, i.e. they must
     /// be both either row- or column-major format, otherwise a
     /// `IncompatibleShape` error is raised.
-    fn least_squares(&self, rhs: &ArrayBase<D2, Ix2>) -> Result<LeastSquaresResult<E, Ix2>> {
+    fn least_squares(&self, rhs: &ArrayRef<E, Ix2>) -> Result<LeastSquaresResult<E, Ix2>> {
         let a = self.to_owned();
         let b = rhs.to_owned();
         a.least_squares_into(b)
@@ -255,11 +249,9 @@ where
 ///
 /// `E` is one of `f32`, `f64`, `c32`, `c64`. `D1`, `D2` can be any
 /// valid representation for `ArrayBase` (over `E`).
-impl<E, D1, D2> LeastSquaresSvdInPlace<D2, E, Ix1> for ArrayBase<D1, Ix2>
+impl<E> LeastSquaresSvdInPlace<E, Ix1> for ArrayRef<E, Ix2>
 where
     E: Scalar + Lapack,
-    D1: DataMut<Elem = E>,
-    D2: DataMut<Elem = E>,
 {
     /// Solve a least squares problem of the form `Ax = rhs`
     /// by calling `A.least_squares(rhs)`, where `rhs` is a
@@ -270,7 +262,7 @@ where
     /// `IncompatibleShape` error is raised.
     fn least_squares_in_place(
         &mut self,
-        rhs: &mut ArrayBase<D2, Ix1>,
+        rhs: &mut ArrayRef<E, Ix1>,
     ) -> Result<LeastSquaresResult<E, Ix1>> {
         if self.shape()[0] != rhs.shape()[0] {
             return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape).into());
@@ -288,14 +280,12 @@ where
     }
 }
 
-fn compute_least_squares_srhs<E, D1, D2>(
-    a: &mut ArrayBase<D1, Ix2>,
-    rhs: &mut ArrayBase<D2, Ix1>,
+fn compute_least_squares_srhs<E>(
+    a: &mut ArrayRef<E, Ix2>,
+    rhs: &mut ArrayRef<E, Ix1>,
 ) -> Result<LeastSquaresResult<E, Ix1>>
 where
     E: Scalar + Lapack,
-    D1: DataMut<Elem = E>,
-    D2: DataMut<Elem = E>,
 {
     let LeastSquaresOwned::<E> {
         singular_values,
@@ -318,11 +308,11 @@ where
     })
 }
 
-fn compute_residual_scalar<E: Scalar, D: Data<Elem = E>>(
+fn compute_residual_scalar<E: Scalar>(
     m: usize,
     n: usize,
     rank: i32,
-    b: &ArrayBase<D, Ix1>,
+    b: &ArrayRef<E, Ix1>,
 ) -> Option<Array<E::Real, Ix0>> {
     if m < n || n != rank as usize {
         return None;
@@ -338,11 +328,9 @@ fn compute_residual_scalar<E: Scalar, D: Data<Elem = E>>(
 ///
 /// `E` is one of `f32`, `f64`, `c32`, `c64`. `D1`, `D2` can be any
 /// valid representation for `ArrayBase` (over `E`).
-impl<E, D1, D2> LeastSquaresSvdInPlace<D2, E, Ix2> for ArrayBase<D1, Ix2>
+impl<E> LeastSquaresSvdInPlace<E, Ix2> for ArrayRef<E, Ix2>
 where
     E: Scalar + Lapack,
-    D1: DataMut<Elem = E>,
-    D2: DataMut<Elem = E>,
 {
     /// Solve a least squares problem of the form `Ax = rhs`
     /// by calling `A.least_squares(rhs)`, where `rhs` is a
@@ -353,7 +341,7 @@ where
     /// `IncompatibleShape` error is raised.
     fn least_squares_in_place(
         &mut self,
-        rhs: &mut ArrayBase<D2, Ix2>,
+        rhs: &mut ArrayRef<E, Ix2>,
     ) -> Result<LeastSquaresResult<E, Ix2>> {
         if self.shape()[0] != rhs.shape()[0] {
             return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape).into());
@@ -375,14 +363,12 @@ where
     }
 }
 
-fn compute_least_squares_nrhs<E, D1, D2>(
-    a: &mut ArrayBase<D1, Ix2>,
-    rhs: &mut ArrayBase<D2, Ix2>,
+fn compute_least_squares_nrhs<E>(
+    a: &mut ArrayRef<E, Ix2>,
+    rhs: &mut ArrayRef<E, Ix2>,
 ) -> Result<LeastSquaresResult<E, Ix2>>
 where
     E: Scalar + Lapack,
-    D1: DataMut<Elem = E>,
-    D2: DataMut<Elem = E>,
 {
     let a_layout = a.layout()?;
     let rhs_layout = rhs.layout()?;
@@ -408,11 +394,11 @@ where
     })
 }
 
-fn compute_residual_array1<E: Scalar, D: Data<Elem = E>>(
+fn compute_residual_array1<E: Scalar>(
     m: usize,
     n: usize,
     rank: i32,
-    b: &ArrayBase<D, Ix2>,
+    b: &ArrayRef<E, Ix2>,
 ) -> Option<Array1<E::Real>> {
     if m < n || n != rank as usize {
         return None;
diff --git a/ndarray-linalg/src/norm.rs b/ndarray-linalg/src/norm.rs
index e0dde1a4..ba28e73f 100644
--- a/ndarray-linalg/src/norm.rs
+++ b/ndarray-linalg/src/norm.rs
@@ -22,10 +22,9 @@ pub trait Norm {
     fn norm_max(&self) -> Self::Output;
 }
 
-impl<A, S, D> Norm for ArrayBase<S, D>
+impl<A, D> Norm for ArrayRef<A, D>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
     D: Dimension,
 {
     type Output = A::Real;
diff --git a/ndarray-linalg/src/operator.rs b/ndarray-linalg/src/operator.rs
index 368650bf..ae6f6cdc 100644
--- a/ndarray-linalg/src/operator.rs
+++ b/ndarray-linalg/src/operator.rs
@@ -10,20 +10,14 @@ pub trait LinearOperator {
     type Elem: Scalar;
 
     /// Apply operator out-place
-    fn apply<S>(&self, a: &ArrayBase<S, Ix1>) -> Array1<S::Elem>
-    where
-        S: Data<Elem = Self::Elem>,
-    {
+    fn apply(&self, a: &ArrayRef<Self::Elem, Ix1>) -> Array1<Self::Elem> {
         let mut a = a.to_owned();
         self.apply_mut(&mut a);
         a
     }
 
     /// Apply operator in-place
-    fn apply_mut<S>(&self, a: &mut ArrayBase<S, Ix1>)
-    where
-        S: DataMut<Elem = Self::Elem>,
-    {
+    fn apply_mut(&self, a: &mut ArrayRef<Self::Elem, Ix1>) {
         let b = self.apply(a);
         azip!((a in a, &b in &b) *a = b);
     }
@@ -38,19 +32,13 @@ pub trait LinearOperator {
     }
 
     /// Apply operator to matrix out-place
-    fn apply2<S>(&self, a: &ArrayBase<S, Ix2>) -> Array2<S::Elem>
-    where
-        S: Data<Elem = Self::Elem>,
-    {
+    fn apply2(&self, a: &ArrayRef<Self::Elem, Ix2>) -> Array2<Self::Elem> {
         let cols: Vec<_> = a.axis_iter(Axis(1)).map(|col| self.apply(&col)).collect();
         hstack(&cols).unwrap()
     }
 
     /// Apply operator to matrix in-place
-    fn apply2_mut<S>(&self, a: &mut ArrayBase<S, Ix2>)
-    where
-        S: DataMut<Elem = Self::Elem>,
-    {
+    fn apply2_mut(&self, a: &mut ArrayRef<Self::Elem, Ix2>) {
         for mut col in a.axis_iter_mut(Axis(1)) {
             self.apply_mut(&mut col)
         }
@@ -73,10 +61,7 @@ where
 {
     type Elem = A;
 
-    fn apply<S>(&self, a: &ArrayBase<S, Ix1>) -> Array1<A>
-    where
-        S: Data<Elem = A>,
-    {
+    fn apply(&self, a: &ArrayRef<A, Ix1>) -> Array1<A> {
         self.dot(a)
     }
 }
diff --git a/ndarray-linalg/src/opnorm.rs b/ndarray-linalg/src/opnorm.rs
index 0a412580..83b46b35 100644
--- a/ndarray-linalg/src/opnorm.rs
+++ b/ndarray-linalg/src/opnorm.rs
@@ -34,10 +34,9 @@ pub trait OperationNorm {
     }
 }
 
-impl<A, S> OperationNorm for ArrayBase<S, Ix2>
+impl<A> OperationNorm for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type Output = A::Real;
 
diff --git a/ndarray-linalg/src/solve.rs b/ndarray-linalg/src/solve.rs
index 79df1f77..e1c043ab 100644
--- a/ndarray-linalg/src/solve.rs
+++ b/ndarray-linalg/src/solve.rs
@@ -72,7 +72,7 @@ pub trait Solve<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of columns
     /// of `A`.
-    fn solve<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
+    fn solve(&self, b: &ArrayRef<A, Ix1>) -> Result<Array1<A>> {
         let mut b = replicate(b);
         self.solve_inplace(&mut b)?;
         Ok(b)
@@ -100,10 +100,7 @@ pub trait Solve<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of columns
     /// of `A`.
-    fn solve_inplace<'a, S: DataMut<Elem = A>>(
-        &self,
-        b: &'a mut ArrayBase<S, Ix1>,
-    ) -> Result<&'a mut ArrayBase<S, Ix1>>;
+    fn solve_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>>;
 
     /// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
     /// is the argument, and `x` is the successful result.
@@ -112,7 +109,7 @@ pub trait Solve<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of rows of
     /// `A`.
-    fn solve_t<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
+    fn solve_t(&self, b: &ArrayRef<A, Ix1>) -> Result<Array1<A>> {
         let mut b = replicate(b);
         self.solve_t_inplace(&mut b)?;
         Ok(b)
@@ -140,10 +137,7 @@ pub trait Solve<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of rows of
     /// `A`.
-    fn solve_t_inplace<'a, S: DataMut<Elem = A>>(
-        &self,
-        b: &'a mut ArrayBase<S, Ix1>,
-    ) -> Result<&'a mut ArrayBase<S, Ix1>>;
+    fn solve_t_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>>;
 
     /// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
     /// is the argument, and `x` is the successful result.
@@ -152,7 +146,7 @@ pub trait Solve<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of rows of
     /// `A`.
-    fn solve_h<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
+    fn solve_h(&self, b: &ArrayRef<A, Ix1>) -> Result<Array1<A>> {
         let mut b = replicate(b);
         self.solve_h_inplace(&mut b)?;
         Ok(b)
@@ -178,10 +172,7 @@ pub trait Solve<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of rows of
     /// `A`.
-    fn solve_h_inplace<'a, S: DataMut<Elem = A>>(
-        &self,
-        b: &'a mut ArrayBase<S, Ix1>,
-    ) -> Result<&'a mut ArrayBase<S, Ix1>>;
+    fn solve_h_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>>;
 }
 
 /// Represents the LU factorization of a matrix `A` as `A = P*L*U`.
@@ -199,13 +190,7 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A> + RawDataClone,
 {
-    fn solve_inplace<'a, Sb>(
-        &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+    fn solve_inplace<'a>(&self, rhs: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>> {
         assert_eq!(
             rhs.len(),
             self.a.len_of(Axis(1)),
@@ -220,13 +205,10 @@ where
         )?;
         Ok(rhs)
     }
-    fn solve_t_inplace<'a, Sb>(
+    fn solve_t_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix1>,
+    ) -> Result<&'a mut ArrayRef<A, Ix1>> {
         assert_eq!(
             rhs.len(),
             self.a.len_of(Axis(0)),
@@ -241,13 +223,10 @@ where
         )?;
         Ok(rhs)
     }
-    fn solve_h_inplace<'a, Sb>(
+    fn solve_h_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix1>,
+    ) -> Result<&'a mut ArrayRef<A, Ix1>> {
         assert_eq!(
             rhs.len(),
             self.a.len_of(Axis(0)),
@@ -264,38 +243,25 @@ where
     }
 }
 
-impl<A, S> Solve<A> for ArrayBase<S, Ix2>
+impl<A> Solve<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
-    fn solve_inplace<'a, Sb>(
-        &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+    fn solve_inplace<'a>(&self, rhs: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>> {
         let f = self.factorize()?;
         f.solve_inplace(rhs)
     }
-    fn solve_t_inplace<'a, Sb>(
+    fn solve_t_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix1>,
+    ) -> Result<&'a mut ArrayRef<A, Ix1>> {
         let f = self.factorize()?;
         f.solve_t_inplace(rhs)
     }
-    fn solve_h_inplace<'a, Sb>(
+    fn solve_h_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix1>,
+    ) -> Result<&'a mut ArrayRef<A, Ix1>> {
         let f = self.factorize()?;
         f.solve_h_inplace(rhs)
     }
@@ -326,10 +292,9 @@ where
     }
 }
 
-impl<A, Si> Factorize<OwnedRepr<A>> for ArrayBase<Si, Ix2>
+impl<A> Factorize<OwnedRepr<A>> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    Si: Data<Elem = A>,
 {
     fn factorize(&self) -> Result<LUFactorized<OwnedRepr<A>>> {
         let mut a: Array2<A> = replicate(self);
@@ -405,10 +370,9 @@ where
     }
 }
 
-impl<A, Si> Inverse for ArrayBase<Si, Ix2>
+impl<A> Inverse for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    Si: Data<Elem = A>,
 {
     type Output = Array2<A>;
 
@@ -520,10 +484,9 @@ where
     }
 }
 
-impl<A, S> Determinant<A> for ArrayBase<S, Ix2>
+impl<A> Determinant<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     fn sln_det(&self) -> Result<(A, A::Real)> {
         self.ensure_square()?;
@@ -612,10 +575,9 @@ where
     }
 }
 
-impl<A, S> ReciprocalConditionNum<A> for ArrayBase<S, Ix2>
+impl<A> ReciprocalConditionNum<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     fn rcond(&self) -> Result<A::Real> {
         self.factorize()?.rcond_into()
diff --git a/ndarray-linalg/src/solveh.rs b/ndarray-linalg/src/solveh.rs
index df0b3a6f..00dd6f92 100644
--- a/ndarray-linalg/src/solveh.rs
+++ b/ndarray-linalg/src/solveh.rs
@@ -65,7 +65,7 @@ pub trait SolveH<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of columns
     /// of `A`.
-    fn solveh<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
+    fn solveh(&self, b: &ArrayRef<A, Ix1>) -> Result<Array1<A>> {
         let mut b = replicate(b);
         self.solveh_inplace(&mut b)?;
         Ok(b)
@@ -96,10 +96,7 @@ pub trait SolveH<A: Scalar> {
     ///
     /// Panics if the length of `b` is not the equal to the number of columns
     /// of `A`.
-    fn solveh_inplace<'a, S: DataMut<Elem = A>>(
-        &self,
-        b: &'a mut ArrayBase<S, Ix1>,
-    ) -> Result<&'a mut ArrayBase<S, Ix1>>;
+    fn solveh_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>>;
 }
 
 /// Represents the Bunch–Kaufman factorization of a Hermitian (or real
@@ -114,13 +111,10 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A>,
 {
-    fn solveh_inplace<'a, Sb>(
+    fn solveh_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix1>,
+    ) -> Result<&'a mut ArrayRef<A, Ix1>> {
         assert_eq!(
             rhs.len(),
             self.a.len_of(Axis(1)),
@@ -137,18 +131,14 @@ where
     }
 }
 
-impl<A, S> SolveH<A> for ArrayBase<S, Ix2>
+impl<A> SolveH<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
-    fn solveh_inplace<'a, Sb>(
+    fn solveh_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix1>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix1>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix1>,
+    ) -> Result<&'a mut ArrayRef<A, Ix1>> {
         let f = self.factorizeh()?;
         f.solveh_inplace(rhs)
     }
@@ -181,10 +171,9 @@ where
     }
 }
 
-impl<A, Si> FactorizeH<OwnedRepr<A>> for ArrayBase<Si, Ix2>
+impl<A> FactorizeH<OwnedRepr<A>> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    Si: Data<Elem = A>,
 {
     fn factorizeh(&self) -> Result<BKFactorized<OwnedRepr<A>>> {
         let mut a: Array2<A> = replicate(self);
@@ -255,10 +244,9 @@ where
     }
 }
 
-impl<A, Si> InverseH for ArrayBase<Si, Ix2>
+impl<A> InverseH for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    Si: Data<Elem = A>,
 {
     type Output = Array2<A>;
 
@@ -317,10 +305,9 @@ pub trait DeterminantHInto {
 }
 
 /// Returns the sign and natural log of the determinant.
-fn bk_sln_det<P, S, A>(uplo: UPLO, ipiv_iter: P, a: &ArrayBase<S, Ix2>) -> (A::Real, A::Real)
+fn bk_sln_det<P, A>(uplo: UPLO, ipiv_iter: P, a: &ArrayRef<A, Ix2>) -> (A::Real, A::Real)
 where
     P: Iterator<Item = i32>,
-    S: Data<Elem = A>,
     A: Scalar + Lapack,
 {
     let layout = a.layout().unwrap();
@@ -424,10 +411,9 @@ where
     }
 }
 
-impl<A, S> DeterminantH for ArrayBase<S, Ix2>
+impl<A> DeterminantH for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type Elem = A;
 
diff --git a/ndarray-linalg/src/svd.rs b/ndarray-linalg/src/svd.rs
index 5c9d59b1..58584958 100644
--- a/ndarray-linalg/src/svd.rs
+++ b/ndarray-linalg/src/svd.rs
@@ -59,10 +59,9 @@ where
     }
 }
 
-impl<A, S> SVD for ArrayBase<S, Ix2>
+impl<A> SVD for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type U = Array2<A>;
     type VT = Array2<A>;
@@ -78,10 +77,9 @@ where
     }
 }
 
-impl<A, S> SVDInplace for ArrayBase<S, Ix2>
+impl<A> SVDInplace for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: DataMut<Elem = A>,
 {
     type U = Array2<A>;
     type VT = Array2<A>;
diff --git a/ndarray-linalg/src/svddc.rs b/ndarray-linalg/src/svddc.rs
index 0b0ae237..3a2aafcb 100644
--- a/ndarray-linalg/src/svddc.rs
+++ b/ndarray-linalg/src/svddc.rs
@@ -35,10 +35,9 @@ pub trait SVDDCInplace {
     ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
 }
 
-impl<A, S> SVDDC for ArrayBase<S, Ix2>
+impl<A> SVDDC for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     type U = Array2<A>;
     type VT = Array2<A>;
@@ -66,10 +65,9 @@ where
     }
 }
 
-impl<A, S> SVDDCInplace for ArrayBase<S, Ix2>
+impl<A> SVDDCInplace for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: DataMut<Elem = A>,
 {
     type U = Array2<A>;
     type VT = Array2<A>;
diff --git a/ndarray-linalg/src/trace.rs b/ndarray-linalg/src/trace.rs
index feb119f2..5692bee3 100644
--- a/ndarray-linalg/src/trace.rs
+++ b/ndarray-linalg/src/trace.rs
@@ -11,10 +11,9 @@ pub trait Trace {
     fn trace(&self) -> Result<Self::Output>;
 }
 
-impl<A, S> Trace for ArrayBase<S, Ix2>
+impl<A> Trace for ArrayRef<A, Ix2>
 where
     A: Scalar + Sum,
-    S: Data<Elem = A>,
 {
     type Output = A;
 
diff --git a/ndarray-linalg/src/triangular.rs b/ndarray-linalg/src/triangular.rs
index 2741d7b0..ded0c801 100644
--- a/ndarray-linalg/src/triangular.rs
+++ b/ndarray-linalg/src/triangular.rs
@@ -146,12 +146,11 @@ pub trait IntoTriangular<T> {
     fn into_triangular(self, uplo: UPLO) -> T;
 }
 
-impl<'a, A, S> IntoTriangular<&'a mut ArrayBase<S, Ix2>> for &'a mut ArrayBase<S, Ix2>
+impl<'a, A> IntoTriangular<&'a mut ArrayRef<A, Ix2>> for &'a mut ArrayRef<A, Ix2>
 where
     A: Zero,
-    S: DataMut<Elem = A>,
 {
-    fn into_triangular(self, uplo: UPLO) -> &'a mut ArrayBase<S, Ix2> {
+    fn into_triangular(self, uplo: UPLO) -> &'a mut ArrayRef<A, Ix2> {
         match uplo {
             UPLO::Upper => {
                 for ((i, j), val) in self.indexed_iter_mut() {
@@ -178,7 +177,7 @@ where
     S: DataMut<Elem = A>,
 {
     fn into_triangular(mut self, uplo: UPLO) -> ArrayBase<S, Ix2> {
-        (&mut self).into_triangular(uplo);
+        (&mut *self).into_triangular(uplo);
         self
     }
 }
diff --git a/ndarray-linalg/src/tridiagonal.rs b/ndarray-linalg/src/tridiagonal.rs
index b5aebbf2..9ff864e6 100644
--- a/ndarray-linalg/src/tridiagonal.rs
+++ b/ndarray-linalg/src/tridiagonal.rs
@@ -23,10 +23,9 @@ pub trait ExtractTridiagonal<A: Scalar> {
     fn extract_tridiagonal(&self) -> Result<Tridiagonal<A>>;
 }
 
-impl<A, S> ExtractTridiagonal<A> for ArrayBase<S, Ix2>
+impl<A> ExtractTridiagonal<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     fn extract_tridiagonal(&self) -> Result<Tridiagonal<A>> {
         let l = self.square_layout()?;
@@ -50,7 +49,7 @@ pub trait SolveTridiagonal<A: Scalar, D: Dimension> {
     /// Solves a system of linear equations `A * x = b` with tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result.
-    fn solve_tridiagonal<S: Data<Elem = A>>(&self, b: &ArrayBase<S, D>) -> Result<Array<A, D>>;
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, D>) -> Result<Array<A, D>>;
     /// Solves a system of linear equations `A * x = b` with tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result.
@@ -61,7 +60,7 @@ pub trait SolveTridiagonal<A: Scalar, D: Dimension> {
     /// Solves a system of linear equations `A^T * x = b` with tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result.
-    fn solve_t_tridiagonal<S: Data<Elem = A>>(&self, b: &ArrayBase<S, D>) -> Result<Array<A, D>>;
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, D>) -> Result<Array<A, D>>;
     /// Solves a system of linear equations `A^T * x = b` with tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result.
@@ -72,7 +71,7 @@ pub trait SolveTridiagonal<A: Scalar, D: Dimension> {
     /// Solves a system of linear equations `A^H * x = b` with tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result.
-    fn solve_h_tridiagonal<S: Data<Elem = A>>(&self, b: &ArrayBase<S, D>) -> Result<Array<A, D>>;
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, D>) -> Result<Array<A, D>>;
     /// Solves a system of linear equations `A^H * x = b` with tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result.
@@ -87,33 +86,33 @@ pub trait SolveTridiagonalInplace<A: Scalar, D: Dimension> {
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result. The value of `x` is also assigned to the
     /// argument.
-    fn solve_tridiagonal_inplace<'a, S: DataMut<Elem = A>>(
+    fn solve_tridiagonal_inplace<'a>(
         &self,
-        b: &'a mut ArrayBase<S, D>,
-    ) -> Result<&'a mut ArrayBase<S, D>>;
+        b: &'a mut ArrayRef<A, D>,
+    ) -> Result<&'a mut ArrayRef<A, D>>;
     /// Solves a system of linear equations `A^T * x = b` tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result. The value of `x` is also assigned to the
     /// argument.
-    fn solve_t_tridiagonal_inplace<'a, S: DataMut<Elem = A>>(
+    fn solve_t_tridiagonal_inplace<'a>(
         &self,
-        b: &'a mut ArrayBase<S, D>,
-    ) -> Result<&'a mut ArrayBase<S, D>>;
+        b: &'a mut ArrayRef<A, D>,
+    ) -> Result<&'a mut ArrayRef<A, D>>;
     /// Solves a system of linear equations `A^H * x = b` tridiagonal
     /// matrix `A`, where `A` is `self`, `b` is the argument, and
     /// `x` is the successful result. The value of `x` is also assigned to the
     /// argument.
-    fn solve_h_tridiagonal_inplace<'a, S: DataMut<Elem = A>>(
+    fn solve_h_tridiagonal_inplace<'a>(
         &self,
-        b: &'a mut ArrayBase<S, D>,
-    ) -> Result<&'a mut ArrayBase<S, D>>;
+        b: &'a mut ArrayRef<A, D>,
+    ) -> Result<&'a mut ArrayRef<A, D>>;
 }
 
 impl<A> SolveTridiagonal<A, Ix2> for LUFactorizedTridiagonal<A>
 where
     A: Scalar + Lapack,
 {
-    fn solve_tridiagonal<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix2>) -> Result<Array<A, Ix2>> {
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -125,10 +124,7 @@ where
         self.solve_tridiagonal_inplace(&mut b)?;
         Ok(b)
     }
-    fn solve_t_tridiagonal<S: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<S, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_t_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -140,10 +136,7 @@ where
         self.solve_t_tridiagonal_inplace(&mut b)?;
         Ok(b)
     }
-    fn solve_h_tridiagonal<S: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<S, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_h_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -161,10 +154,7 @@ impl<A> SolveTridiagonal<A, Ix2> for Tridiagonal<A>
 where
     A: Scalar + Lapack,
 {
-    fn solve_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -176,10 +166,7 @@ where
         self.solve_tridiagonal_inplace(&mut b)?;
         Ok(b)
     }
-    fn solve_t_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_t_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -191,10 +178,7 @@ where
         self.solve_t_tridiagonal_inplace(&mut b)?;
         Ok(b)
     }
-    fn solve_h_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_h_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -208,15 +192,11 @@ where
     }
 }
 
-impl<A, S> SolveTridiagonal<A, Ix2> for ArrayBase<S, Ix2>
+impl<A> SolveTridiagonal<A, Ix2> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
-    fn solve_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -228,10 +208,7 @@ where
         self.solve_tridiagonal_inplace(&mut b)?;
         Ok(b)
     }
-    fn solve_t_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_t_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -243,10 +220,7 @@ where
         self.solve_t_tridiagonal_inplace(&mut b)?;
         Ok(b)
     }
-    fn solve_h_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix2>,
-    ) -> Result<Array<A, Ix2>> {
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, Ix2>) -> Result<Array<A, Ix2>> {
         let mut b = replicate(b);
         self.solve_h_tridiagonal_inplace(&mut b)?;
         Ok(b)
@@ -264,13 +238,10 @@ impl<A> SolveTridiagonalInplace<A, Ix2> for LUFactorizedTridiagonal<A>
 where
     A: Scalar + Lapack,
 {
-    fn solve_tridiagonal_inplace<'a, Sb>(
+    fn solve_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         A::solve_tridiagonal(
             self,
             rhs.layout()?,
@@ -279,13 +250,10 @@ where
         )?;
         Ok(rhs)
     }
-    fn solve_t_tridiagonal_inplace<'a, Sb>(
+    fn solve_t_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         A::solve_tridiagonal(
             self,
             rhs.layout()?,
@@ -294,13 +262,10 @@ where
         )?;
         Ok(rhs)
     }
-    fn solve_h_tridiagonal_inplace<'a, Sb>(
+    fn solve_h_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         A::solve_tridiagonal(
             self,
             rhs.layout()?,
@@ -315,70 +280,51 @@ impl<A> SolveTridiagonalInplace<A, Ix2> for Tridiagonal<A>
 where
     A: Scalar + Lapack,
 {
-    fn solve_tridiagonal_inplace<'a, Sb>(
+    fn solve_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         let f = self.factorize_tridiagonal()?;
         f.solve_tridiagonal_inplace(rhs)
     }
-    fn solve_t_tridiagonal_inplace<'a, Sb>(
+    fn solve_t_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         let f = self.factorize_tridiagonal()?;
         f.solve_t_tridiagonal_inplace(rhs)
     }
-    fn solve_h_tridiagonal_inplace<'a, Sb>(
+    fn solve_h_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         let f = self.factorize_tridiagonal()?;
         f.solve_h_tridiagonal_inplace(rhs)
     }
 }
 
-impl<A, S> SolveTridiagonalInplace<A, Ix2> for ArrayBase<S, Ix2>
+impl<A> SolveTridiagonalInplace<A, Ix2> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
-    fn solve_tridiagonal_inplace<'a, Sb>(
+    fn solve_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         let f = self.factorize_tridiagonal()?;
         f.solve_tridiagonal_inplace(rhs)
     }
-    fn solve_t_tridiagonal_inplace<'a, Sb>(
+    fn solve_t_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         let f = self.factorize_tridiagonal()?;
         f.solve_t_tridiagonal_inplace(rhs)
     }
-    fn solve_h_tridiagonal_inplace<'a, Sb>(
+    fn solve_h_tridiagonal_inplace<'a>(
         &self,
-        rhs: &'a mut ArrayBase<Sb, Ix2>,
-    ) -> Result<&'a mut ArrayBase<Sb, Ix2>>
-    where
-        Sb: DataMut<Elem = A>,
-    {
+        rhs: &'a mut ArrayRef<A, Ix2>,
+    ) -> Result<&'a mut ArrayRef<A, Ix2>> {
         let f = self.factorize_tridiagonal()?;
         f.solve_h_tridiagonal_inplace(rhs)
     }
@@ -388,7 +334,7 @@ impl<A> SolveTridiagonal<A, Ix1> for LUFactorizedTridiagonal<A>
 where
     A: Scalar + Lapack,
 {
-    fn solve_tridiagonal<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array<A, Ix1>> {
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_tridiagonal_into(b)
     }
@@ -400,10 +346,7 @@ where
         let b = self.solve_tridiagonal_into(b)?;
         Ok(flatten(b))
     }
-    fn solve_t_tridiagonal<S: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<S, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_t_tridiagonal_into(b)
     }
@@ -415,10 +358,7 @@ where
         let b = self.solve_t_tridiagonal_into(b)?;
         Ok(flatten(b))
     }
-    fn solve_h_tridiagonal<S: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<S, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_h_tridiagonal_into(b)
     }
@@ -436,10 +376,7 @@ impl<A> SolveTridiagonal<A, Ix1> for Tridiagonal<A>
 where
     A: Scalar + Lapack,
 {
-    fn solve_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_tridiagonal_into(b)
     }
@@ -452,10 +389,7 @@ where
         let b = f.solve_tridiagonal_into(b)?;
         Ok(flatten(b))
     }
-    fn solve_t_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_t_tridiagonal_into(b)
     }
@@ -468,10 +402,7 @@ where
         let b = f.solve_t_tridiagonal_into(b)?;
         Ok(flatten(b))
     }
-    fn solve_h_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_h_tridiagonal_into(b)
     }
@@ -486,15 +417,11 @@ where
     }
 }
 
-impl<A, S> SolveTridiagonal<A, Ix1> for ArrayBase<S, Ix2>
+impl<A> SolveTridiagonal<A, Ix1> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
-    fn solve_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_tridiagonal_into(b)
     }
@@ -507,10 +434,7 @@ where
         let b = f.solve_tridiagonal_into(b)?;
         Ok(flatten(b))
     }
-    fn solve_t_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_t_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_t_tridiagonal_into(b)
     }
@@ -523,10 +447,7 @@ where
         let b = f.solve_t_tridiagonal_into(b)?;
         Ok(flatten(b))
     }
-    fn solve_h_tridiagonal<Sb: Data<Elem = A>>(
-        &self,
-        b: &ArrayBase<Sb, Ix1>,
-    ) -> Result<Array<A, Ix1>> {
+    fn solve_h_tridiagonal(&self, b: &ArrayRef<A, Ix1>) -> Result<Array<A, Ix1>> {
         let b = b.to_owned();
         self.solve_h_tridiagonal_into(b)
     }
@@ -574,10 +495,9 @@ where
     }
 }
 
-impl<A, S> FactorizeTridiagonal<A> for ArrayBase<S, Ix2>
+impl<A> FactorizeTridiagonal<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     fn factorize_tridiagonal(&self) -> Result<LUFactorizedTridiagonal<A>> {
         let a = self.extract_tridiagonal()?;
@@ -626,10 +546,9 @@ where
     }
 }
 
-impl<A, S> DeterminantTridiagonal<A> for ArrayBase<S, Ix2>
+impl<A> DeterminantTridiagonal<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     fn det_tridiagonal(&self) -> Result<A> {
         let tridiag = self.extract_tridiagonal()?;
@@ -684,10 +603,9 @@ where
     }
 }
 
-impl<A, S> ReciprocalConditionNumTridiagonal<A> for ArrayBase<S, Ix2>
+impl<A> ReciprocalConditionNumTridiagonal<A> for ArrayRef<A, Ix2>
 where
     A: Scalar + Lapack,
-    S: Data<Elem = A>,
 {
     fn rcond_tridiagonal(&self) -> Result<A::Real> {
         self.factorize_tridiagonal()?.rcond_tridiagonal_into()