Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor(demo): rework Laplace demo #653

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions examples/demo/Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ authors = ["Stefan Lankes <slankes@eonerc.rwth-aachen.de>"]
edition = "2021"

[dependencies]
ndarray = { version = "0.16", features = ["rayon"] }
rayon = "1.5"

[target.'cfg(target_os = "hermit")'.dependencies]
Expand Down
138 changes: 67 additions & 71 deletions examples/demo/src/laplace.rs
Original file line number Diff line number Diff line change
@@ -1,108 +1,104 @@
//! Jacobi Stencil Iterations
//!
//! This module performs the Jacobi method for solving Laplace's differential equation.

#![allow(clippy::reversed_empty_ranges)]

use std::mem;
use std::time::Instant;
use std::vec;

use ndarray::{s, Array2, ArrayView2, ArrayViewMut2};
use rayon::prelude::*;

const SIZE: usize = if cfg!(debug_assertions) { 16 } else { 64 };
const ITERATIONS: usize = 1000;

pub fn laplace() {
eprintln!();

let matrix = matrix_setup(SIZE, SIZE);
let mut matrix = matrix_setup(SIZE, SIZE);

eprintln!("Laplace iterations");
let now = Instant::now();
let (i, residual) = compute(matrix, SIZE, SIZE);
let residual = compute(matrix.view_mut(), ITERATIONS);
let elapsed = now.elapsed();
eprintln!("{i} iterations: {elapsed:?} (residual: {residual})");
eprintln!("{ITERATIONS} iterations: {elapsed:?} (residual: {residual})");

assert!(residual < 0.001);
}

fn matrix_setup(size_x: usize, size_y: usize) -> vec::Vec<vec::Vec<f64>> {
let mut matrix = vec![vec![0.0; size_x * size_y]; 2];

// top row
for x in 0..size_x {
matrix[0][x] = 1.0;
matrix[1][x] = 1.0;
}
fn matrix_setup(size_x: usize, size_y: usize) -> Array2<f64> {
let mut matrix = Array2::zeros((size_x, size_y));

// bottom row
for x in 0..size_x {
matrix[0][(size_y - 1) * size_x + x] = 1.0;
matrix[1][(size_y - 1) * size_x + x] = 1.0;
}

// left row
for y in 0..size_y {
matrix[0][y * size_x] = 1.0;
matrix[1][y * size_x] = 1.0;
}

// right row
for y in 0..size_y {
matrix[0][y * size_x + size_x - 1] = 1.0;
matrix[1][y * size_x + size_x - 1] = 1.0;
}
matrix.row_mut(0).fill(1.0);
matrix.row_mut(size_x - 1).fill(1.0);
matrix.column_mut(0).fill(1.0);
matrix.column_mut(size_x - 1).fill(1.0);

matrix
}

fn get_residual(matrix: &[f64], size_x: usize, size_y: usize) -> f64 {
(1..size_y - 1)
fn get_residual(matrix: ArrayView2<f64>) -> f64 {
matrix
.slice(s![1..-1, ..])
.outer_iter()
.into_par_iter()
.map(|y| {
let mut local_sum = 0.0;

for x in 1..(size_x - 1) {
let new = (matrix[y * size_x + x - 1]
+ matrix[y * size_x + x + 1]
+ matrix[(y + 1) * size_x + x]
+ matrix[(y - 1) * size_x + x])
* 0.25;

let diff = new - matrix[y * size_x + x];
local_sum += diff * diff;
.enumerate()
.map(|(i, row)| {
let i = i + 1; // To compensate slicing

let up = matrix.row(i - 1);
let here = matrix.row(i);
let down = matrix.row(i + 1);
let len = row.len();
assert_eq!(up.len(), len);
assert_eq!(here.len(), len);
assert_eq!(down.len(), len);

let mut acc = 0.0;
for j in 1..len - 1 {
let sum = up[j] + down[j] + here[j - 1] + here[j + 1];
let new = sum * 0.25;
acc += (new - here[j]).powi(2);
}

local_sum
acc
})
.sum()
}

fn iteration(cur: &[f64], next: &mut [f64], size_x: usize, size_y: usize) {
next.par_chunks_mut(size_y)
.enumerate() // to figure out where this chunk came from
.for_each(|(chunk_index, slice)| {
if chunk_index > 0 && chunk_index < size_y - 1 {
let offset_base = chunk_index * size_x;

for x in 1..size_x - 1 {
slice[x] = (cur[offset_base + x - 1]
+ cur[offset_base + x + 1]
+ cur[offset_base + size_x + x]
+ cur[offset_base - size_x + x])
* 0.25;
}
fn iteration(current: ArrayView2<f64>, mut next: ArrayViewMut2<f64>) {
next.slice_mut(s![1..-1, ..])
.outer_iter_mut()
.into_par_iter()
.enumerate()
.for_each(|(i, mut row)| {
let i = i + 1; // To compensate slicing

let up = current.row(i - 1);
let here = current.row(i);
let down = current.row(i + 1);
let len = row.len();
assert_eq!(up.len(), len);
assert_eq!(here.len(), len);
assert_eq!(down.len(), len);

for j in 1..len - 1 {
let sum = up[j] + down[j] + here[j - 1] + here[j + 1];
row[j] = sum * 0.25;
}
});
}

pub fn compute(mut matrix: vec::Vec<vec::Vec<f64>>, size_x: usize, size_y: usize) -> (usize, f64) {
let mut counter = 0;

while counter < 1000 {
{
// allow a borrow and a reference to the same vector
let (current, next) = matrix.split_at_mut(1);
fn compute(mut matrix: ArrayViewMut2<'_, f64>, iterations: usize) -> f64 {
let mut owned = matrix.to_owned();

iteration(&current[0], &mut next[0], size_x, size_y);
}
matrix.swap(0, 1);
let mut current = matrix.view_mut();
let mut next = owned.view_mut();

counter += 1;
for _ in 0..iterations {
iteration(current.view(), next.view_mut());
mem::swap(&mut current, &mut next);
}

(counter, get_residual(&matrix[0], size_x, size_y))
get_residual(current.view())
}
Loading