Skip to content

Commit 065f81c

Browse files
termoshtttermoshtt
committed
LapackStrict trait for strict memory management API
1 parent 3b37536 commit 065f81c

File tree

4 files changed

+183
-107
lines changed

4 files changed

+183
-107
lines changed

‎lax/src/eigh.rs

Lines changed: 40 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -5,17 +5,18 @@ use crate::{error::*, layout::MatrixLayout};
55
use cauchy::*;
66
use num_traits::{ToPrimitive, Zero};
77

8-
pub trait Eigh: Sized {
9-
type Elem: Scalar;
10-
8+
pub(crate) trait Eigh: Scalar {
119
/// Allocate working memory for eigenvalue problem
12-
fn eigh_work(calc_eigenvec: bool, layout: MatrixLayout, uplo: UPLO) -> Result<Self>;
10+
fn eigh_work(calc_eigenvec: bool, layout: MatrixLayout, uplo: UPLO) -> Result<EighWork<Self>>;
1311

1412
/// Solve eigenvalue problem
15-
fn eigh_calc(&mut self, a: &mut [Self::Elem]) -> Result<&[<Self::Elem as Scalar>::Real]>;
13+
fn eigh_calc<'work>(
14+
work: &'work mut EighWork<Self>,
15+
a: &mut [Self],
16+
) -> Result<&'work [Self::Real]>;
1617
}
1718

