Правильно ли я пишу на Rust?

$ clang-8 --version
clang version 8.0.0-3~ubuntu18.04.1 (tags/RELEASE_800/final)
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin
$ clang++-8 test9.cpp -march=native -Ofast -std=gnu++2a  -stdlib=libc++ -funroll-loops -o test9_cpp
$ ./test9_cpp 3000
13.4062
-287.25
$ clang++-8 test9.cpp -march=native -Ofast -std=gnu++2a  -stdlib=libc++ -funroll-loops -o test9_cpp -mno-fma
$ ./test9_cpp 3000
13.3594
-287.25
$ ./test9_rust 3000
13.0280102s
-287.2500833333321
$ ./test9_rust 3000
12.9767008s
-287.2500833333321
$ ./test9_cpp 3000
13.5312
-287.25

$ rustup default nightly-2019-05-23
info: using existing install for 'nightly-2019-05-23-x86_64-unknown-linux-gnu'
info: default toolchain set to 'nightly-2019-05-23-x86_64-unknown-linux-gnu'

  nightly-2019-05-23-x86_64-unknown-linux-gnu unchanged - rustc 1.36.0-nightly (37ff5d388 2019-05-22)

$ rustc --version --verbose
rustc 1.36.0-nightly (37ff5d388 2019-05-22)
binary: rustc
commit-hash: 37ff5d388f8c004ca248adb635f1cc84d347eda0
commit-date: 2019-05-22
host: x86_64-unknown-linux-gnu
release: 1.36.0-nightly
LLVM version: 8.0
$ rustc src/main.rs -C opt-level=3 -C target-cpu=native -o test9_rust
warning: function is never used: `matrix_print`
  --> src/main.rs:57:1
   |
