1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
use rand::Rng;
pub use crate::vector::Vector;
impl Vector<f64> {
/// Create a linearly spaced vector of f64 with n elements with lower limit a
/// and upper limit b
#[inline]
pub fn linspace( a: f64, b: f64, size: usize ) -> Self {
let mut vec = vec![ 0.0; size ];
let h: f64 = ( b - a ) / ((size as f64) - 1.0);
for i in 0..size {
vec[i] = a + h * (i as f64);
}
Vector{ vec }
}
/// Create a nonuniform vector of f64 with n elements with lower limit a,
/// upper limit b and exponent p (p=1 -> linear)
#[inline]
pub fn powspace( a: f64, b: f64, size: usize, p: f64 ) -> Self {
let mut vec = vec![ 0.0; size ];
for i in 0..size {
vec[i] = a + (b - a) * f64::powf( (i as f64) / ((size as f64) - 1.0), p );
}
Vector{ vec }
}
/// Return the L2 norm: square root of the sum of the squares
#[inline]
pub fn norm_2(&self) -> f64 {
let mut result = 0.0;
for i in 0..self.size() {
result += f64::powf( self.vec[i].abs(), 2.0 );
}
f64::sqrt( result )
}
/// Return the Lp norm: p-th root of the sum of the absolute values
/// raised to the power p
#[inline]
pub fn norm_p(&self, p: f64 ) -> f64 {
let mut result = 0.0;
for i in 0..self.size() {
result += f64::powf( self.vec[i].abs(), p );
}
f64::powf( result, 1.0/p )
}
/// Return the Inf norm: largest absolute value element (p -> infinity)
#[inline]
pub fn norm_inf(&self) -> f64 {
let mut result = self.vec[0].abs();
for i in 1..self.size() {
if result < self.vec[i].abs() {
result = self.vec[i].abs();
}
}
result
}
/// Create a vector containing n random elements between 0 and 1
#[inline]
pub fn random( size: usize ) -> Self {
let mut vec = vec![ 0.0; size ];
let mut rng = rand::thread_rng();
for i in 0..size {
vec[i] = rng.gen::<f64>()
}
Vector{ vec }
}
/// Return the dot product of two vectors v.dot(w) of f64s
#[inline]
pub fn dot_f64(&self, w: &Vector<f64>) -> f64 {
if self.size() != w.size() { panic!( "Vector sizes do not agree dot()." ); }
let num_threads = num_cpus::get();
/*if num_threads < self.size() || num_threads == 1 {
let mut result: f64 = 0.0;
for i in 0..self.size() {
result += self.vec[i] * w.vec[i];
}
return result;
}*/
let chunk_size = self.size() / num_threads;
std::thread::scope(|s| {
let mut threads = Vec::new();
for i in 0..num_threads {
let start = i * chunk_size;
let end = if i == num_threads - 1 { self.size() } else { (i + 1) * chunk_size };
let self_slice = &self.vec[start..end];
let w_slice = &w.vec[start..end];
threads.push( s.spawn(|| {
let mut result: f64 = 0.0;
for i in 0..self_slice.len() {
result += self_slice[i] * w_slice[i];
}
result
}));
}
let mut result: f64 = 0.0;
for thread in threads {
result += thread.join().unwrap();
}
result
})
}
}