diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml
index f8509459..4123d4a6 100644
--- a/.github/workflows/gh-pages.yml
+++ b/.github/workflows/gh-pages.yml
@@ -33,23 +33,21 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
 
       # Generate cargo-doc
-      - uses: actions-rs/cargo@v1
-        with:
-          command: doc
-          args: --no-deps
+      - name: Generate documentation
+        run: cargo doc --no-deps
 
       - name: Setup Pages
-        uses: actions/configure-pages@v2
+        uses: actions/configure-pages@v5
 
       # Upload target/doc directory
       - name: Upload artifact
-        uses: actions/upload-pages-artifact@v1
+        uses: actions/upload-pages-artifact@v3
         with:
           path: 'target/doc'
 
       - name: Deploy to GitHub Pages
         id: deployment
-        uses: actions/deploy-pages@v1
+        uses: actions/deploy-pages@v4
diff --git a/.github/workflows/intel-mkl.yml b/.github/workflows/intel-mkl.yml
index fcb8c845..14a0b443 100644
--- a/.github/workflows/intel-mkl.yml
+++ b/.github/workflows/intel-mkl.yml
@@ -5,44 +5,16 @@ on:
     branches:
       - master
   pull_request: {}
+  workflow_dispatch:
 
 jobs:
-  windows:
-    runs-on: windows-2019
+  intel-mkl:
+    strategy:
+      fail-fast: false
+      matrix:
+        system: [ubuntu-22.04, windows-latest]
+    runs-on: ${{ matrix.system }}
     steps:
-      - uses: actions/checkout@v1
-      - uses: actions-rs/cargo@v1
-        with:
-          command: test
-          args: >
-            --manifest-path=ndarray-linalg/Cargo.toml
-            --no-default-features
-            --features=intel-mkl-static
-
-  linux:
-    runs-on: ubuntu-22.04
-    steps:
-      - uses: actions/checkout@v1
-      - uses: actions-rs/cargo@v1
-        name: cargo test
-        with:
-          command: test
-          args: >
-            --manifest-path=ndarray-linalg/Cargo.toml
-            --no-default-features
-            --features=intel-mkl-static
-
-  linux-container:
-    runs-on: ubuntu-22.04
-    container:
-      image: ghcr.io/rust-math/rust-mkl:1.63.0-2020.1
-    steps:
-      - uses: actions/checkout@v1
-      - uses: actions-rs/cargo@v1
-        name: cargo test
-        with:
-          command: test
-          args: >
-            --manifest-path=ndarray-linalg/Cargo.toml
-            --no-default-features
-            --features=intel-mkl-system
+      - uses: actions/checkout@v4
+      - name: cargo test
+        run: cargo test --manifest-path=ndarray-linalg/Cargo.toml --no-default-features --features=intel-mkl-static  --verbose
diff --git a/.github/workflows/netlib.yml b/.github/workflows/netlib.yml
index d4df49a5..70048503 100644
--- a/.github/workflows/netlib.yml
+++ b/.github/workflows/netlib.yml
@@ -5,6 +5,7 @@ on:
     branches:
       - master
   pull_request: {}
+  workflow_dispatch:
 
 jobs:
   linux:
@@ -15,15 +16,10 @@ jobs:
           - static
     runs-on: ubuntu-22.04
     steps:
-      - uses: actions/checkout@v1
+      - uses: actions/checkout@v4
       - name: apt install gfortran
         run: |
           sudo apt update
           sudo apt install -y gfortran
-      - uses: actions-rs/cargo@v1
-        with:
-          command: test
-          args: >
-            --manifest-path=ndarray-linalg/Cargo.toml
-            --no-default-features
-            --features=netlib-${{ matrix.feature }}
+      - name: cargo test
+        run: cargo test --manifest-path=ndarray-linalg/Cargo.toml --no-default-features --features=netlib-${{ matrix.feature }}
diff --git a/.github/workflows/openblas.yml b/.github/workflows/openblas.yml
index 4e717ba9..5375be70 100644
--- a/.github/workflows/openblas.yml
+++ b/.github/workflows/openblas.yml
@@ -5,12 +5,11 @@ on:
     branches:
       - master
   pull_request: {}
+  workflow_dispatch:
 
 jobs:
   linux:
     runs-on: ubuntu-22.04
-    container:
-      image: rust
     strategy:
       fail-fast: false
       matrix:
@@ -18,20 +17,15 @@ jobs:
           - static
           - system
     steps:
-      - uses: actions/checkout@v1
+      - uses: actions/checkout@v4
       - name: apt install gfortran
         run: |
-          apt update
-          apt install -y gfortran
+          sudo apt update
+          sudo apt install -y gfortran
       - name: Install OpenBLAS by apt
         run: |
-          apt update
-          apt install -y libopenblas-dev
+          sudo apt update
+          sudo apt install -y libopenblas-dev
         if: ${{ contains(matrix.feature, 'system') }}
-      - uses: actions-rs/cargo@v1
-        with:
-          command: test
-          args: >
-            --manifest-path=ndarray-linalg/Cargo.toml
-            --no-default-features
-            --features=openblas-${{ matrix.feature }}
+      - name: cargo test
+        run: cargo test --manifest-path=ndarray-linalg/Cargo.toml --no-default-features --features=openblas-${{ matrix.feature }} 
diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml
index 9cd3919b..66ca1553 100644
--- a/.github/workflows/rust.yml
+++ b/.github/workflows/rust.yml
@@ -5,52 +5,50 @@ on:
     branches:
       - master
   pull_request: {}
+  workflow_dispatch:
 
 jobs:
   check-format:
     runs-on: ubuntu-22.04
     steps:
-    - uses: actions/checkout@v1
-    - uses: actions-rs/cargo@v1
-      with:
-        command: fmt
-        args: -- --check
+    - uses: actions/checkout@v4
+    - name: fmt
+      run: cargo fmt -- --check
 
   check:
     runs-on: ubuntu-22.04
     steps:
-    - uses: actions/checkout@v1
-    - uses: actions-rs/cargo@v1
-      with:
-        command: check
-        args: --all-targets
+    - uses: actions/checkout@v4
+    - name: cargo check
+      run: cargo check --all-targets
 
   check-doc:
     runs-on: ubuntu-22.04
     steps:
-    - uses: actions/checkout@v1
-    - uses: actions-rs/cargo@v1
-      with:
-        command: doc
-        args: --no-deps
+    - uses: actions/checkout@v4
+    - name: cargo doc
+      run: cargo doc --no-deps
 
   clippy:
     runs-on: ubuntu-22.04
     steps:
-    - uses: actions/checkout@v1
-    - uses: actions-rs/cargo@v1
-      with:
-        command: clippy
+    - uses: actions/checkout@v4
+    - name: cargo clippy
+      run: cargo clippy
 
   coverage:
     runs-on: ubuntu-22.04
     container:
-      image: ghcr.io/rust-math/rust-mkl:1.63.0-2020.1
+      image: xd009642/tarpaulin:develop-nightly
       options: --security-opt seccomp=unconfined
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
+      - uses: dtolnay/rust-toolchain@nightly
+      - name: Install Cross
+        uses: taiki-e/install-action@v2
+        with:
+          tool: cargo-tarpaulin
       - name: Generate code coverage