57 | fn matrix_print(n: usize, a: &Vec<f64>) {
   | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   |
   = note: #[warn(dead_code)] on by default

$ ./test9_rust 3000
13.0152259s
-287.2500833333321
$ ./test9_rust 3000
12.9652781s
-287.2500833333321
$ ./test9_cpp 3000
13.3125
-287.25
$ ./test9_cpp 3000
13.3281
-287.25

$ ./cpp 3000
3000
10.873699
-287.250083
$ ./cpp 3000
3000
10.880284
-287.250083
$ ./main_rust 3000
11.063109248s
-287.2500833333321
$ ./main_rust 3000
11.014244154s
-287.2500833333321
$ ./main_rust 3000
11.051293901s
-287.2500833333321

rustc 1.38.0-nightly 
clang version 10.0.0 

#include <vector>
#include <ctime>
#include <valarray>
#include <cstdlib>
#include <cstdio>
#include <numeric>

using matrix_t = std::vector<std::valarray<double>>;
auto matrix = [](size_t n) { return matrix_t{n, std::valarray<double>(n)}; };


void fill(matrix_t & m) {
  auto n = m.size();
  auto tmp = 1. / n / n;
  for(ssize_t i = 0; i < ssize_t(n); ++i)
    for(ssize_t j = 0; j < ssize_t(n); ++j)
      m[i][j] = tmp * (i - j) * (i + j);
}

[[gnu::always_inline]] inline void mult(const matrix_t & a, const matrix_t & b, matrix_t & c) {
  auto n = a.size();
  auto tmp = matrix(n);
  for(size_t i = 0; i < n; ++i) {
    for(size_t j = 0; j < n; ++j)
      tmp[i][j] = b[j][i];
  }
  
  for(size_t i = 0; i < n; ++i) {
    for(size_t j = 0; j < n; ++j) {
//       c[i][j] = [](auto & a, auto & b) {
//         double acc = 0;
//         for(size_t n = 0; n < a.size(); ++n) acc += a[n] * b[n];
//         return acc;
//       }(a[i], tmp[j]);
      c[i][j] = std::inner_product(begin(a[i]), end(a[i]), begin(tmp[j]), 0.);
    }
  }
}


int main(int argc, char * argv[]) {
  size_t N = 3000;
  if(argc > 1) N = std::atol(argv[1]);

  fprintf(stderr, "%lu\n", N);

  auto a = matrix(N), b = matrix(N), c = matrix(N);

  fill(a);
  fill(b);
  
  clock_t t1 = clock();
  mult(a, b, c);
  clock_t t2 = clock();
  fprintf(stderr, "%f\n", (double)(t2 - t1) / CLOCKS_PER_SEC);
  fprintf(stderr, "%f\n", c[N/2][N/2]);
} 

$ clang++ cpp_v3_anon.cpp -march=native -Ofast -std=gnu++2a  -stdlib=libc++ -funroll-loops 
cpp_v3_anon.cpp:1:10: fatal error: 'iostream' file not found
#include <iostream>
         ^~~~~~~~~~
1 error generated.

$ clang++ --version
clang version 8.0.0 (Fedora 8.0.0-1.fc30)
Target: x86_64-unknown-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

$ ./cpp 3000
3000
20.4963
-287.25

$ ./c 3000
19.255953
-287.250083

$ ./rust_v1 3000
18.741719772s
-287.2500833333321

# Cargo.toml
[dependencies]
rayon = "1.1.0"

fn matrix_mult(n: usize, a: &[f64], b: &[f64], c: &mut [f64]) {
    use rayon::prelude::*;

    let mut t = vec![0.0; n * n];

    t.par_chunks_mut(n).enumerate().for_each(|(i, t_row)| {
        t_row
            .iter_mut()
            .enumerate()
            .for_each(|(j, t)| *t = b[j * n + i]);
    });

    c.par_chunks_mut(n)
        .zip(a.par_chunks(n))
        .for_each(|(c_row, a_row)| {
            t.par_chunks(n).zip(c_row).for_each(|(t_row, c)| {
                *c = a_row.iter().zip(t_row).fold(*c, |c, (&a, &t)| c + a * t);
            });
        });
}

fn matrix_mult(n: usize, a: &[f64], b: &[f64], c: &mut [f64]) |	fn matrix_mult_par(n: usize, a: &[f64], b: &[f64], c: &mut [f
    let mut t = vec![0.0; n * n];				    let mut t = vec![0.0; n * n];
    matrix_transpose(n, b, &mut t);				    matrix_transpose(n, b, &mut t);

    let c_rows = c.chunks_exact_mut(n);			      |	    let c_rows = c.par_chunks_mut(n);
    let a_rows = a.chunks_exact(n);			      |	    let a_rows = a.par_chunks(n);

    c_rows.zip(a_rows).for_each(|(c_row, a_row)| {		    c_rows.zip(a_rows).for_each(|(c_row, a_row)| {
        t.chunks_exact(n).zip(c_row).for_each(|(t_col, c)| {  |	        t.par_chunks(n).zip(c_row).for_each(|(t_col, c)| {
            *c = product_sum(*c, a_row, t_col);			            *c = product_sum(*c, a_row, t_col);
        });							        });
    });								    });
}								}

fn matrix_transpose(n: usize, src: &[f64], dst: &mut [f64]) {
    for (i, dst_row) in dst.chunks_exact_mut(n).enumerate() {
        dst_row
            .iter_mut()
            .enumerate()
            .for_each(|(j, t)| *t = src[j * n + i]);
    }
}

fn product_sum(init: f64, x: &[f64], y: &[f64]) -> f64 {
    x.iter()
        .zip(y)
        .fold(init, |acc, (&x, &y)| self::impls::fma(acc, x, y))
}

#[cfg(not(feature = "unstable"))]
mod impls {
    pub fn fma(acc: f64, x: f64, y: f64) -> f64 {
        acc + x * y
    }
}

#[cfg(feature = "unstable")]
mod impls {
    use std::intrinsics::{fadd_fast, fmul_fast};

    pub fn fma(acc: f64, x: f64, y: f64) -> f64 {
        unsafe { fadd_fast(acc, fmul_fast(x, y)) }
    }
}

func multiply(x, y [][]float32) ([][]float32, error) {
    if len(x[0]) != len(y) {
	return nil, errors.New("Can't do matrix multiplication.")
    }
	
    out := make([][]float32, len(x))
    for i := 0; i < len(x); i += 1 {
	for j := 0; j < len(y); j += 1 {
            if len(out[i]) < 1 {
	        out[i] = make([]float32, len(y))      			 
            }
	    out[i][j] += x[i][j] * y[j][i]
	}
    }

    return out, nil
}

class Matrix {
public:
  Matrix(size_t rows, size_t cols)
      : m_rows(rows), m_cols(cols), m_data(rows * cols)
  {
  }

  double& operator()(size_t row, size_t col)
  {
    return m_data[row * m_cols + col];
  };

  const double& operator()(size_t row, size_t col) const
  {
    return m_data[row * m_cols + col];
  };

  size_t rows() const
  {
    return m_rows;
  }

  size_t cols() const
  {
    return m_cols;
  }

private:
  size_t              m_rows;
  size_t              m_cols;
  std::vector<double> m_data;
};

c'' <- let c' = matrixMult n a b c in c' `deepseq` pure c'

type Matrix = [[Double]]

  time = 62ns
     n = 3000
result = -287.25008333333216

use std::{env, time::Instant};

fn main() {
    let args: Vec<_> = env::args().collect();
    let n = args[1].parse().unwrap();

    let a = &matrix_new(n);
    let b = &matrix_new(n);

    let t1 = Instant::now();
    let mut c = (0..n).map(move |_| {
        (0..n).map(move |i| (0..n).map(|j| a[i * n + j] * b[j * n + i]).sum::<f64>())
    });

    let t2 = Instant::now();

    println!("  time = {:?}", t2.duration_since(t1));
    println!("     n = {}", n);
    println!("result = {}", c.nth(n / 2).unwrap().nth(n / 2).unwrap());
}

fn matrix_new(n: usize) -> Vec<f64> {
    let tmp = 1.0 / n as f64 / n as f64;
    (0..n)
        .flat_map(|i| (0..n).map(move |j| tmp * (i as isize - j as isize) as f64 * (i + j) as f64))
        .collect()
}

use std::{env, time::Instant};

fn main() {
    let args: Vec<_> = env::args().collect();
    let n = args[1].parse().unwrap();

    let tmp = 1.0 / n as f64 / n as f64;
    let a = (0..n)
        .map(|i| (0..n).map(move |j| tmp * (i as isize - j as isize) as f64 * (i + j) as f64));
    let b = (0..n)
        .map(|i| (0..n).map(move |j| tmp * (j as isize - i as isize) as f64 * (i + j) as f64));

    let t1 = Instant::now();
    let mut c = a.map(move |row| {
        b.clone().map(move |col| {
            let row = row.clone();
            move || col.zip(row).map(|(b, a)| a * b).sum::<f64>()
        })
    });
    let t2 = Instant::now();

    println!("  time = {:?}", t2.duration_since(t1));
    println!("     n = {}", n);
    println!("result = {}", c.nth(n / 2).unwrap().nth(n / 2).unwrap()());
}

func multiply(x, y [][]float32) ([][]float32, error) {
    if len(x[0]) != len(y) {
	return nil, errors.New("Can't do matrix multiplication.")
    }
	
    out := make([][]float32, len(x))
    for i := range x {
	out[i] = make([]float32, len(y))    
	for j := range y {
	    out[i][j] += x[i][j] * y[j][i]
	}
    }

    return out, nil
}