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
|
#![feature(portable_simd)]
use core_simd::simd::*;
fn a(i: usize, j: usize) -> f64 {
((i + j) * (i + j + 1) / 2 + i + 1) as f64
}
fn mult_av(v: &[f64], out: &mut [f64]) {
assert!(v.len() == out.len());
assert!(v.len() % 2 == 0);
for (i, out) in out.iter_mut().enumerate() {
let mut sum = f64x2::splat(0.0);
let mut j = 0;
while j < v.len() {
let b = f64x2::from_slice(&v[j..]);
let a = f64x2::from_array([a(i, j), a(i, j + 1)]);
sum += b / a;
j += 2
}
*out = sum.reduce_sum();
}
}
fn mult_atv(v: &[f64], out: &mut [f64]) {
assert!(v.len() == out.len());
assert!(v.len() % 2 == 0);
for (i, out) in out.iter_mut().enumerate() {
let mut sum = f64x2::splat(0.0);
let mut j = 0;
while j < v.len() {
let b = f64x2::from_slice(&v[j..]);
let a = f64x2::from_array([a(j, i), a(j + 1, i)]);
sum += b / a;
j += 2
}
*out = sum.reduce_sum();
}
}
fn mult_atav(v: &[f64], out: &mut [f64], tmp: &mut [f64]) {
mult_av(v, tmp);
mult_atv(tmp, out);
}
pub fn spectral_norm(n: usize) -> f64 {
assert!(n % 2 == 0, "only even lengths are accepted");
let mut u = vec![1.0; n];
let mut v = u.clone();
let mut tmp = u.clone();
for _ in 0..10 {
mult_atav(&u, &mut v, &mut tmp);
mult_atav(&v, &mut u, &mut tmp);
}
(dot(&u, &v) / dot(&v, &v)).sqrt()
}
fn dot(x: &[f64], y: &[f64]) -> f64 {
// This is auto-vectorized:
x.iter().zip(y).map(|(&x, &y)| x * y).sum()
}
#[cfg(test)]
#[test]
fn test() {
assert_eq!(format!("{:.9}", spectral_norm(100)), "1.274219991");
}
fn main() {
// Empty main to make cargo happy
}
|
