nalgebra-lapack: unify API of LU.solve and Cholesky.solve with nalgebra.

This commit is contained in:
Sébastien Crozet 2017-08-13 19:52:51 +02:00 committed by Sébastien Crozet
parent 00039c0a76
commit 053de0576f
4 changed files with 90 additions and 22 deletions

View File

@ -1,8 +1,9 @@
use num::Zero; use num::Zero;
use num_complex::Complex; use num_complex::Complex;
use na::{Scalar, DefaultAllocator, MatrixN, MatrixMN}; use na::{Scalar, DefaultAllocator, Matrix, MatrixN, MatrixMN};
use na::dimension::Dim; use na::dimension::Dim;
use na::storage::Storage;
use na::allocator::Allocator; use na::allocator::Allocator;
use lapack::fortran as interface; use lapack::fortran as interface;
@ -48,8 +49,22 @@ impl<N: CholeskyScalar + Zero, D: Dim> Cholesky<N, D>
/// Solves the symmetric-definite-positive linear system `self * x = b`, where `x` is the /// Solves the symmetric-definite-positive linear system `self * x = b`, where `x` is the
/// unknown to be determined. /// unknown to be determined.
pub fn solve<R2: Dim, C2: Dim>(&self, mut b: MatrixMN<N, R2, C2>) pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>>
-> Option<MatrixMN<N, R2, C2>> where S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2> {
let mut res = b.clone_owned();
if self.solve_mut(&mut res) {
Some(res)
}
else {
None
}
}
/// Solves in-place the symmetric-definite-positive linear system `self * x = b`, where `x` is
/// the unknown to be determined.
pub fn solve_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
where DefaultAllocator: Allocator<N, R2, C2> { where DefaultAllocator: Allocator<N, R2, C2> {
let dim = self.l.nrows(); let dim = self.l.nrows();
@ -62,9 +77,7 @@ impl<N: CholeskyScalar + Zero, D: Dim> Cholesky<N, D>
let mut info = 0; let mut info = 0;
N::xpotrs(b'L', dim as i32, nrhs, self.l.as_slice(), lda, b.as_mut_slice(), ldb, &mut info); N::xpotrs(b'L', dim as i32, nrhs, self.l.as_slice(), lda, b.as_mut_slice(), ldb, &mut info);
lapack_check!(info); lapack_test!(info)
Some(b)
} }
/// Computes the inverse of the decomposed matrix. /// Computes the inverse of the decomposed matrix.

View File

@ -19,3 +19,9 @@ macro_rules! lapack_panic(
assert!($info == 0, "Lapack error."); assert!($info == 0, "Lapack error.");
); );
); );
macro_rules! lapack_test(
($info: expr) => (
$info == 0
);
);

View File

@ -105,8 +105,7 @@ impl<N: LUScalar, R: Dim, C: Dim> LU<N, R, C>
1, self.p.len() as i32, self.p.as_slice(), -1); 1, self.p.len() as i32, self.p.as_slice(), -1);
} }
fn generic_solve<R2: Dim, C2: Dim>(&self, trans: u8, mut b: MatrixMN<N, R2, C2>) fn generic_solve_mut<R2: Dim, C2: Dim>(&self, trans: u8, b: &mut MatrixMN<N, R2, C2>) -> bool
-> Option<MatrixMN<N, R2, C2>>
where DefaultAllocator: Allocator<N, R2, C2> + where DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> { Allocator<i32, R2> {
@ -121,39 +120,89 @@ impl<N: LUScalar, R: Dim, C: Dim> LU<N, R, C>
let mut info = 0; let mut info = 0;
N::xgetrs(trans, dim as i32, nrhs, self.lu.as_slice(), lda, self.p.as_slice(), N::xgetrs(trans, dim as i32, nrhs, self.lu.as_slice(), lda, self.p.as_slice(),
b.as_mut_slice(), ldb, &mut info); b.as_mut_slice(), ldb, &mut info);
lapack_check!(info); lapack_test!(info)
Some(b)
} }
/// Solves the linear system `self * x = b`, where `x` is the unknown to be determined. /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
pub fn solve<R2: Dim, C2: Dim>(&self, b: MatrixMN<N, R2, C2>) pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>>
-> Option<MatrixMN<N, R2, C2>> where S2: Storage<N, R2, C2>,
where DefaultAllocator: Allocator<N, R2, C2> + DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> { Allocator<i32, R2> {
self.generic_solve(b'N', b) let mut res = b.clone_owned();
if self.generic_solve_mut(b'N', &mut res) {
Some(res)
}
else {
None
}
} }
/// Solves the linear system `self.transpose() * x = b`, where `x` is the unknown to be /// Solves the linear system `self.transpose() * x = b`, where `x` is the unknown to be
/// determined. /// determined.
pub fn solve_transpose<R2: Dim, C2: Dim>(&self, b: MatrixMN<N, R2, C2>) pub fn solve_transpose<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>)
-> Option<MatrixMN<N, R2, C2>> -> Option<MatrixMN<N, R2, C2>>
where DefaultAllocator: Allocator<N, R2, C2> + where S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> { Allocator<i32, R2> {
self.generic_solve(b'T', b) let mut res = b.clone_owned();
if self.generic_solve_mut(b'T', &mut res) {
Some(res)
}
else {
None
}
} }
/// Solves the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to /// Solves the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to
/// be determined. /// be determined.
pub fn solve_conjugate_transpose<R2: Dim, C2: Dim>(&self, b: MatrixMN<N, R2, C2>) pub fn solve_conjugate_transpose<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>)
-> Option<MatrixMN<N, R2, C2>> -> Option<MatrixMN<N, R2, C2>>
where S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> {
let mut res = b.clone_owned();
if self.generic_solve_mut(b'T', &mut res) {
Some(res)
}
else {
None
}
}
/// Solves in-place the linear system `self * x = b`, where `x` is the unknown to be determined.
///
/// Retuns `false` if no solution was found (the decomposed matrix is singular).
pub fn solve_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
where DefaultAllocator: Allocator<N, R2, C2> + where DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> { Allocator<i32, R2> {
self.generic_solve(b'T', b) self.generic_solve_mut(b'N', b)
}
/// Solves in-place the linear system `self.transpose() * x = b`, where `x` is the unknown to be
/// determined.
///
/// Retuns `false` if no solution was found (the decomposed matrix is singular).
pub fn solve_transpose_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
where DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> {
self.generic_solve_mut(b'T', b)
}
/// Solves in-place the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to
/// be determined.
///
/// Retuns `false` if no solution was found (the decomposed matrix is singular).
pub fn solve_conjugate_transpose_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
where DefaultAllocator: Allocator<N, R2, C2> +
Allocator<i32, R2> {
self.generic_solve_mut(b'T', b)
} }
} }

View File

@ -169,7 +169,7 @@ impl<N: Real, D: DimMin<D, Output = D>> LU<N, D, D>
/// ///
/// Returns `None` if `self` is not invertible. /// Returns `None` if `self` is not invertible.
pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>> pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>>
where S2: StorageMut<N, R2, C2>, where S2: Storage<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>, ShapeConstraint: SameNumberOfRows<R2, D>,
DefaultAllocator: Allocator<N, R2, C2> { DefaultAllocator: Allocator<N, R2, C2> {
let mut res = b.clone_owned(); let mut res = b.clone_owned();