From 433db1a3896a2835476476c1e99796501c92bf1c Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Sat, 24 Sep 2022 23:56:02 +0900
Subject: [PATCH] WIP: EighWork<c64>

---
 lax/src/eigh.rs | 68 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 68 insertions(+)

diff --git a/lax/src/eigh.rs b/lax/src/eigh.rs
index cd83f2b8..3c255d11 100644
--- a/lax/src/eigh.rs
+++ b/lax/src/eigh.rs
@@ -23,6 +23,74 @@ pub trait Eigh_: Scalar {
     ) -> Result<Vec<Self::Real>>;
 }
 
+pub struct EighWork<T: Scalar> {
+    pub jobz: JobEv,
+    pub eigs: Vec<MaybeUninit<T::Real>>,
+    pub work: Vec<MaybeUninit<T>>,
+    pub rwork: Option<Vec<MaybeUninit<T::Real>>>,
+}
+
+pub trait EighWorkImpl: Sized {
+    type Elem: Scalar;
+    fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result<Self>;
+    fn calc(&mut self, a: &mut [Self::Elem]) -> Result<&[<Self::Elem as Scalar>::Real]>;
+    fn eval(self, a: &mut [Self::Elem]) -> Result<Vec<<Self::Elem as Scalar>::Real>>;
+}
+
+impl EighWorkImpl for EighWork<c64> {
+    type Elem = c64;
+
+    fn new(calc_eigenvectors: bool, layout: MatrixLayout) -> Result<Self> {
+        assert_eq!(layout.len(), layout.lda());
+        let n = layout.len();
+        let jobz = if calc_eigenvectors {
+            JobEv::All
+        } else {
+            JobEv::None
+        };
+        let uplo = UPLO::Upper; // dummy, working memory is not affected by UPLO
+        let mut eigs = vec_uninit(n as usize);
+        let mut rwork = vec_uninit(3 * n as usize - 2 as usize);
+
+        // calc work size
+        let mut info = 0;
+        let mut work_size = [c64::zero()];
+        unsafe {
+            lapack_sys::zheev_(
+                jobz.as_ptr(),
+                uplo.as_ptr(),
+                &n,
+                std::ptr::null_mut(),
+                &n,
+                AsPtr::as_mut_ptr(&mut eigs),
+                AsPtr::as_mut_ptr(&mut work_size),
+                &(-1),
+                AsPtr::as_mut_ptr(&mut rwork),
+                &mut info,
+            );
+        }
+        info.as_lapack_result()?;
+
+        // actual ev
+        let lwork = work_size[0].to_usize().unwrap();
+        let work = vec_uninit(lwork);
+
+        Ok(EighWork {
+            eigs,
+            jobz,
+            work,
+            rwork: Some(rwork),
+        })
+    }
+
+    fn calc(&mut self, _a: &mut [Self::Elem]) -> Result<&[<Self::Elem as Scalar>::Real]> {
+        todo!()
+    }
+    fn eval(self, _a: &mut [Self::Elem]) -> Result<Vec<<Self::Elem as Scalar>::Real>> {
+        todo!()
+    }
+}
+
 macro_rules! impl_eigh {
     (@real, $scalar:ty, $ev:path) => {
         impl_eigh!(@body, $scalar, $ev, );