Skip to content

Commit 419ce95

Browse files
termoshtttermoshtt
committed
Merge branch 'master' into fix-eig
2 parents d944821 + 685fc3c commit 419ce95

File tree

10 files changed

+426
-46
lines changed

10 files changed

+426
-46
lines changed

‎CHANGELOG.md

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,16 @@ Updated dependencies
66
- ndarray 0.15 https://github.com/rust-ndarray/ndarray-linalg/pull/273
77
- cauchy 0.4 (num-complex 0.4, rand 0.8), lapack 0.18 https://github.com/rust-ndarray/ndarray-linalg/pull/276
88

9+
Fixed
10+
-----
11+
- Fix Solve::solve_h_* for complex inputs with standard layout https://github.com/rust-ndarray/ndarray-linalg/pull/296
12+
- Add checks for matching shapes in Solve, SolveH, and EighInplace https://github.com/rust-ndarray/ndarray-linalg/pull/290
13+
914
Changed
1015
--------
16+
- Avoid unnecessary calculation of eigenvectors in EigVals::eigvals and relax the DataMut bound https://github.com/rust-ndarray/ndarray-linalg/pull/286
17+
- Add basic documentation for eigh module https://github.com/rust-ndarray/ndarray-linalg/pull/283
18+
- Port usage of uninitialized to build_uninit https://github.com/rust-ndarray/ndarray-linalg/pull/287
1119
- Relax type bounds for LeastSquaresSvd family https://github.com/rust-ndarray/ndarray-linalg/pull/272
1220

1321
0.13.1 - 13 March 2021

‎lax/src/solve.rs

