Skip to content

Commit 156d920

Browse files
bangconghuynhbangconghuynh
committed
Introduced an enum to handle finite and indeterminate generalized eigenvalues.
1 parent f1ed32a commit 156d920

File tree

4 files changed

+66
-37
lines changed

4 files changed

+66
-37
lines changed

‎lax/src/eig_generalized.rs

Lines changed: 25 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,19 @@
99
//!
1010
use std::mem::MaybeUninit;
1111

12-
use crate::{error::*, layout::MatrixLayout, *};
1312
use crate::eig::reconstruct_eigenvectors;
13+
use crate::{error::*, layout::MatrixLayout, *};
1414
use cauchy::*;
1515
use num_traits::{ToPrimitive, Zero};
1616

17+
pub enum GeneralizedEigenvalue<T: Scalar> {
18+
/// Finite generalized eigenvalue: `Finite(α/β, (α, β))`
19+
Finite(T, (T, T)),
20+
21+
/// Indeterminate generalized eigenvalue: `Indeterminate((α, β))`
22+
Indeterminate((T, T)),
23+
}
24+
1725
#[non_exhaustive]
1826
pub struct EigGeneralizedWork<T: Scalar> {
1927
/// Problem size
@@ -254,22 +262,28 @@ macro_rules! impl_eig_generalized_work_c {
254262
}
255263

256264
impl EigGeneralizedOwned<$c> {
257-
pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<Option<$c>> {
265+
pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<GeneralizedEigenvalue<$c>> {
258266
self.alpha
259267
.iter()
260268
.zip(self.beta.iter())
261269
.map(|(alpha, beta)| {
262270
if let Some(thresh) = thresh_opt {
263271
if beta.abs() < thresh {
264-
None
272+
GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
265273
} else {
266-
Some(alpha / beta)
274+
GeneralizedEigenvalue::Finite(
275+
alpha / beta,
276+
(alpha.clone(), beta.clone()),
277+
)
267278
}
268279
} else {
269280
if beta.is_zero() {
270-
None
281+
GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
271282
} else {
272-
Some(alpha / beta)
283+
GeneralizedEigenvalue::Finite(
284+
alpha / beta,
285+
(alpha.clone(), beta.clone()),
286+
)
273287
}
274288
}
275289
})
@@ -441,22 +455,22 @@ macro_rules! impl_eig_generalized_work_r {
441455
}
442456

443457
impl EigGeneralizedOwned<$f> {
444-
pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<Option<$c>> {
458+
pub fn calc_eigs(&self, thresh_opt: Option<$f>) -> Vec<GeneralizedEigenvalue<$c>> {
445459
self.alpha
446460
.iter()
447461
.zip(self.beta.iter())
448462
.map(|(alpha, beta)| {
449463
if let Some(thresh) = thresh_opt {
450464
if beta.abs() < thresh {
451-
None
465+
GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
452466
} else {
453-
Some(alpha / beta)
467+
GeneralizedEigenvalue::Finite(alpha / beta, (alpha.clone(), beta.clone()))
454468
}
455469
} else {
456470
if beta.is_zero() {
457-
None
471+
GeneralizedEigenvalue::Indeterminate((alpha.clone(), beta.clone()))
458472
} else {
459-
Some(alpha / beta)
473+
GeneralizedEigenvalue::Finite(alpha / beta, (alpha.clone(), beta.clone()))
460474
}
461475
}
462476
})

‎lax/src/lib.rs

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,8 @@ pub mod svddc;
105105
pub mod triangular;
106106
pub mod tridiagonal;
107107

108+
pub use crate::eig_generalized::GeneralizedEigenvalue;
109+
108110
pub use self::flags::*;
109111
pub use self::least_squares::LeastSquaresOwned;
110112
pub use self::svd::{SvdOwned, SvdRef};
@@ -133,7 +135,7 @@ pub trait Lapack: Scalar {
133135
a: &mut [Self],
134136
b: &mut [Self],
135137
thresh_opt: Option<Self::Real>,
136-
) -> Result<(Vec<Option<Self::Complex>>, Vec<Self::Complex>)>;
138+
) -> Result<(Vec<GeneralizedEigenvalue<Self::Complex>>, Vec<Self::Complex>)>;
137139

138140
/// Compute right eigenvalue and eigenvectors for a symmetric or Hermitian matrix
139141
fn eigh(
@@ -348,7 +350,7 @@ macro_rules! impl_lapack {
348350
a: &mut [Self],
349351
b: &mut [Self],
350352
thresh_opt: Option<Self::Real>,
351-
) -> Result<(Vec<Option<Self::Complex>>, Vec<Self::Complex>)> {
353+
) -> Result<(Vec<GeneralizedEigenvalue<Self::Complex>>, Vec<Self::Complex>)> {
352354
use eig_generalized::*;
353355
let work = EigGeneralizedWork::<$s>::new(calc_v, l)?;
354356
let eig_generalized_owned = work.eval(a, b)?;

‎ndarray-linalg/src/eig.rs

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
use crate::error::*;
44
use crate::layout::*;
55
use crate::types::*;
6+
pub use lax::GeneralizedEigenvalue;
67
use ndarray::*;
78

89
#[cfg_attr(doc, katexit::katexit)]
@@ -107,12 +108,12 @@ pub trait EigGeneralized {
107108
/// [ 4.44, -7.77, 0.00, 1.11, 5.55],
108109
/// [-8.88, 6.66, -3.33, 2.22, -9.99],
109110
/// ];
110-
/// let (eigs, vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
111+
/// let (geneigs, vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
111112
///
112113
/// let a = a.map(|v| v.as_c());
113114
/// let b = b.map(|v| v.as_c());
114-
/// for (e_opt, vec) in eigs.iter().zip(vecs.axis_iter(Axis(1))) {
115-
/// if let Some(e) = e_opt.as_ref() {
115+
/// for (ge, vec) in geneigs.iter().zip(vecs.axis_iter(Axis(1))) {
116+
/// if let GeneralizedEigenvalue::Finite(e, _) = ge {
116117
/// let ebv = b.dot(&vec).map(|v| v * e);
117118
/// let av = a.dot(&vec);
118119
/// assert_close_l2!(&av, &ebv, 1e-5);
@@ -136,7 +137,7 @@ where
136137
A: Scalar + Lapack,
137138
S: Data<Elem = A>,
138139
{
139-
type EigVal = Array1<Option<A::Complex>>;
140+
type EigVal = Array1<GeneralizedEigenvalue<A::Complex>>;
140141
type EigVec = Array2<A::Complex>;
141142
type Real = A::Real;
142143

‎ndarray-linalg/tests/eig_generalized.rs

Lines changed: 32 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -15,21 +15,24 @@ fn real_a_real_b_3x3_full_rank() {
1515
[-3.0, 1.0, 6.0],
1616
[ 4.0, -5.0, 1.0],
1717
];
18-
let (eigvals_opt, eigvecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
18+
let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
1919

2020
let a = a.map(|v| v.as_c());
2121
let b = b.map(|v| v.as_c());
22-
for (e_opt, vec) in eigvals_opt.iter().zip(eigvecs.columns()) {
23-
if let Some(e) = e_opt.as_ref() {
22+
for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
23+
if let GeneralizedEigenvalue::Finite(e, _) = ge {
2424
let ebv = b.dot(&vec).map(|v| v * e);
2525
let av = a.dot(&vec);
2626
assert_close_l2!(&av, &ebv, 1e-7);
2727
}
2828
}
2929

30-
let mut eigvals = eigvals_opt
30+
let mut eigvals = geneigvals
3131
.iter()
32-
.filter_map(|&e_opt: &Option<c64>| e_opt)
32+
.filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
33+
GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
34+
GeneralizedEigenvalue::Indeterminate(_) => None,
35+
})
3336
.collect::<Vec<_>>();
3437
eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
3538
let eigvals = Array1::from_vec(eigvals);
@@ -55,21 +58,24 @@ fn real_a_real_b_3x3_nullity_1() {
5558
[0.0, 1.0, 1.0],
5659
[1.0, -1.0, 0.0],
5760
];
58-
let (eigvals_opt, eigvecs) = (a.clone(), b.clone()).eig_generalized(Some(1e-4)).unwrap();
61+
let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(Some(1e-4)).unwrap();
5962

6063
let a = a.map(|v| v.as_c());
6164
let b = b.map(|v| v.as_c());
62-
for (e_opt, vec) in eigvals_opt.iter().zip(eigvecs.columns()) {
63-
if let Some(e) = e_opt.as_ref() {
65+
for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
66+
if let GeneralizedEigenvalue::Finite(e, _) = ge {
6467
let ebv = b.dot(&vec).map(|v| v * e);
6568
let av = a.dot(&vec);
6669
assert_close_l2!(&av, &ebv, 1e-7);
6770
}
6871
}
6972

70-
let mut eigvals = eigvals_opt
73+
let mut eigvals = geneigvals
7174
.iter()
72-
.filter_map(|&e_opt: &Option<c64>| e_opt)
75+
.filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
76+
GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
77+
GeneralizedEigenvalue::Indeterminate(_) => None,
78+
})
7379
.collect::<Vec<_>>();
7480
eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
7581
let eigvals = Array1::from_vec(eigvals);
@@ -95,21 +101,24 @@ fn complex_a_complex_b_3x3_full_rank() {
95101
[c64::new( 0.0, -3.0), c64::new( 2.0, 2.0), c64::new(-4.0, 0.0)],
96102
[c64::new( 5.0, 5.0), c64::new(-1.5, 1.5), c64::new( 0.0, -2.0)],
97103
];
98-
let (eigvals_opt, eigvecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
104+
let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
99105

100106
let a = a.map(|v| v.as_c());
101107
let b = b.map(|v| v.as_c());
102-
for (e_opt, vec) in eigvals_opt.iter().zip(eigvecs.columns()) {
103-
if let Some(e) = e_opt.as_ref() {
108+
for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
109+
if let GeneralizedEigenvalue::Finite(e, _) = ge {
104110
let ebv = b.dot(&vec).map(|v| v * e);
105111
let av = a.dot(&vec);
106112
assert_close_l2!(&av, &ebv, 1e-7);
107113
}
108114
}
109115

110-
let mut eigvals = eigvals_opt
116+
let mut eigvals = geneigvals
111117
.iter()
112-
.filter_map(|&e_opt: &Option<c64>| e_opt)
118+
.filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
119+
GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
120+
GeneralizedEigenvalue::Indeterminate(_) => None,
121+
})
113122
.collect::<Vec<_>>();
114123
eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
115124
let eigvals = Array1::from_vec(eigvals);
@@ -139,21 +148,24 @@ fn complex_a_complex_b_3x3_nullity_1() {
139148
[c64::new( 7.85029, 7.02144), c64::new(9.23225, -0.479451), c64::new(13.9507, -16.5402)],
140149
[c64::new(-4.47803, 3.98981), c64::new(9.44434, -4.519970), c64::new(40.9006, -23.5060)],
141150
];
142-
let (eigvals_opt, eigvecs) = (a.clone(), b.clone()).eig_generalized(Some(1e-4)).unwrap();
151+
let (geneigvals, eigvecs) = (a.clone(), b.clone()).eig_generalized(Some(1e-4)).unwrap();
143152

144153
let a = a.map(|v| v.as_c());
145154
let b = b.map(|v| v.as_c());
146-
for (e_opt, vec) in eigvals_opt.iter().zip(eigvecs.columns()) {
147-
if let Some(e) = e_opt.as_ref() {
155+
for (ge, vec) in geneigvals.iter().zip(eigvecs.columns()) {
156+
if let GeneralizedEigenvalue::Finite(e, _) = ge {
148157
let ebv = b.dot(&vec).map(|v| v * e);
149158
let av = a.dot(&vec);
150159
assert_close_l2!(&av, &ebv, 1e-7);
151160
}
152161
}
153162

154-
let mut eigvals = eigvals_opt
163+
let mut eigvals = geneigvals
155164
.iter()
156-
.filter_map(|&e_opt: &Option<c64>| e_opt)
165+
.filter_map(|ge: &GeneralizedEigenvalue<c64>| match ge {
166+
GeneralizedEigenvalue::Finite(e, _) => Some(e.clone()),
167+
GeneralizedEigenvalue::Indeterminate(_) => None,
168+
})
157169
.collect::<Vec<_>>();
158170
eigvals.sort_by(|a, b| a.re().partial_cmp(&b.re()).unwrap());
159171
let eigvals = Array1::from_vec(eigvals);

0 commit comments

Comments
 (0)