diff --git a/Cargo.toml b/Cargo.toml index b477bf401..34150ca00 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -37,11 +37,11 @@ optional = true [dependencies] blas-sys = { version = "0.5", optional = true, default-features = false } -# rblas is deprecatedish +# deprecated! use ndarray-rblas instead rblas = { version = "0.0.13", optional = true } [features] - +# deprecated! Will be default assign_ops = [] blas = ["blas-sys"] diff --git a/README.rst b/README.rst index e2da4bd26..5b8cdc864 100644 --- a/README.rst +++ b/README.rst @@ -57,8 +57,7 @@ Status and Lookout + ``.fold()`` and ``.zip_mut_with()`` are the most efficient ways to perform single traversal and lock step traversal respectively. + ``.iter()`` and ``.iter_mut()`` are efficient for contiguous arrays. - -- There is experimental bridging to the linear algebra package ``rblas``. + + Can use BLAS in some operations (``dot`` and ``mat_mul``). Crate Feature Flags ------------------- @@ -78,7 +77,7 @@ your `Cargo.toml`. - ``rblas`` - - Optional, stable + - **Deprecated:** replaced by separate crate ``ndarray-rblas`` - Enables ``rblas`` integration How to use with cargo:: diff --git a/ndarray-rblas/Cargo.toml b/ndarray-rblas/Cargo.toml new file mode 100644 index 000000000..4eb5d1373 --- /dev/null +++ b/ndarray-rblas/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "ndarray-rblas" +version = "0.1.0" +authors = ["bluss"] +license = "MIT/Apache-2.0" + +repository = "https://github.com/bluss/rust-ndarray" + +description = "`rblas` bindings for `ndarray`." + +keywords = ["multidimensional", "matrix", "blas"] + +[dependencies] +rblas = "0.0.13" +ndarray = { version = "0.4.1" } + +[dev-dependencies] +num = { version = "0.1", default-features = false } diff --git a/ndarray-rblas/LICENSE-APACHE b/ndarray-rblas/LICENSE-APACHE new file mode 100644 index 000000000..16fe87b06 --- /dev/null +++ b/ndarray-rblas/LICENSE-APACHE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + +TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + +1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + +2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + +3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + +4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + +5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + +6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + +7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + +8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + +9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + +END OF TERMS AND CONDITIONS + +APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + +Copyright [yyyy] [name of copyright owner] + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/ndarray-rblas/LICENSE-MIT b/ndarray-rblas/LICENSE-MIT new file mode 100644 index 000000000..9203baa05 --- /dev/null +++ b/ndarray-rblas/LICENSE-MIT @@ -0,0 +1,25 @@ +Copyright (c) 2015 + +Permission is hereby granted, free of charge, to any +person obtaining a copy of this software and associated +documentation files (the "Software"), to deal in the +Software without restriction, including without +limitation the rights to use, copy, modify, merge, +publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software +is furnished to do so, subject to the following +conditions: + +The above copyright notice and this permission notice +shall be included in all copies or substantial portions +of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF +ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED +TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A +PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT +SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR +IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/ndarray-rblas/src/lib.rs b/ndarray-rblas/src/lib.rs new file mode 100644 index 000000000..00ffe2016 --- /dev/null +++ b/ndarray-rblas/src/lib.rs @@ -0,0 +1,397 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +//! Experimental BLAS (Basic Linear Algebra Subprograms) integration +//! +//! Depends on crate [`rblas`], ([docs]). +//! +//! [`rblas`]: https://crates.io/crates/rblas/ +//! [docs]: http://mikkyang.github.io/rust-blas/doc/rblas/ +//! +//! ``` +//! extern crate ndarray; +//! extern crate ndarray_rblas; +//! extern crate rblas; +//! +//! use rblas::Gemv; +//! use rblas::attribute::Transpose; +//! +//! use ndarray::{arr1, arr2}; +//! use ndarray_rblas::AsBlas; +//! +//! fn main() { +//! // Gemv is the operation y = α a x + β y +//! let alpha = 1.; +//! let mut a = arr2(&[[1., 2., 3.], +//! [4., 5., 6.], +//! [7., 8., 9.]]); +//! let x = [1., 0., 1.]; +//! let beta = 1.; +//! let mut y = arr1(&[0., 0., 0.]); +//! +//! Gemv::gemv(Transpose::NoTrans, &alpha, &a.blas(), &x[..], +//! &beta, &mut y.blas()); +//! +//! assert_eq!(y, arr1(&[4., 10., 16.])); +//! } +//! +//! ``` +//! +//! Use the methods in trait `AsBlas` to convert an array into a view that +//! implements rblas’ `Vector` or `Matrix` traits. +//! +//! Blas supports strided vectors and matrices; Matrices need to be contiguous +//! in their lowest dimension, so they will be copied into c-contiguous layout +//! automatically if needed. You should be able to use blocks sliced out +//! from a larger matrix without copying. Use the transpose flags in blas +//! instead of transposing with `ndarray`. +//! +//! Blas has its own error reporting system and will not panic on errors (that +//! I know), instead output its own error conditions, for example on dimension +//! mismatch in a matrix multiplication. +//! + +extern crate rblas; +extern crate ndarray; + +use std::os::raw::{c_int}; + +use rblas::{ + Matrix, + Vector, +}; +use ndarray::{ + ShapeError, + ErrorKind, + ArrayView, + ArrayViewMut, + Data, + DataOwned, + DataMut, + Dimension, + ArrayBase, + Ix, Ixs, +}; + + +/// ***Requires crate feature `"rblas"`*** +pub struct BlasArrayView<'a, A: 'a, D>(ArrayView<'a, A, D>); +impl<'a, A, D: Copy> Copy for BlasArrayView<'a, A, D> { } +impl<'a, A, D: Clone> Clone for BlasArrayView<'a, A, D> { + fn clone(&self) -> Self { + BlasArrayView(self.0.clone()) + } +} + +/// ***Requires crate feature `"rblas"`*** +pub struct BlasArrayViewMut<'a, A: 'a, D>(ArrayViewMut<'a, A, D>); + +struct Priv(T); + +/// Return `true` if the innermost dimension is contiguous (includes +/// the special cases of 0 or 1 length in that axis). +fn is_inner_contiguous(a: &ArrayBase) -> bool + where S: Data, + D: Dimension, +{ + let ndim = a.ndim(); + if ndim == 0 { + return true; + } + a.shape()[ndim - 1] <= 1 || a.strides()[ndim - 1] == 1 +} + +/// If the array is not in the standard layout, copy all elements +/// into the standard layout so that the array is C-contiguous. +fn ensure_standard_layout(a: &mut ArrayBase) + where S: DataOwned, + D: Dimension, + A: Clone +{ + if !a.is_standard_layout() { + let d = a.dim(); + let v: Vec = a.iter().cloned().collect(); + *a = ArrayBase::from_vec_dim(d, v).unwrap(); + } +} + + +impl Priv> + where S: Data, + D: Dimension +{ + fn size_check(&self) -> Result<(), ShapeError> { + let max = c_int::max_value(); + let self_ = &self.0; + for (&dim, &stride) in self_.shape().iter().zip(self_.strides()) { + if dim > max as Ix || stride > max as Ixs { + return Err(ShapeError::from_kind(ErrorKind::RangeLimited)); + } + } + Ok(()) + } + + fn contiguous_check(&self) -> Result<(), ShapeError> { + // FIXME: handle transposed? + if is_inner_contiguous(&self.0) { + Ok(()) + } else { + Err(ShapeError::from_kind(ErrorKind::IncompatibleLayout)) + } + } +} + +impl<'a, A, D> Priv> + where D: Dimension +{ + pub fn into_blas_view(self) -> Result, ShapeError> { + if self.0.ndim() > 1 { + try!(self.contiguous_check()); + } + try!(self.size_check()); + Ok(BlasArrayView(self.0)) + } +} + +impl<'a, A, D> Priv> + where D: Dimension +{ + fn into_blas_view_mut(self) -> Result, ShapeError> { + if self.0.ndim() > 1 { + try!(self.contiguous_check()); + } + try!(self.size_check()); + Ok(BlasArrayViewMut(self.0)) + } +} +/* +*/ + +/// Convert an array into a blas friendly wrapper. +/// +/// Note that `blas` suppors four different element types: `f32`, `f64`, +/// `Complex`, and `Complex`. +/// +/// ***Requires crate feature `"rblas"`*** +pub trait AsBlas { + /// Return an array view implementing Vector (1D) or Matrix (2D) + /// traits. + /// + /// Elements are copied if needed to produce a contiguous matrix.
+ /// The result is always mutable, due to the requirement of having write + /// access to update the layout either way. Breaks sharing if the array is + /// an `RcArray`. + /// + /// **Errors** if any dimension is larger than `c_int::MAX`. + fn blas_checked(&mut self) -> Result, ShapeError> + where S: DataOwned + DataMut, + A: Clone; + + /// Equivalent to `.blas_checked().unwrap()` + /// + /// **Panics** if there was a an error in `.blas_checked()`. + fn blas(&mut self) -> BlasArrayViewMut + where S: DataOwned + DataMut, + A: Clone + { + self.blas_checked().unwrap() + } + + /// Return a read-only array view implementing Vector (1D) or Matrix (2D) + /// traits. + /// + /// The array must already be in a blas compatible layout: its innermost + /// dimension must be contiguous. + /// + /// **Errors** if any dimension is larger than `c_int::MAX`.
+ /// **Errors** if the inner dimension is not c-contiguous. + /// + /// Layout requirements may be loosened in the future. + fn blas_view_checked(&self) -> Result, ShapeError> + where S: Data; + + /// `bv` stands for **b**las **v**iew. + /// + /// Equivalent to `.blas_view_checked().unwrap()` + /// + /// **Panics** if there was a an error in `.blas_view_checked()`. + fn bv(&self) -> BlasArrayView + where S: Data, + { + self.blas_view_checked().unwrap() + } + + /// Return a read-write array view implementing Vector (1D) or Matrix (2D) + /// traits. + /// + /// The array must already be in a blas compatible layout: its innermost + /// dimension must be contiguous. + /// + /// **Errors** if any dimension is larger than `c_int::MAX`.
+ /// **Errors** if the inner dimension is not c-contiguous. + /// + /// Layout requirements may be loosened in the future. + fn blas_view_mut_checked(&mut self) -> Result, ShapeError> + where S: DataMut; + + /// `bvm` stands for **b**las **v**iew **m**ut. + /// + /// Equivalent to `.blas_view_mut_checked().unwrap()` + /// + /// **Panics** if there was a an error in `.blas_view_mut_checked()`. + fn bvm(&mut self) -> BlasArrayViewMut + where S: DataMut, + { + self.blas_view_mut_checked().unwrap() + } + /* + + /// Equivalent to `.blas_checked().unwrap()`, except elements + /// are not copied to make the array contiguous: instead just + /// dimensions and strides are adjusted, and elements end up in + /// arbitrary location. Useful if the content of the array doesn't matter. + /// + /// **Panics** if there was a an error in `blas_checked`. + fn blas_overwrite(&mut self) -> BlasArrayViewMut + where S: DataMut; + */ +} + +/// ***Requires crate feature `"rblas"`*** +impl AsBlas for ArrayBase + where S: Data, + D: Dimension, +{ + fn blas_checked(&mut self) -> Result, ShapeError> + where S: DataOwned + DataMut, + A: Clone, + { + try!(Priv(self.view()).size_check()); + match self.ndim() { + 0 | 1 => { } + 2 => { + if !is_inner_contiguous(self) { + ensure_standard_layout(self); + } + } + _n => ensure_standard_layout(self), + } + Priv(self.view_mut()).into_blas_view_mut() + } + + fn blas_view_checked(&self) -> Result, ShapeError> + where S: Data + { + Priv(self.view()).into_blas_view() + } + + fn blas_view_mut_checked(&mut self) -> Result, ShapeError> + where S: DataMut, + { + Priv(self.view_mut()).into_blas_view_mut() + } + + /* + fn blas_overwrite(&mut self) -> BlasArrayViewMut + where S: DataMut, + { + self.size_check().unwrap(); + if self.dim.ndim() > 1 { + self.force_standard_layout(); + } + BlasArrayViewMut(self.view_mut()) + } + */ +} + +/// **Panics** if `as_mut_ptr` is called on a read-only view. +impl<'a, A> Vector for BlasArrayView<'a, A, Ix> { + fn len(&self) -> c_int { + self.0.len() as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + panic!("ndarray: as_mut_ptr called on BlasArrayView (not mutable)"); + } + + // increment: stride + fn inc(&self) -> c_int { + self.0.strides()[0] as c_int + } +} + +impl<'a, A> Vector for BlasArrayViewMut<'a, A, Ix> { + fn len(&self) -> c_int { + self.0.len() as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + self.0.as_mut_ptr() + } + + // increment: stride + fn inc(&self) -> c_int { + self.0.strides()[0] as c_int + } +} + +/// **Panics** if `as_mut_ptr` is called on a read-only view. +impl<'a, A> Matrix for BlasArrayView<'a, A, (Ix, Ix)> { + fn rows(&self) -> c_int { + self.0.dim().0 as c_int + } + + fn cols(&self) -> c_int { + self.0.dim().1 as c_int + } + + // leading dimension == stride between each row + fn lead_dim(&self) -> c_int { + debug_assert!(self.cols() <= 1 || self.0.strides()[1] == 1); + self.0.strides()[0] as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + panic!("ndarray: as_mut_ptr called on BlasArrayView (not mutable)"); + } +} + +impl<'a, A> Matrix for BlasArrayViewMut<'a, A, (Ix, Ix)> { + fn rows(&self) -> c_int { + self.0.dim().0 as c_int + } + + fn cols(&self) -> c_int { + self.0.dim().1 as c_int + } + + // leading dimension == stride between each row + fn lead_dim(&self) -> c_int { + debug_assert!(self.cols() <= 1 || self.0.strides()[1] == 1); + self.0.strides()[0] as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + self.0.as_mut_ptr() + } +} diff --git a/ndarray-rblas/tests/blas.rs b/ndarray-rblas/tests/blas.rs new file mode 100644 index 000000000..5825e0835 --- /dev/null +++ b/ndarray-rblas/tests/blas.rs @@ -0,0 +1,164 @@ +extern crate rblas; +extern crate num; +#[macro_use] extern crate ndarray; +extern crate ndarray_rblas; + +use rblas::Gemm; +use rblas::attribute::Transpose; + +use num::Float; + +use ndarray::{ + OwnedArray, + ArrayView, + ArrayViewMut, + arr2, + Ix, + ShapeError +}; + +use ndarray_rblas::AsBlas; + +#[test] +fn strided_matrix() { + // smoke test, a matrix multiplication of uneven size + let (n, m) = (45, 33); + let mut a = OwnedArray::linspace(0., ((n * m) - 1) as f32, n as usize * m as usize ).into_shape((n, m)).unwrap(); + let mut b = OwnedArray::eye(m); + let mut res = OwnedArray::zeros(a.dim()); + Gemm::gemm(&1., Transpose::NoTrans, &a.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, a); + + // matrix multiplication, strided + let mut aprim = a.to_shared(); + aprim.islice(s![0..12, 0..11]); + println!("{:?}", aprim.shape()); + println!("{:?}", aprim.strides()); + let mut b = OwnedArray::eye(aprim.shape()[1]); + let mut res = OwnedArray::zeros(aprim.dim()); + Gemm::gemm(&1., Transpose::NoTrans, &aprim.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, aprim); + + // Transposed matrix multiply + let (np, mp) = aprim.dim(); + let mut res = OwnedArray::zeros((mp, np)); + let mut b = OwnedArray::eye(np); + Gemm::gemm(&1., Transpose::Trans, &aprim.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + let mut at = aprim.clone(); + at.swap_axes(0, 1); + assert_eq!(at, res); + + // strided, needs copy + let mut abis = a.to_shared(); + abis.islice(s![0..12, ..;2]); + println!("{:?}", abis.shape()); + println!("{:?}", abis.strides()); + let mut b = OwnedArray::eye(abis.shape()[1]); + let mut res = OwnedArray::zeros(abis.dim()); + Gemm::gemm(&1., Transpose::NoTrans, &abis.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, abis); +} + +#[test] +fn strided_view() { + // smoke test, a matrix multiplication of uneven size + let (n, m) = (45, 33); + let mut a = OwnedArray::linspace(0., ((n * m) - 1) as f32, n as usize * m as usize ).into_shape((n, m)).unwrap(); + let mut b = OwnedArray::eye(m); + let mut res = OwnedArray::zeros(a.dim()); + Gemm::gemm(&1., + Transpose::NoTrans, &a.blas_view_mut_checked().unwrap(), + Transpose::NoTrans, &b.blas_view_mut_checked().unwrap(), + &0., &mut res.blas_view_mut_checked().unwrap()); + assert_eq!(res, a); + + // matrix multiplication, strided + let aprim = a.slice(s![0..12, 0..11]); + let mut b = OwnedArray::eye(aprim.shape()[1]); + let mut res = OwnedArray::zeros(aprim.dim()); + Gemm::gemm(&1., + Transpose::NoTrans, &aprim.bv(), + Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, aprim); + + // test out with matrices where lower axis is strided but has length 1 + let mut a3 = arr2(&[[1., 2., 3.]]); + a3.swap_axes(0, 1); + let mut b = OwnedArray::eye(a3.shape()[1]); + let mut res = arr2(&[[0., 0., 0.]]); + res.swap_axes(0, 1); + Gemm::gemm(&1., + Transpose::NoTrans, &a3.bvm(), + Transpose::NoTrans, &b.blas(), + &0., &mut res.bvm()); + assert_eq!(res, a3); +} + +#[test] +fn as_blas() { + let mut a = OwnedArray::::zeros((4, 4)); + assert!(a.blas_view_mut_checked().is_ok()); + a.swap_axes(0, 1); + assert!(a.blas_view_mut_checked().is_err()); + a.swap_axes(0, 1); + + { + // increased row stride + let mut b = a.slice_mut(s![..;2, ..]); + assert!(b.blas_view_mut_checked().is_ok()); + b.bvm(); // no panic + } + { + // inner dimension is not contig + let mut b = a.slice_mut(s![.., ..;2]); + assert!(b.blas_view_mut_checked().is_err()); + } + { + // inner dimension is length 1, is ok again + let mut b = a.slice_mut(s![.., ..;4]); + assert!(b.blas_view_mut_checked().is_ok()); + b.bvm(); + } +} + +type Ix2 = (Ix, Ix); +fn dot(a: ArrayView, b: ArrayView, + c: &mut ArrayViewMut) -> Result<(), ShapeError> + where F: Gemm + Float, +{ + let at = Transpose::NoTrans; + let bt = Transpose::NoTrans; + + let ba = try!(a.blas_view_checked()); + let bb = try!(b.blas_view_checked()); + let mut bc = try!(c.blas_view_mut_checked()); + F::gemm(&F::one(), + at, &ba, + bt, &bb, + &F::zero(), &mut bc); + Ok(()) +} + +#[test] +fn test_dot() { + let mut a = arr2(&[[1., 2.], + [0., 3.]]); + + let b = a.clone(); + let mut res = a.clone(); + res.assign_scalar(&0.); + + dot(a.view(), b.view(), &mut res.view_mut()).unwrap(); + println!("{:?}", res); + + a.swap_axes(0, 1); + res.assign_scalar(&0.); + + let result = dot(a.view(), b.view(), &mut res.view_mut()); + assert!(result.is_err()); +} diff --git a/src/blas.rs b/src/blas.rs index 25dd5f0d0..4ffabd15f 100644 --- a/src/blas.rs +++ b/src/blas.rs @@ -55,7 +55,8 @@ //! I know), instead output its own error conditions, for example on dimension //! mismatch in a matrix multiplication. //! -#![cfg_attr(has_deprecated, deprecated(note="`rblas` integration will move to its own crate."))] +#![cfg_attr(has_deprecated, deprecated(note="`rblas` integration has moved to crate `ndarray-rblas`, use it instead."))] +#![allow(deprecated)] use std::os::raw::{c_int}; diff --git a/src/error.rs b/src/error.rs index ed39d93f4..e5ff22b92 100644 --- a/src/error.rs +++ b/src/error.rs @@ -24,6 +24,11 @@ impl ShapeError { pub fn kind(&self) -> ErrorKind { self.repr } + + /// Create a new `ShapeError` + pub fn from_kind(error: ErrorKind) -> Self { + from_kind(error) + } } /// Error code for an error related to array shape or layout. diff --git a/src/impl_methods.rs b/src/impl_methods.rs index f09232f7d..fa9e00227 100644 --- a/src/impl_methods.rs +++ b/src/impl_methods.rs @@ -561,6 +561,19 @@ impl ArrayBase where S: Data, D: Dimension true } + #[inline(always)] + pub fn as_ptr(&self) -> *const A { + self.ptr + } + + #[inline(always)] + pub fn as_mut_ptr(&mut self) -> *mut A + where S: DataMut + { + self.ensure_unique(); // for RcArray + self.ptr + } + /// Return the array’s data as a slice, if it is contiguous and /// the element order corresponds to the memory order. Return `None` otherwise. pub fn as_slice(&self) -> Option<&[A]> { diff --git a/src/lib.rs b/src/lib.rs index 1b0be291f..33f97b07d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -47,7 +47,7 @@ //! + `.fold()` and `.zip_mut_with()` are the most efficient ways to //! perform single traversal and lock step traversal respectively. //! + `.iter()` and `.iter_mut()` are efficient for contiguous arrays. -//! - There is experimental bridging to the linear algebra package `rblas`. +//! + Can use BLAS in some operations (`dot` and `mat_mul`). //! //! ## Crate Feature Flags //! @@ -61,7 +61,7 @@ //! - Optional, stable //! - Enables serialization support //! - `rblas` -//! - Optional, stable +//! - ***Deprecated:*** replaced by separate crate `ndarray-rblas` //! - Enables `rblas` integration //! #![cfg_attr(all(feature = "assign_ops", not(has_assign)), diff --git a/tests/blas.rs b/tests/blas.rs index 7779e7190..7336e61b0 100644 --- a/tests/blas.rs +++ b/tests/blas.rs @@ -1,4 +1,5 @@ #![cfg(feature = "rblas")] +#![allow(deprecated)] extern crate rblas; extern crate num;