Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,6 @@
# Remove Cargo.lock from gitignore if creating an executable, leave it for libraries
# More information here http://doc.crates.io/guide.html#cargotoml-vs-cargolock
Cargo.lock

# cargo fmt
*.bk
4 changes: 2 additions & 2 deletions Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,9 @@ lapack = "0.11.1"
num-traits = "0.1.36"

[dependencies.ndarray]
version = "0.6.9"
version = "0.7"
features = ["blas"]

[dev-dependencies]
rand = "0.3.14"
ndarray-rand = "0.2.0"
ndarray-rand = "0.3"
1 change: 1 addition & 0 deletions rustfmt.toml
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
max_width = 120
use_try_shorthand = true
18 changes: 4 additions & 14 deletions src/eigh.rs
Original file line number Diff line number Diff line change
@@ -1,30 +1,20 @@
//! Implement eigenvalue decomposition of Hermite matrix

use lapack::fortran::*;
use lapack::c::*;
use num_traits::Zero;

use error::LapackError;

pub trait ImplEigh: Sized {
fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
fn eigh(layout: Layout, n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
}

macro_rules! impl_eigh {
($scalar:ty, $syev:path) => {
impl ImplEigh for $scalar {
fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
fn eigh(layout: Layout, n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
let mut w = vec![Self::zero(); n];
let mut work = vec![Self::zero(); 4 * n];
let mut info = 0;
$syev(b'V',
b'U',
n as i32,
&mut a,
n as i32,
&mut w,
&mut work,
4 * n as i32,
&mut info);
let info = $syev(layout, b'V', b'U', n as i32, &mut a, n as i32, &mut w);
if info == 0 {
Ok((w, a))
} else {
Expand Down
28 changes: 28 additions & 0 deletions src/error.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

use std::error;
use std::fmt;
use ndarray::Ixs;

#[derive(Debug)]
pub struct LapackError {
Expand Down Expand Up @@ -44,17 +45,37 @@ impl error::Error for NotSquareError {
}
}

#[derive(Debug)]
pub struct StrideError {
pub s0: Ixs,
pub s1: Ixs,
}

impl fmt::Display for StrideError {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "invalid stride: s0={}, s1={}", self.s0, self.s1)
}
}

impl error::Error for StrideError {
fn description(&self) -> &str {
"invalid stride"
}
}

#[derive(Debug)]
pub enum LinalgError {
NotSquare(NotSquareError),
Lapack(LapackError),
Stride(StrideError),
}

impl fmt::Display for LinalgError {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
match *self {
LinalgError::NotSquare(ref err) => err.fmt(f),
LinalgError::Lapack(ref err) => err.fmt(f),
LinalgError::Stride(ref err) => err.fmt(f),
}
}
}
Expand All @@ -64,6 +85,7 @@ impl error::Error for LinalgError {
match *self {
LinalgError::NotSquare(ref err) => err.description(),
LinalgError::Lapack(ref err) => err.description(),
LinalgError::Stride(ref err) => err.description(),
}
}
}
Expand All @@ -79,3 +101,9 @@ impl From<LapackError> for LinalgError {
LinalgError::Lapack(err)
}
}

impl From<StrideError> for LinalgError {
fn from(err: StrideError) -> LinalgError {
LinalgError::Stride(err)
}
}
24 changes: 13 additions & 11 deletions src/hermite.rs
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
//! Define trait for Hermite matrices

use ndarray::{Ix2, Array, LinalgScalar};
use std::fmt::Debug;
use ndarray::prelude::*;
use ndarray::LinalgScalar;
use num_traits::float::Float;
use lapack::c::Layout;

Expand All @@ -26,20 +25,24 @@ pub trait HermiteMatrix: SquareMatrix + Matrix {
fn cholesky(self) -> Result<Self, LinalgError>;
}

impl<A> HermiteMatrix for Array<A, (Ix, Ix)>
impl<A> HermiteMatrix for Array<A, Ix2>
where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug
{
fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
try!(self.check_square());
self.check_square()?;
let layout = self.layout()?;
let (rows, cols) = self.size();
let (w, a) = try!(ImplEigh::eigh(rows, self.into_raw_vec()));
let (w, a) = ImplEigh::eigh(layout, rows, self.into_raw_vec())?;
let ea = Array::from_vec(w);
let va = Array::from_vec(a).into_shape((rows, cols)).unwrap().reversed_axes();
let va = match layout {
Layout::ColumnMajor => Array::from_vec(a).into_shape((rows, cols)).unwrap().reversed_axes(),
Layout::RowMajor => Array::from_vec(a).into_shape((rows, cols)).unwrap(),
};
Ok((ea, va))
}
fn ssqrt(self) -> Result<Self, LinalgError> {
let (n, _) = self.size();
let (e, v) = try!(self.eigh());
let (e, v) = self.eigh()?;
let mut res = Array::zeros((n, n));
for i in 0..n {
for j in 0..n {
Expand All @@ -49,11 +52,10 @@ impl<A> HermiteMatrix for Array<A, (Ix, Ix)>
Ok(v.dot(&res))
}
fn cholesky(self) -> Result<Self, LinalgError> {
try!(self.check_square());
println!("layout = {:?}", self.layout());
self.check_square()?;
let (n, _) = self.size();
let layout = self.layout();
let a = try!(ImplCholesky::cholesky(layout, n, self.into_raw_vec()));
let layout = self.layout()?;
let a = ImplCholesky::cholesky(layout, n, self.into_raw_vec())?;
let mut c = match layout {
Layout::RowMajor => Array::from_vec(a).into_shape((n, n)).unwrap(),
Layout::ColumnMajor => Array::from_vec(a).into_shape((n, n)).unwrap().reversed_axes(),
Expand Down
81 changes: 33 additions & 48 deletions src/matrix.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ use ndarray::prelude::*;
use ndarray::LinalgScalar;
use lapack::c::Layout;

use error::LapackError;
use error::{LinalgError, StrideError};
use qr::ImplQR;
use svd::ImplSVD;
use norm::ImplNorm;
Expand All @@ -20,19 +20,19 @@ pub trait Matrix: Sized {
/// number of (rows, columns)
fn size(&self) -> (usize, usize);
/// Layout (C/Fortran) of matrix
fn layout(&self) -> Layout;
fn layout(&self) -> Result<Layout, StrideError>;
/// Operator norm for L-1 norm
fn norm_1(&self) -> Self::Scalar;
/// Operator norm for L-inf norm
fn norm_i(&self) -> Self::Scalar;
/// Frobenius norm
fn norm_f(&self) -> Self::Scalar;
/// singular-value decomposition (SVD)
fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError>;
fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError>;
/// QR decomposition
fn qr(self) -> Result<(Self, Self), LapackError>;
fn qr(self) -> Result<(Self, Self), LinalgError>;
/// LU decomposition
fn lu(self) -> Result<(Self::Permutator, Self, Self), LapackError>;
fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError>;
/// permutate matrix (inplace)
fn permutate(&mut self, p: &Self::Permutator);
/// permutate matrix (outplace)
Expand All @@ -42,22 +42,28 @@ pub trait Matrix: Sized {
}
}

impl<A> Matrix for Array<A, (Ix, Ix)>
impl<A> Matrix for Array<A, Ix2>
where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + LinalgScalar + Debug
{
type Scalar = A;
type Vector = Array<A, Ix>;
type Vector = Array<A, Ix1>;
type Permutator = Vec<i32>;

fn size(&self) -> (usize, usize) {
(self.rows(), self.cols())
}
fn layout(&self) -> Layout {
fn layout(&self) -> Result<Layout, StrideError> {
let strides = self.strides();
if min(strides[0], strides[1]) != 1 {
return Err(StrideError {
s0: strides[0],
s1: strides[1],
});;
}
if strides[0] < strides[1] {
Layout::ColumnMajor
Ok(Layout::ColumnMajor)
} else {
Layout::RowMajor
Ok(Layout::RowMajor)
}
}
fn norm_1(&self) -> Self::Scalar {
Expand All @@ -82,35 +88,24 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
let (m, n) = self.size();
ImplNorm::norm_f(m, n, self.clone().into_raw_vec())
}
fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError> {
let strides = self.strides();
let (m, n) = if strides[0] > strides[1] {
self.size()
} else {
let (n, m) = self.size();
(m, n)
};
let (u, s, vt) = try!(ImplSVD::svd(m, n, self.clone().into_raw_vec()));
fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
let (n, m) = self.size();
let layout = self.layout()?;
let (u, s, vt) = ImplSVD::svd(layout, m, n, self.clone().into_raw_vec())?;
let sv = Array::from_vec(s);
if strides[0] > strides[1] {
let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
Ok((va, sv, ua))
} else {
let ua = Array::from_vec(u).into_shape((n, n)).unwrap().reversed_axes();
let va = Array::from_vec(vt).into_shape((m, m)).unwrap().reversed_axes();
Ok((ua, sv, va))
let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
match layout {
Layout::RowMajor => Ok((ua, sv, va)),
Layout::ColumnMajor => Ok((ua.reversed_axes(), sv, va.reversed_axes())),
}
}
fn qr(self) -> Result<(Self, Self), LapackError> {
fn qr(self) -> Result<(Self, Self), LinalgError> {
let (n, m) = self.size();
let strides = self.strides();
let k = min(n, m);
let (q, r) = if strides[0] < strides[1] {
try!(ImplQR::qr(m, n, self.clone().into_raw_vec()))
} else {
try!(ImplQR::lq(n, m, self.clone().into_raw_vec()))
};
let layout = self.layout()?;
let (q, r) = ImplQR::qr(layout, m, n, self.clone().into_raw_vec())?;
let (qa, ra) = if strides[0] < strides[1] {
(Array::from_vec(q).into_shape((m, n)).unwrap().reversed_axes(),
Array::from_vec(r).into_shape((m, n)).unwrap().reversed_axes())
Expand All @@ -136,23 +131,14 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
}
Ok((qm, rm))
}
fn lu(self) -> Result<(Self::Permutator, Self, Self), LapackError> {
fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError> {
let (n, m) = self.size();
println!("n={}, m={}", n, m);
let k = min(n, m);
let (p, mut a) = match self.layout() {
Layout::ColumnMajor => {
println!("ColumnMajor");
let (p, l) = ImplSolve::lu(self.layout(), n, m, self.clone().into_raw_vec())?;
(p, Array::from_vec(l).into_shape((m, n)).unwrap().reversed_axes())
}
Layout::RowMajor => {
println!("RowMajor");
let (p, l) = ImplSolve::lu(self.layout(), n, m, self.clone().into_raw_vec())?;
(p, Array::from_vec(l).into_shape((n, m)).unwrap())
}
let (p, l) = ImplSolve::lu(self.layout()?, n, m, self.clone().into_raw_vec())?;
let mut a = match self.layout()? {
Layout::ColumnMajor => Array::from_vec(l).into_shape((m, n)).unwrap().reversed_axes(),
Layout::RowMajor => Array::from_vec(l).into_shape((n, m)).unwrap(),
};
println!("a (after LU) = \n{:?}", &a);
let mut lm = Array::zeros((n, k));
for ((i, j), val) in lm.indexed_iter_mut() {
if i > j {
Expand All @@ -171,7 +157,6 @@ impl<A> Matrix for Array<A, (Ix, Ix)>
} else {
a
};
println!("am = \n{:?}", am);
Ok((p, lm, am))
}
fn permutate(&mut self, ipiv: &Self::Permutator) {
Expand Down
12 changes: 4 additions & 8 deletions src/norm.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
//! Implement Norms for matrices

use lapack::fortran::*;
use num_traits::Zero;
use lapack::c::*;

pub trait ImplNorm: Sized {
fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self;
Expand All @@ -13,16 +12,13 @@ macro_rules! impl_norm {
($scalar:ty, $lange:path) => {
impl ImplNorm for $scalar {
fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
let mut work = Vec::<Self>::new();
$lange(b'o', m as i32, n as i32, &mut a, m as i32, &mut work)
$lange(Layout::ColumnMajor, b'o', m as i32, n as i32, &mut a, m as i32)
}
fn norm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
let mut work = vec![Self::zero(); m];
$lange(b'i', m as i32, n as i32, &mut a, m as i32, &mut work)
$lange(Layout::ColumnMajor, b'i', m as i32, n as i32, &mut a, m as i32)
}
fn norm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
let mut work = Vec::<Self>::new();
$lange(b'f', m as i32, n as i32, &mut a, m as i32, &mut work)
$lange(Layout::ColumnMajor, b'f', m as i32, n as i32, &mut a, m as i32)
}
}
}} // end macro_rules
Expand Down
Loading