Update blas and lapack

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"
@@ -21,5 +21,4 @@ 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 }

src/cholesky.rs

@@ -34,7 +34,7 @@ where
     S: DataMut<Elem = A>,
 {
     fn cholesky_into(mut self, uplo: UPLO) -> Result<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))
     }
 }
@@ -45,7 +45,7 @@ where
     S: DataMut<Elem = A>,
 {
     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<Self::Output> {
         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))
     }
 }

src/eigh.rs

@@ -63,7 +63,7 @@ where
     type EigVal = Array1<A::Real>;
 
     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<A::Real>;
 
     fn eigvalsh_mut(&mut self, uplo: UPLO) -> Result<Self::EigVal> {
-        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))
     }
 }

src/lapack_traits/cholesky.rs

@@ -9,13 +9,13 @@ 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 = $potrf(l.lapacke_layout(), uplo as u8, n, &mut a, n);
         into_result(info, ())

src/lapack_traits/eigh.rs

@@ -11,13 +11,13 @@ 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<Vec<Self::Real>>;
+    unsafe fn eigh(calc_eigenvec: bool, MatrixLayout, UPLO, a: &mut [Self]) -> Result<Vec<Self::Real>>;
 }
 
 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<Vec<Self::Real>> {
+    unsafe fn eigh(calc_v: bool, l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result<Vec<Self::Real>> {
         let (n, _) = l.size();
         let jobz = if calc_v { b'V' } else { b'N' };
         let mut w = vec![Self::Real::zero(); n as usize];

src/lapack_traits/opnorm.rs

@@ -24,13 +24,13 @@ 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)) => $lange(cm, t as u8, lda, col, a, lda),
             MatrixLayout::C((row, lda)) => $lange(cm, t.transpose() as u8, lda, row, a, lda),

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<Vec<Self>>;
-    fn q(MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>;
-    fn qr(MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>>;
+    unsafe fn householder(MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>>;
+    unsafe fn q(MatrixLayout, a: &mut [Self], tau: &[Self]) -> Result<()>;
+    unsafe fn qr(MatrixLayout, a: &mut [Self]) -> Result<Vec<Self>>;
 }
 
 macro_rules! impl_qr {
     ($scalar:ty, $qrf:path, $gqr:path) => {
 impl QR_ for $scalar {
-    fn householder(l: MatrixLayout, mut a: &mut [Self]) -> Result<Vec<Self>> {
+    unsafe fn householder(l: MatrixLayout, mut a: &mut [Self]) -> Result<Vec<Self>> {
         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);
         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 = $gqr(l.lapacke_layout(), row, k, k, &mut a, l.lda(), &tau);
         into_result(info, ())
     }
 
-    fn qr(l: MatrixLayout, mut a: &mut [Self]) -> Result<Vec<Self>> {
+    unsafe fn qr(l: MatrixLayout, mut a: &mut [Self]) -> Result<Vec<Self>> {
         let tau = Self::householder(l, a)?;
         let r = Vec::from(&*a);
         Self::q(l, a, &tau)?;

src/lapack_traits/solve.rs

@@ -12,30 +12,30 @@ pub type Pivot = Vec<i32>;
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`
 pub trait Solve_: Sized {
-    fn lu(MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
-    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<Pivot>;
+    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<Pivot> {
+    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)
     }
 
-    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 = $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;

src/lapack_traits/svd.rs

@@ -29,14 +29,14 @@ pub struct SVDOutput<A: AssociatedReal> {
 
 /// Wraps `*gesvd`
 pub trait SVD_: AssociatedReal {
-    fn svd(MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result<SVDOutput<Self>>;
+    unsafe fn svd(MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result<SVDOutput<Self>>;
 }
 
 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<SVDOutput<Self>> {
+    unsafe fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, mut a: &mut [Self]) -> Result<SVDOutput<Self>> {
         let (m, n) = l.size();
         let k = ::std::cmp::min(n, m);
         let lda = l.lda();

src/lapack_traits/triangular.rs

@@ -16,32 +16,33 @@ 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 = $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 = $trtrs(al.lapacke_layout(), uplo as u8, Transpose::No as u8, diag as u8, n, nrhs, a, lda, &mut b, ldb);
         into_result(info, ())
     }

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;

src/opnorm.rs

@@ -43,6 +43,6 @@ where
     fn opnorm(&self, t: NormType) -> Result<Self::Output> {
         let l = self.layout()?;
         let a = self.as_allocated()?;
-        Ok(A::opnorm(t, l, a))
+        Ok(unsafe { A::opnorm(t, l, a) })
     }
 }

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)))

src/solve.rs

@@ -23,13 +23,15 @@ where
     where
         Sb: DataMut<Elem = A>,
     {
-        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<Elem = A>,
 {
     pub fn into_inverse(mut self) -> Result<ArrayBase<S, Ix2>> {
-        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<Elem = A>,
 {
     fn factorize_into(mut self) -> Result<Factorized<S>> {
-        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<Factorized<So>> {
         let mut a: ArrayBase<So, Ix2> = 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 })
     }
 }

src/svd.rs

@@ -73,7 +73,7 @@ where
 
     fn svd_mut(&mut self, calc_u: bool, calc_vt: bool) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
         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(

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)
     }
 }