Merge remote-tracking branch 'origin/master'

Author: termoshtt

src/cholesky.rs

@@ -1,6 +1,36 @@
-//! Cholesky decomposition
+//! Cholesky decomposition of Hermitian (or real symmetric) positive definite matrices
 //!
-//! https://en.wikipedia.org/wiki/Cholesky_decomposition
+//! See the [Wikipedia page about Cholesky
+//! decomposition](https://en.wikipedia.org/wiki/Cholesky_decomposition) for
+//! more information.
+//!
+//! # Example
+//!
+//! Calculate `L` in the Cholesky decomposition `A = L * L^H`, where `A` is a
+//! Hermitian (or real symmetric) positive definite matrix:
+//!
+//! ```
+//! #[macro_use]
+//! extern crate ndarray;
+//! extern crate ndarray_linalg;
+//!
+//! use ndarray::prelude::*;
+//! use ndarray_linalg::{CholeskyInto, UPLO};
+//! # fn main() {
+//!
+//! let a: Array2<f64> = array![
+//!     [  4.,  12., -16.],
+//!     [ 12.,  37., -43.],
+//!     [-16., -43.,  98.]
+//! ];
+//! let lower = a.cholesky_into(UPLO::Lower).unwrap();
+//! assert!(lower.all_close(&array![
+//!     [ 2., 0., 0.],
+//!     [ 6., 1., 0.],
+//!     [-8., 5., 3.]
+//! ], 1e-9));
+//! # }
+//! ```
 
 use ndarray::*;
 
@@ -12,19 +42,44 @@ use super::types::*;
 
 pub use lapack_traits::UPLO;
 
-/// Cholesky decomposition of matrix reference
+/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix reference
 pub trait Cholesky {
     type Output;
+    /// Computes the Cholesky decomposition of the Hermitian (or real
+    /// symmetric) positive definite matrix.
+    ///
+    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
+    /// U^H * U` using the upper triangular portion of `A` and returns `U`.
+    /// Otherwise, if the argument is `UPLO::Lower`, computes the decomposition
+    /// `A = L * L^H` using the lower triangular portion of `A` and returns
+    /// `L`.
     fn cholesky(&self, UPLO) -> Result<Self::Output>;
 }
 
-/// Cholesky decomposition
+/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix
 pub trait CholeskyInto: Sized {
+    /// Computes the Cholesky decomposition of the Hermitian (or real
+    /// symmetric) positive definite matrix.
+    ///
+    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
+    /// U^H * U` using the upper triangular portion of `A` and returns `U`.
+    /// Otherwise, if the argument is `UPLO::Lower`, computes the decomposition
+    /// `A = L * L^H` using the lower triangular portion of `A` and returns
+    /// `L`.
     fn cholesky_into(self, UPLO) -> Result<Self>;
 }
 
-/// Cholesky decomposition of mutable reference of matrix
+/// Cholesky decomposition of Hermitian (or real symmetric) positive definite mutable reference of matrix
 pub trait CholeskyMut {
+    /// Computes the Cholesky decomposition of the Hermitian (or real
+    /// symmetric) positive definite matrix, storing the result in `self` and
+    /// returning it.
+    ///
+    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
+    /// U^H * U` using the upper triangular portion of `A` and returns `U`.
+    /// Otherwise, if the argument is `UPLO::Lower`, computes the decomposition
+    /// `A = L * L^H` using the lower triangular portion of `A` and returns
+    /// `L`.
     fn cholesky_mut(&mut self, UPLO) -> Result<&mut Self>;
 }
 

src/solveh.rs

@@ -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::*;
 
@@ -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^H * P^T`.
 pub struct FactorizedH<S: Data> {
     pub a: ArrayBase<S, Ix2>,
     pub ipiv: Pivot,
@@ -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(
@@ -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>>;
 }
 
@@ -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>;
 }