use lapacke;
use error::*;
use layout::MatrixLayout;
use num_traits::Zero;
use types::*;
use super::{Pivot, Transpose, into_result};
use super::NormType;
pub trait Solve_: AssociatedReal + Sized {
unsafe fn lu(MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
unsafe fn inv(MatrixLayout, a: &mut [Self], &Pivot) -> Result<()>;
unsafe fn rcond(MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real>;
unsafe fn solve(MatrixLayout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;
}
macro_rules! impl_solve {
($scalar:ty, $getrf:path, $getri:path, $gecon:path, $getrs:path) => {
impl Solve_ for $scalar {
unsafe fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot> {
let (row, col) = l.size();
let k = ::std::cmp::min(row, col);
let mut ipiv = vec![0; k as usize];
let info = $getrf(l.lapacke_layout(), row, col, a, l.lda(), &mut ipiv);
into_result(info, ipiv)
}
unsafe fn inv(l: MatrixLayout, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
let (n, _) = l.size();
let info = $getri(l.lapacke_layout(), n, a, l.lda(), ipiv);
into_result(info, ())
}
unsafe fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
let (n, _) = l.size();
let mut rcond = Self::Real::zero();
let info = $gecon(l.lapacke_layout(), NormType::One as u8, n, a, l.lda(), anorm, &mut rcond);
into_result(info, rcond)
}
unsafe fn solve(l: MatrixLayout, t: Transpose, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()> {
let (n, _) = l.size();
let nrhs = 1;
let ldb = 1;
let info = $getrs(l.lapacke_layout(), t as u8, n, nrhs, a, l.lda(), ipiv, b, ldb);
into_result(info, ())
}
}
}}
impl_solve!(
f64,
lapacke::dgetrf,
lapacke::dgetri,
lapacke::dgecon,
lapacke::dgetrs
);
impl_solve!(
f32,
lapacke::sgetrf,
lapacke::sgetri,
lapacke::sgecon,
lapacke::sgetrs
);
impl_solve!(
c64,
lapacke::zgetrf,
lapacke::zgetri,
lapacke::zgecon,
lapacke::zgetrs
);
impl_solve!(
c32,
lapacke::cgetrf,
lapacke::cgetri,
lapacke::cgecon,
lapacke::cgetrs
);