Skip to content

Commit 6ecf638

Browse files
termoshtttermoshtt
committed
Merge branch 'master' into rcarray
2 parents 1e75543 + 7601cb5 commit 6ecf638

File tree

9 files changed

+214
-126
lines changed

9 files changed

+214
-126
lines changed

‎src/lib.rs

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ pub mod triangular;
4848
pub mod qr;
4949
pub mod svd;
5050
pub mod eigh;
51-
pub mod norm;
51+
pub mod opnorm;
5252
pub mod solve;
5353
pub mod cholesky;
54+
55+
pub mod util;

‎src/matrix.rs

Lines changed: 20 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -8,11 +8,11 @@ use lapack::c::Layout;
88
use error::{LinalgError, StrideError};
99
use qr::ImplQR;
1010
use svd::ImplSVD;
11-
use norm::ImplNorm;
11+
use opnorm::ImplOpNorm;
1212
use solve::ImplSolve;
1313

14-
pub trait MFloat: ImplQR + ImplSVD + ImplNorm + ImplSolve + NdFloat {}
15-
impl<A: ImplQR + ImplSVD + ImplNorm + ImplSolve + NdFloat> MFloat for A {}
14+
pub trait MFloat: ImplQR + ImplSVD + ImplOpNorm + ImplSolve + NdFloat {}
15+
impl<A: ImplQR + ImplSVD + ImplOpNorm + ImplSolve + NdFloat> MFloat for A {}
1616