Lines changed: 35 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -75,18 +75,49 @@ macro_rules! impl_solve {
7575
ipiv: &Pivot,
7676
b: &mut [Self],
7777
) -> Result<()> {
78-
let t = match l {
78+
// If the array has C layout, then it needs to be handled
79+
// specially, since LAPACK expects a Fortran-layout array.
80+
// Reinterpreting a C layout array as Fortran layout is
81+
// equivalent to transposing it. So, we can handle the "no
82+
// transpose" and "transpose" cases by swapping to "transpose"
83+
// or "no transpose", respectively. For the "Hermite" case, we
84+
// can take advantage of the following:
85+
//
86+
// ```text
87+
// A^H x = b
88+
// ⟺ conj(A^T) x = b
89+
// ⟺ conj(conj(A^T) x) = conj(b)
90+
// ⟺ conj(conj(A^T)) conj(x) = conj(b)
91+
// ⟺ A^T conj(x) = conj(b)
92+
// ```
93+
//
94+
// So, we can handle this case by switching to "no transpose"
95+
// (which is equivalent to transposing the array since it will
96+
// be reinterpreted as Fortran layout) and applying the
97+
// elementwise conjugate to `x` and `b`.
98+
let (t, conj) = match l {
7999
MatrixLayout::C { .. } => match t {
80-
Transpose::No => Transpose::Transpose,
81-
Transpose::Transpose | Transpose::Hermite => Transpose::No,
100+
Transpose::No => (Transpose::Transpose, false),
101+
Transpose::Transpose => (Transpose::No, false),
102+
Transpose::Hermite => (Transpose::No, true),
82103
},
83-
_ => t,
104+
MatrixLayout::F { .. } => (t, false),
84105
};
85106
let (n, _) = l.size();
86107
let nrhs = 1;
87108
let ldb = l.lda();
88109
let mut info = 0;
110+
if conj {
111+
for b_elem in &mut *b {
112+
*b_elem = b_elem.conj();
113+
}
114+
}
89115
unsafe { $getrs(t as u8, n, nrhs, a, l.lda(), ipiv, b, ldb, &mut info) };
116+
if conj {
117+
for b_elem in &mut *b {
118+
*b_elem = b_elem.conj();
119+
}
120+
}
90121
info.as_lapack_result()?;
91122
Ok(())
92123
}

‎ndarray-linalg/src/eig.rs

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -66,13 +66,13 @@ pub trait EigVals {
6666
impl<A, S> EigVals for ArrayBase<S, Ix2>
6767
where
6868
A: Scalar + Lapack,
69-
S: DataMut<Elem = A>,
69+
S: Data<Elem = A>,
7070
{
7171
type EigVal = Array1<A::Complex>;
7272

7373
fn eigvals(&self) -> Result<Self::EigVal> {
7474
let mut a = self.to_owned();
75-
let (s, _) = A::eig(true, a.square_layout()?, a.as_allocated_mut()?)?;
75+
let (s, _) = A::eig(false, a.square_layout()?, a.as_allocated_mut()?)?;
7676
Ok(ArrayBase::from(s))
7777
}
7878
}

‎ndarray-linalg/src/eigh.rs

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,17 @@ where
144144
{
145145
type EigVal = Array1<A::Real>;
146146

147+
/// Solves the generalized eigenvalue problem.
148+
///
149+
/// # Panics
150+
///
151+
/// Panics if the shapes of the matrices are different.
147152
fn eigh_inplace(&mut self, uplo: UPLO) -> Result<(Self::EigVal, &mut Self)> {
153+
assert_eq!(
154+
self.0.shape(),
155+
self.1.shape(),
156+
"The shapes of the matrices must be identical.",
157+
);
148158
let layout = self.0.square_layout()?;
149159
// XXX Force layout to be Fortran (see #146)
150160
match layout {

‎ndarray-linalg/src/solve.rs

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,59 +77,103 @@ pub use lax::{Pivot, Transpose};
7777
pub trait Solve<A: Scalar> {
7878
/// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
7979
/// is the argument, and `x` is the successful result.
80+
///
81+
/// # Panics
82+
///
83+
/// Panics if the length of `b` is not the equal to the number of columns
84+
/// of `A`.
8085
fn solve<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
8186
let mut b = replicate(b);
8287
self.solve_inplace(&mut b)?;
8388
Ok(b)
8489
}
90+
8591
/// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
8692
/// is the argument, and `x` is the successful result.
93+
///
94+
/// # Panics
95+
///
96+
/// Panics if the length of `b` is not the equal to the number of columns
97+
/// of `A`.
8798
fn solve_into<S: DataMut<Elem = A>>(
8899
&self,
89100
mut b: ArrayBase<S, Ix1>,
90101
) -> Result<ArrayBase<S, Ix1>> {
91102
self.solve_inplace(&mut b)?;
92103
Ok(b)
93104
}
105+
94106
/// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
95107
/// is the argument, and `x` is the successful result.
108+
///
109+
/// # Panics
110+
///
111+
/// Panics if the length of `b` is not the equal to the number of columns
112+
/// of `A`.
96113
fn solve_inplace<'a, S: DataMut<Elem = A>>(
97114
&self,
98115
b: &'a mut ArrayBase<S, Ix1>,
99116
) -> Result<&'a mut ArrayBase<S, Ix1>>;
100117

101118
/// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
102119
/// is the argument, and `x` is the successful result.
120+
///
121+
/// # Panics
122+
///
123+
/// Panics if the length of `b` is not the equal to the number of rows of
124+
/// `A`.
103125
fn solve_t<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
104126
let mut b = replicate(b);
105127
self.solve_t_inplace(&mut b)?;
106128
Ok(b)
107129
}
130+
108131
/// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
109132
/// is the argument, and `x` is the successful result.
133+
///
134+
/// # Panics
135+
///
136+
/// Panics if the length of `b` is not the equal to the number of rows of
137+
/// `A`.
110138
fn solve_t_into<S: DataMut<Elem = A>>(
111139
&self,
112140
mut b: ArrayBase<S, Ix1>,
113141
) -> Result<ArrayBase<S, Ix1>> {
114142
self.solve_t_inplace(&mut b)?;
115143
Ok(b)
116144
}
145+
117146
/// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
118147
/// is the argument, and `x` is the successful result.
148+
///
149+
/// # Panics
150+
///
151+
/// Panics if the length of `b` is not the equal to the number of rows of
152+
/// `A`.
119153
fn solve_t_inplace<'a, S: DataMut<Elem = A>>(
120154
&self,
121155
b: &'a mut ArrayBase<S, Ix1>,
122156
) -> Result<&'a mut ArrayBase<S, Ix1>>;
123157

124158
/// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
125159
/// is the argument, and `x` is the successful result.
160+
///
161+
/// # Panics
162+
///
163+
/// Panics if the length of `b` is not the equal to the number of rows of
164+
/// `A`.
126165
fn solve_h<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
127166
let mut b = replicate(b);
128167
self.solve_h_inplace(&mut b)?;
129168
Ok(b)
130169
}
131170
/// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
132171
/// is the argument, and `x` is the successful result.
172+
///
173+
/// # Panics
174+
///
175+
/// Panics if the length of `b` is not the equal to the number of rows of
176+
/// `A`.
133177
fn solve_h_into<S: DataMut<Elem = A>>(
134178
&self,
135179
mut b: ArrayBase<S, Ix1>,
@@ -139,6 +183,11 @@ pub trait Solve<A: Scalar> {
139183
}
140184
/// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
141185
/// is the argument, and `x` is the successful result.
186+
///
187+
/// # Panics
188+
///
189+
/// Panics if the length of `b` is not the equal to the number of rows of
190+
/// `A`.
142191
fn solve_h_inplace<'a, S: DataMut<Elem = A>>(
143192
&self,
144193
b: &'a mut ArrayBase<S, Ix1>,
@@ -167,6 +216,11 @@ where
167216
where
168217
Sb: DataMut<Elem = A>,
169218
{
219+
assert_eq!(
220+
rhs.len(),
221+
self.a.len_of(Axis(1)),
222+
"The length of `rhs` must be compatible with the shape of the factored matrix.",
223+
);
170224
A::solve(
171225
self.a.square_layout()?,
172226
Transpose::No,
@@ -183,6 +237,11 @@ where
183237
where
184238
Sb: DataMut<Elem = A>,
185239
{
240+
assert_eq!(
241+
rhs.len(),
242+
self.a.len_of(Axis(0)),
243+
"The length of `rhs` must be compatible with the shape of the factored matrix.",
244+
);
186245
A::solve(
187246
self.a.square_layout()?,
188247
Transpose::Transpose,
@@ -199,6 +258,11 @@ where
199258
where
200259
Sb: DataMut<Elem = A>,
201260
{
261+
assert_eq!(
262+
rhs.len(),
263+
self.a.len_of(Axis(0)),
264+
"The length of `rhs` must be compatible with the shape of the factored matrix.",
265+
);
202266
A::solve(
203267
self.a.square_layout()?,
204268
Transpose::Hermite,

‎ndarray-linalg/src/solveh.rs

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -69,25 +69,42 @@ pub trait SolveH<A: Scalar> {
6969
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
7070
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
7171
/// `x` is the successful result.
72+
///
73+
/// # Panics
74+
///
75+
/// Panics if the length of `b` is not the equal to the number of columns
76+
/// of `A`.
7277
fn solveh<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
7378
let mut b = replicate(b);
7479
self.solveh_inplace(&mut b)?;
7580
Ok(b)
7681
}
82+
7783
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
7884
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
7985
/// `x` is the successful result.
86+
///
87+
/// # Panics
88+
///
89+
/// Panics if the length of `b` is not the equal to the number of columns
90+
/// of `A`.
8091
fn solveh_into<S: DataMut<Elem = A>>(
8192
&self,
8293
mut b: ArrayBase<S, Ix1>,
8394
) -> Result<ArrayBase<S, Ix1>> {
8495
self.solveh_inplace(&mut b)?;
8596
Ok(b)
8697
}
98+
8799
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
88100
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
89101
/// `x` is the successful result. The value of `x` is also assigned to the
90102
/// argument.
103+
///
104+
/// # Panics
105+
///
106+
/// Panics if the length of `b` is not the equal to the number of columns
107+
/// of `A`.
91108
fn solveh_inplace<'a, S: DataMut<Elem = A>>(
92109
&self,
93110
b: &'a mut ArrayBase<S, Ix1>,
@@ -113,6 +130,11 @@ where
113130
where
114131
Sb: DataMut<Elem = A>,
115132
{
133+
assert_eq!(
134+
rhs.len(),
135+
self.a.len_of(Axis(1)),
136+
"The length of `rhs` must be compatible with the shape of the factored matrix.",
137+
);
116138
A::solveh(
117139
self.a.square_layout()?,
118140
UPLO::Upper,

0 commit comments

Comments
 (0)