From 37544642a390906607d8bf258607ed755e6cab86 Mon Sep 17 00:00:00 2001 From: Ivan Ukhov Date: Sun, 9 Jul 2017 07:50:04 +0200 Subject: [PATCH 1/3] Update blas and lapack --- Cargo.toml | 4 ++-- src/lapack_traits/cholesky.rs | 2 +- src/lapack_traits/eigh.rs | 2 +- src/lapack_traits/opnorm.rs | 4 ++-- src/lapack_traits/qr.rs | 4 ++-- src/lapack_traits/solve.rs | 6 +++--- src/lapack_traits/svd.rs | 2 +- src/lapack_traits/triangular.rs | 4 ++-- 8 files changed, 14 insertions(+), 14 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index dbe37f12..08ca5966 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,5 +21,5 @@ enum-error-derive = "0.1" num-traits = "0.1" num-complex = "0.1" ndarray = { version = "0.9", default-features = false, features = ["blas"] } -lapack = { version = "0.11", default-features = false } -blas = { version = "0.15", default-features = false } +lapack = { version = "0.13", default-features = false } +blas = { version = "0.16", default-features = false } diff --git a/src/lapack_traits/cholesky.rs b/src/lapack_traits/cholesky.rs index 94767d97..66dcf1aa 100644 --- a/src/lapack_traits/cholesky.rs +++ b/src/lapack_traits/cholesky.rs @@ -17,7 +17,7 @@ macro_rules! impl_cholesky { impl Cholesky_ for $scalar { fn cholesky(l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result<()> { let (n, _) = l.size(); - let info = $potrf(l.lapacke_layout(), uplo as u8, n, &mut a, n); + let info = unsafe { $potrf(l.lapacke_layout(), uplo as u8, n, &mut a, n) }; into_result(info, ()) } } diff --git a/src/lapack_traits/eigh.rs b/src/lapack_traits/eigh.rs index 490660f5..d66aa09f 100644 --- a/src/lapack_traits/eigh.rs +++ b/src/lapack_traits/eigh.rs @@ -21,7 +21,7 @@ impl Eigh_ for $scalar { let (n, _) = l.size(); let jobz = if calc_v { b'V' } else { b'N' }; let mut w = vec![Self::Real::zero(); n as usize]; - let info = $ev(l.lapacke_layout(), jobz, uplo as u8, n, &mut a, n, &mut w); + let info = unsafe { $ev(l.lapacke_layout(), jobz, uplo as u8, n, &mut a, n, &mut w) }; into_result(info, w) } } diff --git a/src/lapack_traits/opnorm.rs b/src/lapack_traits/opnorm.rs index 50c5cde1..b11736ea 100644 --- a/src/lapack_traits/opnorm.rs +++ b/src/lapack_traits/opnorm.rs @@ -32,8 +32,8 @@ macro_rules! impl_opnorm { impl OperatorNorm_ for $scalar { fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real { match l { - MatrixLayout::F((col, lda)) => $lange(cm, t as u8, lda, col, a, lda), - MatrixLayout::C((row, lda)) => $lange(cm, t.transpose() as u8, lda, row, a, lda), + MatrixLayout::F((col, lda)) => unsafe { $lange(cm, t as u8, lda, col, a, lda) }, + MatrixLayout::C((row, lda)) => unsafe { $lange(cm, t.transpose() as u8, lda, row, a, lda) }, } } } diff --git a/src/lapack_traits/qr.rs b/src/lapack_traits/qr.rs index 1f06eb1d..d35d31f1 100644 --- a/src/lapack_traits/qr.rs +++ b/src/lapack_traits/qr.rs @@ -24,14 +24,14 @@ impl QR_ for $scalar { let (row, col) = l.size(); let k = min(row, col); let mut tau = vec![Self::zero(); k as usize]; - let info = $qrf(l.lapacke_layout(), row, col, &mut a, l.lda(), &mut tau); + let info = unsafe { $qrf(l.lapacke_layout(), row, col, &mut a, l.lda(), &mut tau) }; into_result(info, tau) } fn q(l: MatrixLayout, mut a: &mut [Self], tau: &[Self]) -> Result<()> { let (row, col) = l.size(); let k = min(row, col); - let info = $gqr(l.lapacke_layout(), row, k, k, &mut a, l.lda(), &tau); + let info = unsafe { $gqr(l.lapacke_layout(), row, k, k, &mut a, l.lda(), &tau) }; into_result(info, ()) } diff --git a/src/lapack_traits/solve.rs b/src/lapack_traits/solve.rs index 4bda4374..d1deccad 100644 --- a/src/lapack_traits/solve.rs +++ b/src/lapack_traits/solve.rs @@ -25,13 +25,13 @@ impl Solve_ for $scalar { 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); + let info = unsafe { $getrf(l.lapacke_layout(), row, col, a, l.lda(), &mut ipiv) }; into_result(info, ipiv) } 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); + let info = unsafe { $getri(l.lapacke_layout(), n, a, l.lda(), ipiv) }; into_result(info, ()) } @@ -39,7 +39,7 @@ impl Solve_ for $scalar { 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); + let info = unsafe { $getrs(l.lapacke_layout(), t as u8, n, nrhs, a, l.lda(), ipiv, b, ldb) }; into_result(info, ()) } } diff --git a/src/lapack_traits/svd.rs b/src/lapack_traits/svd.rs index 75233910..8db47b03 100644 --- a/src/lapack_traits/svd.rs +++ b/src/lapack_traits/svd.rs @@ -52,7 +52,7 @@ impl SVD_ for $scalar { }; let mut s = vec![Self::Real::zero(); k as usize]; let mut superb = vec![Self::Real::zero(); (k-2) as usize]; - let info = $gesvd(l.lapacke_layout(), ju as u8, jvt as u8, m, n, &mut a, lda, &mut s, &mut u, ldu, &mut vt, ldvt, &mut superb); + let info = unsafe { $gesvd(l.lapacke_layout(), ju as u8, jvt as u8, m, n, &mut a, lda, &mut s, &mut u, ldu, &mut vt, ldvt, &mut superb) }; into_result(info, SVDOutput { s: s, u: if ldu > 0 { Some(u) } else { None }, diff --git a/src/lapack_traits/triangular.rs b/src/lapack_traits/triangular.rs index e86622e0..a853c489 100644 --- a/src/lapack_traits/triangular.rs +++ b/src/lapack_traits/triangular.rs @@ -27,7 +27,7 @@ impl Triangular_ for $scalar { fn inv_triangular(l: MatrixLayout, uplo: UPLO, diag: Diag, a: &mut [Self]) -> Result<()> { let (n, _) = l.size(); let lda = l.lda(); - let info = $trtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda); + let info = unsafe { $trtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda) }; into_result(info, ()) } @@ -42,7 +42,7 @@ impl Triangular_ for $scalar { println!("lda = {}", lda); println!("nrhs = {}", nrhs); println!("ldb = {}", ldb); - let info = $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb); + let info = unsafe { $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb) }; into_result(info, ()) } } From 22cc4059b8a3efc4b6c3af5f2e9d034d10fa81bc Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 10 Jul 2017 17:15:10 +0900 Subject: [PATCH 2/3] Drop blas --- Cargo.toml | 5 ++--- src/lib.rs | 1 - 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 08ca5966..0a49ad1f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,8 +11,8 @@ license = "MIT" [features] default = ["openblas"] -openblas = ["blas/openblas", "lapack/openblas"] -netlib = ["blas/netlib", "lapack/netlib"] +openblas = ["lapack/openblas"] +netlib = ["lapack/netlib"] [dependencies] rand = "0.3" @@ -22,4 +22,3 @@ num-traits = "0.1" num-complex = "0.1" ndarray = { version = "0.9", default-features = false, features = ["blas"] } lapack = { version = "0.13", default-features = false } -blas = { version = "0.16", default-features = false } diff --git a/src/lib.rs b/src/lib.rs index f1e52fd5..fab1aadf 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -17,7 +17,6 @@ //! - [generator functions](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) -extern crate blas; extern crate lapack; extern crate num_traits; extern crate num_complex; From 00b7e814628967e01d45eef6c5584270815d934c Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 10 Jul 2017 17:20:26 +0900 Subject: [PATCH 3/3] Make lapack_traits functions unsafe --- src/cholesky.rs | 6 +++--- src/eigh.rs | 4 ++-- src/lapack_traits/cholesky.rs | 6 +++--- src/lapack_traits/eigh.rs | 6 +++--- src/lapack_traits/opnorm.rs | 8 ++++---- src/lapack_traits/qr.rs | 16 ++++++++-------- src/lapack_traits/solve.rs | 18 +++++++++--------- src/lapack_traits/svd.rs | 6 +++--- src/lapack_traits/triangular.rs | 25 +++++++++++++------------ src/opnorm.rs | 2 +- src/qr.rs | 4 ++-- src/solve.rs | 32 ++++++++++++++++++-------------- src/svd.rs | 2 +- src/triangular.rs | 2 +- 14 files changed, 71 insertions(+), 66 deletions(-) diff --git a/src/cholesky.rs b/src/cholesky.rs index ba2fb36b..26bb2f83 100644 --- a/src/cholesky.rs +++ b/src/cholesky.rs @@ -34,7 +34,7 @@ where S: DataMut, { fn cholesky_into(mut self, uplo: UPLO) -> Result { - A::cholesky(self.square_layout()?, uplo, self.as_allocated_mut()?)?; + unsafe { A::cholesky(self.square_layout()?, uplo, self.as_allocated_mut()?)? }; Ok(self.into_triangular(uplo)) } } @@ -45,7 +45,7 @@ where S: DataMut, { fn cholesky_mut(&mut self, uplo: UPLO) -> Result<&mut Self> { - A::cholesky(self.square_layout()?, uplo, self.as_allocated_mut()?)?; + unsafe { A::cholesky(self.square_layout()?, uplo, self.as_allocated_mut()?)? }; Ok(self.into_triangular(uplo)) } } @@ -59,7 +59,7 @@ where fn cholesky(&self, uplo: UPLO) -> Result { let mut a = replicate(self); - A::cholesky(a.square_layout()?, uplo, a.as_allocated_mut()?)?; + unsafe { A::cholesky(a.square_layout()?, uplo, a.as_allocated_mut()?)? }; Ok(a.into_triangular(uplo)) } } diff --git a/src/eigh.rs b/src/eigh.rs index 6e9cb2fb..332c5462 100644 --- a/src/eigh.rs +++ b/src/eigh.rs @@ -63,7 +63,7 @@ where type EigVal = Array1; fn eigh_mut(&mut self, uplo: UPLO) -> Result<(Self::EigVal, &mut Self)> { - let s = A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)?; + let s = unsafe { A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)? }; Ok((ArrayBase::from_vec(s), self)) } } @@ -119,7 +119,7 @@ where type EigVal = Array1; fn eigvalsh_mut(&mut self, uplo: UPLO) -> Result { - let s = A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)?; + let s = unsafe { A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)? }; Ok(ArrayBase::from_vec(s)) } } diff --git a/src/lapack_traits/cholesky.rs b/src/lapack_traits/cholesky.rs index 66dcf1aa..b26be99d 100644 --- a/src/lapack_traits/cholesky.rs +++ b/src/lapack_traits/cholesky.rs @@ -9,15 +9,15 @@ use types::*; use super::{UPLO, into_result}; pub trait Cholesky_: Sized { - fn cholesky(MatrixLayout, UPLO, a: &mut [Self]) -> Result<()>; + unsafe fn cholesky(MatrixLayout, UPLO, a: &mut [Self]) -> Result<()>; } macro_rules! impl_cholesky { ($scalar:ty, $potrf:path) => { impl Cholesky_ for $scalar { - fn cholesky(l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result<()> { + unsafe fn cholesky(l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result<()> { let (n, _) = l.size(); - let info = unsafe { $potrf(l.lapacke_layout(), uplo as u8, n, &mut a, n) }; + let info = $potrf(l.lapacke_layout(), uplo as u8, n, &mut a, n); into_result(info, ()) } } diff --git a/src/lapack_traits/eigh.rs b/src/lapack_traits/eigh.rs index d66aa09f..de6ae285 100644 --- a/src/lapack_traits/eigh.rs +++ b/src/lapack_traits/eigh.rs @@ -11,17 +11,17 @@ use super::{UPLO, into_result}; /// Wraps `*syev` for real and `*heev` for complex pub trait Eigh_: AssociatedReal { - fn eigh(calc_eigenvec: bool, MatrixLayout, UPLO, a: &mut [Self]) -> Result>; + unsafe fn eigh(calc_eigenvec: bool, MatrixLayout, UPLO, a: &mut [Self]) -> Result>; } macro_rules! impl_eigh { ($scalar:ty, $ev:path) => { impl Eigh_ for $scalar { - fn eigh(calc_v: bool, l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result> { + unsafe fn eigh(calc_v: bool, l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result> { let (n, _) = l.size(); let jobz = if calc_v { b'V' } else { b'N' }; let mut w = vec![Self::Real::zero(); n as usize]; - let info = unsafe { $ev(l.lapacke_layout(), jobz, uplo as u8, n, &mut a, n, &mut w) }; + let info = $ev(l.lapacke_layout(), jobz, uplo as u8, n, &mut a, n, &mut w); into_result(info, w) } } diff --git a/src/lapack_traits/opnorm.rs b/src/lapack_traits/opnorm.rs index b11736ea..81163cff 100644 --- a/src/lapack_traits/opnorm.rs +++ b/src/lapack_traits/opnorm.rs @@ -24,16 +24,16 @@ impl NormType { } pub trait OperatorNorm_: AssociatedReal { - fn opnorm(NormType, MatrixLayout, &[Self]) -> Self::Real; + unsafe fn opnorm(NormType, MatrixLayout, &[Self]) -> Self::Real; } macro_rules! impl_opnorm { ($scalar:ty, $lange:path) => { impl OperatorNorm_ for $scalar { - fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real { + unsafe fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real { match l { - MatrixLayout::F((col, lda)) => unsafe { $lange(cm, t as u8, lda, col, a, lda) }, - MatrixLayout::C((row, lda)) => unsafe { $lange(cm, t.transpose() as u8, lda, row, a, lda) }, + MatrixLayout::F((col, lda)) => $lange(cm, t as u8, lda, col, a, lda), + MatrixLayout::C((row, lda)) => $lange(cm, t.transpose() as u8, lda, row, a, lda), } } } diff --git a/src/lapack_traits/qr.rs b/src/lapack_traits/qr.rs index d35d31f1..290a0cf7 100644 --- a/src/lapack_traits/qr.rs +++ b/src/lapack_traits/qr.rs @@ -12,30 +12,30 @@ use super::into_result; /// Wraps `*geqrf` and `*orgqr` (`*ungqr` for complex numbers) pub trait QR_: Sized { - fn householder(MatrixLayout, a: &mut [Self]) -> Result>; - fn q(MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>; - fn qr(MatrixLayout, a: &mut [Self]) -> Result>; + unsafe fn householder(MatrixLayout, a: &mut [Self]) -> Result>; + unsafe fn q(MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>; + unsafe fn qr(MatrixLayout, a: &mut [Self]) -> Result>; } macro_rules! impl_qr { ($scalar:ty, $qrf:path, $gqr:path) => { impl QR_ for $scalar { - fn householder(l: MatrixLayout, mut a: &mut [Self]) -> Result> { + unsafe fn householder(l: MatrixLayout, mut a: &mut [Self]) -> Result> { let (row, col) = l.size(); let k = min(row, col); let mut tau = vec![Self::zero(); k as usize]; - let info = unsafe { $qrf(l.lapacke_layout(), row, col, &mut a, l.lda(), &mut tau) }; + let info = $qrf(l.lapacke_layout(), row, col, &mut a, l.lda(), &mut tau); into_result(info, tau) } - fn q(l: MatrixLayout, mut a: &mut [Self], tau: &[Self]) -> Result<()> { + unsafe fn q(l: MatrixLayout, mut a: &mut [Self], tau: &[Self]) -> Result<()> { let (row, col) = l.size(); let k = min(row, col); - let info = unsafe { $gqr(l.lapacke_layout(), row, k, k, &mut a, l.lda(), &tau) }; + let info = $gqr(l.lapacke_layout(), row, k, k, &mut a, l.lda(), &tau); into_result(info, ()) } - fn qr(l: MatrixLayout, mut a: &mut [Self]) -> Result> { + unsafe fn qr(l: MatrixLayout, mut a: &mut [Self]) -> Result> { let tau = Self::householder(l, a)?; let r = Vec::from(&*a); Self::q(l, a, &tau)?; diff --git a/src/lapack_traits/solve.rs b/src/lapack_traits/solve.rs index d1deccad..1b7a2b1a 100644 --- a/src/lapack_traits/solve.rs +++ b/src/lapack_traits/solve.rs @@ -12,34 +12,34 @@ pub type Pivot = Vec; /// Wraps `*getrf`, `*getri`, and `*getrs` pub trait Solve_: Sized { - fn lu(MatrixLayout, a: &mut [Self]) -> Result; - fn inv(MatrixLayout, a: &mut [Self], &Pivot) -> Result<()>; - fn solve(MatrixLayout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>; + unsafe fn lu(MatrixLayout, a: &mut [Self]) -> Result; + unsafe fn inv(MatrixLayout, a: &mut [Self], &Pivot) -> Result<()>; + unsafe fn solve(MatrixLayout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>; } macro_rules! impl_solve { ($scalar:ty, $getrf:path, $getri:path, $getrs:path) => { impl Solve_ for $scalar { - fn lu(l: MatrixLayout, a: &mut [Self]) -> Result { + unsafe fn lu(l: MatrixLayout, a: &mut [Self]) -> Result { let (row, col) = l.size(); let k = ::std::cmp::min(row, col); let mut ipiv = vec![0; k as usize]; - let info = unsafe { $getrf(l.lapacke_layout(), row, col, a, l.lda(), &mut ipiv) }; + let info = $getrf(l.lapacke_layout(), row, col, a, l.lda(), &mut ipiv); into_result(info, ipiv) } - fn inv(l: MatrixLayout, a: &mut [Self], ipiv: &Pivot) -> Result<()> { + unsafe fn inv(l: MatrixLayout, a: &mut [Self], ipiv: &Pivot) -> Result<()> { let (n, _) = l.size(); - let info = unsafe { $getri(l.lapacke_layout(), n, a, l.lda(), ipiv) }; + let info = $getri(l.lapacke_layout(), n, a, l.lda(), ipiv); into_result(info, ()) } - fn solve(l: MatrixLayout, t: Transpose, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()> { + 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 = unsafe { $getrs(l.lapacke_layout(), t as u8, n, nrhs, a, l.lda(), ipiv, b, ldb) }; + let info = $getrs(l.lapacke_layout(), t as u8, n, nrhs, a, l.lda(), ipiv, b, ldb); into_result(info, ()) } } diff --git a/src/lapack_traits/svd.rs b/src/lapack_traits/svd.rs index 8db47b03..f1784f21 100644 --- a/src/lapack_traits/svd.rs +++ b/src/lapack_traits/svd.rs @@ -29,14 +29,14 @@ pub struct SVDOutput { /// Wraps `*gesvd` pub trait SVD_: AssociatedReal { - fn svd(MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result>; + unsafe fn svd(MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result>; } macro_rules! impl_svd { ($scalar:ty, $gesvd:path) => { impl SVD_ for $scalar { - fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, mut a: &mut [Self]) -> Result> { + unsafe fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, mut a: &mut [Self]) -> Result> { let (m, n) = l.size(); let k = ::std::cmp::min(n, m); let lda = l.lda(); @@ -52,7 +52,7 @@ impl SVD_ for $scalar { }; let mut s = vec![Self::Real::zero(); k as usize]; let mut superb = vec![Self::Real::zero(); (k-2) as usize]; - let info = unsafe { $gesvd(l.lapacke_layout(), ju as u8, jvt as u8, m, n, &mut a, lda, &mut s, &mut u, ldu, &mut vt, ldvt, &mut superb) }; + let info = $gesvd(l.lapacke_layout(), ju as u8, jvt as u8, m, n, &mut a, lda, &mut s, &mut u, ldu, &mut vt, ldvt, &mut superb); into_result(info, SVDOutput { s: s, u: if ldu > 0 { Some(u) } else { None }, diff --git a/src/lapack_traits/triangular.rs b/src/lapack_traits/triangular.rs index a853c489..de885f3a 100644 --- a/src/lapack_traits/triangular.rs +++ b/src/lapack_traits/triangular.rs @@ -16,33 +16,34 @@ pub enum Diag { /// Wraps `*trtri` and `*trtrs` pub trait Triangular_: Sized { - fn inv_triangular(l: MatrixLayout, UPLO, Diag, a: &mut [Self]) -> Result<()>; - fn solve_triangular(al: MatrixLayout, bl: MatrixLayout, UPLO, Diag, a: &[Self], b: &mut [Self]) -> Result<()>; + unsafe fn inv_triangular(l: MatrixLayout, UPLO, Diag, a: &mut [Self]) -> Result<()>; + unsafe fn solve_triangular( + al: MatrixLayout, + bl: MatrixLayout, + UPLO, + Diag, + a: &[Self], + b: &mut [Self], + ) -> Result<()>; } macro_rules! impl_triangular { ($scalar:ty, $trtri:path, $trtrs:path) => { impl Triangular_ for $scalar { - fn inv_triangular(l: MatrixLayout, uplo: UPLO, diag: Diag, a: &mut [Self]) -> Result<()> { + unsafe fn inv_triangular(l: MatrixLayout, uplo: UPLO, diag: Diag, a: &mut [Self]) -> Result<()> { let (n, _) = l.size(); let lda = l.lda(); - let info = unsafe { $trtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda) }; + let info = $trtri(l.lapacke_layout(), uplo as u8, diag as u8, n, a, lda); into_result(info, ()) } - fn solve_triangular(al: MatrixLayout, bl: MatrixLayout, uplo: UPLO, diag: Diag, a: &[Self], mut b: &mut [Self]) -> Result<()> { + unsafe fn solve_triangular(al: MatrixLayout, bl: MatrixLayout, uplo: UPLO, diag: Diag, a: &[Self], mut b: &mut [Self]) -> Result<()> { let (n, _) = al.size(); let lda = al.lda(); let (_, nrhs) = bl.size(); let ldb = bl.lda(); - println!("al = {:?}", al); - println!("bl = {:?}", bl); - println!("n = {}", n); - println!("lda = {}", lda); - println!("nrhs = {}", nrhs); - println!("ldb = {}", ldb); - let info = unsafe { $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb) }; + let info = $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb); into_result(info, ()) } } diff --git a/src/opnorm.rs b/src/opnorm.rs index d5a19239..2ba7590b 100644 --- a/src/opnorm.rs +++ b/src/opnorm.rs @@ -43,6 +43,6 @@ where fn opnorm(&self, t: NormType) -> Result { let l = self.layout()?; let a = self.as_allocated()?; - Ok(A::opnorm(t, l, a)) + Ok(unsafe { A::opnorm(t, l, a) }) } } diff --git a/src/qr.rs b/src/qr.rs index b847406b..a8d941fb 100644 --- a/src/qr.rs +++ b/src/qr.rs @@ -61,7 +61,7 @@ where fn qr_square_mut<'a>(&'a mut self) -> Result<(&'a mut Self, Self::R)> { let l = self.square_layout()?; - let r = A::qr(l, self.as_allocated_mut()?)?; + let r = unsafe { A::qr(l, self.as_allocated_mut()?)? }; let r: Array2<_> = into_matrix(l, r)?; Ok((self, r.into_triangular(UPLO::Upper))) } @@ -108,7 +108,7 @@ where let m = self.cols(); let k = ::std::cmp::min(n, m); let l = self.layout()?; - let r = A::qr(l, self.as_allocated_mut()?)?; + let r = unsafe { A::qr(l, self.as_allocated_mut()?)? }; let r: Array2<_> = into_matrix(l, r)?; let q = self; Ok((take_slice(&q, n, k), take_slice_upper(&r, k, m))) diff --git a/src/solve.rs b/src/solve.rs index b526c0cb..f6cd3ef5 100644 --- a/src/solve.rs +++ b/src/solve.rs @@ -23,13 +23,15 @@ where where Sb: DataMut, { - A::solve( - self.a.square_layout()?, - t, - self.a.as_allocated()?, - &self.ipiv, - rhs.as_slice_mut().unwrap(), - )?; + unsafe { + A::solve( + self.a.square_layout()?, + t, + self.a.as_allocated()?, + &self.ipiv, + rhs.as_slice_mut().unwrap(), + )? + }; Ok(rhs) } } @@ -40,11 +42,13 @@ where S: DataMut, { pub fn into_inverse(mut self) -> Result> { - A::inv( - self.a.square_layout()?, - self.a.as_allocated_mut()?, - &self.ipiv, - )?; + unsafe { + A::inv( + self.a.square_layout()?, + self.a.as_allocated_mut()?, + &self.ipiv, + )? + }; Ok(self.a) } } @@ -63,7 +67,7 @@ where S: DataMut, { fn factorize_into(mut self) -> Result> { - let ipiv = A::lu(self.layout()?, self.as_allocated_mut()?)?; + let ipiv = unsafe { A::lu(self.layout()?, self.as_allocated_mut()?)? }; Ok(Factorized { a: self, ipiv: ipiv, @@ -79,7 +83,7 @@ where { fn factorize(&self) -> Result> { let mut a: ArrayBase = replicate(self); - let ipiv = A::lu(a.layout()?, a.as_allocated_mut()?)?; + let ipiv = unsafe { A::lu(a.layout()?, a.as_allocated_mut()?)? }; Ok(Factorized { a: a, ipiv: ipiv }) } } diff --git a/src/svd.rs b/src/svd.rs index 324fee82..c6054000 100644 --- a/src/svd.rs +++ b/src/svd.rs @@ -73,7 +73,7 @@ where fn svd_mut(&mut self, calc_u: bool, calc_vt: bool) -> Result<(Option, Self::Sigma, Option)> { let l = self.layout()?; - let svd_res = A::svd(l, calc_u, calc_vt, self.as_allocated_mut()?)?; + let svd_res = unsafe { A::svd(l, calc_u, calc_vt, self.as_allocated_mut()?)? }; let (n, m) = l.size(); let u = svd_res.u.map(|u| into_matrix(l.resized(n, n), u).unwrap()); let vt = svd_res.vt.map( diff --git a/src/triangular.rs b/src/triangular.rs index a9307310..1bb16a50 100644 --- a/src/triangular.rs +++ b/src/triangular.rs @@ -71,7 +71,7 @@ where transpose_data(b)?; } let lb = b.layout()?; - A::solve_triangular(la, lb, uplo, diag, a_, b.as_allocated_mut()?)?; + unsafe { A::solve_triangular(la, lb, uplo, diag, a_, b.as_allocated_mut()?)? }; Ok(b) } }