Merge pull request #55 from termoshtt/convert

Split Array conversion functions

Author: termoshtt

src/cholesky.rs

@@ -3,6 +3,7 @@
 use ndarray::*;
 use num_traits::Zero;
 
+use super::convert::*;
 use super::error::*;
 use super::layout::*;
 use super::triangular::IntoTriangular;

src/convert.rs

@@ -0,0 +1,98 @@
+use ndarray::*;
+
+use super::error::*;
+use super::layout::*;
+
+pub fn into_col<S>(a: ArrayBase<S, Ix1>) -> ArrayBase<S, Ix2>
+where
+    S: Data,
+{
+    let n = a.len();
+    a.into_shape((n, 1)).unwrap()
+}
+
+pub fn into_row<S>(a: ArrayBase<S, Ix1>) -> ArrayBase<S, Ix2>
+where
+    S: Data,
+{
+    let n = a.len();
+    a.into_shape((1, n)).unwrap()
+}
+
+pub fn flatten<S>(a: ArrayBase<S, Ix2>) -> ArrayBase<S, Ix1>
+where
+    S: Data,
+{
+    let n = a.len();
+    a.into_shape((n)).unwrap()
+}
+
+pub fn into_matrix<A, S>(l: MatrixLayout, a: Vec<A>) -> Result<ArrayBase<S, Ix2>>
+where
+    S: DataOwned<Elem = A>,
+{
+    Ok(ArrayBase::from_shape_vec(l.as_shape(), a)?)
+}
+
+fn uninitialized<A, S>(l: MatrixLayout) -> ArrayBase<S, Ix2>
+where
+    A: Copy,
+    S: DataOwned<Elem = A>,
+{
+    unsafe { ArrayBase::uninitialized(l.as_shape()) }
+}
+
+pub fn replicate<A, Sv, So, D>(a: &ArrayBase<Sv, D>) -> ArrayBase<So, D>
+where
+    A: Copy,
+    Sv: Data<Elem = A>,
+    So: DataOwned<Elem = A> + DataMut,
+    D: Dimension,
+{
+    let mut b = unsafe { ArrayBase::uninitialized(a.dim()) };
+    b.assign(a);
+    b
+}
+
+fn clone_with_layout<A, Si, So>(l: MatrixLayout, a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
+where
+    A: Copy,
+    Si: Data<Elem = A>,
+    So: DataOwned<Elem = A> + DataMut,
+{
+    let mut b = uninitialized(l);
+    b.assign(a);
+    b
+}
+
+pub fn transpose_data<A, S>(a: &mut ArrayBase<S, Ix2>) -> Result<&mut ArrayBase<S, Ix2>>
+where
+    A: Copy,
+    S: DataOwned<Elem = A> + DataMut,
+{
+    let l = a.layout()?.toggle_order();
+    let new = clone_with_layout(l, a);
+    ::std::mem::replace(a, new);
+    Ok(a)
+}
+
+pub fn generalize<A, S, D>(a: Array<A, D>) -> ArrayBase<S, D>
+where
+    S: DataOwned<Elem = A>,
+    D: Dimension,
+{
+    // FIXME
+    // https://github.com/bluss/rust-ndarray/issues/325
+    let strides: Vec<isize> = a.strides().to_vec();
+    let new = if a.is_standard_layout() {
+        ArrayBase::from_shape_vec(a.dim(), a.into_raw_vec()).unwrap()
+    } else {
+        ArrayBase::from_shape_vec(a.dim().f(), a.into_raw_vec()).unwrap()
+    };
+    assert_eq!(
+        new.strides(),
+        strides.as_slice(),
+        "Custom stride is not supported"
+    );
+    new
+}

src/eigh.rs

@@ -2,6 +2,7 @@
 
 use ndarray::*;
 
+use super::convert::*;
 use super::error::*;
 use super::layout::*;
 

src/generate.rs

@@ -4,8 +4,8 @@ use ndarray::*;
 use rand::*;
 use std::ops::*;
 
+use super::convert::*;
 use super::error::*;
-use super::layout::*;
 use super::types::*;
 
 /// Hermite conjugate matrix

src/lapack_traits/cholesky.rs

@@ -3,19 +3,19 @@
 use lapack::c;
 
 use error::*;
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 use super::{UPLO, into_result};
 
 pub trait Cholesky_: Sized {
-    fn cholesky(Layout, UPLO, a: &mut [Self]) -> Result<()>;
+    fn cholesky(MatrixLayout, UPLO, a: &mut [Self]) -> Result<()>;
 }
 
 macro_rules! impl_cholesky {
     ($scalar:ty, $potrf:path) => {
 impl Cholesky_ for $scalar {
-    fn cholesky(l: Layout, uplo: UPLO, mut a: &mut [Self]) -> Result<()> {
+    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

@@ -4,20 +4,20 @@ use lapack::c;
 use num_traits::Zero;
 
 use error::*;
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 use super::{UPLO, into_result};
 
 /// Wraps `*syev` for real and `*heev` for complex
 pub trait Eigh_: AssociatedReal {
-    fn eigh(calc_eigenvec: bool, Layout, UPLO, a: &mut [Self]) -> Result<Vec<Self::Real>>;
+    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: Layout, uplo: UPLO, mut a: &mut [Self]) -> Result<Vec<Self::Real>> {
+    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

@@ -3,7 +3,7 @@
 use lapack::c;
 use lapack::c::Layout::ColumnMajor as cm;
 
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 #[repr(u8)]
@@ -24,16 +24,16 @@ impl NormType {
 }
 
 pub trait OperatorNorm_: AssociatedReal {
-    fn opnorm(NormType, Layout, &[Self]) -> Self::Real;
+    fn opnorm(NormType, MatrixLayout, &[Self]) -> Self::Real;
 }
 
 macro_rules! impl_opnorm {
     ($scalar:ty, $lange:path) => {
 impl OperatorNorm_ for $scalar {
-    fn opnorm(t: NormType, l: Layout, a: &[Self]) -> Self::Real {
+    fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real {
         match l {
-            Layout::F((col, lda)) => $lange(cm, t as u8, lda, col, a, lda),
-            Layout::C((row, lda)) => $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),
         }
     }
 }

src/lapack_traits/qr.rs

@@ -5,37 +5,37 @@ use num_traits::Zero;
 use std::cmp::min;
 
 use error::*;
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 use super::into_result;
 
 /// Wraps `*geqrf` and `*orgqr` (`*ungqr` for complex numbers)
 pub trait QR_: Sized {
-    fn householder(Layout, a: &mut [Self]) -> Result<Vec<Self>>;
-    fn q(Layout, a: &mut [Self], tau: &[Self]) -> Result<()>;
-    fn qr(Layout, a: &mut [Self]) -> Result<Vec<Self>>;
+    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>>;
 }
 
 macro_rules! impl_qr {
     ($scalar:ty, $qrf:path, $gqr:path) => {
 impl QR_ for $scalar {
-    fn householder(l: Layout, mut a: &mut [Self]) -> Result<Vec<Self>> {
+    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: Layout, mut a: &mut [Self], tau: &[Self]) -> Result<()> {
+    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: Layout, mut a: &mut [Self]) -> Result<Vec<Self>> {
+    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

@@ -3,7 +3,7 @@
 use lapack::c;
 
 use error::*;
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 use super::{Transpose, into_result};
@@ -12,30 +12,30 @@ pub type Pivot = Vec<i32>;
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`
 pub trait Solve_: Sized {
-    fn lu(Layout, a: &mut [Self]) -> Result<Pivot>;
-    fn inv(Layout, a: &mut [Self], &Pivot) -> Result<()>;
-    fn solve(Layout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;
+    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<()>;
 }
 
 macro_rules! impl_solve {
     ($scalar:ty, $getrf:path, $getri:path, $getrs:path) => {
 
 impl Solve_ for $scalar {
-    fn lu(l: Layout, a: &mut [Self]) -> Result<Pivot> {
+    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: Layout, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
+    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: Layout, t: Transpose, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()> {
+    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

@@ -4,7 +4,7 @@ use lapack::c;
 use num_traits::Zero;
 
 use error::*;
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 use super::into_result;
@@ -29,14 +29,14 @@ pub struct SVDOutput<A: AssociatedReal> {
 
 /// Wraps `*gesvd`
 pub trait SVD_: AssociatedReal {
-    fn svd(Layout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result<SVDOutput<Self>>;
+    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: Layout, calc_u: bool, calc_vt: bool, mut a: &mut [Self]) -> Result<SVDOutput<Self>> {
+    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

@@ -4,7 +4,7 @@ use lapack::c;
 
 use super::{Transpose, UPLO, into_result};
 use error::*;
-use layout::Layout;
+use layout::MatrixLayout;
 use types::*;
 
 #[derive(Debug, Clone, Copy)]
@@ -16,22 +16,22 @@ pub enum Diag {
 
 /// Wraps `*trtri` and `*trtrs`
 pub trait Triangular_: Sized {
-    fn inv_triangular(l: Layout, UPLO, Diag, a: &mut [Self]) -> Result<()>;
-    fn solve_triangular(al: Layout, bl: Layout, UPLO, Diag, a: &[Self], b: &mut [Self]) -> Result<()>;
+    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<()>;
 }
 
 macro_rules! impl_triangular {
     ($scalar:ty, $trtri:path, $trtrs:path) => {
 
 impl Triangular_ for $scalar {
-    fn inv_triangular(l: Layout, uplo: UPLO, diag: Diag, a: &mut [Self]) -> Result<()> {
+    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: Layout, bl: Layout, uplo: UPLO, diag: Diag, a: &[Self], mut b: &mut [Self]) -> Result<()> {
+    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();