lax/
eig_generalized.rs

1//! Generalized eigenvalue problem for general matrices
2//!
3//! LAPACK correspondance
4//! ----------------------
5//!
6//! | f32   | f64   | c32   | c64   |
7//! |:------|:------|:------|:------|
8//! | sggev | dggev | cggev | zggev |
9//!
10use std::mem::MaybeUninit;
11
12use crate::eig::reconstruct_eigenvectors;
13use crate::{error::*, layout::MatrixLayout, *};
14use cauchy::*;
15use num_traits::{ToPrimitive, Zero};
16
17#[derive(Clone, PartialEq, Eq)]
18pub enum GeneralizedEigenvalue<T: Scalar> {
19    /// Finite generalized eigenvalue: `Finite(α/β, (α, β))`
20    Finite(T, (T, T)),
21
22    /// Indeterminate generalized eigenvalue: `Indeterminate((α, β))`
23    Indeterminate((T, T)),
24}
25
26impl<T: Scalar> std::fmt::Display for GeneralizedEigenvalue<T> {
27    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
28        match self {
29            Self::Finite(e, (a, b)) => write!(f, "{e:.3e} ({a:.3e}/{b:.3e})"),
30            Self::Indeterminate((a, b)) => write!(f, "∞ ({a:.3e}/{b:.3e})"),
31        }
32    }
33}
34
35#[non_exhaustive]
36pub struct EigGeneralizedWork<T: Scalar> {
37    /// Problem size
38    pub n: i32,
39    /// Compute right eigenvectors or not
40    pub jobvr: JobEv,
41    /// Compute left eigenvectors or not
42    pub jobvl: JobEv,
43
44    /// Eigenvalues: alpha (numerators)
45    pub alpha: Vec<MaybeUninit<T::Complex>>,
46    /// Eigenvalues: beta (denominators)
47    pub beta: Vec<MaybeUninit<T::Complex>>,
48    /// Real part of alpha (eigenvalue numerators) used in real routines
49    pub alpha_re: Option<Vec<MaybeUninit<T::Real>>>,
50    /// Imaginary part of alpha (eigenvalue numerators) used in real routines
51    pub alpha_im: Option<Vec<MaybeUninit<T::Real>>>,
52    /// Real part of beta (eigenvalue denominators) used in real routines
53    pub beta_re: Option<Vec<MaybeUninit<T::Real>>>,
54    /// Imaginary part of beta (eigenvalue denominators) used in real routines
55    pub beta_im: Option<Vec<MaybeUninit<T::Real>>>,
56
57    /// Left eigenvectors
58    pub vc_l: Option<Vec<MaybeUninit<T::Complex>>>,
59    /// Left eigenvectors used in real routines
60    pub vr_l: Option<Vec<MaybeUninit<T::Real>>>,
61    /// Right eigenvectors
62    pub vc_r: Option<Vec<MaybeUninit<T::Complex>>>,
63    /// Right eigenvectors used in real routines
64    pub vr_r: Option<Vec<MaybeUninit<T::Real>>>,
65
66    /// Working memory
67    pub work: Vec<MaybeUninit<T>>,
68    /// Working memory with `T::Real`
69    pub rwork: Option<Vec<MaybeUninit<T::Real>>>,
70}
71
72impl<T> EigGeneralizedWork<T>
73where
74    T: Scalar,
75    EigGeneralizedWork<T>: EigGeneralizedWorkImpl<Elem = T>,
76{
77    /// Create new working memory for eigenvalues compution.
78    pub fn new(calc_v: bool, l: MatrixLayout) -> Result<Self> {
79        EigGeneralizedWorkImpl::new(calc_v, l)
80    }
81
82    /// Compute eigenvalues and vectors on this working memory.
83    pub fn calc(&mut self, a: &mut [T], b: &mut [T]) -> Result<EigGeneralizedRef<T>> {
84        EigGeneralizedWorkImpl::calc(self, a, b)
85    }
86
87    /// Compute eigenvalues and vectors by consuming this working memory.
88    pub fn eval(self, a: &mut [T], b: &mut [T]) -> Result<EigGeneralizedOwned<T>> {
89        EigGeneralizedWorkImpl::eval(self, a, b)
90    }
91}
92
93/// Owned result of eigenvalue problem by [EigGeneralizedWork::eval]
94#[derive(Debug, Clone, PartialEq)]
95pub struct EigGeneralizedOwned<T: Scalar> {
96    /// Eigenvalues
97    pub alpha: Vec<T::Complex>,
98
99    pub beta: Vec<T::Complex>,
100
101    /// Right eigenvectors
102    pub vr: Option<Vec<T::Complex>>,
103
104    /// Left eigenvectors
105    pub vl: Option<Vec<T::Complex>>,
106}
107
108/// Reference result of eigenvalue problem by [EigGeneralizedWork::calc]
109#[derive(Debug, Clone, PartialEq)]
110pub struct EigGeneralizedRef<'work, T: Scalar> {
111    /// Eigenvalues
112    pub alpha: &'work [T::Complex],
113
114    pub beta: &'work [T::Complex],
115
116    /// Right eigenvectors
117    pub vr: Option<&'work [T::Complex]>,
118
119    /// Left eigenvectors
120    pub vl: Option<&'work [T::Complex]>,
121}
122
123/// Helper trait for implementing [EigGeneralizedWork] methods
124pub trait EigGeneralizedWorkImpl: Sized {
125    type Elem: Scalar;
126    fn new(calc_v: bool, l: MatrixLayout) -> Result<Self>;
127    fn calc<'work>(
128        &'work mut self,
129        a: &mut [Self::Elem],
130        b: &mut [Self::Elem],
131    ) -> Result<EigGeneralizedRef<'work, Self::Elem>>;
132    fn eval(
133        self,
134        a: &mut [Self::Elem],
135        b: &mut [Self::Elem],
136    ) -> Result<EigGeneralizedOwned<Self::Elem>>;
137}
138
139macro_rules! impl_eig_generalized_work_c {
140    ($f:ty, $c:ty, $ev:path) => {
141        impl EigGeneralizedWorkImpl for EigGeneralizedWork<$c> {
142            type Elem = $c;
143
144            fn new(calc_v: bool, l: MatrixLayout) -> Result<Self> {
145                let (n, _) = l.size();
146                let (jobvl, jobvr) = if calc_v {
147                    match l {
148                        MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
149                        MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
150                    }
151                } else {
152                    (JobEv::None, JobEv::None)
153                };
154                let mut rwork = vec_uninit(8 * n as usize);
155
156                let mut alpha = vec_uninit(n as usize);
157                let mut beta = vec_uninit(n as usize);
158
159                let mut vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
160                let mut vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
161
162                // calc work size
163                let mut info = 0;
164                let mut work_size = [<$c>::zero()];
165                unsafe {
166                    $ev(
167                        jobvl.as_ptr(),
168                        jobvr.as_ptr(),
169                        &n,
170                        std::ptr::null_mut(),
171                        &n,
172                        std::ptr::null_mut(),
173                        &n,
174                        AsPtr::as_mut_ptr(&mut alpha),
175                        AsPtr::as_mut_ptr(&mut beta),
176                        AsPtr::as_mut_ptr(vc_l.as_deref_mut().unwrap_or(&mut [])),
177                        &n,
178                        AsPtr::as_mut_ptr(vc_r.as_deref_mut().unwrap_or(&mut [])),
179                        &n,
180                        AsPtr::as_mut_ptr(&mut work_size),
181                        &(-1),
182                        AsPtr::as_mut_ptr(&mut rwork),
183                        &mut info,
184                    )
185                };
186                info.as_lapack_result()?;
187
188                let lwork = work_size[0].to_usize().unwrap();
189                let work: Vec<MaybeUninit<$c>> = vec_uninit(lwork);
190                Ok(Self {
191                    n,
192                    jobvl,
193                    jobvr,
194                    alpha,
195                    beta,
196                    alpha_re: None,
197                    alpha_im: None,
198                    beta_re: None,
199                    beta_im: None,
200                    rwork: Some(rwork),
201                    vc_l,
202                    vc_r,
203                    vr_l: None,
204                    vr_r: None,
205                    work,
206                })
207            }
208
209            fn calc<'work>(
210                &'work mut self,
211                a: &mut [Self::Elem],
212                b: &mut [Self::Elem],
213            ) -> Result<EigGeneralizedRef<'work, Self::Elem>> {
214                let lwork = self.work.len().to_i32().unwrap();
215                let mut info = 0;
216                unsafe {
217                    $ev(
218                        self.jobvl.as_ptr(),
219                        self.jobvr.as_ptr(),
220                        &self.n,
221                        AsPtr::as_mut_ptr(a),
222                        &self.n,
223                        AsPtr::as_mut_ptr(b),
224                        &self.n,
225                        AsPtr::as_mut_ptr(&mut self.alpha),
226                        AsPtr::as_mut_ptr(&mut self.beta),
227                        AsPtr::as_mut_ptr(self.vc_l.as_deref_mut().unwrap_or(&mut [])),
228                        &self.n,
229                        AsPtr::as_mut_ptr(self.vc_r.as_deref_mut().unwrap_or(&mut [])),
230                        &self.n,
231                        AsPtr::as_mut_ptr(&mut self.work),
232                        &lwork,
233                        AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
234                        &mut info,
235                    )
236                };
237                info.as_lapack_result()?;
238                // Hermite conjugate
239                if let Some(vl) = self.vc_l.as_mut() {
240                    for value in vl {
241                        let value = unsafe { value.assume_init_mut() };
242                        value.im = -value.im;
243                    }
244                }
245                Ok(EigGeneralizedRef {
246                    alpha: unsafe { self.alpha.slice_assume_init_ref() },
247                    beta: unsafe { self.beta.slice_assume_init_ref() },
248                    vl: self
249                        .vc_l
250                        .as_ref()
251                        .map(|v| unsafe { v.slice_assume_init_ref() }),
252                    vr: self
253                        .vc_r
254                        .as_ref()
255                        .map(|v| unsafe { v.slice_assume_init_ref() }),
256                })
257            }
258
259            fn eval(
260                mut self,
261                a: &mut [Self::Elem],
262                b: &mut [Self::Elem],
263            ) -> Result<EigGeneralizedOwned<Self::Elem>> {
264                let _eig_generalized_ref = self.calc(a, b)?;
265                Ok(EigGeneralizedOwned {
266                    alpha: unsafe { self.alpha.assume_init() },
267                    beta: unsafe { self.beta.assume_init() },
268                    vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
269                    vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
270                })
271            }
272        }
273
274        impl EigGeneralizedOwned<$c> {
275            pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<GeneralizedEigenvalue<$c>> {
276                self.alpha
277                    .iter()
278                    .zip(self.beta.iter())
279                    .map(|(alpha, beta)| {
280                        if let Some(thresh) = thresh_opt {
281                            if beta.abs() < thresh {
282                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
283                            } else {
284                                GeneralizedEigenvalue::Finite(
285                                    alpha / beta,
286                                    (alpha.clone(), beta.clone()),
287                                )
288                            }
289                        } else {
290                            if beta.is_zero() {
291                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
292                            } else {
293                                GeneralizedEigenvalue::Finite(
294                                    alpha / beta,
295                                    (alpha.clone(), beta.clone()),
296                                )
297                            }
298                        }
299                    })
300                    .collect::<Vec<_>>()
301            }
302        }
303    };
304}
305
306impl_eig_generalized_work_c!(f32, c32, lapack_sys::cggev_);
307impl_eig_generalized_work_c!(f64, c64, lapack_sys::zggev_);
308
309macro_rules! impl_eig_generalized_work_r {
310    ($f:ty, $c:ty, $ev:path) => {
311        impl EigGeneralizedWorkImpl for EigGeneralizedWork<$f> {
312            type Elem = $f;
313
314            fn new(calc_v: bool, l: MatrixLayout) -> Result<Self> {
315                let (n, _) = l.size();
316                let (jobvl, jobvr) = if calc_v {
317                    match l {
318                        MatrixLayout::C { .. } => (JobEv::All, JobEv::None),
319                        MatrixLayout::F { .. } => (JobEv::None, JobEv::All),
320                    }
321                } else {
322                    (JobEv::None, JobEv::None)
323                };
324                let mut alpha_re = vec_uninit(n as usize);
325                let mut alpha_im = vec_uninit(n as usize);
326                let mut beta_re = vec_uninit(n as usize);
327                let mut vr_l = jobvl.then(|| vec_uninit((n * n) as usize));
328                let mut vr_r = jobvr.then(|| vec_uninit((n * n) as usize));
329                let vc_l = jobvl.then(|| vec_uninit((n * n) as usize));
330                let vc_r = jobvr.then(|| vec_uninit((n * n) as usize));
331
332                // calc work size
333                let mut info = 0;
334                let mut work_size: [$f; 1] = [0.0];
335                unsafe {
336                    $ev(
337                        jobvl.as_ptr(),
338                        jobvr.as_ptr(),
339                        &n,
340                        std::ptr::null_mut(),
341                        &n,
342                        std::ptr::null_mut(),
343                        &n,
344                        AsPtr::as_mut_ptr(&mut alpha_re),
345                        AsPtr::as_mut_ptr(&mut alpha_im),
346                        AsPtr::as_mut_ptr(&mut beta_re),
347                        AsPtr::as_mut_ptr(vr_l.as_deref_mut().unwrap_or(&mut [])),
348                        &n,
349                        AsPtr::as_mut_ptr(vr_r.as_deref_mut().unwrap_or(&mut [])),
350                        &n,
351                        AsPtr::as_mut_ptr(&mut work_size),
352                        &(-1),
353                        &mut info,
354                    )
355                };
356                info.as_lapack_result()?;
357
358                // actual ev
359                let lwork = work_size[0].to_usize().unwrap();
360                let work = vec_uninit(lwork);
361
362                Ok(Self {
363                    n,
364                    jobvr,
365                    jobvl,
366                    alpha: vec_uninit(n as usize),
367                    beta: vec_uninit(n as usize),
368                    alpha_re: Some(alpha_re),
369                    alpha_im: Some(alpha_im),
370                    beta_re: Some(beta_re),
371                    beta_im: None,
372                    rwork: None,
373                    vr_l,
374                    vr_r,
375                    vc_l,
376                    vc_r,
377                    work,
378                })
379            }
380
381            fn calc<'work>(
382                &'work mut self,
383                a: &mut [Self::Elem],
384                b: &mut [Self::Elem],
385            ) -> Result<EigGeneralizedRef<'work, Self::Elem>> {
386                let lwork = self.work.len().to_i32().unwrap();
387                let mut info = 0;
388                unsafe {
389                    $ev(
390                        self.jobvl.as_ptr(),
391                        self.jobvr.as_ptr(),
392                        &self.n,
393                        AsPtr::as_mut_ptr(a),
394                        &self.n,
395                        AsPtr::as_mut_ptr(b),
396                        &self.n,
397                        AsPtr::as_mut_ptr(self.alpha_re.as_mut().unwrap()),
398                        AsPtr::as_mut_ptr(self.alpha_im.as_mut().unwrap()),
399                        AsPtr::as_mut_ptr(self.beta_re.as_mut().unwrap()),
400                        AsPtr::as_mut_ptr(self.vr_l.as_deref_mut().unwrap_or(&mut [])),
401                        &self.n,
402                        AsPtr::as_mut_ptr(self.vr_r.as_deref_mut().unwrap_or(&mut [])),
403                        &self.n,
404                        AsPtr::as_mut_ptr(&mut self.work),
405                        &lwork,
406                        &mut info,
407                    )
408                };
409                info.as_lapack_result()?;
410
411                let alpha_re = self
412                    .alpha_re
413                    .as_ref()
414                    .map(|e| unsafe { e.slice_assume_init_ref() })
415                    .unwrap();
416                let alpha_im = self
417                    .alpha_im
418                    .as_ref()
419                    .map(|e| unsafe { e.slice_assume_init_ref() })
420                    .unwrap();
421                let beta_re = self
422                    .beta_re
423                    .as_ref()
424                    .map(|e| unsafe { e.slice_assume_init_ref() })
425                    .unwrap();
426                reconstruct_eigs_optional_im(alpha_re, Some(alpha_im), &mut self.alpha);
427                reconstruct_eigs_optional_im(beta_re, None, &mut self.beta);
428
429                if let Some(v) = self.vr_l.as_ref() {
430                    let v = unsafe { v.slice_assume_init_ref() };
431                    reconstruct_eigenvectors(true, alpha_im, v, self.vc_l.as_mut().unwrap());
432                }
433                if let Some(v) = self.vr_r.as_ref() {
434                    let v = unsafe { v.slice_assume_init_ref() };
435                    reconstruct_eigenvectors(false, alpha_im, v, self.vc_r.as_mut().unwrap());
436                }
437
438                Ok(EigGeneralizedRef {
439                    alpha: unsafe { self.alpha.slice_assume_init_ref() },
440                    beta: unsafe { self.beta.slice_assume_init_ref() },
441                    vl: self
442                        .vc_l
443                        .as_ref()
444                        .map(|v| unsafe { v.slice_assume_init_ref() }),
445                    vr: self
446                        .vc_r
447                        .as_ref()
448                        .map(|v| unsafe { v.slice_assume_init_ref() }),
449                })
450            }
451
452            fn eval(
453                mut self,
454                a: &mut [Self::Elem],
455                b: &mut [Self::Elem],
456            ) -> Result<EigGeneralizedOwned<Self::Elem>> {
457                let _eig_generalized_ref = self.calc(a, b)?;
458                Ok(EigGeneralizedOwned {
459                    alpha: unsafe { self.alpha.assume_init() },
460                    beta: unsafe { self.beta.assume_init() },
461                    vl: self.vc_l.map(|v| unsafe { v.assume_init() }),
462                    vr: self.vc_r.map(|v| unsafe { v.assume_init() }),
463                })
464            }
465        }
466
467        impl EigGeneralizedOwned<$f> {
468            pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<GeneralizedEigenvalue<$c>> {
469                self.alpha
470                    .iter()
471                    .zip(self.beta.iter())
472                    .map(|(alpha, beta)| {
473                        if let Some(thresh) = thresh_opt {
474                            if beta.abs() < thresh {
475                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
476                            } else {
477                                GeneralizedEigenvalue::Finite(
478                                    alpha / beta,
479                                    (alpha.clone(), beta.clone()),
480                                )
481                            }
482                        } else {
483                            if beta.is_zero() {
484                                GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
485                            } else {
486                                GeneralizedEigenvalue::Finite(
487                                    alpha / beta,
488                                    (alpha.clone(), beta.clone()),
489                                )
490                            }
491                        }
492                    })
493                    .collect::<Vec<_>>()
494            }
495        }
496    };
497}
498impl_eig_generalized_work_r!(f32, c32, lapack_sys::sggev_);
499impl_eig_generalized_work_r!(f64, c64, lapack_sys::dggev_);
500
501/// Create complex eigenvalues from real and optional imaginary parts.
502fn reconstruct_eigs_optional_im<T: Scalar>(
503    re: &[T],
504    im_opt: Option<&[T]>,
505    eigs: &mut [MaybeUninit<T::Complex>],
506) {
507    let n = eigs.len();
508    assert_eq!(re.len(), n);
509
510    if let Some(im) = im_opt {
511        assert_eq!(im.len(), n);
512        for i in 0..n {
513            eigs[i].write(T::complex(re[i], im[i]));
514        }
515    } else {
516        for i in 0..n {
517            eigs[i].write(T::complex(re[i], T::zero()));
518        }
519    }
520}
lax/eig_generalized.rs

lax/
eig_generalized.rs