-        run: |
-          cargo tarpaulin --features=intel-mkl-system --out Xml
+        run: cargo +nightly tarpaulin --features=intel-mkl-static --out xml
       - name: Upload to codecov.io
-        uses: codecov/codecov-action@v1
+        uses: codecov/codecov-action@v5
diff --git a/README.md b/README.md
index 530cc589..85b1428c 100644
--- a/README.md
+++ b/README.md
@@ -62,7 +62,7 @@ Supported features are following:
 
 ### For library developer
 
-If you creating a library depending on this crate, we encourage you not to link any backend:
+If you are creating a library depending on this crate, we encourage you not to link any backend:
 
 ```toml
 [dependencies]
diff --git a/lax/Cargo.toml b/lax/Cargo.toml
index 2f095ea8..076d5cf2 100644
--- a/lax/Cargo.toml
+++ b/lax/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "lax"
-version = "0.16.0"
+version = "0.17.0"
 authors = ["Toshiki Teramura <toshiki.teramura@gmail.com>"]
 edition = "2018"
 
@@ -29,10 +29,10 @@ intel-mkl-static = ["intel-mkl-src/mkl-static-lp64-seq"]
 intel-mkl-system = ["intel-mkl-src/mkl-dynamic-lp64-seq"]
 
 [dependencies]
-thiserror = "1.0.24"
+thiserror = "2.0.0"
 cauchy = "0.4.0"
 num-traits = "0.2.14"
-lapack-sys = "0.14.0"
+lapack-sys = "0.15.0"
 katexit = "0.1.2"
 
 [dependencies.intel-mkl-src]
@@ -41,7 +41,7 @@ default-features = false
 optional = true
 
 [dependencies.netlib-src]
-version = "0.8.0"
+version = "0.9.0"
 optional = true
 features = ["cblas"]
 default-features = false
diff --git a/lax/src/eig.rs b/lax/src/eig.rs
index 710beb9c..f02035bb 100644
--- a/lax/src/eig.rs
+++ b/lax/src/eig.rs
@@ -404,7 +404,7 @@ impl_eig_work_r!(f64, lapack_sys::dgeev_);
 ///
 /// In the C-layout case, we need the conjugates of the left
 /// eigenvectors, so the signs should be reversed.
-fn reconstruct_eigenvectors<T: Scalar>(
+pub(crate) fn reconstruct_eigenvectors<T: Scalar>(
     take_hermite_conjugate: bool,
     eig_im: &[T],
     vr: &[T],
diff --git a/lax/src/eig_generalized.rs b/lax/src/eig_generalized.rs
new file mode 100644
index 00000000..ea99dbdb
--- /dev/null
+++ b/lax/src/eig_generalized.rs
@@ -0,0 +1,520 @@
+//! Generalized eigenvalue problem for general matrices
+//!
+//! LAPACK correspondance
+//! ----------------------
+//!
+//! | f32   | f64   | c32   | c64   |
+//! |:------|:------|:------|:------|
+//! | sggev | dggev | cggev | zggev |
+//!
+use std::mem::MaybeUninit;
+
+use crate::eig::reconstruct_eigenvectors;
+use crate::{error::*, layout::MatrixLayout, *};
+use cauchy::*;
+use num_traits::{ToPrimitive, Zero};
+
+#[derive(Clone, PartialEq, Eq)]
+pub enum GeneralizedEigenvalue<T: Scalar> {
+    /// Finite generalized eigenvalue: `Finite(α/β, (α, β))`
+    Finite(T, (T, T)),
+
+    /// Indeterminate generalized eigenvalue: `Indeterminate((α, β))`
+    Indeterminate((T, T)),
+}
+
+impl<T: Scalar> std::fmt::Display for GeneralizedEigenvalue<T> {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::Finite(e, (a, b)) => write!(f, "{e:.3e} ({a:.3e}/{b:.3e})"),
+            Self::Indeterminate((a, b)) => write!(f, "∞ ({a:.3e}/{b:.3e})"),
+        }
+    }
+}
+
+#[non_exhaustive]
+pub struct EigGeneralizedWork<T: Scalar> {
+    /// Problem size
+    pub n: i32,
+    /// Compute right eigenvectors or not
+    pub jobvr: JobEv,
+    /// Compute left eigenvectors or not
+    pub jobvl: JobEv,
+
+    /// Eigenvalues: alpha (numerators)
+    pub alpha: Vec<MaybeUninit<T::Complex>>,
+    /// Eigenvalues: beta (denominators)
+    pub beta: Vec<MaybeUninit<T::Complex>>,
+    /// Real part of alpha (eigenvalue numerators) used in real routines
+    pub alpha_re: Option<Vec<MaybeUninit<T::Real>>>,
+    /// Imaginary part of alpha (eigenvalue numerators) used in real routines
+    pub alpha_im: Option<Vec<MaybeUninit<T::Real>>>,
+    /// Real part of beta (eigenvalue denominators) used in real routines
+    pub beta_re: Option<Vec<MaybeUninit<T::Real>>>,
+    /// Imaginary part of beta (eigenvalue denominators) used in real routines
+    pub beta_im: Option<Vec<MaybeUninit<T::Real>>>,
+
+    /// Left eigenvectors
+    pub vc_l: Option<Vec<MaybeUninit<T::Complex>>>,
+    /// Left eigenvectors used in real routines
+    pub vr_l: Option<Vec<MaybeUninit<T::Real>>>,
+    /// Right eigenvectors
+    pub vc_r: Option<Vec<MaybeUninit<T::Complex>>>,
+    /// Right eigenvectors used in real routines
+    pub vr_r: Option<Vec<MaybeUninit<T::Real>>>,
+
+    /// Working memory
+    pub work: Vec<MaybeUninit<T>>,
+    /// Working memory with `T::Real`
+    pub rwork: Option<Vec<MaybeUninit<T::Real>>>,
+}
+
+impl<T> EigGeneralizedWork<T>
+where
+    T: Scalar,
+    EigGeneralizedWork<T>: EigGeneralizedWorkImpl<Elem = T>,
+{
+    /// Create new working memory for eigenvalues compution.
+    pub fn new(calc_v: bool, l: MatrixLayout) -> Result<Self> {
+        EigGeneralizedWorkImpl::new(calc_v, l)
+    }
+
+    /// Compute eigenvalues and vectors on this working memory.
+    pub fn calc(&mut self, a: &mut [T], b: &mut [T]) -> Result<EigGeneralizedRef<T>> {
+        EigGeneralizedWorkImpl::calc(self, a, b)
+    }
+
+    /// Compute eigenvalues and vectors by consuming this working memory.
+    pub fn eval(self, a: &mut [T], b: &mut [T]) -> Result<EigGeneralizedOwned<T>> {
+        EigGeneralizedWorkImpl::eval(self, a, b)
+    }
+}
+
+/// Owned result of eigenvalue problem by [EigGeneralizedWork::eval]
+#[derive(Debug, Clone, PartialEq)]
+pub struct EigGeneralizedOwned<T: Scalar> {
+    /// Eigenvalues
+    pub alpha: Vec<T::Complex>,
+
+    pub beta: Vec<T::Complex>,
+
+    /// Right eigenvectors
+    pub vr: Option<Vec<T::Complex>>,
+
+    /// Left eigenvectors
+    pub vl: Option<Vec<T::Complex>>,
+}
+
+/// Reference result of eigenvalue problem by [EigGeneralizedWork::calc]
+#[derive(Debug, Clone, PartialEq)]
+pub struct EigGeneralizedRef<'work, T: Scalar> {
+    /// Eigenvalues
+    pub alpha: &'work [T::Complex],
+
+    pub beta: &'work [T::Complex],
+
+    /// Right eigenvectors
+    pub vr: Option<&'work [T::Complex]>,
+
+    /// Left eigenvectors
+    pub vl: Option<&'work [T::Complex]>,
+}
+
+/// Helper trait for implementing [EigGeneralizedWork] methods
+pub trait EigGeneralizedWorkImpl: Sized {
+    type Elem: Scalar;
+    fn new(calc_v: bool, l: MatrixLayout) -> Result<Self>;
+    fn calc<'work>(
+        &'work mut self,
+        a: &mut [Self::Elem],
+        b: &mut [Self::Elem],
+    ) -> Result<EigGeneralizedRef<'work, Self::Elem>>;
+    fn eval(
+        self,
+        a: &mut [Self::Elem],
+        b: &mut [Self::Elem],
+    ) -> Result<EigGeneralizedOwned<Self::Elem>>;
+}
+
+macro_rules! impl_eig_generalized_work_c {
+    ($f:ty, $c:ty, $ev:path) => {
+        impl EigGeneralizedWorkImpl for EigGeneralizedWork<$c> {
+            type Elem = $c;
+
+            fn new(calc_v: bool, l: MatrixLayout) -> Result<Self> {
+                let (n, _) = l.size();
+                let (jobvl, jobvr) = if calc_v {
+                    match l {
+                        MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
+                        MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
+                    }
+                } else {
+                    (JobEv::None, JobEv::None)
+                };
+                let mut rwork = vec_uninit(8 * n as usize);
+
+                let mut alpha = vec_uninit(n as usize);
+                let mut beta = vec_uninit(n as usize);
+
+                let mut vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
+                let mut vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
+
+                // calc work size
+                let mut info = 0;
+                let mut work_size = [<$c>::zero()];
+                unsafe {
+                    $ev(
+                        jobvl.as_ptr(),
+                        jobvr.as_ptr(),
+                        &n,
+                        std::ptr::null_mut(),
+                        &n,
+                        std::ptr::null_mut(),
+                        &n,
+                        AsPtr::as_mut_ptr(&mut alpha),
+                        AsPtr::as_mut_ptr(&mut beta),
+                        AsPtr::as_mut_ptr(vc_l.as_deref_mut().unwrap_or(&mut [])),
+                        &n,
+                        AsPtr::as_mut_ptr(vc_r.as_deref_mut().unwrap_or(&mut [])),
+                        &n,
+                        AsPtr::as_mut_ptr(&mut work_size),
+                        &(-1),
+                        AsPtr::as_mut_ptr(&mut rwork),
+                        &mut info,
+                    )
+                };
+                info.as_lapack_result()?;
+
+                let lwork = work_size[0].to_usize().unwrap();
+                let work: Vec<MaybeUninit<$c>> = vec_uninit(lwork);
+                Ok(Self {
+                    n,
+                    jobvl,
+                    jobvr,
+                    alpha,
+                    beta,
+                    alpha_re: None,
+                    alpha_im: None,
+                    beta_re: None,
+                    beta_im: None,
+                    rwork: Some(rwork),
+                    vc_l,
+                    vc_r,
+                    vr_l: None,
+                    vr_r: None,
+                    work,
+                })
+            }
+
+            fn calc<'work>(
+                &'work mut self,
+                a: &mut [Self::Elem],
+                b: &mut [Self::Elem],
+            ) -> Result<EigGeneralizedRef<'work, Self::Elem>> {
+                let lwork = self.work.len().to_i32().unwrap();
+                let mut info = 0;
+                unsafe {
+                    $ev(
+                        self.jobvl.as_ptr(),
+                        self.jobvr.as_ptr(),
+                        &self.n,
+                        AsPtr::as_mut_ptr(a),
+                        &self.n,
+                        AsPtr::as_mut_ptr(b),
+                        &self.n,
+                        AsPtr::as_mut_ptr(&mut self.alpha),
+                        AsPtr::as_mut_ptr(&mut self.beta),
+                        AsPtr::as_mut_ptr(self.vc_l.as_deref_mut().unwrap_or(&mut [])),
+                        &self.n,
+                        AsPtr::as_mut_ptr(self.vc_r.as_deref_mut().unwrap_or(&mut [])),
+                        &self.n,
+                        AsPtr::as_mut_ptr(&mut self.work),
+                        &lwork,
+                        AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
+                        &mut info,
+                    )
+                };
+                info.as_lapack_result()?;
+                // Hermite conjugate
+                if let Some(vl) = self.vc_l.as_mut() {
+                    for value in vl {
+                        let value = unsafe { value.assume_init_mut() };
+                        value.im = -value.im;
+                    }
+                }
+                Ok(EigGeneralizedRef {
+                    alpha: unsafe { self.alpha.slice_assume_init_ref() },
+                    beta: unsafe { self.beta.slice_assume_init_ref() },
+                    vl: self
+                        .vc_l
+                        .as_ref()
+                        .map(|v| unsafe { v.slice_assume_init_ref() }),
+                    vr: self
+                        .vc_r
+                        .as_ref()
+                        .map(|v| unsafe { v.slice_assume_init_ref() }),
+                })
+            }
+
+            fn eval(
+                mut self,
+                a: &mut [Self::Elem],
+                b: &mut [Self::Elem],
+            ) -> Result<EigGeneralizedOwned<Self::Elem>> {
+                let _eig_generalized_ref = self.calc(a, b)?;
+                Ok(EigGeneralizedOwned {
+                    alpha: unsafe { self.alpha.assume_init() },
+                    beta: unsafe { self.beta.assume_init() },
+                    vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
+                    vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
+                })
+            }
+        }
+
+        impl EigGeneralizedOwned<$c> {
+            pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<GeneralizedEigenvalue<$c>> {
+                self.alpha
+                    .iter()
+                    .zip(self.beta.iter())
+                    .map(|(alpha, beta)| {
+                        if let Some(thresh) = thresh_opt {
+                            if beta.abs() < thresh {
+                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+                            } else {
+                                GeneralizedEigenvalue::Finite(
+                                    alpha / beta,
+                                    (alpha.clone(), beta.clone()),
+                                )
+                            }
+                        } else {
+                            if beta.is_zero() {
+                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+                            } else {
+                                GeneralizedEigenvalue::Finite(
+                                    alpha / beta,
+                                    (alpha.clone(), beta.clone()),
+                                )
+                            }
+                        }
+                    })
+                    .collect::<Vec<_>>()
+            }
+        }
+    };
+}
+
+impl_eig_generalized_work_c!(f32, c32, lapack_sys::cggev_);
+impl_eig_generalized_work_c!(f64, c64, lapack_sys::zggev_);
+
+macro_rules! impl_eig_generalized_work_r {
+    ($f:ty, $c:ty, $ev:path) => {
+        impl EigGeneralizedWorkImpl for EigGeneralizedWork<$f> {
+            type Elem = $f;
+
+            fn new(calc_v: bool, l: MatrixLayout) -> Result<Self> {
+                let (n, _) = l.size();
+                let (jobvl, jobvr) = if calc_v {
+                    match l {
+                        MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
+                        MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
+                    }
+                } else {
+                    (JobEv::None, JobEv::None)
+                };
+                let mut alpha_re = vec_uninit(n as usize);
+                let mut alpha_im = vec_uninit(n as usize);
+                let mut beta_re = vec_uninit(n as usize);
+                let mut vr_l = jobvl.then(|| vec_uninit((n * n) as usize));
+                let mut vr_r = jobvr.then(|| vec_uninit((n * n) as usize));
+                let vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
+                let vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
+
+                // calc work size
+                let mut info = 0;
+                let mut work_size: [$f; 1] = [0.0];
+                unsafe {
+                    $ev(
+                        jobvl.as_ptr(),
+                        jobvr.as_ptr(),
+                        &n,
+                        std::ptr::null_mut(),
+                        &n,
+                        std::ptr::null_mut(),
+                        &n,
+                        AsPtr::as_mut_ptr(&mut alpha_re),
+                        AsPtr::as_mut_ptr(&mut alpha_im),
+                        AsPtr::as_mut_ptr(&mut beta_re),
+                        AsPtr::as_mut_ptr(vr_l.as_deref_mut().unwrap_or(&mut [])),
+                        &n,
+                        AsPtr::as_mut_ptr(vr_r.as_deref_mut().unwrap_or(&mut [])),
+                        &n,
+                        AsPtr::as_mut_ptr(&mut work_size),
+                        &(-1),
+                        &mut info,
+                    )
+                };
+                info.as_lapack_result()?;
+
+                // actual ev
+                let lwork = work_size[0].to_usize().unwrap();
+                let work = vec_uninit(lwork);
+
+                Ok(Self {
+                    n,
+                    jobvr,
+                    jobvl,
+                    alpha: vec_uninit(n as usize),
+                    beta: vec_uninit(n as usize),
+                    alpha_re: Some(alpha_re),
+                    alpha_im: Some(alpha_im),
+                    beta_re: Some(beta_re),
+                    beta_im: None,
+                    rwork: None,
+                    vr_l,
+                    vr_r,
+                    vc_l,
+                    vc_r,
+                    work,
+                })
+            }
+
+            fn calc<'work>(
+                &'work mut self,
+                a: &mut [Self::Elem],
+                b: &mut [Self::Elem],
+            ) -> Result<EigGeneralizedRef<'work, Self::Elem>> {
+                let lwork = self.work.len().to_i32().unwrap();
+                let mut info = 0;
+                unsafe {
+                    $ev(
+                        self.jobvl.as_ptr(),
+                        self.jobvr.as_ptr(),
+                        &self.n,
+                        AsPtr::as_mut_ptr(a),
+                        &self.n,
+                        AsPtr::as_mut_ptr(b),
+                        &self.n,
+                        AsPtr::as_mut_ptr(self.alpha_re.as_mut().unwrap()),
+                        AsPtr::as_mut_ptr(self.alpha_im.as_mut().unwrap()),
+                        AsPtr::as_mut_ptr(self.beta_re.as_mut().unwrap()),
+                        AsPtr::as_mut_ptr(self.vr_l.as_deref_mut().unwrap_or(&mut [])),
+                        &self.n,
+                        AsPtr::as_mut_ptr(self.vr_r.as_deref_mut().unwrap_or(&mut [])),
+                        &self.n,
+                        AsPtr::as_mut_ptr(&mut self.work),
+                        &lwork,
+                        &mut info,
+                    )
+                };
+                info.as_lapack_result()?;
+
+                let alpha_re = self
+                    .alpha_re
+                    .as_ref()
+                    .map(|e| unsafe { e.slice_assume_init_ref() })
+                    .unwrap();
+                let alpha_im = self
+                    .alpha_im
+                    .as_ref()
+                    .map(|e| unsafe { e.slice_assume_init_ref() })
+                    .unwrap();
+                let beta_re = self
+                    .beta_re
+                    .as_ref()
+                    .map(|e| unsafe { e.slice_assume_init_ref() })
+                    .unwrap();
+                reconstruct_eigs_optional_im(alpha_re, Some(alpha_im), &mut self.alpha);
+                reconstruct_eigs_optional_im(beta_re, None, &mut self.beta);
+
+                if let Some(v) = self.vr_l.as_ref() {
+                    let v = unsafe { v.slice_assume_init_ref() };
+                    reconstruct_eigenvectors(true, alpha_im, v, self.vc_l.as_mut().unwrap());
+                }
+                if let Some(v) = self.vr_r.as_ref() {
+                    let v = unsafe { v.slice_assume_init_ref() };
+                    reconstruct_eigenvectors(false, alpha_im, v, self.vc_r.as_mut().unwrap());
+                }
+
+                Ok(EigGeneralizedRef {
+                    alpha: unsafe { self.alpha.slice_assume_init_ref() },
+                    beta: unsafe { self.beta.slice_assume_init_ref() },
+                    vl: self
+                        .vc_l
+                        .as_ref()
+                        .map(|v| unsafe { v.slice_assume_init_ref() }),
+                    vr: self
+                        .vc_r
+                        .as_ref()
+                        .map(|v| unsafe { v.slice_assume_init_ref() }),
+                })
+            }
+
+            fn eval(
+                mut self,
+                a: &mut [Self::Elem],
+                b: &mut [Self::Elem],
+            ) -> Result<EigGeneralizedOwned<Self::Elem>> {
+                let _eig_generalized_ref = self.calc(a, b)?;
+                Ok(EigGeneralizedOwned {
+                    alpha: unsafe { self.alpha.assume_init() },
+                    beta: unsafe { self.beta.assume_init() },
+                    vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
+                    vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
+                })
+            }
+        }
+
+        impl EigGeneralizedOwned<$f> {
+            pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<GeneralizedEigenvalue<$c>> {
+                self.alpha
+                    .iter()
+                    .zip(self.beta.iter())
+                    .map(|(alpha, beta)| {
+                        if let Some(thresh) = thresh_opt {
+                            if beta.abs() < thresh {
+                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+                            } else {
+                                GeneralizedEigenvalue::Finite(
+                                    alpha / beta,
+                                    (alpha.clone(), beta.clone()),
+                                )
+                            }
+                        } else {
+                            if beta.is_zero() {
+                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
+                            } else {
+                                GeneralizedEigenvalue::Finite(
+                                    alpha / beta,
+                                    (alpha.clone(), beta.clone()),
+                                )
+                            }
+                        }
+                    })
+                    .collect::<Vec<_>>()
+            }
+        }
+    };
+}
+impl_eig_generalized_work_r!(f32, c32, lapack_sys::sggev_);
+impl_eig_generalized_work_r!(f64, c64, lapack_sys::dggev_);
+
+/// Create complex eigenvalues from real and optional imaginary parts.
+fn reconstruct_eigs_optional_im<T: Scalar>(
+    re: &[T],
+    im_opt: Option<&[T]>,
+    eigs: &mut [MaybeUninit<T::Complex>],
+) {
+    let n = eigs.len();
+    assert_eq!(re.len(), n);
+
+    if let Some(im) = im_opt {
+        assert_eq!(im.len(), n);
+        for i in 0..n {
+            eigs[i].write(T::complex(re[i], im[i]));
+        }
+    } else {
+        for i in 0..n {
+            eigs[i].write(T::complex(re[i], T::zero()));
+        }
+    }
+}
diff --git a/lax/src/flags.rs b/lax/src/flags.rs
index dd9dde3d..f9dea20d 100644
--- a/lax/src/flags.rs
+++ b/lax/src/flags.rs
@@ -1,4 +1,5 @@
 //! Charactor flags, e.g. `'T'`, used in LAPACK API
+use core::ffi::c_char;
 
 /// Upper/Lower specification for seveal usages
 #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
@@ -17,8 +18,8 @@ impl UPLO {
     }
 
     /// To use Fortran LAPACK API in lapack-sys crate
-    pub fn as_ptr(&self) -> *const i8 {
-        self as *const UPLO as *const i8
+    pub fn as_ptr(&self) -> *const c_char {
+        self as *const UPLO as *const c_char
     }
 }
 
@@ -32,8 +33,8 @@ pub enum Transpose {
 
 impl Transpose {
     /// To use Fortran LAPACK API in lapack-sys crate
-    pub fn as_ptr(&self) -> *const i8 {
-        self as *const Transpose as *const i8
+    pub fn as_ptr(&self) -> *const c_char {
+        self as *const Transpose as *const c_char
     }
 }
 
@@ -55,8 +56,8 @@ impl NormType {
     }
 
     /// To use Fortran LAPACK API in lapack-sys crate
-    pub fn as_ptr(&self) -> *const i8 {
-        self as *const NormType as *const i8
+    pub fn as_ptr(&self) -> *const c_char {
+        self as *const NormType as *const c_char
     }
 }
 
@@ -87,8 +88,8 @@ impl JobEv {
     }
 
     /// To use Fortran LAPACK API in lapack-sys crate
-    pub fn as_ptr(&self) -> *const i8 {
-        self as *const JobEv as *const i8
+    pub fn as_ptr(&self) -> *const c_char {
+        self as *const JobEv as *const c_char
     }
 }
 
@@ -117,8 +118,8 @@ impl JobSvd {
         }
     }
 
-    pub fn as_ptr(&self) -> *const i8 {
-        self as *const JobSvd as *const i8
+    pub fn as_ptr(&self) -> *const c_char {
+        self as *const JobSvd as *const c_char
     }
 }
 
@@ -133,7 +134,7 @@ pub enum Diag {
 }
 
 impl Diag {
-    pub fn as_ptr(&self) -> *const i8 {
-        self as *const Diag as *const i8
+    pub fn as_ptr(&self) -> *const c_char {
+        self as *const Diag as *const c_char
     }
 }
diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index af48f257..680ff0db 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -62,6 +62,7 @@
 //! there are several types of eigenvalue problem API
 //!
 //! - [eig] module for eigenvalue problem for general matrix.
+//! - [eig_generalized] module for generalized eigenvalue problem for general matrix.
 //! - [eigh] module for eigenvalue problem for symmetric/Hermitian matrix.
 //! - [eigh_generalized] module for generalized eigenvalue problem for symmetric/Hermitian matrix.
 //!
@@ -87,6 +88,7 @@ extern crate netlib_src as _src;
 pub mod alloc;
 pub mod cholesky;
 pub mod eig;
+pub mod eig_generalized;
 pub mod eigh;
 pub mod eigh_generalized;
 pub mod error;
@@ -103,6 +105,8 @@ pub mod svddc;
 pub mod triangular;
 pub mod tridiagonal;
 
+pub use crate::eig_generalized::GeneralizedEigenvalue;
+
 pub use self::flags::*;
 pub use self::least_squares::LeastSquaresOwned;
 pub use self::svd::{SvdOwned, SvdRef};
@@ -124,6 +128,18 @@ pub trait Lapack: Scalar {
         a: &mut [Self],
     ) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)>;
 
+    /// Compute right eigenvalue and eigenvectors for a general matrix
+    fn eig_generalized(
+        calc_v: bool,
+        l: MatrixLayout,
+        a: &mut [Self],
+        b: &mut [Self],
+        thresh_opt: Option<Self::Real>,
+    ) -> Result<(
+        Vec<GeneralizedEigenvalue<Self::Complex>>,
+        Vec<Self::Complex>,
+    )>;
+
     /// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix
     fn eigh(
         calc_eigenvec: bool,
@@ -331,6 +347,25 @@ macro_rules! impl_lapack {
                 Ok((eigs, vr.or(vl).unwrap_or_default()))
             }
 
+            fn eig_generalized(
+                calc_v: bool,
+                l: MatrixLayout,
+                a: &mut [Self],
+                b: &mut [Self],
+                thresh_opt: Option<Self::Real>,
+            ) -> Result<(
+                Vec<GeneralizedEigenvalue<Self::Complex>>,
+                Vec<Self::Complex>,
+            )> {
+                use eig_generalized::*;
+                let work = EigGeneralizedWork::<$s>::new(calc_v, l)?;
+                let eig_generalized_owned = work.eval(a, b)?;
+                let eigs = eig_generalized_owned.calc_eigs(thresh_opt);
+                let vr = eig_generalized_owned.vr;
+                let vl = eig_generalized_owned.vl;
+                Ok((eigs, vr.or(vl).unwrap_or_default()))
+            }
+
             fn eigh(
                 calc_eigenvec: bool,
                 layout: MatrixLayout,
diff --git a/ndarray-linalg/Cargo.toml b/ndarray-linalg/Cargo.toml
index e06fa675..8c5c83a8 100644
--- a/ndarray-linalg/Cargo.toml
+++ b/ndarray-linalg/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name    = "ndarray-linalg"
-version = "0.16.0"
+version = "0.17.0"
 authors = ["Toshiki Teramura <toshiki.teramura@gmail.com>"]
 edition = "2018"
 
@@ -13,7 +13,8 @@ readme        = "../README.md"
 categories    = ["algorithms", "science"]
 
 [features]
-default   = []
+default   = ["blas"]
+blas      = ["ndarray/blas"]
 
 netlib    = ["lax/netlib"]
 openblas  = ["lax/openblas"]
@@ -34,23 +35,23 @@ katexit = "0.1.2"
 num-complex = "0.4.0"
 num-traits  = "0.2.14"
 rand = "0.8.3"
-thiserror = "1.0.24"
+thiserror = "2.0.0"
 
 [dependencies.ndarray]
-version = "0.15.2"
-features = ["blas", "approx", "std"]
+version = "0.16.0"
+features = ["approx", "std"]
 default-features = false
 
 [dependencies.lax]
-version = "0.16.0-rc.0"
+version = "0.17.0"
 path = "../lax"
 default-features = false
 
 [dev-dependencies]
 paste = "1.0.5"
-criterion = "0.3.4"
+criterion = "0.5.1"
 # Keep the same version as ndarray's dependency!
-approx = { version = "0.4.0", features = ["num-complex"] }
+approx = { version = "0.5", features = ["num-complex"] }
 rand_pcg = "0.3.1"
 
 [[bench]]
diff --git a/ndarray-linalg/benches/eig_generalized.rs b/ndarray-linalg/benches/eig_generalized.rs
new file mode 100644
index 00000000..d1f5621b
--- /dev/null
+++ b/ndarray-linalg/benches/eig_generalized.rs
@@ -0,0 +1,40 @@
+use criterion::*;
+use ndarray::*;
+use ndarray_linalg::*;
+
+fn eig_generalized_small(c: &mut Criterion) {
+    let mut group = c.benchmark_group("eig");
+    for &n in &[4, 8, 16, 32, 64, 128] {
+        group.bench_with_input(BenchmarkId::new("vecs/C/r", n), &n, |c, n| {
+            let a: Array2<f64> = random((*n, *n));
+            let b: Array2<f64> = random((*n, *n));
+            c.iter(|| {
+                let (_e, _vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+            })
+        });
+        group.bench_with_input(BenchmarkId::new("vecs/F/r", n), &n, |c, n| {
+            let a: Array2<f64> = random((*n, *n).f());
+            let b: Array2<f64> = random((*n, *n).f());
+            c.iter(|| {
+                let (_e, _vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+            })
+        });
+        group.bench_with_input(BenchmarkId::new("vecs/C/c", n), &n, |c, n| {
+            let a: Array2<c64> = random((*n, *n));
+            let b: Array2<c64> = random((*n, *n));
+            c.iter(|| {
+                let (_e, _vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+            })
+        });
+        group.bench_with_input(BenchmarkId::new("vecs/F/c", n), &n, |c, n| {
+            let a: Array2<c64> = random((*n, *n).f());
+            let b: Array2<c64> = random((*n, *n).f());
+            c.iter(|| {
+                let (_e, _vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+            })
+        });
+    }
+}
+
+criterion_group!(eig, eig_generalized_small);
+criterion_main!(eig);
diff --git a/ndarray-linalg/src/assert.rs b/ndarray-linalg/src/assert.rs
index 670aabc8..74c1618b 100644
--- a/ndarray-linalg/src/assert.rs
+++ b/ndarray-linalg/src/assert.rs
@@ -86,7 +86,7 @@ where
 }
 
 macro_rules! generate_assert {
-    ($assert:ident, $close:path) => {
+    ($assert:ident, $close:tt) => {
         #[macro_export]
         macro_rules! $assert {
             ($test: expr,$truth: expr,$tol: expr) => {
diff --git a/ndarray-linalg/src/convert.rs b/ndarray-linalg/src/convert.rs
index 43002966..c808211e 100644
--- a/ndarray-linalg/src/convert.rs
+++ b/ndarray-linalg/src/convert.rs
@@ -12,7 +12,7 @@ where
     S: Data,
 {
     let n = a.len();
-    a.into_shape((n, 1)).unwrap()
+    a.into_shape_with_order((n, 1)).unwrap()
 }
 
 pub fn into_row<S>(a: ArrayBase<S, Ix1>) -> ArrayBase<S, Ix2>
@@ -20,7 +20,7 @@ where
     S: Data,
 {
     let n = a.len();
-    a.into_shape((1, n)).unwrap()
+    a.into_shape_with_order((1, n)).unwrap()
 }
 
 pub fn flatten<S>(a: ArrayBase<S, Ix2>) -> ArrayBase<S, Ix1>
@@ -28,7 +28,7 @@ where
     S: Data,
 {
     let n = a.len();
-    a.into_shape(n).unwrap()
+    a.into_shape_with_order(n).unwrap()
 }
 
 pub fn into_matrix<A, S>(l: MatrixLayout, a: Vec<A>) -> Result<ArrayBase<S, Ix2>>
@@ -99,9 +99,9 @@ where
     // https://github.com/bluss/rust-ndarray/issues/325
     let strides: Vec<isize> = a.strides().to_vec();
     let new = if a.is_standard_layout() {
-        ArrayBase::from_shape_vec(a.dim(), a.into_raw_vec()).unwrap()
+        ArrayBase::from_shape_vec(a.dim(), a.into_raw_vec_and_offset().0).unwrap()
     } else {
-        ArrayBase::from_shape_vec(a.dim().f(), a.into_raw_vec()).unwrap()
+        ArrayBase::from_shape_vec(a.dim().f(), a.into_raw_vec_and_offset().0).unwrap()
     };
     assert_eq!(
         new.strides(),
diff --git a/ndarray-linalg/src/eig.rs b/ndarray-linalg/src/eig.rs
index bde5bd7d..03e3ee03 100644
--- a/ndarray-linalg/src/eig.rs
+++ b/ndarray-linalg/src/eig.rs
@@ -3,6 +3,7 @@
 use crate::error::*;
 use crate::layout::*;
 use crate::types::*;
+pub use lax::GeneralizedEigenvalue;
 use ndarray::*;
 
 #[cfg_attr(doc, katexit::katexit)]
@@ -77,3 +78,86 @@ where
         Ok(ArrayBase::from(s))
     }
 }
+
+#[cfg_attr(doc, katexit::katexit)]
+/// Eigenvalue decomposition of general matrix reference
+pub trait EigGeneralized {
+    /// EigVec is the right eivenvector
+    type EigVal;
+    type EigVec;
+    type Real;
+    /// Calculate eigenvalues with the right eigenvector
+    ///
+    /// $$ A u_i = \lambda_i B u_i $$
+    ///
+    /// ```
+    /// use ndarray::*;
+    /// use ndarray_linalg::*;
+    ///
+    /// let a: Array2<f64> = array![
+    ///     [-1.01,  0.86, -4.60,  3.31, -4.81],
+    ///     [ 3.98,  0.53, -7.04,  5.29,  3.55],
+    ///     [ 3.30,  8.26, -3.89,  8.20, -1.51],
+    ///     [ 4.43,  4.96, -7.66, -7.33,  6.18],
+    ///     [ 7.31, -6.43, -6.16,  2.47,  5.58],
+    /// ];
+    /// let b: Array2<f64> = array![
+    ///     [ 1.23, -4.56,  7.89,  0.12, -3.45],
+    ///     [ 6.78, -9.01,  2.34, -5.67,  8.90],
+    ///     [-1.11,  3.33, -6.66,  9.99, -2.22],
+    ///     [ 4.44, -7.77,  0.00,  1.11,  5.55],
+    ///     [-8.88,  6.66, -3.33,  2.22, -9.99],
+    /// ];
+    /// let (geneigs, vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+    ///
+    /// let a = a.map(|v| v.as_c());
+    /// let b = b.map(|v| v.as_c());
+    /// for (ge, vec) in geneigs.iter().zip(vecs.axis_iter(Axis(1))) {
+    ///     if let GeneralizedEigenvalue::Finite(e, _) = ge {
+    ///         let ebv = b.dot(&vec).map(|v| v * e);
+    ///         let av = a.dot(&vec);
+    ///         assert_close_l2!(&av, &ebv, 1e-5);
+    ///     }
+    /// }
+    /// ```
+    ///
+    /// # Arguments
+    ///
+    /// * `thresh_opt` - An optional threshold for determining approximate zero |β| values when
+    /// computing the eigenvalues as α/β. If `None`, no approximate comparisons to zero will be
+    /// made.
+    fn eig_generalized(
+        &self,
+        thresh_opt: Option<Self::Real>,
+    ) -> Result<(Self::EigVal, Self::EigVec)>;
+}
+
+impl<A, S> EigGeneralized for (ArrayBase<S, Ix2>, ArrayBase<S, Ix2>)
+where
+    A: Scalar + Lapack,
+    S: Data<Elem = A>,
+{
+    type EigVal = Array1<GeneralizedEigenvalue<A::Complex>>;
+    type EigVec = Array2<A::Complex>;
+    type Real = A::Real;
+
+    fn eig_generalized(
+        &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 layout = a.square_layout()?;
+        let (s, t) = A::eig_generalized(
+            true,
+            layout,
+            a.as_allocated_mut()?,
+            b.as_allocated_mut()?,
+            thresh_opt,
+        )?;
+        let n = layout.len() as usize;
+        Ok((
+            ArrayBase::from(s),
+            Array2::from_shape_vec((n, n).f(), t).unwrap(),
+        ))
+    }
+}
diff --git a/ndarray-linalg/tests/det.rs b/ndarray-linalg/tests/det.rs
index d14fc8d0..40dafd57 100644
--- a/ndarray-linalg/tests/det.rs
+++ b/ndarray-linalg/tests/det.rs
@@ -94,7 +94,7 @@ fn det_zero_nonsquare() {
             assert!(a.sln_det_into().is_err());
         };
     }
