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( +pub(crate) fn reconstruct_eigenvectors( 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 { + /// Finite generalized eigenvalue: `Finite(α/β, (α, β))` + Finite(T, (T, T)), + + /// Indeterminate generalized eigenvalue: `Indeterminate((α, β))` + Indeterminate((T, T)), +} + +impl std::fmt::Display for GeneralizedEigenvalue { + 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 { + /// 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>, + /// Eigenvalues: beta (denominators) + pub beta: Vec>, + /// Real part of alpha (eigenvalue numerators) used in real routines + pub alpha_re: Option>>, + /// Imaginary part of alpha (eigenvalue numerators) used in real routines + pub alpha_im: Option>>, + /// Real part of beta (eigenvalue denominators) used in real routines + pub beta_re: Option>>, + /// Imaginary part of beta (eigenvalue denominators) used in real routines + pub beta_im: Option>>, + + /// Left eigenvectors + pub vc_l: Option>>, + /// Left eigenvectors used in real routines + pub vr_l: Option>>, + /// Right eigenvectors + pub vc_r: Option>>, + /// Right eigenvectors used in real routines + pub vr_r: Option>>, + + /// Working memory + pub work: Vec>, + /// Working memory with `T::Real` + pub rwork: Option>>, +} + +impl EigGeneralizedWork +where + T: Scalar, + EigGeneralizedWork: EigGeneralizedWorkImpl, +{ + /// Create new working memory for eigenvalues compution. + pub fn new(calc_v: bool, l: MatrixLayout) -> Result { + 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> { + 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> { + EigGeneralizedWorkImpl::eval(self, a, b) + } +} + +/// Owned result of eigenvalue problem by [EigGeneralizedWork::eval] +#[derive(Debug, Clone, PartialEq)] +pub struct EigGeneralizedOwned { + /// Eigenvalues + pub alpha: Vec, + + pub beta: Vec, + + /// Right eigenvectors + pub vr: Option>, + + /// Left eigenvectors + pub vl: Option>, +} + +/// 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; + fn calc<'work>( + &'work mut self, + a: &mut [Self::Elem], + b: &mut [Self::Elem], + ) -> Result>; + fn eval( + self, + a: &mut [Self::Elem], + b: &mut [Self::Elem], + ) -> Result>; +} + +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 { + 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> = 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> { + 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> { + 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> { + 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::>() + } + } + }; +} + +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 { + 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> { + 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> { + 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> { + 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::>() + } + } + }; +} +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( + re: &[T], + im_opt: Option<&[T]>, + eigs: &mut [MaybeUninit], +) { + 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/lib.rs b/lax/src/lib.rs index af48f257..fdc6d2f2 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,15 @@ pub trait Lapack: Scalar { a: &mut [Self], ) -> Result<(Vec, Vec)>; + /// 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, + ) -> Result<(Vec>, Vec)>; + /// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix fn eigh( calc_eigenvec: bool, @@ -331,6 +344,22 @@ 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, + ) -> Result<(Vec>, Vec)> { + 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/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 = random((*n, *n)); + let b: Array2 = 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 = random((*n, *n).f()); + let b: Array2 = 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 = random((*n, *n)); + let b: Array2 = 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 = random((*n, *n).f()); + let b: Array2 = 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/eig.rs b/ndarray-linalg/src/eig.rs index bde5bd7d..a664d1bb 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,81 @@ 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 = 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 = 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, + ) -> Result<(Self::EigVal, Self::EigVec)>; +} + +impl EigGeneralized for (ArrayBase, ArrayBase) +where + A: Scalar + Lapack, + S: Data, +{ + type EigVal = Array1>; + type EigVec = Array2; + type Real = A::Real; + + fn eig_generalized( + &self, + thresh_opt: Option, + ) -> 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/eig_generalized.rs b/ndarray-linalg/tests/eig_generalized.rs new file mode 100644 index 00000000..bb41abfe --- /dev/null +++ b/ndarray-linalg/tests/eig_generalized.rs @@ -0,0 +1,181 @@ +use ndarray::*; +use ndarray_linalg::*; + +#[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| match ge { + GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()), + GeneralizedEigenvalue::Indeterminate(_) => None, + }) + .collect::>(); + 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| match ge { + GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()), + GeneralizedEigenvalue::Indeterminate(_) => None, + }) + .collect::>(); + 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| match ge { + GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()), + GeneralizedEigenvalue::Indeterminate(_) => None, + }) + .collect::>(); + 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| match ge { + GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()), + GeneralizedEigenvalue::Indeterminate(_) => None, + }) + .collect::>(); + 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 + ); +}