Skip to content

Refactoring Factorized{H} #81

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 4 commits into from
Sep 23, 2017
Merged
Show file tree
Hide file tree
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
74 changes: 46 additions & 28 deletions src/solve.rs
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ pub use lapack_traits::{Pivot, Transpose};
///
/// If you plan to solve many equations with the same `A` matrix but different
/// `b` vectors, it's faster to factor the `A` matrix once using the
/// `Factorize` trait, and then solve using the `Factorized` struct.
/// `Factorize` trait, and then solve using the `LUFactorized` struct.
pub trait Solve<A: Scalar> {
/// Solves a system of linear equations `A * x = b` where `A` is `self`, `b`
/// is the argument, and `x` is the successful result.
Expand Down Expand Up @@ -125,15 +125,15 @@ pub trait Solve<A: Scalar> {
}

/// Represents the LU factorization of a matrix `A` as `A = P*L*U`.
pub struct Factorized<S: Data> {
pub struct LUFactorized<S: Data> {
/// The factors `L` and `U`; the unit diagonal elements of `L` are not
/// stored.
pub a: ArrayBase<S, Ix2>,
/// The pivot indices that define the permutation matrix `P`.
pub ipiv: Pivot,
}

impl<A, S> Solve<A> for Factorized<S>
impl<A, S> Solve<A> for LUFactorized<S>
where
A: Scalar,
S: Data<Elem = A>,
Expand Down Expand Up @@ -213,46 +213,29 @@ where
}
}

impl<A, S> Factorized<S>
where
A: Scalar,
S: DataMut<Elem = A>,
{
/// Computes the inverse of the factorized matrix.
pub fn into_inverse(mut self) -> Result<ArrayBase<S, Ix2>> {
unsafe {
A::inv(
self.a.square_layout()?,
self.a.as_allocated_mut()?,
&self.ipiv,
)?
};
Ok(self.a)
}
}

/// An interface for computing LU factorizations of matrix refs.
pub trait Factorize<S: Data> {
/// Computes the LU factorization `A = P*L*U`, where `P` is a permutation
/// matrix.
fn factorize(&self) -> Result<Factorized<S>>;
fn factorize(&self) -> Result<LUFactorized<S>>;
}

/// An interface for computing LU factorizations of matrices.
pub trait FactorizeInto<S: Data> {
/// Computes the LU factorization `A = P*L*U`, where `P` is a permutation
/// matrix.
fn factorize_into(self) -> Result<Factorized<S>>;
fn factorize_into(self) -> Result<LUFactorized<S>>;
}

impl<A, S> FactorizeInto<S> for ArrayBase<S, Ix2>
where
A: Scalar,
S: DataMut<Elem = A>,
{
fn factorize_into(mut self) -> Result<Factorized<S>> {
fn factorize_into(mut self) -> Result<LUFactorized<S>> {
let ipiv = unsafe { A::lu(self.layout()?, self.as_allocated_mut()?)? };
Ok(Factorized {
Ok(LUFactorized {
a: self,
ipiv: ipiv,
})
Expand All @@ -264,10 +247,10 @@ where
A: Scalar,
Si: Data<Elem = A>,
{
fn factorize(&self) -> Result<Factorized<OwnedRepr<A>>> {
fn factorize(&self) -> Result<LUFactorized<OwnedRepr<A>>> {
let mut a: Array2<A> = replicate(self);
let ipiv = unsafe { A::lu(a.layout()?, a.as_allocated_mut()?)? };
Ok(Factorized { a: a, ipiv: ipiv })
Ok(LUFactorized { a: a, ipiv: ipiv })
}
}

Expand All @@ -285,6 +268,41 @@ pub trait InverseInto {
fn inv_into(self) -> Result<Self::Output>;
}

impl<A, S> InverseInto for LUFactorized<S>
where
A: Scalar,
S: DataMut<Elem = A>,
{
type Output = ArrayBase<S, Ix2>;

fn inv_into(mut self) -> Result<ArrayBase<S, Ix2>> {
unsafe {
A::inv(
self.a.square_layout()?,
self.a.as_allocated_mut()?,
&self.ipiv,
)?
};
Ok(self.a)
}
}

impl<A, S> Inverse for LUFactorized<S>
where
A: Scalar,
S: Data<Elem = A>,
{
type Output = Array2<A>;

fn inv(&self) -> Result<Array2<A>> {
let f = LUFactorized {
a: replicate(&self.a),
ipiv: self.ipiv.clone(),
};
f.inv_into()
}
}

impl<A, S> InverseInto for ArrayBase<S, Ix2>
where
A: Scalar,
Expand All @@ -294,7 +312,7 @@ where

fn inv_into(self) -> Result<Self::Output> {
let f = self.factorize_into()?;
f.into_inverse()
f.inv_into()
}
}

Expand All @@ -307,6 +325,6 @@ where

fn inv(&self) -> Result<Self::Output> {
let f = self.factorize()?;
f.into_inverse()
f.inv_into()
}
}
82 changes: 47 additions & 35 deletions src/solveh.rs
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ pub use lapack_traits::{Pivot, UPLO};
/// 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.
/// using the `BKFactorized` struct.
pub trait SolveH<A: Scalar> {
/// 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
Expand All @@ -89,12 +89,12 @@ pub trait SolveH<A: Scalar> {

/// 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 struct BKFactorized<S: Data> {
pub a: ArrayBase<S, Ix2>,
pub ipiv: Pivot,
}

impl<A, S> SolveH<A> for FactorizedH<S>
impl<A, S> SolveH<A> for BKFactorized<S>
where
A: Scalar,
S: Data<Elem = A>,
Expand Down Expand Up @@ -131,54 +131,30 @@ where
}


impl<A, S> FactorizedH<S>
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(
self.a.square_layout()?,
UPLO::Upper,
self.a.as_allocated_mut()?,
&self.ipiv,
)?
};
Ok(self.a)
}
}

/// 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>>;
fn factorizeh(&self) -> Result<BKFactorized<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>>;
fn factorizeh_into(self) -> Result<BKFactorized<S>>;
}

