1- //! Solve linear problems
1+ //! Solve systems of linear equations and invert matrices
2+ //!
3+ //! # Examples
4+ //!
5+ //! Solve `A * x = b`:
6+ //!
7+ //! ```
8+ //! #[macro_use]
9+ //! extern crate ndarray;
10+ //! extern crate ndarray_linalg;
11+ //!
12+ //! use ndarray::prelude::*;
13+ //! use ndarray_linalg::Solve;
14+ //! # fn main() {
15+ //!
16+ //! let a: Array2<f64> = array![[3., 2., -1.], [2., -2., 4.], [-2., 1., -2.]];
17+ //! let b: Array1<f64> = array![1., -2., 0.];
18+ //! let x = a.solve_into(b).unwrap();
19+ //! assert!(x.all_close(&array![1., -2., -2.], 1e-9));
20+ //!
21+ //! # }
22+ //! ```
23+ //!
24+ //! There are also special functions for solving `A^T * x = b` and
25+ //! `A^H * x = b`.
26+ //!
27+ //! If you are solving multiple systems of linear equations with the same
28+ //! coefficient matrix `A`, it's faster to compute the LU factorization once at
29+ //! the beginning than solving directly using `A`:
30+ //!
31+ //! ```
32+ //! # extern crate ndarray;
33+ //! # extern crate ndarray_linalg;
34+ //! use ndarray::prelude::*;
35+ //! use ndarray_linalg::*;
36+ //! # fn main() {
37+ //!
38+ //! let a: Array2<f64> = random((3, 3));
39+ //! let f = a.factorize_into().unwrap(); // LU factorize A (A is consumed)
40+ //! for _ in 0..10 {
41+ //! let b: Array1<f64> = random(3);
42+ //! let x = f.solve_into(b).unwrap(); // Solve A * x = b using factorized L, U
43+ //! }
44+ //!
45+ //! # }
46+ //! ```
247
348use ndarray:: * ;
449
@@ -9,43 +54,82 @@ use super::types::*;
954
1055pub use lapack_traits:: { Pivot , Transpose } ;
1156
57+ /// An interface for solving systems of linear equations.
58+ ///
59+ /// There are three groups of methods:
60+ ///
61+ /// * `solve*` (normal) methods solve `A * x = b` for `x`.
62+ /// * `solve_t*` (transpose) methods solve `A^T * x = b` for `x`.
63+ /// * `solve_h*` (Hermitian conjugate) methods solve `A^H * x = b` for `x`.
64+ ///
65+ /// Within each group, there are three methods that handle ownership differently:
66+ ///
67+ /// * `*` methods take a reference to `b` and return `x` as a new array.
68+ /// * `*_into` methods take ownership of `b`, store the result in it, and return it.
69+ /// * `*_mut` methods take a mutable reference to `b` and store the result in that array.
70+ ///
71+ /// If you plan to solve many equations with the same `A` matrix but different
72+ /// `b` vectors, it's faster to factor the `A` matrix once using the
73+ /// `Factorize` trait, and then solve using the `Factorized` struct.
1274pub trait Solve < A : Scalar > {
13- fn solve < S : Data < Elem = A > > ( & self , a : & ArrayBase < S , Ix1 > ) -> Result < Array1 < A > > {
14- let mut a = replicate ( a) ;
15- self . solve_mut ( & mut a) ?;
16- Ok ( a)
75+ /// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
76+ /// is the argument, and `x` is the successful result.
77+ fn solve < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix1 > ) -> Result < Array1 < A > > {
78+ let mut b = replicate ( b) ;
79+ self . solve_mut ( & mut b) ?;
80+ Ok ( b)
1781 }
18- fn solve_into < S : DataMut < Elem = A > > ( & self , mut a : ArrayBase < S , Ix1 > ) -> Result < ArrayBase < S , Ix1 > > {
19- self . solve_mut ( & mut a) ?;
20- Ok ( a)
82+ /// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
83+ /// is the argument, and `x` is the successful result.
84+ fn solve_into < S : DataMut < Elem = A > > ( & self , mut b : ArrayBase < S , Ix1 > ) -> Result < ArrayBase < S , Ix1 > > {
85+ self . solve_mut ( & mut b) ?;
86+ Ok ( b)
2187 }
88+ /// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
89+ /// is the argument, and `x` is the successful result.
2290 fn solve_mut < ' a , S : DataMut < Elem = A > > ( & self , & ' a mut ArrayBase < S , Ix1 > ) -> Result < & ' a mut ArrayBase < S , Ix1 > > ;
2391
24- fn solve_t < S : Data < Elem = A > > ( & self , a : & ArrayBase < S , Ix1 > ) -> Result < Array1 < A > > {
25- let mut a = replicate ( a) ;
26- self . solve_t_mut ( & mut a) ?;
27- Ok ( a)
92+ /// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
93+ /// is the argument, and `x` is the successful result.
94+ fn solve_t < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix1 > ) -> Result < Array1 < A > > {
95+ let mut b = replicate ( b) ;
96+ self . solve_t_mut ( & mut b) ?;
97+ Ok ( b)
2898 }
29- fn solve_t_into < S : DataMut < Elem = A > > ( & self , mut a : ArrayBase < S , Ix1 > ) -> Result < ArrayBase < S , Ix1 > > {
30- self . solve_t_mut ( & mut a) ?;
31- Ok ( a)
99+ /// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
100+ /// is the argument, and `x` is the successful result.
101+ fn solve_t_into < S : DataMut < Elem = A > > ( & self , mut b : ArrayBase < S , Ix1 > ) -> Result < ArrayBase < S , Ix1 > > {
102+ self . solve_t_mut ( & mut b) ?;
103+ Ok ( b)
32104 }
105+ /// Solves a system of linear equations `A^T * x = b` where `A` is `self`, `b`
106+ /// is the argument, and `x` is the successful result.
33107 fn solve_t_mut < ' a , S : DataMut < Elem = A > > ( & self , & ' a mut ArrayBase < S , Ix1 > ) -> Result < & ' a mut ArrayBase < S , Ix1 > > ;
34108
35- fn solve_h < S : Data < Elem = A > > ( & self , a : & ArrayBase < S , Ix1 > ) -> Result < Array1 < A > > {
36- let mut a = replicate ( a) ;
37- self . solve_h_mut ( & mut a) ?;
38- Ok ( a)
109+ /// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
110+ /// is the argument, and `x` is the successful result.
111+ fn solve_h < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix1 > ) -> Result < Array1 < A > > {
112+ let mut b = replicate ( b) ;
113+ self . solve_h_mut ( & mut b) ?;
114+ Ok ( b)
39115 }
40- fn solve_h_into < S : DataMut < Elem = A > > ( & self , mut a : ArrayBase < S , Ix1 > ) -> Result < ArrayBase < S , Ix1 > > {
41- self . solve_h_mut ( & mut a) ?;
42- Ok ( a)
116+ /// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
117+ /// is the argument, and `x` is the successful result.
118+ fn solve_h_into < S : DataMut < Elem = A > > ( & self , mut b : ArrayBase < S , Ix1 > ) -> Result < ArrayBase < S , Ix1 > > {
119+ self . solve_h_mut ( & mut b) ?;
120+ Ok ( b)
43121 }
122+ /// Solves a system of linear equations `A^H * x = b` where `A` is `self`, `b`
123+ /// is the argument, and `x` is the successful result.
44124 fn solve_h_mut < ' a , S : DataMut < Elem = A > > ( & self , & ' a mut ArrayBase < S , Ix1 > ) -> Result < & ' a mut ArrayBase < S , Ix1 > > ;
45125}
46126
127+ /// Represents the LU factorization of a matrix `A` as `A = P*L*U`.
47128pub struct Factorized < S : Data > {
129+ /// The factors `L` and `U`; the unit diagonal elements of `L` are not
130+ /// stored.
48131 pub a : ArrayBase < S , Ix2 > ,
132+ /// The pivot indices that define the permutation matrix `P`.
49133 pub ipiv : Pivot ,
50134}
51135
@@ -134,6 +218,7 @@ where
134218 A : Scalar ,
135219 S : DataMut < Elem = A > ,
136220{
221+ /// Computes the inverse of the factorized matrix.
137222 pub fn into_inverse ( mut self ) -> Result < ArrayBase < S , Ix2 > > {
138223 unsafe {
139224 A :: inv (
@@ -146,11 +231,17 @@ where
146231 }
147232}
148233
234+ /// An interface for computing LU factorizations of matrix refs.
149235pub trait Factorize < S : Data > {
236+ /// Computes the LU factorization `A = P*L*U`, where `P` is a permutation
237+ /// matrix.
150238 fn factorize ( & self ) -> Result < Factorized < S > > ;
151239}
152240
241+ /// An interface for computing LU factorizations of matrices.
153242pub trait FactorizeInto < S : Data > {
243+ /// Computes the LU factorization `A = P*L*U`, where `P` is a permutation
244+ /// matrix.
154245 fn factorize_into ( self ) -> Result < Factorized < S > > ;
155246}
156247
@@ -180,13 +271,17 @@ where
180271 }
181272}
182273
274+ /// An interface for inverting matrix refs.
183275pub trait Inverse {
184276 type Output ;
277+ /// Computes the inverse of the matrix.
185278 fn inv ( & self ) -> Result < Self :: Output > ;
186279}
187280
281+ /// An interface for inverting matrices.
188282pub trait InverseInto {
189283 type Output ;
284+ /// Computes the inverse of the matrix.
190285 fn inv_into ( self ) -> Result < Self :: Output > ;
191286}
192287
0 commit comments