Merge Eigh_ into Lapack trait

Author: termoshtt

lax/src/eigh.rs

@@ -1,28 +1,17 @@
+//! Eigenvalue problem for symmetric/Hermite matricies
+//!
+//! LAPACK correspondance
+//! ----------------------
+//!
+//! | f32   | f64   | c32   | c64   |
+//! |:------|:------|:------|:------|
+//! | ssyev | dsyev | cheev | zheev |
+
 use super::*;
 use crate::{error::*, layout::MatrixLayout};
 use cauchy::*;
 use num_traits::{ToPrimitive, Zero};
 
-#[cfg_attr(doc, katexit::katexit)]
-/// Eigenvalue problem for symmetric/hermite matrix
-pub trait Eigh_: Scalar {
-    /// Compute right eigenvalue and eigenvectors $Ax = \lambda x$
-    ///
-    /// LAPACK correspondance
-    /// ----------------------
-    ///
-    /// | f32   | f64   | c32   | c64   |
-    /// |:------|:------|:------|:------|
-    /// | ssyev | dsyev | cheev | zheev |
-    ///
-    fn eigh(
-        calc_eigenvec: bool,
-        layout: MatrixLayout,
-        uplo: UPLO,
-        a: &mut [Self],
-    ) -> Result<Vec<Self::Real>>;
-}
-
 pub struct EighWork<T: Scalar> {
     pub n: i32,
     pub jobz: JobEv,
@@ -199,78 +188,3 @@ macro_rules! impl_eigh_work_r {
 }
 impl_eigh_work_r!(f64, lapack_sys::dsyev_);
 impl_eigh_work_r!(f32, lapack_sys::ssyev_);
-
-macro_rules! impl_eigh {
-    (@real, $scalar:ty, $ev:path) => {
-        impl_eigh!(@body, $scalar, $ev, );
-    };
-    (@complex, $scalar:ty, $ev:path) => {
-        impl_eigh!(@body, $scalar, $ev, rwork);
-    };
-    (@body, $scalar:ty, $ev:path, $($rwork_ident:ident),*) => {
-        impl Eigh_ for $scalar {
-            fn eigh(
-                calc_v: bool,
-                layout: MatrixLayout,
-                uplo: UPLO,
-                a: &mut [Self],
-            ) -> Result<Vec<Self::Real>> {
-                assert_eq!(layout.len(), layout.lda());
-                let n = layout.len();
-                let jobz = if calc_v { JobEv::All } else { JobEv::None };
-                let mut eigs: Vec<MaybeUninit<Self::Real>> = vec_uninit(n as usize);
-
-                $(
-                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = vec_uninit(3 * n as usize - 2 as usize);
-                )*
-
-                // calc work size
-                let mut info = 0;
-                let mut work_size = [Self::zero()];
-                unsafe {
-                    $ev(
-                        jobz.as_ptr() ,
-                        uplo.as_ptr(),
-                        &n,
-                        AsPtr::as_mut_ptr(a),
-                        &n,
-                        AsPtr::as_mut_ptr(&mut eigs),
-                        AsPtr::as_mut_ptr(&mut work_size),
-                        &(-1),
-                        $(AsPtr::as_mut_ptr(&mut $rwork_ident),)*
-                        &mut info,
-                    );
-                }
-                info.as_lapack_result()?;
-
-                // actual ev
-                let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = vec_uninit(lwork);
-                let lwork = lwork as i32;
-                unsafe {
-                    $ev(
-                        jobz.as_ptr(),
-                        uplo.as_ptr(),
-                        &n,
-                        AsPtr::as_mut_ptr(a),
-                        &n,
-                        AsPtr::as_mut_ptr(&mut eigs),
-                        AsPtr::as_mut_ptr(&mut work),
-                        &lwork,
-                        $(AsPtr::as_mut_ptr(&mut $rwork_ident),)*
-                        &mut info,
-                    );
-                }
-                info.as_lapack_result()?;
-
-                let eigs = unsafe { eigs.assume_init() };
-                Ok(eigs)
-            }
-        }
-    };
-} // impl_eigh!
-
-impl_eigh!(@real, f64, lapack_sys::dsyev_);
-impl_eigh!(@real, f32, lapack_sys::ssyev_);
-impl_eigh!(@complex, c64, lapack_sys::zheev_);
-impl_eigh!(@complex, c32, lapack_sys::cheev_);

lax/src/lib.rs

@@ -131,25 +131,31 @@ pub trait Lapack:
     + Solve_
     + Solveh_
     + Cholesky_
-    + Eigh_
     + EighGeneralized_
     + Triangular_
     + Tridiagonal_
     + Rcond_
     + LeastSquaresSvdDivideConquer_
 {
-    /// Compute right eigenvalue and eigenvectors
+    /// Compute right eigenvalue and eigenvectors for a general matrix
     fn eig(
         calc_v: bool,
         l: MatrixLayout,
         a: &mut [Self],
     ) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)>;
+
+    /// Compute right eigenvalue and eigenvectors for a symmetric or hermite matrix
+    fn eigh(
+        calc_eigenvec: bool,
+        layout: MatrixLayout,
+        uplo: UPLO,
+        a: &mut [Self],
+    ) -> Result<Vec<Self::Real>>;
 }
 
 macro_rules! impl_lapack {
     ($s:ty) => {
         impl Lapack for $s {
-            /// Compute right eigenvalue and eigenvectors
             fn eig(
                 calc_v: bool,
                 l: MatrixLayout,
@@ -160,6 +166,17 @@ macro_rules! impl_lapack {
                 let EigOwned { eigs, vr, vl } = work.eval(a)?;
                 Ok((eigs, vr.or(vl).unwrap_or_default()))
             }
+
+            fn eigh(
+                calc_eigenvec: bool,
+                layout: MatrixLayout,
+                uplo: UPLO,
+                a: &mut [Self],
+            ) -> Result<Vec<Self::Real>> {
+                use eigh::*;
+                let work = EighWork::<$s>::new(calc_eigenvec, layout)?;
+                work.eval(uplo, a)
+            }
         }
     };
 }