Merge pull request #206 from rust-ndarray/lapack

Rewrite LAPACKE by Rust, call LAPACK directly

Author: termoshtt

lax/Cargo.toml

@@ -11,10 +11,10 @@ netlib     = ["lapack-src/netlib", "blas-src/netlib"]
 openblas   = ["lapack-src/openblas", "blas-src/openblas"]
 
 [dependencies]
-thiserror = "1"
-cauchy = "0.2"
-lapacke = "0.2.0"
+thiserror = "1.0"
+cauchy = "0.2.0"
 num-traits = "0.2"
+lapack = "0.16.0"
 
 [dependencies.blas-src]
 version = "0.6.1"

lax/src/cholesky.rs

@@ -1,56 +1,90 @@
 //! Cholesky decomposition
 
 use super::*;
-use crate::{error::*, layout::MatrixLayout};
+use crate::{error::*, layout::*};
 use cauchy::*;
 
 pub trait Cholesky_: Sized {
     /// Cholesky: wrapper of `*potrf`
     ///
     /// **Warning: Only the portion of `a` corresponding to `UPLO` is written.**
-    unsafe fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+    fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+
     /// Wrapper of `*potri`
     ///
     /// **Warning: Only the portion of `a` corresponding to `UPLO` is written.**
-    unsafe fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+    fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()>;
+
     /// Wrapper of `*potrs`
-    unsafe fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self])
-        -> Result<()>;
+    fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>;
 }
 
 macro_rules! impl_cholesky {
     ($scalar:ty, $trf:path, $tri:path, $trs:path) => {
         impl Cholesky_ for $scalar {
-            unsafe fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
+            fn cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
                 let (n, _) = l.size();
-                $trf(l.lapacke_layout(), uplo as u8, n, a, n).as_lapack_result()?;
+                if matches!(l, MatrixLayout::C { .. }) {
+                    square_transpose(l, a);
+                }
+                let mut info = 0;
+                unsafe {
+                    $trf(uplo as u8, n, a, n, &mut info);
+                }
+                info.as_lapack_result()?;
+                if matches!(l, MatrixLayout::C { .. }) {
+                    square_transpose(l, a);
+                }
                 Ok(())
             }
 
-            unsafe fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
+            fn inv_cholesky(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<()> {
                 let (n, _) = l.size();
-                $tri(l.lapacke_layout(), uplo as u8, n, a, l.lda()).as_lapack_result()?;
+                if matches!(l, MatrixLayout::C { .. }) {
+                    square_transpose(l, a);
+                }
+                let mut info = 0;
+                unsafe {
+                    $tri(uplo as u8, n, a, l.lda(), &mut info);
+                }
+                info.as_lapack_result()?;
+                if matches!(l, MatrixLayout::C { .. }) {
+                    square_transpose(l, a);
+                }
                 Ok(())
             }
 
-            unsafe fn solve_cholesky(
+            fn solve_cholesky(
                 l: MatrixLayout,
-                uplo: UPLO,
+                mut uplo: UPLO,
                 a: &[Self],
                 b: &mut [Self],
             ) -> Result<()> {
                 let (n, _) = l.size();
                 let nrhs = 1;
-                let ldb = 1;
-                $trs(l.lapacke_layout(), uplo as u8, n, nrhs, a, l.lda(), b, ldb)
-                    .as_lapack_result()?;
+                let mut info = 0;
+                if matches!(l, MatrixLayout::C { .. }) {
+                    uplo = uplo.t();
+                    for val in b.iter_mut() {
+                        *val = val.conj();
+                    }
+                }
+                unsafe {
+                    $trs(uplo as u8, n, nrhs, a, l.lda(), b, n, &mut info);
+                }
+                info.as_lapack_result()?;
+                if matches!(l, MatrixLayout::C { .. }) {
+                    for val in b.iter_mut() {
+                        *val = val.conj();
+                    }
+                }
                 Ok(())
             }
         }
     };
 } // end macro_rules
 
-impl_cholesky!(f64, lapacke::dpotrf, lapacke::dpotri, lapacke::dpotrs);
-impl_cholesky!(f32, lapacke::spotrf, lapacke::spotri, lapacke::spotrs);
-impl_cholesky!(c64, lapacke::zpotrf, lapacke::zpotri, lapacke::zpotrs);
-impl_cholesky!(c32, lapacke::cpotrf, lapacke::cpotri, lapacke::cpotrs);
+impl_cholesky!(f64, lapack::dpotrf, lapack::dpotri, lapack::dpotrs);
+impl_cholesky!(f32, lapack::spotrf, lapack::spotri, lapack::spotrs);
+impl_cholesky!(c64, lapack::zpotrf, lapack::zpotri, lapack::zpotrs);
+impl_cholesky!(c32, lapack::cpotrf, lapack::cpotri, lapack::cpotrs);