Skip to content

Add documentation for solveh module #76

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Sep 21, 2017
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
111 changes: 103 additions & 8 deletions src/solveh.rs
Original file line number Diff line number Diff line change
@@ -1,4 +1,53 @@
//! Solve Hermite/Symmetric linear problems
//! Solve Hermitian (or real symmetric) linear problems and invert Hermitian
//! (or real symmetric) matrices
//!
//! **Note that only the upper triangular portion of the matrix is used.**
//!
//! # Examples
//!
//! Solve `A * x = b`, where `A` is a Hermitian (or real symmetric) matrix:
//!
//! ```
//! #[macro_use]
//! extern crate ndarray;
//! extern crate ndarray_linalg;
//!
//! use ndarray::prelude::*;
//! use ndarray_linalg::SolveH;
//! # fn main() {
//!
//! let a: Array2<f64> = array![
//! [3., 2., -1.],
//! [2., -2., 4.],
//! [-1., 4., 5.]
//! ];
//! let b: Array1<f64> = array![11., -12., 1.];
//! let x = a.solveh_into(b).unwrap();
//! assert!(x.all_close(&array![1., 3., -2.], 1e-9));
//!
//! # }
//! ```
//!
//! If you are solving multiple systems of linear equations with the same
//! Hermitian or real symmetric coefficient matrix `A`, it's faster to compute
//! the factorization once at the beginning than solving directly using `A`:
//!
//! ```
//! # extern crate ndarray;
//! # extern crate ndarray_linalg;
//! use ndarray::prelude::*;
//! use ndarray_linalg::*;
//! # fn main() {
//!
//! let a: Array2<f64> = random((3, 3));
//! let f = a.factorizeh_into().unwrap(); // Factorize A (A is consumed)
//! for _ in 0..10 {
//! let b: Array1<f64> = random(3);
//! let x = f.solveh_into(b).unwrap(); // Solve A * x = b using the factorization
//! }
//!
//! # }
//! ```

use ndarray::*;

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

pub use lapack_traits::{Pivot, UPLO};

/// An interface for solving systems of Hermitian (or real symmetric) linear equations.
///
/// If you plan to solve many equations with the same Hermitian (or real
/// symmetric) coefficient matrix `A` but different `b` vectors, it's faster to
/// factor the `A` matrix once using the `FactorizeH` trait, and then solve
/// using the `FactorizedH` struct.
pub trait SolveH<A: Scalar> {
fn solveh<S: Data<Elem = A>>(&self, a: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
let mut a = replicate(a);
self.solveh_mut(&mut a)?;
Ok(a)
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
/// `x` is the successful result.
fn solveh<S: Data<Elem = A>>(&self, b: &ArrayBase<S, Ix1>) -> Result<Array1<A>> {
let mut b = replicate(b);
self.solveh_mut(&mut b)?;
Ok(b)
}
fn solveh_into<S: DataMut<Elem = A>>(&self, mut a: ArrayBase<S, Ix1>) -> Result<ArrayBase<S, Ix1>> {
self.solveh_mut(&mut a)?;
Ok(a)
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
/// `x` is the successful result.
fn solveh_into<S: DataMut<Elem = A>>(&self, mut b: ArrayBase<S, Ix1>) -> Result<ArrayBase<S, Ix1>> {
self.solveh_mut(&mut b)?;
Ok(b)
}
/// Solves a system of linear equations `A * x = b` with Hermitian (or real
/// symmetric) matrix `A`, where `A` is `self`, `b` is the argument, and
/// `x` is the successful result. The value of `x` is also assigned to the
/// argument.
fn solveh_mut<'a, S: DataMut<Elem = A>>(&self, &'a mut ArrayBase<S, Ix1>) -> Result<&'a mut ArrayBase<S, Ix1>>;
}

/// Represents the Bunch–Kaufman factorization of a Hermitian (or real
/// symmetric) matrix as `A = P * U * D * U^T * P^T`.
pub struct FactorizedH<S: Data> {
pub a: ArrayBase<S, Ix2>,
pub ipiv: Pivot,
Expand Down Expand Up @@ -69,6 +136,12 @@ where
A: Scalar,
S: DataMut<Elem = A>,
{
/// Computes the inverse of the factorized matrix.
///
/// **Warning: The inverse is stored only in the upper triangular portion
/// of the result matrix!** If you want the lower triangular portion to be
/// correct, you must fill it in according to the results in the upper
/// triangular portion.
pub fn into_inverseh(mut self) -> Result<ArrayBase<S, Ix2>> {
unsafe {
A::invh(
Expand All @@ -82,11 +155,19 @@ where
}
}

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

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

Expand Down Expand Up @@ -116,13 +197,27 @@ where
}
}

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

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

Expand Down