Merge branch 'outer'

Author: termoshtt

Cargo.toml

@@ -10,7 +10,8 @@ keywords = ["ndarray", "lapack", "matrix"]
 license = "MIT"
 
 [dependencies]
-lapack = "0.11.1"
+lapack = "0.11"
+blas = "0.15"
 num-traits = "0.1.36"
 
 [dependencies.ndarray]

src/lib.rs

@@ -32,6 +32,7 @@
 //! - [symmetric square root](hermite/trait.HermiteMatrix.html#tymethod.ssqrt)
 //! - [Cholesky factorization](hermite/trait.HermiteMatrix.html#tymethod.cholesky)
 
+extern crate blas;
 extern crate lapack;
 extern crate num_traits;
 #[macro_use(s)]
@@ -45,6 +46,7 @@ pub mod square;
 pub mod hermite;
 pub mod triangular;
 
+pub mod outer;
 pub mod qr;
 pub mod svd;
 pub mod eigh;

src/outer.rs

@@ -0,0 +1,18 @@
+
+use blas::c::{Layout, dger, sger};
+
+pub trait ImplOuter: Sized {
+    fn outer(m: usize, n: usize, a: &[Self], b: &[Self], ab: &mut [Self]);
+}
+
+macro_rules! impl_cholesky {
+    ($scalar:ty, $ger:path) => {
+impl ImplOuter for $scalar {
+    fn outer(m: usize, n: usize, a: &[Self], b: &[Self], mut ab: &mut [Self]) {
+        $ger(Layout::ColumnMajor, m as i32, n as i32, 1.0, a, 1, b, 1, ab, m as i32);
+    }
+}
+}} // end macro_rules
+
+impl_cholesky!(f64, dger);
+impl_cholesky!(f32, sger);

src/vector.rs

@@ -1,8 +1,9 @@
 //! Define trait for vectors
 
 use std::iter::Sum;
-use ndarray::{LinalgScalar, ArrayBase, Data, Dimension};
+use ndarray::{Array, NdFloat, Ix1, Ix2, LinalgScalar, ArrayBase, Data, Dimension};
 use num_traits::float::Float;
+use super::outer::ImplOuter;
 
 /// Methods for vectors
 pub trait Vector {
@@ -55,3 +56,19 @@ impl<A: Float> Squared for A {
         self.abs()
     }
 }
+
+pub fn outer<A, S1, S2>(a: &ArrayBase<S1, Ix1>, b: &ArrayBase<S2, Ix1>) -> Array<A, Ix2>
+    where A: NdFloat + ImplOuter,
+          S1: Data<Elem = A>,
+          S2: Data<Elem = A>
+{
+    let m = a.len();
+    let n = b.len();
+    let mut ab = Array::zeros((n, m));
+    ImplOuter::outer(m,
+                     n,
+                     a.as_slice_memory_order().unwrap(),
+                     b.as_slice_memory_order().unwrap(),
+                     ab.as_slice_memory_order_mut().unwrap());
+    ab.reversed_axes()
+}

tests/outer.rs

@@ -0,0 +1,20 @@
+include!("header.rs");
+use ndarray_linalg::vector::outer;
+
+#[test]
+fn outer_() {
+    let dist = RealNormal::<f64>::new(0.0, 1.0);
+    let m = 2;
+    let n = 3;
+    let a = Array::random(m, dist);
+    let b = Array::random(n, dist);
+    println!("a = \n{:?}", &a);
+    println!("b = \n{:?}", &b);
+    let ab = outer(&a, &b);
+    println!("ab = \n{:?}", &ab);
+    for i in 0..m {
+        for j in 0..n {
+            ab[(i, j)].assert_close(a[i] * b[j], 1e-7);
+        }
+    }
+}