diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index 51efffd1..8f7d2d7e 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -68,21 +68,22 @@ extern crate openblas_src as _src;
 #[cfg(any(feature = "netlib-system", feature = "netlib-static"))]
 extern crate netlib_src as _src;
 
-pub mod cholesky;
-pub mod eig;
-pub mod eigh;
 pub mod error;
 pub mod layout;
-pub mod least_squares;
-pub mod opnorm;
-pub mod qr;
-pub mod rcond;
-pub mod solve;
-pub mod solveh;
-pub mod svd;
-pub mod svddc;
-pub mod triangular;
-pub mod tridiagonal;
+
+mod cholesky;
+mod eig;
+mod eigh;
+mod least_squares;
+mod opnorm;
+mod qr;
+mod rcond;
+mod solve;
+mod solveh;
+mod svd;
+mod svddc;
+mod triangular;
+mod tridiagonal;
 
 pub use self::cholesky::*;
 pub use self::eig::*;
@@ -159,7 +160,7 @@ pub enum NormType {
 }
 
 impl NormType {
-    pub(crate) fn transpose(self) -> Self {
+    pub fn transpose(self) -> Self {
         match self {
             NormType::One => NormType::Infinity,
             NormType::Infinity => NormType::One,
diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs
index 0c594d92..32bec9fa 100644
--- a/lax/src/opnorm.rs
+++ b/lax/src/opnorm.rs
@@ -1,11 +1,10 @@
 //! Operator norms of matrices
 
+use super::NormType;
 use crate::layout::MatrixLayout;
 use cauchy::*;
 use num_traits::Zero;
 
-pub use super::NormType;
-
 pub trait OperatorNorm_: Scalar {
     fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real;
 }
diff --git a/ndarray-linalg/src/cholesky.rs b/ndarray-linalg/src/cholesky.rs
index 3f445305..5dd40cae 100644
--- a/ndarray-linalg/src/cholesky.rs
+++ b/ndarray-linalg/src/cholesky.rs
@@ -52,7 +52,7 @@ use crate::layout::*;
 use crate::triangular::IntoTriangular;
 use crate::types::*;
 
-pub use crate::lapack::UPLO;
+pub use lax::UPLO;
 
 /// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix
 pub struct CholeskyFactorized<S: Data> {
diff --git a/ndarray-linalg/src/convert.rs b/ndarray-linalg/src/convert.rs
index 1742ac4d..e1446e96 100644
--- a/ndarray-linalg/src/convert.rs
+++ b/ndarray-linalg/src/convert.rs
@@ -1,9 +1,9 @@
 //! utilities for convert array
 
+use lax::UPLO;
 use ndarray::*;
 
 use super::error::*;
-use super::lapack::UPLO;
 use super::layout::*;
 use super::types::*;
 
diff --git a/ndarray-linalg/src/error.rs b/ndarray-linalg/src/error.rs
index 893a7a85..c073b985 100644
--- a/ndarray-linalg/src/error.rs
+++ b/ndarray-linalg/src/error.rs
@@ -14,7 +14,7 @@ pub enum LinalgError {
 
     /// LAPACK subroutine returns non-zero code
     #[error(transparent)]
-    Lapack(#[from] lapack::error::Error),
+    Lapack(#[from] lax::error::Error),
 
     /// Strides of the array is not supported
     #[error("invalid stride: s0={}, s1={}", s0, s1)]
diff --git a/ndarray-linalg/src/layout.rs b/ndarray-linalg/src/layout.rs
index cadbcf2d..9ca772bb 100644
--- a/ndarray-linalg/src/layout.rs
+++ b/ndarray-linalg/src/layout.rs
@@ -3,7 +3,7 @@
 use super::error::*;
 use ndarray::*;
 
-pub use lapack::layout::MatrixLayout;
+pub use lax::layout::MatrixLayout;
 
 pub trait AllocatedArray {
     type Elem;
diff --git a/ndarray-linalg/src/least_squares.rs b/ndarray-linalg/src/least_squares.rs
index 0ff518ad..03583a25 100644
--- a/ndarray-linalg/src/least_squares.rs
+++ b/ndarray-linalg/src/least_squares.rs
@@ -60,10 +60,10 @@
 //! // `a` and `b` have been moved, no longer valid
 //! ```
 
+use lax::*;
 use ndarray::*;
 
 use crate::error::*;
-use crate::lapack::least_squares::*;
 use crate::layout::*;
 use crate::types::*;
 
@@ -253,7 +253,7 @@ where
 /// valid representation for `ArrayBase`.
 impl<E, D> LeastSquaresSvdInPlace<D, E, Ix1> for ArrayBase<D, Ix2>
 where
-    E: Scalar + Lapack + LeastSquaresSvdDivideConquer_,
+    E: Scalar + Lapack,
     D: DataMut<Elem = E>,
 {
     /// Solve a least squares problem of the form `Ax = rhs`
diff --git a/ndarray-linalg/src/lib.rs b/ndarray-linalg/src/lib.rs
index 33fbacef..ba9d10c5 100644
--- a/ndarray-linalg/src/lib.rs
+++ b/ndarray-linalg/src/lib.rs
@@ -47,7 +47,6 @@
 
 #[macro_use]
 extern crate ndarray;
-extern crate lax as lapack;
 
 pub mod assert;
 pub mod cholesky;
diff --git a/ndarray-linalg/src/lobpcg/eig.rs b/ndarray-linalg/src/lobpcg/eig.rs
index f5f3aa3f..b0b45009 100644
--- a/ndarray-linalg/src/lobpcg/eig.rs
+++ b/ndarray-linalg/src/lobpcg/eig.rs
@@ -1,5 +1,7 @@
 use super::lobpcg::{lobpcg, LobpcgResult, Order};
-use crate::{generate, Lapack, Scalar};
+use crate::{generate, Scalar};
+use lax::Lapack;
+
 ///! Implements truncated eigenvalue decomposition
 ///
 use ndarray::prelude::*;
diff --git a/ndarray-linalg/src/lobpcg/lobpcg.rs b/ndarray-linalg/src/lobpcg/lobpcg.rs
index c942a743..91e6f5b4 100644
--- a/ndarray-linalg/src/lobpcg/lobpcg.rs
+++ b/ndarray-linalg/src/lobpcg/lobpcg.rs
@@ -4,7 +4,8 @@
 ///which can be used as a solver for large symmetric positive definite eigenproblems.
 use crate::error::{LinalgError, Result};
 use crate::{cholesky::*, close_l2, eigh::*, norm::*, triangular::*};
-use crate::{Lapack, Scalar};
+use cauchy::Scalar;
+use lax::Lapack;
 use ndarray::prelude::*;
 use ndarray::{Data, OwnedRepr, ScalarOperand};
 use num_traits::{Float, NumCast};
@@ -338,7 +339,7 @@ pub fn lobpcg<
         let result = p_ap
             .as_ref()
             .ok_or(LinalgError::Lapack(
-                lapack::error::Error::LapackComputationalFailure { return_code: 1 },
+                lax::error::Error::LapackComputationalFailure { return_code: 1 },
             ))
             .and_then(|(active_p, active_ap)| {
                 let xap = x.t().dot(active_ap);
diff --git a/ndarray-linalg/src/lobpcg/svd.rs b/ndarray-linalg/src/lobpcg/svd.rs
index d857386f..a796364d 100644
--- a/ndarray-linalg/src/lobpcg/svd.rs
+++ b/ndarray-linalg/src/lobpcg/svd.rs
@@ -3,7 +3,9 @@
 ///! This module computes the k largest/smallest singular values/vectors for a dense matrix.
 use super::lobpcg::{lobpcg, LobpcgResult, Order};
 use crate::error::Result;
-use crate::{generate, Lapack, Scalar};
+use crate::generate;
+use cauchy::Scalar;
+use lax::Lapack;
 use ndarray::prelude::*;
 use ndarray::ScalarOperand;
 use num_traits::{Float, NumCast};
diff --git a/ndarray-linalg/src/opnorm.rs b/ndarray-linalg/src/opnorm.rs
index b4956867..c0621aa3 100644
--- a/ndarray-linalg/src/opnorm.rs
+++ b/ndarray-linalg/src/opnorm.rs
@@ -1,14 +1,14 @@
 //! Operator norm
 
+use lax::Tridiagonal;
 use ndarray::*;
 
 use crate::convert::*;
 use crate::error::*;
 use crate::layout::*;
-use crate::tridiagonal::Tridiagonal;
 use crate::types::*;
 
-pub use crate::lapack::NormType;
+pub use lax::NormType;
 
 /// Operator norm using `*lange` LAPACK routines
 ///
diff --git a/ndarray-linalg/src/qr.rs b/ndarray-linalg/src/qr.rs
index ae7b2c25..441e7cef 100644
--- a/ndarray-linalg/src/qr.rs
+++ b/ndarray-linalg/src/qr.rs
@@ -11,7 +11,7 @@ use crate::layout::*;
 use crate::triangular::*;
 use crate::types::*;
 
-pub use crate::lapack::UPLO;
+pub use lax::UPLO;
 
 /// QR decomposition for matrix reference
 ///
diff --git a/ndarray-linalg/src/solve.rs b/ndarray-linalg/src/solve.rs
index fd4b3017..e58b27a0 100644
--- a/ndarray-linalg/src/solve.rs
+++ b/ndarray-linalg/src/solve.rs
@@ -55,7 +55,7 @@ use crate::layout::*;
 use crate::opnorm::OperationNorm;
 use crate::types::*;
 
-pub use crate::lapack::{Pivot, Transpose};
+pub use lax::{Pivot, Transpose};
 
 /// An interface for solving systems of linear equations.
 ///
@@ -468,7 +468,7 @@ where
         self.ensure_square()?;
         match self.factorize() {
             Ok(fac) => fac.sln_det(),
-            Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure {..}) =>
+            Err(LinalgError::Lapack(e)) if matches!(e, lax::error::Error::LapackComputationalFailure {..}) =>
             {
                 // The determinant is zero.
                 Ok((A::zero(), A::Real::neg_infinity()))
@@ -487,7 +487,7 @@ where
         self.ensure_square()?;
         match self.factorize_into() {
             Ok(fac) => fac.sln_det_into(),
-            Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure { .. }) =>
+            Err(LinalgError::Lapack(e)) if matches!(e, lax::error::Error::LapackComputationalFailure { .. }) =>
             {
                 // The determinant is zero.
                 Ok((A::zero(), A::Real::neg_infinity()))
diff --git a/ndarray-linalg/src/solveh.rs b/ndarray-linalg/src/solveh.rs
index 102158aa..04c16fec 100644
--- a/ndarray-linalg/src/solveh.rs
+++ b/ndarray-linalg/src/solveh.rs
@@ -57,7 +57,7 @@ use crate::error::*;
 use crate::layout::*;
 use crate::types::*;
 
-pub use crate::lapack::{Pivot, UPLO};
+pub use lax::{Pivot, UPLO};
 
 /// An interface for solving systems of Hermitian (or real symmetric) linear equations.
 ///
@@ -426,7 +426,7 @@ where
     fn sln_deth(&self) -> Result<(A::Real, A::Real)> {
         match self.factorizeh() {
             Ok(fac) => Ok(fac.sln_deth()),
-            Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure {..}) =>
+            Err(LinalgError::Lapack(e)) if matches!(e, lax::error::Error::LapackComputationalFailure {..}) =>
             {
                 // Determinant is zero.
                 Ok((A::Real::zero(), A::Real::neg_infinity()))
@@ -451,7 +451,7 @@ where
     fn sln_deth_into(self) -> Result<(A::Real, A::Real)> {
         match self.factorizeh_into() {
             Ok(fac) => Ok(fac.sln_deth_into()),
-            Err(LinalgError::Lapack(e)) if matches!(e, lapack::error::Error::LapackComputationalFailure {..}) =>
+            Err(LinalgError::Lapack(e)) if matches!(e, lax::error::Error::LapackComputationalFailure {..}) =>
             {
                 // Determinant is zero.
                 Ok((A::Real::zero(), A::Real::neg_infinity()))
diff --git a/ndarray-linalg/src/svddc.rs b/ndarray-linalg/src/svddc.rs
index c79fa0e2..017354c2 100644
--- a/ndarray-linalg/src/svddc.rs
+++ b/ndarray-linalg/src/svddc.rs
@@ -3,7 +3,7 @@
 use super::{convert::*, error::*, layout::*, types::*};
 use ndarray::*;
 
-pub use lapack::svddc::UVTFlag;
+pub use lax::UVTFlag;
 
 /// Singular-value decomposition of matrix (copying) by divide-and-conquer
 pub trait SVDDC {
diff --git a/ndarray-linalg/src/triangular.rs b/ndarray-linalg/src/triangular.rs
index b378dcd2..2741d7b0 100644
--- a/ndarray-linalg/src/triangular.rs
+++ b/ndarray-linalg/src/triangular.rs
@@ -1,15 +1,15 @@
 //! Methods for triangular matrices
 
+use lax::*;
 use ndarray::*;
 use num_traits::Zero;
 
 use super::convert::*;
 use super::error::*;
-use super::lapack::*;
 use super::layout::*;
 use super::types::*;
 
-pub use super::lapack::Diag;
+pub use lax::Diag;
 
 /// solve a triangular system with upper triangular matrix
 pub trait SolveTriangular<A, S, D>
diff --git a/ndarray-linalg/src/tridiagonal.rs b/ndarray-linalg/src/tridiagonal.rs
index 92eae05b..b603a5ee 100644
--- a/ndarray-linalg/src/tridiagonal.rs
+++ b/ndarray-linalg/src/tridiagonal.rs
@@ -4,13 +4,13 @@
 
 use super::convert::*;
 use super::error::*;
-use super::lapack::*;
 use super::layout::*;
 use cauchy::Scalar;
+use lax::*;
 use ndarray::*;
 use num_traits::One;
 
-pub use lapack::tridiagonal::{LUFactorizedTridiagonal, Tridiagonal};
+pub use lax::{LUFactorizedTridiagonal, Tridiagonal};
 
 /// An interface for making a Tridiagonal struct.
 pub trait ExtractTridiagonal<A: Scalar> {
diff --git a/ndarray-linalg/src/types.rs b/ndarray-linalg/src/types.rs
index d42f6e24..57d947bc 100644
--- a/ndarray-linalg/src/types.rs
+++ b/ndarray-linalg/src/types.rs
@@ -1,7 +1,7 @@
 //! Basic types and their methods for linear algebra
 
-pub use super::lapack::Lapack;
 pub use cauchy::Scalar;
+pub use lax::Lapack;
 
 pub use num_complex::Complex32 as c32;
 pub use num_complex::Complex64 as c64;