18-
/// Working memory for symmetric/Hermitian eigenvalue problem. See [Eigh trait](trait.Eigh.html)
19+
/// Working memory for symmetric/Hermitian eigenvalue problem. See [LapackStrict trait](trait.LapackStrict.html)
1920
pub struct EighWork<T: Scalar> {
2021
jobz: u8,
2122
uplo: UPLO,
@@ -29,17 +30,15 @@ pub struct EighWork<T: Scalar> {
2930

3031
macro_rules! impl_eigh_work_real {
3132
($scalar:ty, $ev:path) => {
32-
impl Eigh for EighWork<$scalar> {
33-
type Elem = $scalar;
34-
35-
fn eigh_work(calc_v: bool, layout: MatrixLayout, uplo: UPLO) -> Result<Self> {
33+
impl Eigh for $scalar {
34+
fn eigh_work(calc_v: bool, layout: MatrixLayout, uplo: UPLO) -> Result<EighWork<Self>> {
3635
assert_eq!(layout.len(), layout.lda());
3736
let n = layout.len();
3837
let jobz = if calc_v { b'V' } else { b'N' };
3938
let mut eigs = unsafe { vec_uninit(n as usize) };
4039

4140
let mut info = 0;
42-
let mut work_size = [Self::Elem::zero()];
41+
let mut work_size = [Self::zero()];
4342
unsafe {
4443
$ev(
4544
jobz,
@@ -66,28 +65,28 @@ macro_rules! impl_eigh_work_real {
6665
})
6766
}
6867

69-
fn eigh_calc(
70-
&mut self,
71-
a: &mut [Self::Elem],
72-
) -> Result<&[<Self::Elem as Scalar>::Real]> {
73-
assert_eq!(a.len(), (self.n * self.n) as usize);
68+
fn eigh_calc<'work>(
69+
work: &'work mut EighWork<Self>,
70+
a: &mut [Self],
71+
) -> Result<&'work [Self::Real]> {
72+
assert_eq!(a.len(), (work.n * work.n) as usize);
7473
let mut info = 0;
75-
let lwork = self.work.len() as i32;
74+
let lwork = work.work.len() as i32;
7675
unsafe {
7776
$ev(
78-
self.jobz,
79-
self.uplo as u8,
80-
self.n,
77+
work.jobz,
78+
work.uplo as u8,
79+
work.n,
8180
a,
82-
self.n,
83-
&mut self.eigs,
84-
&mut self.work,
81+
work.n,
82+
&mut work.eigs,
83+
&mut work.work,
8584
lwork,
8685
&mut info,
8786
);
8887
}
8988
info.as_lapack_result()?;
90-
Ok(&self.eigs)
89+
Ok(&work.eigs)
9190
}
9291
}
9392
};
@@ -98,18 +97,16 @@ impl_eigh_work_real!(f64, lapack::dsyev);
9897

9998
macro_rules! impl_eigh_work_complex {
10099
($scalar:ty, $ev:path) => {
101-
impl Eigh for EighWork<$scalar> {
102-
type Elem = $scalar;
103-
104-
fn eigh_work(calc_v: bool, layout: MatrixLayout, uplo: UPLO) -> Result<Self> {
100+
impl Eigh for $scalar {
101+
fn eigh_work(calc_v: bool, layout: MatrixLayout, uplo: UPLO) -> Result<EighWork<Self>> {
105102
assert_eq!(layout.len(), layout.lda());
106103
let n = layout.len();
107104
let jobz = if calc_v { b'V' } else { b'N' };
108105
let mut eigs = unsafe { vec_uninit(n as usize) };
109106

110107
let mut a = [];
111108
let mut info = 0;
112-
let mut work_size = [Self::Elem::zero()];
109+
let mut work_size = [Self::zero()];
113110
let mut rwork = Vec::with_capacity(3 * n as usize - 2);
114111
unsafe {
115112
$ev(
@@ -138,29 +135,29 @@ macro_rules! impl_eigh_work_complex {
138135
})
139136
}
140137

141-
fn eigh_calc(
142-
&mut self,
143-
a: &mut [Self::Elem],
144-
) -> Result<&[<Self::Elem as Scalar>::Real]> {
145-
assert_eq!(a.len(), (self.n * self.n) as usize);
138+
fn eigh_calc<'work>(
139+
work: &'work mut EighWork<Self>,
140+
a: &mut [Self],
141+
) -> Result<&'work [Self::Real]> {
142+
assert_eq!(a.len(), (work.n * work.n) as usize);
146143
let mut info = 0;
147-
let lwork = self.work.len() as i32;
144+
let lwork = work.work.len() as i32;
148145
unsafe {
149146
$ev(
150-
self.jobz,
151-
self.uplo as u8,
152-
self.n,
147+
work.jobz,
148+
work.uplo as u8,
149+
work.n,
153150
a,
154-
self.n,
155-
&mut self.eigs,
156-
&mut self.work,
151+
work.n,
152+
&mut work.eigs,
153+
&mut work.work,
157154
lwork,
158-
self.rwork.as_mut().unwrap(),
155+
work.rwork.as_mut().unwrap(),
159156
&mut info,
160157
);
161158
}
162159
info.as_lapack_result()?;
163-
Ok(&self.eigs)
160+
Ok(&work.eigs)
164161
}
165162
}
166163
};

‎lax/src/eigh_generalized.rs

Lines changed: 49 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -3,36 +3,25 @@ use crate::{error::*, layout::MatrixLayout};
33
use cauchy::*;
44
use num_traits::{ToPrimitive, Zero};
55

6-
/// Types of generalized eigenvalue problem
7-
#[allow(dead_code)] // FIXME create interface to use ABxlx and BAxlx
8-
#[repr(i32)]
9-
pub enum ITYPE {
10-
/// Solve $ A x = \lambda B x $
11-
AxlBx = 1,
12-
/// Solve $ A B x = \lambda x $
13-
ABxlx = 2,
14-
/// Solve $ B A x = \lambda x $
15-
BAxlx = 3,
16-
}
17-
186
/// Generalized eigenvalue problem for Symmetric/Hermite matrices
19-
pub trait EighGeneralized: Sized {
20-
type Elem: Scalar;
21-
7+
pub(crate) trait EighGeneralized: Scalar {
228
/// Allocate working memory
23-
fn eigh_generalized_work(calc_eigenvec: bool, layout: MatrixLayout, uplo: UPLO)
24-
-> Result<Self>;
9+
fn eigh_generalized_work(
10+
calc_eigenvec: bool,
11+
layout: MatrixLayout,
12+
uplo: UPLO,
13+
) -> Result<EighGeneralizedWork<Self>>;
2514

2615
/// Solve generalized eigenvalue problem
27-
fn eigh_generalized_calc(
28-
&mut self,
29-
a: &mut [Self::Elem],
30-
b: &mut [Self::Elem],
31-
) -> Result<&[<Self::Elem as Scalar>::Real]>;
16+
fn eigh_generalized_calc<'work>(
17+
work: &'work mut EighGeneralizedWork<Self>,
18+
a: &mut [Self],
19+
b: &mut [Self],
20+
) -> Result<&'work [Self::Real]>;
3221
}
3322

3423
/// Working memory for symmetric/Hermitian generalized eigenvalue problem.
35-
/// See [EighGeneralized trait](trait.EighGeneralized.html)
24+
/// See [LapackStrict trait](trait.LapackStrict.html)
3625
pub struct EighGeneralizedWork<T: Scalar> {
3726
jobz: u8,
3827
uplo: UPLO,
@@ -46,21 +35,19 @@ pub struct EighGeneralizedWork<T: Scalar> {
4635

4736
macro_rules! impl_eigh_work_real {
4837
($scalar:ty, $ev:path) => {
49-
impl EighGeneralized for EighGeneralizedWork<$scalar> {
50-
type Elem = $scalar;
51-
38+
impl EighGeneralized for $scalar {
5239
fn eigh_generalized_work(
5340
calc_v: bool,
5441
layout: MatrixLayout,
5542
uplo: UPLO,
56-
) -> Result<Self> {
43+
) -> Result<EighGeneralizedWork<Self>> {
5744
assert_eq!(layout.len(), layout.lda());
5845
let n = layout.len();
5946
let jobz = if calc_v { b'V' } else { b'N' };
6047
let mut eigs = unsafe { vec_uninit(n as usize) };
6148

6249
let mut info = 0;
63-
let mut work_size = [Self::Elem::zero()];
50+
let mut work_size = [Self::zero()];
6451
unsafe {
6552
$ev(
6653
&[ITYPE::AxlBx as i32],
@@ -90,32 +77,32 @@ macro_rules! impl_eigh_work_real {
9077
})
9178
}
9279

93-
fn eigh_generalized_calc(
94-
&mut self,
95-
a: &mut [Self::Elem],
96-
b: &mut [Self::Elem],
97-
) -> Result<&[<Self::Elem as Scalar>::Real]> {
98-
assert_eq!(a.len(), (self.n * self.n) as usize);
80+
fn eigh_generalized_calc<'work>(
81+
work: &'work mut EighGeneralizedWork<Self>,
82+
a: &mut [Self],
83+
b: &mut [Self],
84+
) -> Result<&'work [Self::Real]> {
85+
assert_eq!(a.len(), (work.n * work.n) as usize);
9986
let mut info = 0;
100-
let lwork = self.work.len() as i32;
87+
let lwork = work.work.len() as i32;
10188
unsafe {
10289
$ev(
10390
&[ITYPE::AxlBx as i32],
104-
self.jobz,
105-
self.uplo as u8,
106-
self.n,
91+
work.jobz,
92+
work.uplo as u8,
93+
work.n,
10794
a,
108-
self.n,
95+
work.n,
10996
b,
110-
self.n,
111-
&mut self.eigs,
112-
&mut self.work,
97+
work.n,
98+
&mut work.eigs,
99+
&mut work.work,
113100
lwork,
114101
&mut info,
115102
);
116103
}
117104
info.as_lapack_result()?;
118-
Ok(&self.eigs)
105+
Ok(&work.eigs)
119106
}
120107
}
121108
};
@@ -126,14 +113,12 @@ impl_eigh_work_real!(f64, lapack::dsygv);
126113

127114
macro_rules! impl_eigh_work_complex {
128115
($scalar:ty, $ev:path) => {
129-
impl EighGeneralized for EighGeneralizedWork<$scalar> {
130-
type Elem = $scalar;
131-
116+
impl EighGeneralized for $scalar {
132117
fn eigh_generalized_work(
133118
calc_v: bool,
134119
layout: MatrixLayout,
135120
uplo: UPLO,
136-
) -> Result<Self> {
121+
) -> Result<EighGeneralizedWork<Self>> {
137122
assert_eq!(layout.len(), layout.lda());
138123
let n = layout.len();
139124
let jobz = if calc_v { b'V' } else { b'N' };
@@ -142,7 +127,7 @@ macro_rules! impl_eigh_work_complex {
142127
let mut eigs = unsafe { vec_uninit(n as usize) };
143128

144129
let mut info = 0;
145-
let mut work_size = [Self::Elem::zero()];
130+
let mut work_size = [Self::zero()];
146131
let mut rwork = unsafe { vec_uninit(3 * n as usize - 2) };
147132
unsafe {
148133
$ev(
@@ -174,33 +159,33 @@ macro_rules! impl_eigh_work_complex {
174159
})
175160
}
176161

177-
fn eigh_generalized_calc(
178-
&mut self,
179-
a: &mut [Self::Elem],
180-
b: &mut [Self::Elem],
181-
) -> Result<&[<Self::Elem as Scalar>::Real]> {
182-
assert_eq!(a.len(), (self.n * self.n) as usize);
162+
fn eigh_generalized_calc<'work>(
163+
work: &'work mut EighGeneralizedWork<Self>,
164+
a: &mut [Self],
165+
b: &mut [Self],
166+
) -> Result<&'work [Self::Real]> {
167+
assert_eq!(a.len(), (work.n * work.n) as usize);
183168
let mut info = 0;
184-
let lwork = self.work.len() as i32;
169+
let lwork = work.work.len() as i32;
185170
unsafe {
186171
$ev(
187172
&[ITYPE::AxlBx as i32],
188-
self.jobz,
189-
self.uplo as u8,
190-
self.n,
173+
work.jobz,
174+
work.uplo as u8,
175+
work.n,
191176
a,
192-
self.n,
177+
work.n,
193178
b,
194-
self.n,
195-
&mut self.eigs,
196-
&mut self.work,
179+
work.n,
180+
&mut work.eigs,
181+
&mut work.work,
197182
lwork,
198-
self.rwork.as_mut().unwrap(),
183+
work.rwork.as_mut().unwrap(),
199184
&mut info,
200185
);
201186
}
202187
info.as_lapack_result()?;
203-
Ok(&self.eigs)
188+
Ok(&work.eigs)
204189
}
205190
}
206191
};

‎lax/src/lib.rs

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,7 @@ pub mod qr;
8080
pub mod rcond;
8181
pub mod solve;
8282
pub mod solveh;
83+
mod strict;
8384
pub mod svd;
8485
pub mod svddc;
8586
mod traits;
@@ -95,6 +96,7 @@ pub use self::qr::*;
9596
pub use self::rcond::*;
9697
pub use self::solve::*;
9798
pub use self::solveh::*;
99+
pub use self::strict::*;
98100
pub use self::svd::*;
99101
pub use self::svddc::*;
100102
pub use self::traits::*;
@@ -146,6 +148,18 @@ impl NormType {
146148
}
147149
}
148150

151+
/// Types of generalized eigenvalue problem
152+
#[allow(dead_code)] // FIXME create interface to use ABxlx and BAxlx
153+
#[repr(i32)]
154+
pub enum ITYPE {
155+
/// Solve $ A x = \lambda B x $
156+
AxlBx = 1,
157+
/// Solve $ A B x = \lambda x $
158+
ABxlx = 2,
159+
/// Solve $ B A x = \lambda x $
160+
BAxlx = 3,
161+
}
162+
149163
/// Create a vector without initialization
150164
///
151165
/// Safety

0 commit comments

Comments
 (0)