Replace Array2D with rlst matrix (#123)

* Use faster kernel from bempp_kernel

* cargo fmt

* Add assembly benchmark

* directly transpose points

* fmt

* benchmark singular and nonsingular parts separately

* clippy

* remove timing example (as its now replaced by benchmark)

* const generics

* working on moving to using rlst storage formats

* fix tests

* only tabulate once

* get element working without Array2D

* make grid use rlst formats

* Remove Array2D, use rlst instead

* adjust dependency test

* &[] not &vec![]

* fmt

* more vec!s

* more vec!s

* update parallel grid

* Laplace

bem/Cargo.toml

@@ -24,9 +24,22 @@ bempp-tools = { path = "../tools"}
 bempp-traits = { path = "../traits"}
 bempp-element = { path = "../element"}
 bempp-grid = { path = "../grid"}
+bempp-kernel = { path = "../kernel"}
 bempp-quadrature = { path = "../quadrature"}
 approx = "0.5"
 itertools = "0.10"
 mpi = { version = "0.6.*", optional = true }
 num = "0.4"
 rayon = "1.7"
+rlst = { git = "https://github.com/linalg-rs/rlst.git" }
+rlst-blis-src = { git = "https://github.com/linalg-rs/rlst.git" }
+rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" }
+rlst-algorithms = { git = "https://github.com/linalg-rs/rlst.git" }
+
+[dev-dependencies]
+criterion = { version = "0.3", features = ["html_reports"]}
+
+[[bench]]
+name = "assembly_benchmark"
+harness = false
+

bem/benches/assembly_benchmark.rs

@@ -0,0 +1,121 @@
+use bempp_bem::assembly::{assemble_batched, batched, BoundaryOperator, PDEType};
+use bempp_bem::function_space::SerialFunctionSpace;
+use bempp_element::element::create_element;
+use bempp_grid::shapes::regular_sphere;
+use bempp_kernel::laplace_3d;
+use bempp_tools::arrays::zero_matrix;
+use bempp_traits::bem::DofMap;
+use bempp_traits::bem::FunctionSpace;
+use bempp_traits::cell::ReferenceCellType;
+use bempp_traits::element::{Continuity, ElementFamily};
+use criterion::{criterion_group, criterion_main, Criterion};
+
+pub fn full_assembly_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("assembly");
+    group.sample_size(20);
+
+    for i in 3..5 {
+        let grid = regular_sphere(i);
+        let element = create_element(
+            ElementFamily::Lagrange,
+            ReferenceCellType::Triangle,
+            0,
+            Continuity::Discontinuous,
+        );
+
+        let space = SerialFunctionSpace::new(&grid, &element);
+        let mut matrix = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size()));
+
+        group.bench_function(
+            &format!(
+                "Full assembly of {}x{} matrix",
+                space.dofmap().global_size(),
+                space.dofmap().global_size()
+            ),
+            |b| {
+                b.iter(|| {
+                    assemble_batched(
+                        &mut matrix,
+                        BoundaryOperator::SingleLayer,
+                        PDEType::Laplace,
+                        &space,
+                        &space,
+                    )
+                })
+            },
+        );
+    }
+    group.finish();
+}
+
+pub fn assembly_parts_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("assembly");
+    group.sample_size(20);
+
+    for i in 3..5 {
+        let grid = regular_sphere(i);
+        let element = create_element(
+            ElementFamily::Lagrange,
+            ReferenceCellType::Triangle,
+            0,
+            Continuity::Discontinuous,
+        );
+
+        let space = SerialFunctionSpace::new(&grid, &element);
+        let mut matrix = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size()));
+
+        let colouring = space.compute_cell_colouring();
+
+        /*
+        group.bench_function(
+            &format!(
+                "Assembly of singular terms of {}x{} matrix",
+                space.dofmap().global_size(),
+                space.dofmap().global_size()
+            ),
+            |b| {
+                b.iter(|| {
+                    batched::assemble_singular(
+                        &mut matrix,
+                        &laplace_3d::Laplace3dKernel::new(),
+                        false,
+                        false,
+                        &space,
+                        &space,
+                        &colouring,
+                        &colouring,
+                        128,
+                    )
+                })
+            },
+        );
+        */
+        group.bench_function(
+            &format!(
+                "Assembly of non-singular terms of {}x{} matrix",
+                space.dofmap().global_size(),
+                space.dofmap().global_size()
+            ),
+            |b| {
+                b.iter(|| {
+                    batched::assemble_nonsingular::<16, 16>(
+                        &mut matrix,
+                        &laplace_3d::Laplace3dKernel::new(),
+                        false,
+                        false,
+                        &space,
+                        &space,
+                        &colouring,
+                        &colouring,
+                        128,
+                    )
+                })
+            },
+        );
+    }
+    group.finish();
+}
+
+// criterion_group!(benches, full_assembly_benchmark, assembly_parts_benchmark);
+criterion_group!(benches, assembly_parts_benchmark);
+criterion_main!(benches);

bem/examples/assemble.rs

@@ -1,13 +1,12 @@
-use bempp_bem::assembly::{assemble_dense, BoundaryOperator, PDEType};
+use bempp_bem::assembly::{assemble_batched, BoundaryOperator, PDEType};
 use bempp_bem::function_space::SerialFunctionSpace;
 use bempp_element::element::create_element;
 use bempp_grid::shapes::regular_sphere;
-use bempp_tools::arrays::Array2D;
+use bempp_tools::arrays::zero_matrix;
 use bempp_traits::bem::DofMap;
 use bempp_traits::bem::FunctionSpace;
 use bempp_traits::cell::ReferenceCellType;
 use bempp_traits::element::{Continuity, ElementFamily};
-use num::complex::Complex;
 
 fn main() {
     println!("Creating grid");
@@ -34,16 +33,14 @@ fn main() {
         space1.dofmap().global_size(),
         space0.dofmap().global_size()
     );
-    let mut matrix = Array2D::<Complex<f64>>::new((
-        space1.dofmap().global_size(),
-        space0.dofmap().global_size(),
-    ));
+    let mut matrix =
+        zero_matrix::<f64>((space1.dofmap().global_size(), space0.dofmap().global_size()));
 
     println!("Assembling dense matrix (complex)");
-    assemble_dense(
+    assemble_batched(
         &mut matrix,
         BoundaryOperator::SingleLayer,
-        PDEType::Helmholtz(5.0),
+        PDEType::Laplace,
         &space0,
         &space1,
     );