Skip to content

Commit d0d9186

Browse files
committed
Add documentation for solveh module
1 parent b74ad3d commit d0d9186

File tree

1 file changed

+103
-8
lines changed

1 file changed

+103
-8
lines changed

src/solveh.rs

Lines changed: 103 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,53 @@
1-
//! Solve Hermite/Symmetric linear problems
1+
//! Solve Hermitian (or real symmetric) linear problems and invert Hermitian
2+
//! (or real symmetric) matrices
3+
//!
4+
//! **Note that only the upper triangular portion of the matrix is used.**
5+
//!
6+
//! # Examples
7+
//!
8+
//! Solve `A * x = b`, where `A` is a Hermitian (or real symmetric) matrix:
9+
//!
10+
//! ```
11+
//! #[macro_use]
12+
//! extern crate ndarray;
13+
//! extern crate ndarray_linalg;
14+
//!
15+
//! use ndarray::prelude::*;
16+
//! use ndarray_linalg::SolveH;
17+
//! # fn main() {
18+
//!
19+
//! let a: Array2<f64> = array![
20+
//! [3., 2., -1.],
21+
//! [2., -2., 4.],
22+
//! [-1., 4., 5.]
23+
//! ];
24+
//! let b: Array1<f64> = array![11., -12., 1.];
25+
//! let x = a.solveh_into(b).unwrap();
26+
//! assert!(x.all_close(&array![1., 3., -2.], 1e-9));
27+
//!
28+
//! # }
29+
//! ```
30+
//!
31+
//! If you are solving multiple systems of linear equations with the same
32+
//! Hermitian or real symmetric coefficient matrix `A`, it's faster to compute
33+
//! the factorization once at the beginning than solving directly using `A`:
34+
//!
35+
//! ```
36+
//! # extern crate ndarray;
37+
//! # extern crate ndarray_linalg;
38+
//! use ndarray::prelude::*;
39+
//! use ndarray_linalg::*;
40+
//! # fn main() {
41+
//!
42+
//! let a: Array2<f64> = random((3, 3));
43+
//! let f = a.factorizeh_into().unwrap(); // Factorize A (A is consumed)
44+
//! for _ in 0..10 {
45+
//! let b: Array1<f64> = random(3);
46+
//! let x = f.solveh_into(b).unwrap(); // Solve A * x = b using the factorization
47+
//! }
48+
//!
49+
//! # }
50+
//! ```
251
352
use ndarray::*;
453

@@ -9,19 +58,37 @@ use super::types::*;
958

1059
pub use lapack_traits::{Pivot, UPLO};
1160

61+
/// An interface for solving systems of Hermitian (or real symmetric) linear equations.
62+
///
63+
/// If you plan to solve many equations with the same Hermitian (or real
64+
/// symmetric) coefficient matrix `A` but different `b` vectors, it's faster to
65+
/// factor the `A` matrix once using the `FactorizeH` trait, and then solve
66+
/// using the `FactorizedH` struct.
1267
pub trait SolveH<A: Scalar> {
13-
fn solveh<S: Data<Elem = A>>(&self, a: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
14-
let mut a = replicate(a);
15-
self.solveh_mut(&mut a)?;
16-
Ok(a)
68+
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
69+
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
70+
/// `x` is the successful result.
71+
fn solveh<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
72+
let mut b = replicate(b);
73+
self.solveh_mut(&mut b)?;
74+
Ok(b)
1775
}
18-
fn solveh_into<S: DataMut<Elem = A>>(&self, mut a: ArrayBase<S, Ix1>) -> Result<ArrayBase<S, Ix1>> {
19-
self.solveh_mut(&mut a)?;
20-
Ok(a)
76+
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
77+
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
78+
/// `x` is the successful result.
79+
fn solveh_into<S: DataMut<Elem = A>>(&self, mut b: ArrayBase<S, Ix1>) -> Result<ArrayBase<S, Ix1>> {
80+
self.solveh_mut(&mut b)?;
81+
Ok(b)
2182
}
83+
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
84+
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
85+
/// `x` is the successful result. The value of `x` is also assigned to the
86+
/// argument.
2287
fn solveh_mut<'a, S: DataMut<Elem = A>>(&self, &'a mut ArrayBase<S, Ix1>) -> Result<&'a mut ArrayBase<S, Ix1>>;
2388
}
2489

90+
/// Represents the Bunch–Kaufman factorization of a Hermitian (or real
91+
/// symmetric) matrix as `A = P * U * D * U^T * P^T`.
2592
pub struct FactorizedH<S: Data> {
2693
pub a: ArrayBase<S, Ix2>,
2794
pub ipiv: Pivot,
@@ -69,6 +136,12 @@ where
69136
A: Scalar,
70137
S: DataMut<Elem = A>,
71138
{
139+
/// Computes the inverse of the factorized matrix.
140+
///
141+
/// **Warning: The inverse is stored only in the upper triangular portion
142+
/// of the result matrix!** If you want the lower triangular portion to be
143+
/// correct, you must fill it in according to the results in the upper
144+
/// triangular portion.
72145
pub fn into_inverseh(mut self) -> Result<ArrayBase<S, Ix2>> {
73146
unsafe {
74147
A::invh(
@@ -82,11 +155,19 @@ where
82155
}
83156
}
84157

158+
/// An interface for computing the Bunch–Kaufman factorization of Hermitian (or
159+
/// real symmetric) matrix refs.
85160
pub trait FactorizeH<S: Data> {
161+
/// Computes the Bunch–Kaufman factorization of a Hermitian (or real
162+
/// symmetric) matrix.
86163
fn factorizeh(&self) -> Result<FactorizedH<S>>;
87164
}
88165

166+
/// An interface for computing the Bunch–Kaufman factorization of Hermitian (or
167+
/// real symmetric) matrices.
89168
pub trait FactorizeHInto<S: Data> {
169+
/// Computes the Bunch–Kaufman factorization of a Hermitian (or real
170+
/// symmetric) matrix.
90171
fn factorizeh_into(self) -> Result<FactorizedH<S>>;
91172
}
92173

@@ -116,13 +197,27 @@ where
116197
}
117198
}
118199

200+
/// An interface for inverting Hermitian (or real symmetric) matrix refs.
119201
pub trait InverseH {
120202
type Output;
203+
/// Computes the inverse of the Hermitian (or real symmetric) matrix.
204+
///
205+
/// **Warning: The inverse is stored only in the upper triangular portion
206+
/// of the result matrix!** If you want the lower triangular portion to be
207+
/// correct, you must fill it in according to the results in the upper
208+
/// triangular portion.
121209
fn invh(&self) -> Result<Self::Output>;
122210
}
123211

212+
/// An interface for inverting Hermitian (or real symmetric) matrices.
124213
pub trait InverseHInto {
125214
type Output;
215+
/// Computes the inverse of the Hermitian (or real symmetric) matrix.
216+
///
217+
/// **Warning: The inverse is stored only in the upper triangular portion
218+
/// of the result matrix!** If you want the lower triangular portion to be
219+
/// correct, you must fill it in according to the results in the upper
220+
/// triangular portion.
126221
fn invh_into(self) -> Result<Self::Output>;
127222
}
128223

0 commit comments

Comments
 (0)