-    for &shape in &[(1, 2).into_shape(), (1, 2).f()] {
+    for &shape in &[(1, 2).into_shape_with_order(), (1, 2).f()] {
         det_zero_nonsquare!(f64, shape);
         det_zero_nonsquare!(f32, shape);
         det_zero_nonsquare!(c64, shape);
@@ -186,7 +186,7 @@ fn det_nonsquare() {
         };
     }
     for &dims in &[(1, 0), (1, 2), (2, 1), (2, 3)] {
-        for &shape in &[dims.into_shape(), dims.f()] {
+        for &shape in &[dims.into_shape_with_order(), dims.f()] {
             det_nonsquare!(f64, shape);
             det_nonsquare!(f32, shape);
             det_nonsquare!(c64, shape);
diff --git a/ndarray-linalg/tests/deth.rs b/ndarray-linalg/tests/deth.rs
index 4785aadc..9079c342 100644
--- a/ndarray-linalg/tests/deth.rs
+++ b/ndarray-linalg/tests/deth.rs
@@ -60,7 +60,7 @@ fn deth_zero_nonsquare() {
             assert!(a.sln_deth_into().is_err());
         };
     }
-    for &shape in &[(1, 2).into_shape(), (1, 2).f()] {
+    for &shape in &[(1, 2).into_shape_with_order(), (1, 2).f()] {
         deth_zero_nonsquare!(f64, shape);
         deth_zero_nonsquare!(f32, shape);
         deth_zero_nonsquare!(c64, shape);
@@ -138,7 +138,7 @@ fn deth_nonsquare() {
         };
     }
     for &dims in &[(1, 0), (1, 2), (2, 1), (2, 3)] {
-        for &shape in &[dims.into_shape(), dims.f()] {
+        for &shape in &[dims.into_shape_with_order(), dims.f()] {
             deth_nonsquare!(f64, shape);
             deth_nonsquare!(f32, shape);
             deth_nonsquare!(c64, shape);
diff --git a/ndarray-linalg/tests/eig_generalized.rs b/ndarray-linalg/tests/eig_generalized.rs
new file mode 100644
index 00000000..06df81ec
--- /dev/null
+++ b/ndarray-linalg/tests/eig_generalized.rs
@@ -0,0 +1,190 @@
+use ndarray::*;
+use ndarray_linalg::*;
+
+#[test]
+fn generalized_eigenvalue_fmt() {
+    let ge0 = GeneralizedEigenvalue::Finite(0.1, (1.0, 10.0));
+    assert_eq!(ge0.to_string(), "1.000e-1 (1.000e0/1.000e1)".to_string());
+
+    let ge1 = GeneralizedEigenvalue::Indeterminate((1.0, 0.0));
+    assert_eq!(ge1.to_string(), "∞ (1.000e0/0.000e0)".to_string());
+}
+
+#[test]
+fn real_a_real_b_3x3_full_rank() {
+    #[rustfmt::skip]
+    let a = array![
+        [ 2.0, 1.0, 8.0],
+        [-2.0, 0.0, 3.0],
+        [ 7.0, 6.0, 5.0],
+    ];
+    #[rustfmt::skip]
+    let b = array![
+        [ 1.0,  2.0, -7.0],
+        [-3.0,  1.0,  6.0],
+        [ 4.0, -5.0,  1.0],
+    ];
+    let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+
+    let a = a.map(|v| v.as_c());
+    let b = b.map(|v| v.as_c());
+    for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
+        if let GeneralizedEigenvalue::Finite(e, _) = ge {
+            let ebv = b.dot(&vec).map(|v| v * e);
+            let av = a.dot(&vec);
+            assert_close_l2!(&av, &ebv, 1e-7);
+        }
+    }
+
+    let mut eigvals = geneigvals
+        .iter()
+        .filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
+            GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
+            GeneralizedEigenvalue::Indeterminate(_) => None,
+        })
+        .collect::<Vec<_>>();
+    eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
+    let eigvals = Array1::from_vec(eigvals);
+    // Reference eigenvalues from Mathematica
+    assert_close_l2!(
+        &eigvals,
+        &array![-0.4415795111, 0.5619249537, 50.87965456].map(c64::from),
+        1e-7
+    );
+}
+
+#[test]
+fn real_a_real_b_3x3_nullity_1() {
+    #[rustfmt::skip]
+    let a = array![
+        [ 2.0, 1.0, 8.0],
+        [-2.0, 0.0, 3.0],
+        [ 7.0, 6.0, 5.0],
+    ];
+    #[rustfmt::skip]
+    let b = array![
+        [1.0,  2.0, 3.0],
+        [0.0,  1.0, 1.0],
+        [1.0, -1.0, 0.0],
+    ];
+    let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(Some(1e-4)).unwrap();
+
+    let a = a.map(|v| v.as_c());
+    let b = b.map(|v| v.as_c());
+    for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
+        if let GeneralizedEigenvalue::Finite(e, _) = ge {
+            let ebv = b.dot(&vec).map(|v| v * e);
+            let av = a.dot(&vec);
+            assert_close_l2!(&av, &ebv, 1e-7);
+        }
+    }
+
+    let mut eigvals = geneigvals
+        .iter()
+        .filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
+            GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
+            GeneralizedEigenvalue::Indeterminate(_) => None,
+        })
+        .collect::<Vec<_>>();
+    eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
+    let eigvals = Array1::from_vec(eigvals);
+    // Reference eigenvalues from Mathematica
+    assert_close_l2!(
+        &eigvals,
+        &array![-12.91130192, 3.911301921].map(c64::from),
+        1e-7
+    );
+}
+
+#[test]
+fn complex_a_complex_b_3x3_full_rank() {
+    #[rustfmt::skip]
+    let a = array![
+        [c64::new(1.0,  2.0), c64::new(-3.0,  0.5), c64::new( 0.0, -1.0)],
+        [c64::new(2.5, -4.0), c64::new( 1.0,  1.0), c64::new(-1.5,  2.5)],
+        [c64::new(0.0,  0.0), c64::new( 3.0, -2.0), c64::new( 4.0,  4.0)],
+    ];
+    #[rustfmt::skip]
+    let b = array![
+        [c64::new(-2.0,  1.0), c64::new( 3.5, -1.0), c64::new( 1.0,  1.0)],
+        [c64::new( 0.0, -3.0), c64::new( 2.0,  2.0), c64::new(-4.0,  0.0)],
+        [c64::new( 5.0,  5.0), c64::new(-1.5,  1.5), c64::new( 0.0, -2.0)],
+    ];
+    let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
+
+    let a = a.map(|v| v.as_c());
+    let b = b.map(|v| v.as_c());
+    for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
+        if let GeneralizedEigenvalue::Finite(e, _) = ge {
+            let ebv = b.dot(&vec).map(|v| v * e);
+            let av = a.dot(&vec);
+            assert_close_l2!(&av, &ebv, 1e-7);
+        }
+    }
+
+    let mut eigvals = geneigvals
+        .iter()
+        .filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
+            GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
+            GeneralizedEigenvalue::Indeterminate(_) => None,
+        })
+        .collect::<Vec<_>>();
+    eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
+    let eigvals = Array1::from_vec(eigvals);
+    // Reference eigenvalues from Mathematica
+    assert_close_l2!(
+        &eigvals,
+        &array![
+            c64::new(-0.701598, -1.71262),
+            c64::new(-0.67899, -0.0172468),
+            c64::new(0.59059, 0.276034)
+        ],
+        1e-5
+    );
+}
+
+#[test]
+fn complex_a_complex_b_3x3_nullity_1() {
+    #[rustfmt::skip]
+    let a = array![
+        [c64::new(1.0,  2.0), c64::new(-3.0,  0.5), c64::new( 0.0, -1.0)],
+        [c64::new(2.5, -4.0), c64::new( 1.0,  1.0), c64::new(-1.5,  2.5)],
+        [c64::new(0.0,  0.0), c64::new( 3.0, -2.0), c64::new( 4.0,  4.0)],
+    ];
+    #[rustfmt::skip]
+    let b = array![
+        [c64::new(-2.55604, -4.10176), c64::new(9.03944,  3.745000), c64::new(35.4641,  21.1704)],
+        [c64::new( 7.85029,  7.02144), c64::new(9.23225, -0.479451), c64::new(13.9507, -16.5402)],
+        [c64::new(-4.47803,  3.98981), c64::new(9.44434, -4.519970), c64::new(40.9006, -23.5060)],
+    ];
+    let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(Some(1e-4)).unwrap();
+
+    let a = a.map(|v| v.as_c());
+    let b = b.map(|v| v.as_c());
+    for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
+        if let GeneralizedEigenvalue::Finite(e, _) = ge {
+            let ebv = b.dot(&vec).map(|v| v * e);
+            let av = a.dot(&vec);
+            assert_close_l2!(&av, &ebv, 1e-7);
+        }
+    }
+
+    let mut eigvals = geneigvals
+        .iter()
+        .filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
+            GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
+            GeneralizedEigenvalue::Indeterminate(_) => None,
+        })
+        .collect::<Vec<_>>();
+    eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
+    let eigvals = Array1::from_vec(eigvals);
+    // Reference eigenvalues from Mathematica
+    assert_close_l2!(
+        &eigvals,
+        &array![
+            c64::new(-0.0620674, -0.270016),
+            c64::new(0.0218236, 0.0602709),
+        ],
+        1e-5
+    );
+}
diff --git a/ndarray-linalg/tests/inv.rs b/ndarray-linalg/tests/inv.rs
index 030773c1..93ff1aad 100644
--- a/ndarray-linalg/tests/inv.rs
+++ b/ndarray-linalg/tests/inv.rs
@@ -106,7 +106,9 @@ fn inv_into_random_complex() {
 #[should_panic]
 fn inv_error() {
     // do not have inverse
-    let a = Array::<f64, _>::zeros(9).into_shape((3, 3)).unwrap();
+    let a = Array::<f64, _>::zeros(9)
+        .into_shape_with_order((3, 3))
+        .unwrap();
     let a_inv = a.inv().unwrap();
     println!("{:?}", a_inv);
 }
diff --git a/ndarray-linalg/tests/opnorm.rs b/ndarray-linalg/tests/opnorm.rs
index abd41748..cd45d258 100644
--- a/ndarray-linalg/tests/opnorm.rs
+++ b/ndarray-linalg/tests/opnorm.rs
@@ -14,11 +14,13 @@ fn gen(i: usize, j: usize, rev: bool) -> Array2<f64> {
     let n = (i * j + 1) as f64;
     if rev {
         Array::range(1., n, 1.)
-            .into_shape((j, i))
+            .into_shape_with_order((j, i))
             .unwrap()
             .reversed_axes()
     } else {
-        Array::range(1., n, 1.).into_shape((i, j)).unwrap()
+        Array::range(1., n, 1.)
+            .into_shape_with_order((i, j))
+            .unwrap()
     }
 }
 
diff --git a/ndarray-linalg/tests/solve.rs b/ndarray-linalg/tests/solve.rs
index 57c29b67..074350ce 100644
--- a/ndarray-linalg/tests/solve.rs
+++ b/ndarray-linalg/tests/solve.rs
@@ -254,7 +254,7 @@ fn rcond_hilbert() {
     macro_rules! rcond_hilbert {
         ($elem:ty, $rows:expr, $atol:expr) => {
             let a = Array2::<$elem>::from_shape_fn(($rows, $rows), |(i, j)| {
-                1. / (i as $elem + j as $elem - 1.)
+                1. / (i as $elem + j as $elem + 1.)
             });
             assert_aclose!(a.rcond().unwrap(), 0., $atol);
             assert_aclose!(a.rcond_into().unwrap(), 0., $atol);