1 #![feature(portable_simd)]
3 use core_simd::simd::*;
5 fn a(i: usize, j: usize) -> f64 {
6 ((i + j) * (i + j + 1) / 2 + i + 1) as f64
9 fn mult_av(v: &[f64], out: &mut [f64]) {
10 assert!(v.len() == out.len());
11 assert!(v.len() % 2 == 0);
13 for (i, out) in out.iter_mut().enumerate() {
14 let mut sum = f64x2::splat(0.0);
18 let b = f64x2::from_slice(&v[j..]);
19 let a = f64x2::from_array([a(i, j), a(i, j + 1)]);
23 *out = sum.reduce_sum();
27 fn mult_atv(v: &[f64], out: &mut [f64]) {
28 assert!(v.len() == out.len());
29 assert!(v.len() % 2 == 0);
31 for (i, out) in out.iter_mut().enumerate() {
32 let mut sum = f64x2::splat(0.0);
36 let b = f64x2::from_slice(&v[j..]);
37 let a = f64x2::from_array([a(j, i), a(j + 1, i)]);
41 *out = sum.reduce_sum();
45 fn mult_atav(v: &[f64], out: &mut [f64], tmp: &mut [f64]) {
50 pub fn spectral_norm(n: usize) -> f64 {
51 assert!(n % 2 == 0, "only even lengths are accepted");
53 let mut u = vec![1.0; n];
54 let mut v = u.clone();
55 let mut tmp = u.clone();
58 mult_atav(&u, &mut v, &mut tmp);
59 mult_atav(&v, &mut u, &mut tmp);
61 (dot(&u, &v) / dot(&v, &v)).sqrt()
64 fn dot(x: &[f64], y: &[f64]) -> f64 {
65 // This is auto-vectorized:
66 x.iter().zip(y).map(|(&x, &y)| x * y).sum()
72 assert_eq!(format!("{:.9}", spectral_norm(100)), "1.274219991");
76 // Empty main to make cargo happy