diff --git a/sbp/Cargo.toml b/sbp/Cargo.toml
index e3ce98f..a1ba176 100644
--- a/sbp/Cargo.toml
+++ b/sbp/Cargo.toml
@@ -20,7 +20,7 @@ sparse = ["sprs"]
 serde1 = ["serde", "ndarray/serde"]
 
 [dev-dependencies]
-criterion = "0.3.2"
+iai = "0.1.1"
 
 [[bench]]
 name = "sbpoperators"
diff --git a/sbp/benches/sbpoperators.rs b/sbp/benches/sbpoperators.rs
index 4da115c..9f1a468 100644
--- a/sbp/benches/sbpoperators.rs
+++ b/sbp/benches/sbpoperators.rs
@@ -1,118 +1,72 @@
-use criterion::{criterion_group, criterion_main, Criterion};
+use ndarray::Array2;
 use sbp::operators::{self, SbpOperator2d};
 use sbp::Float;
 
-fn performance_benchmark(c: &mut Criterion) {
-    let mut group = c.benchmark_group("SBPoperators");
-    group.sample_size(25);
+const W: usize = 64;
+const H: usize = 64;
 
-    let w = 64;
-    let h = 64;
-
-    let x = ndarray::Array2::from_shape_fn((w, h), |(j, i)| (j * h + i) as Float);
-
-    group.bench_function("upwind4 diffxi", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::Upwind4.diffxi(x.view(), res.view_mut()))
-    });
-    group.bench_function("upwind9 diffxi", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::Upwind9.diffxi(x.view(), res.view_mut()))
-    });
-    group.bench_function("trad4 diffxi", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::SBP4.diffxi(x.view(), res.view_mut()))
-    });
-    group.bench_function("trad8 diffxi", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::SBP8.diffxi(x.view(), res.view_mut()))
-    });
-
-    group.bench_function("upwind4 diffeta", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::Upwind4.diffeta(x.view(), res.view_mut()))
-    });
-    group.bench_function("upwind9 diffeta", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::Upwind9.diffeta(x.view(), res.view_mut()))
-    });
-    group.bench_function("trad4 diffeta", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::SBP4.diffeta(x.view(), res.view_mut()))
-    });
-    group.bench_function("trad8 diffeta", |b| {
-        let mut res = x.clone();
-        b.iter(|| operators::SBP8.diffeta(x.view(), res.view_mut()))
-    });
-
-    #[cfg(feature = "sparse")]
-    {
-        use sbp::utils::kronecker_product;
-
-        fn make_operators(
-            op: impl SbpOperator2d,
-            nx: usize,
-            ny: usize,
-        ) -> (sprs::CsMat<Float>, sprs::CsMat<Float>) {
-            let dx = op.op_xi().diff_matrix(nx);
-            let dx = kronecker_product(sprs::CsMat::eye(ny).view(), dx.view());
-
-            let dy = op.op_eta().diff_matrix(ny);
-            let dy = kronecker_product(dy.view(), sprs::CsMat::eye(nx).view());
-
-            (dy, dx)
-        }
-
-        {
-            let (dy, dx) = make_operators(operators::Upwind4, w, h);
-            group.bench_function("upwind4 diffxi matrix", |b| {
-                let mut res = ndarray::Array2::zeros(x.raw_dim());
-                b.iter(|| {
-                    sprs::prod::mul_acc_mat_vec_csr(
-                        dx.view(),
-                        x.as_slice().unwrap(),
-                        res.as_slice_mut().unwrap(),
-                    )
-                })
-            });
-            group.bench_function("upwind4 diffeta matrix", |b| {
-                let mut res = ndarray::Array2::zeros(x.raw_dim());
-                b.iter(|| {
-                    sprs::prod::mul_acc_mat_vec_csr(
-                        dy.view(),
-                        x.as_slice().unwrap(),
-                        res.as_slice_mut().unwrap(),
-                    )
-                })
-            });
-        }
-        {
-            let (dy, dx) = make_operators(operators::Upwind9, w, h);
-            group.bench_function("upwind9 diffxi matrix", |b| {
-                let mut res = ndarray::Array2::zeros(x.raw_dim());
-                b.iter(|| {
-                    sprs::prod::mul_acc_mat_vec_csr(
-                        dx.view(),
-                        x.as_slice().unwrap(),
-                        res.as_slice_mut().unwrap(),
-                    )
-                })
-            });
-            group.bench_function("upwind9 diffeta matrix", |b| {
-                let mut res = ndarray::Array2::zeros(x.raw_dim());
-                b.iter(|| {
-                    sprs::prod::mul_acc_mat_vec_csr(
-                        dy.view(),
-                        x.as_slice().unwrap(),
-                        res.as_slice_mut().unwrap(),
-                    )
-                })
-            });
-        }
-    }
-
-    group.finish();
+fn baseline() {
+    let _x = Array2::<Float>::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let _res = Array2::<Float>::zeros((W, H));
 }
 
-criterion_group!(benches, performance_benchmark);
-criterion_main!(benches);
+fn upwind4_diffxi() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::Upwind4.diffxi(x.view(), res.view_mut());
+}
+fn upwind9_diffxi() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::Upwind9.diffxi(x.view(), res.view_mut());
+}
+fn trad4_diffxi() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::SBP4.diffxi(x.view(), res.view_mut());
+}
+fn trad8_diffxi() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::SBP8.diffxi(x.view(), res.view_mut());
+}
+fn upwind4_diffeta() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::Upwind4.diffeta(x.view(), res.view_mut());
+}
+fn upwind9_diffeta() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::Upwind9.diffeta(x.view(), res.view_mut());
+}
+fn trad4_diffeta() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::SBP4.diffeta(x.view(), res.view_mut());
+}
+fn trad8_diffeta() {
+    let x = Array2::from_shape_fn((W, H), |(j, i)| (j * W + i) as Float);
+    let mut res = Array2::zeros((W, H));
+
+    operators::SBP8.diffeta(x.view(), res.view_mut());
+}
+
+iai::main!(
+    baseline,
+    upwind4_diffxi,
+    upwind9_diffxi,
+    trad4_diffxi,
+    trad8_diffxi,
+    upwind4_diffeta,
+    upwind9_diffeta,
+    trad4_diffeta,
+    trad8_diffeta
+);