1717
/// Methods for general matrices
1818
pub trait Matrix: Sized {
@@ -24,11 +24,11 @@ pub trait Matrix: Sized {
2424
/// Layout (C/Fortran) of matrix
2525
fn layout(&self) -> Result<Layout, StrideError>;
2626
/// Operator norm for L-1 norm
27-
fn norm_1(&self) -> Self::Scalar;
27+
fn opnorm_1(&self) -> Self::Scalar;
2828
/// Operator norm for L-inf norm
29-
fn norm_i(&self) -> Self::Scalar;
29+
fn opnorm_i(&self) -> Self::Scalar;
3030
/// Frobenius norm
31-
fn norm_f(&self) -> Self::Scalar;
31+
fn opnorm_f(&self) -> Self::Scalar;
3232
/// singular-value decomposition (SVD)
3333
fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError>;
3434
/// QR decomposition
@@ -84,27 +84,27 @@ impl<A: MFloat> Matrix for Array<A, Ix2> {
8484
fn layout(&self) -> Result<Layout, StrideError> {
8585
check_layout(self.strides())
8686
}
87-
fn norm_1(&self) -> Self::Scalar {
87+
fn opnorm_1(&self) -> Self::Scalar {
8888
let (m, n) = self.size();
8989
let strides = self.strides();
9090
if strides[0] > strides[1] {
91-
ImplNorm::norm_i(n, m, self.clone().into_raw_vec())
91+
ImplOpNorm::opnorm_i(n, m, self.clone().into_raw_vec())
9292
} else {
93-
ImplNorm::norm_1(m, n, self.clone().into_raw_vec())
93+
ImplOpNorm::opnorm_1(m, n, self.clone().into_raw_vec())
9494
}
9595
}
96-
fn norm_i(&self) -> Self::Scalar {
96+
fn opnorm_i(&self) -> Self::Scalar {
9797
let (m, n) = self.size();
9898
let strides = self.strides();
9999
if strides[0] > strides[1] {
100-
ImplNorm::norm_1(n, m, self.clone().into_raw_vec())
100+
ImplOpNorm::opnorm_1(n, m, self.clone().into_raw_vec())
101101
} else {
102-
ImplNorm::norm_i(m, n, self.clone().into_raw_vec())
102+
ImplOpNorm::opnorm_i(m, n, self.clone().into_raw_vec())
103103
}
104104
}
105-
fn norm_f(&self) -> Self::Scalar {
105+
fn opnorm_f(&self) -> Self::Scalar {
106106
let (m, n) = self.size();
107-
ImplNorm::norm_f(m, n, self.clone().into_raw_vec())
107+
ImplOpNorm::opnorm_f(m, n, self.clone().into_raw_vec())
108108
}
109109
fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
110110
let (n, m) = self.size();
@@ -192,17 +192,17 @@ impl<A: MFloat> Matrix for RcArray<A, Ix2> {
192192
fn layout(&self) -> Result<Layout, StrideError> {
193193
check_layout(self.strides())
194194
}
195-
fn norm_1(&self) -> Self::Scalar {
195+
fn opnorm_1(&self) -> Self::Scalar {
196196
// XXX unnecessary clone
197-
self.to_owned().norm_1()
197+
self.to_owned().opnorm_1()
198198
}
199-
fn norm_i(&self) -> Self::Scalar {
199+
fn opnorm_i(&self) -> Self::Scalar {
200200
// XXX unnecessary clone
201-
self.to_owned().norm_i()
201+
self.to_owned().opnorm_i()
202202
}
203-
fn norm_f(&self) -> Self::Scalar {
203+
fn opnorm_f(&self) -> Self::Scalar {
204204
// XXX unnecessary clone
205-
self.to_owned().norm_f()
205+
self.to_owned().opnorm_f()
206206
}
207207
fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
208208
// XXX unnecessary clone (should use into_owned())

‎src/norm.rs

Lines changed: 0 additions & 27 deletions
This file was deleted.

‎src/opnorm.rs

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
//! Implement Norms for matrices
2+
3+
use lapack::c::*;
4+
5+
pub trait ImplOpNorm: Sized {
6+
fn opnorm_1(m: usize, n: usize, a: Vec<Self>) -> Self;
7+
fn opnorm_i(m: usize, n: usize, a: Vec<Self>) -> Self;
8+
fn opnorm_f(m: usize, n: usize, a: Vec<Self>) -> Self;
9+
}
10+
11+
macro_rules! impl_opnorm {
12+
($scalar:ty, $lange:path) => {
13+
impl ImplOpNorm for $scalar {
14+
fn opnorm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
15+
$lange(Layout::ColumnMajor, b'o', m as i32, n as i32, &mut a, m as i32)
16+
}
17+
fn opnorm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
18+
$lange(Layout::ColumnMajor, b'i', m as i32, n as i32, &mut a, m as i32)
19+
}
20+
fn opnorm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
21+
$lange(Layout::ColumnMajor, b'f', m as i32, n as i32, &mut a, m as i32)
22+
}
23+
}
24+
}} // end macro_rules
25+
26+
impl_opnorm!(f64, dlange);
27+
impl_opnorm!(f32, slange);

‎src/util.rs

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
//! module for topologcal vector space
2+
//!
3+
4+
use std::iter::Sum;
5+
use ndarray::{ArrayBase, Data, Dimension, LinalgScalar};
6+
use num_traits::Float;
7+
use super::vector::*;
8+
9+
pub fn all_close_max<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>,
10+
truth: &ArrayBase<S2, D>,
11+
atol: Tol)
12+
-> Result<Tol, Tol>
13+
where A: LinalgScalar + Squared<Output = Tol>,
14+
Tol: Float + Sum,
15+
S1: Data<Elem = A>,
16+
S2: Data<Elem = A>,
17+
D: Dimension
18+
{
19+
let tol = (test - truth).norm_max();
20+
if tol < atol { Ok(tol) } else { Err(tol) }
21+
}
22+
23+
pub fn all_close_l1<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: Tol) -> Result<Tol, Tol>
24+
where A: LinalgScalar + Squared<Output = Tol>,
25+
Tol: Float + Sum,
26+
S1: Data<Elem = A>,
27+
S2: Data<Elem = A>,
28+
D: Dimension
29+
{
30+
let tol = (test - truth).norm_l1() / truth.norm_l1();
31+
if tol < rtol { Ok(tol) } else { Err(tol) }
32+
}
33+
34+
pub fn all_close_l2<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: Tol) -> Result<Tol, Tol>
35+
where A: LinalgScalar + Squared<Output = Tol>,
36+
Tol: Float + Sum,
37+
S1: Data<Elem = A>,
38+
S2: Data<Elem = A>,
39+
D: Dimension
40+
{
41+
let tol = (test - truth).norm_l2() / truth.norm_l2();
42+
if tol < rtol { Ok(tol) } else { Err(tol) }
43+
}

‎src/vector.rs

Lines changed: 43 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,24 +1,57 @@
11
//! Define trait for vectors
22
3-
use ndarray::{NdFloat, Array, RcArray, Ix1};
3+
use std::iter::Sum;
4+
use ndarray::{LinalgScalar, ArrayBase, Data, Dimension};
5+
use num_traits::float::Float;
46

57
/// Methods for vectors
68
pub trait Vector {
79
type Scalar;
10+
/// rename of norm_l2
11+
fn norm(&self) -> Self::Scalar {
12+
self.norm_l2()
13+
}
14+
/// L-1 norm
15+
fn norm_l1(&self) -> Self::Scalar;
816
/// L-2 norm
9-
fn norm(&self) -> Self::Scalar;
17+
fn norm_l2(&self) -> Self::Scalar;
18+
/// maximum norm
19+
fn norm_max(&self) -> Self::Scalar;
1020
}
1121

12-
impl<A: NdFloat> Vector for Array<A, Ix1> {
13-
type Scalar = A;
14-
fn norm(&self) -> Self::Scalar {
15-
self.dot(&self).sqrt()
22+
impl<A, S, D, T> Vector for ArrayBase<S, D>
23+
where A: LinalgScalar + Squared<Output = T>,
24+
T: Float + Sum,
25+
S: Data<Elem = A>,
26+
D: Dimension
27+
{
28+
type Scalar = T;
29+
fn norm_l1(&self) -> Self::Scalar {
30+
self.iter().map(|x| x.sq_abs()).sum()
31+
}
32+
fn norm_l2(&self) -> Self::Scalar {
33+
self.iter().map(|x| x.squared()).sum::<T>().sqrt()
34+
}
35+
fn norm_max(&self) -> Self::Scalar {
36+
self.iter().fold(T::zero(), |f, &val| {
37+
let v = val.sq_abs();
38+
if f > v { f } else { v }
39+
})
1640
}
1741
}
1842

19-
impl<A: NdFloat> Vector for RcArray<A, Ix1> {
20-
type Scalar = A;
21-
fn norm(&self) -> Self::Scalar {
22-
self.dot(&self).sqrt()
43+
pub trait Squared {
44+
type Output;
45+
fn squared(&self) -> Self::Output;
46+
fn sq_abs(&self) -> Self::Output;
47+
}
48+
49+
impl<A: Float> Squared for A {
50+
type Output = A;
51+
fn squared(&self) -> A {
52+
*self * *self
53+
}
54+
fn sq_abs(&self) -> A {
55+
self.abs()
2356
}
2457
}

‎tests/norm.rs

Lines changed: 0 additions & 68 deletions
This file was deleted.

‎tests/opnorm.rs

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
include!("header.rs");
2+
3+
#[test]
4+
fn matrix_opnorm_square() {
5+
let a = Array::range(1., 10., 1.).into_shape((3, 3)).unwrap();
6+
a.opnorm_1().assert_close(18.0, 1e-7);
7+
a.opnorm_i().assert_close(24.0, 1e-7);
8+
a.opnorm_f().assert_close(285.0.sqrt(), 1e-7);
9+
}
10+
11+
#[test]
12+
fn matrix_opnorm_square_t() {
13+
let a = Array::range(1., 10., 1.).into_shape((3, 3)).unwrap().reversed_axes();
14+
a.opnorm_1().assert_close(24.0, 1e-7);
15+
a.opnorm_i().assert_close(18.0, 1e-7);
16+
a.opnorm_f().assert_close(285.0.sqrt(), 1e-7);
17+
}
18+
19+
#[test]
20+
fn matrix_opnorm_3x4() {
21+
let a = Array::range(1., 13., 1.).into_shape((3, 4)).unwrap();
22+
a.opnorm_1().assert_close(24.0, 1e-7);
23+
a.opnorm_i().assert_close(42.0, 1e-7);
24+
a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
25+
}
26+
27+
#[test]
28+
fn matrix_opnorm_3x4_t() {
29+
let a = Array::range(1., 13., 1.)
30+
.into_shape((3, 4))
31+
.unwrap()
32+
.reversed_axes();
33+
a.opnorm_1().assert_close(42.0, 1e-7);
34+
a.opnorm_i().assert_close(24.0, 1e-7);
35+
a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
36+
}
37+
38+
#[test]
39+
fn matrix_opnorm_4x3() {
40+
let a = Array::range(1., 13., 1.).into_shape((4, 3)).unwrap();
41+
a.opnorm_1().assert_close(30.0, 1e-7);
42+
a.opnorm_i().assert_close(33.0, 1e-7);
43+
a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
44+
}
45+
46+
#[test]
47+
fn matrix_opnorm_4x3_t() {
48+
let a = Array::range(1., 13., 1.)
49+
.into_shape((4, 3))
50+
.unwrap()
51+
.reversed_axes();
52+
a.opnorm_1().assert_close(33.0, 1e-7);
53+
a.opnorm_i().assert_close(30.0, 1e-7);
54+
a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
55+
}

0 commit comments

Comments
 (0)