impl<A, S> FactorizeHInto<S> for ArrayBase<S, Ix2>
where
A: Scalar,
S: DataMut<Elem = A>,
{
fn factorizeh_into(mut self) -> Result<FactorizedH<S>> {
fn factorizeh_into(mut self) -> Result<BKFactorized<S>> {
let ipiv = unsafe { A::bk(self.layout()?, UPLO::Upper, self.as_allocated_mut()?)? };
Ok(FactorizedH {
Ok(BKFactorized {
a: self,
ipiv: ipiv,
})
Expand All @@ -190,10 +166,10 @@ where
A: Scalar,
Si: Data<Elem = A>,
{
fn factorizeh(&self) -> Result<FactorizedH<OwnedRepr<A>>> {
fn factorizeh(&self) -> Result<BKFactorized<OwnedRepr<A>>> {
let mut a: Array2<A> = replicate(self);
let ipiv = unsafe { A::bk(a.layout()?, UPLO::Upper, a.as_allocated_mut()?)? };
Ok(FactorizedH { a: a, ipiv: ipiv })
Ok(BKFactorized { a: a, ipiv: ipiv })
}
}

Expand Down Expand Up @@ -221,6 +197,42 @@ pub trait InverseHInto {
fn invh_into(self) -> Result<Self::Output>;
}

impl<A, S> InverseHInto for BKFactorized<S>
where
A: Scalar,
S: DataMut<Elem = A>,
{
type Output = ArrayBase<S, Ix2>;

fn invh_into(mut self) -> Result<ArrayBase<S, Ix2>> {
unsafe {
A::invh(
self.a.square_layout()?,
UPLO::Upper,
self.a.as_allocated_mut()?,
&self.ipiv,
)?
};
Ok(self.a)
}
}

impl<A, S> InverseH for BKFactorized<S>
where
A: Scalar,
S: Data<Elem = A>,
{
type Output = Array2<A>;

fn invh(&self) -> Result<Self::Output> {
let f = BKFactorized {
a: replicate(&self.a),
ipiv: self.ipiv.clone(),
};
f.invh_into()
}
}

impl<A, S> InverseHInto for ArrayBase<S, Ix2>
where
A: Scalar,
Expand All @@ -230,7 +242,7 @@ where

fn invh_into(self) -> Result<Self::Output> {
let f = self.factorizeh_into()?;
f.into_inverseh()
f.invh_into()
}
}

Expand All @@ -243,6 +255,6 @@ where

fn invh(&self) -> Result<Self::Output> {
let f = self.factorizeh()?;
f.into_inverseh()
f.invh_into()
}
}