diff --git a/crates/network/apsis-isl/src/batch.rs b/crates/network/apsis-isl/src/batch.rs
new file mode 100644
index 0000000..4e9500f
--- /dev/null
+++ b/crates/network/apsis-isl/src/batch.rs
@@ -0,0 +1,50 @@
+//! Deciding many links at once, optionally across several threads.
+//!
+//! Building a constellation's topology asks the same question for every pair
+//! of satellites, which is a quadratic number of independent geometry tests.
+//! Each one is small, but there are so many that the loop is the dominant
+//! cost of building the graph, and the tests do not depend on each other, so
+//! they are the natural place to use several cores.
+//!
+//! The parallel entry point returns exactly what the sequential one does, in
+//! the same order, so a caller can switch between them without changing the
+//! graph it builds.
+
+use crate::geometry::mutual_visibility;
+use apsis_math::{Vec3, parallel::map_index};
+
+/// Returns the visibility of every ordered pair of positions, in order.
+///
+/// The result is row major: the entry at `row * positions.len() + column` is
+/// the visibility from `positions[row]` to `positions[column]`. A satellite is
+/// always visible to itself, which is the caller's cue to skip the diagonal.
+///
+/// A pair whose test is undefined, because a position is not finite or the
+/// body radius is not positive, yields `false`.
+pub fn visibility_matrix(positions: &[Vec3], body_radius: f64) -> Vec<bool> {
+    let len = positions.len();
+    let mut matrix = vec![false; len * len];
+    for (row, first) in positions.iter().enumerate() {
+        for (column, second) in positions.iter().enumerate() {
+            matrix[row * len + column] =
+                mutual_visibility(*first, *second, body_radius).unwrap_or(false);
+        }
+    }
+    matrix
+}
+
+/// Returns the visibility of every ordered pair of positions, across several
+/// threads.
+///
+/// The result matches [`visibility_matrix`] element for element and in the
+/// same order. The pairs are independent, so splitting the flat index space
+/// across threads needs no coordination and produces the same matrix however
+/// the work is divided.
+pub fn visibility_matrix_parallel(positions: &[Vec3], body_radius: f64) -> Vec<bool> {
+    let len = positions.len();
+    map_index(len * len, |index| {
+        let row = index / len;
+        let column = index - row * len;
+        mutual_visibility(positions[row], positions[column], body_radius).unwrap_or(false)
+    })
+}
diff --git a/crates/network/apsis-isl/src/lib.rs b/crates/network/apsis-isl/src/lib.rs
index c82bccf..3d20526 100644
--- a/crates/network/apsis-isl/src/lib.rs
+++ b/crates/network/apsis-isl/src/lib.rs
@@ -36,9 +36,11 @@
 #![forbid(unsafe_code)]
 #![deny(missing_docs)]
 
+mod batch;
 mod geometry;
 mod link;
 
+pub use batch::{visibility_matrix, visibility_matrix_parallel};
 pub use geometry::{closest_approach_to_centre, link_geometry, mutual_visibility, sun_blocks_link};
 pub use link::{LinkGeometry, LinkState};
 
diff --git a/crates/network/apsis-isl/tests/batch.rs b/crates/network/apsis-isl/tests/batch.rs
new file mode 100644
index 0000000..67f0ad5
--- /dev/null
+++ b/crates/network/apsis-isl/tests/batch.rs
@@ -0,0 +1,123 @@
+//! Tests for the batch visibility matrix.
+//!
+//! The parallel form must produce exactly the matrix the sequential form
+//! does, in the same row major order. Every test compares the two rather than
+//! against a constant, so a content or ordering difference is caught.
+
+use apsis_isl::{visibility_matrix, visibility_matrix_parallel};
+use apsis_math::Vec3;
+
+/// The body radius used throughout, in kilometers.
+const RADIUS: f64 = 6_378.137;
+
+/// Returns a set of positions spread around a shell.
+fn shell(count: usize) -> Vec<Vec3> {
+    (0..count)
+        .map(|index| {
+            let angle = index as f64 * 0.3;
+            let radius = 6_900.0;
+            Vec3::new(
+                radius * angle.cos(),
+                radius * angle.sin(),
+                radius * 0.2 * (index as f64 * 0.1).sin(),
+            )
+        })
+        .collect()
+}
+
+#[test]
+fn a_single_position_gives_a_one_by_one_matrix() {
+    let positions = shell(1);
+    assert_eq!(
+        visibility_matrix_parallel(&positions, RADIUS),
+        visibility_matrix(&positions, RADIUS)
+    );
+    // A satellite is visible to itself over a body radius below its own.
+    assert!(visibility_matrix(&positions, RADIUS)[0]);
+}
+
+#[test]
+fn a_shell_matches_the_sequential_matrix() {
+    // Long enough that the parallel path takes the multi thread branch.
+    for count in [2_usize, 16, 64, 200] {
+        let positions = shell(count);
+        assert_eq!(
+            visibility_matrix_parallel(&positions, RADIUS),
+            visibility_matrix(&positions, RADIUS),
+            "count {count} differs"
+        );
+    }
+}
+
+#[test]
+fn an_empty_set_gives_an_empty_matrix() {
+    assert!(visibility_matrix_parallel(&[], RADIUS).is_empty());
+}
+
+#[test]
+fn the_matrix_is_row_major_and_symmetric() {
+    // Row major: the entry at row * len + column is the visibility from the
+    // row satellite to the column one. For a spherical body the geometry is
+    // symmetric, so the matrix is symmetric, which is what a wrong ordering
+    // of a row major layout would break.
+    let positions = shell(12);
+    let matrix = visibility_matrix_parallel(&positions, RADIUS);
+    let len = positions.len();
+    assert_eq!(matrix.len(), len * len);
+    for row in 0..len {
+        for column in 0..len {
+            assert_eq!(
+                matrix[row * len + column],
+                matrix[column * len + row],
+                "not symmetric at ({row}, {column})"
+            );
+        }
+    }
+}
+
+#[test]
+fn the_matrix_matches_the_scalar_test() {
+    // Each entry must agree with calling the geometry directly, so a wrong
+    // pairing of row and column is caught even when the matrix is symmetric.
+    let positions = shell(10);
+    let matrix = visibility_matrix_parallel(&positions, RADIUS);
+    let len = positions.len();
+    for row in 0..len {
+        for column in 0..len {
+            let expected = apsis_isl::mutual_visibility(positions[row], positions[column], RADIUS)
+                .unwrap_or(false);
+            assert_eq!(matrix[row * len + column], expected, "at ({row}, {column})");
+        }
+    }
+}
+
+#[test]
+fn an_unusable_body_radius_yields_all_false() {
+    // A non positive radius makes every test undefined, which the sequential
+    // form reports as false; the parallel form must do the same.
+    let positions = shell(8);
+    for radius in [0.0, -1.0] {
+        let matrix = visibility_matrix_parallel(&positions, radius);
+        assert!(matrix.iter().all(|&visible| !visible));
+    }
+}
+
+#[test]
+fn a_position_that_is_not_finite_yields_false_in_its_row() {
+    // A satellite with a non finite position cannot see anything. The rest
+    // of the matrix is unaffected.
+    let mut positions = shell(6);
+    positions[2] = Vec3::new(f64::NAN, 0.0, 0.0);
+    let matrix = visibility_matrix_parallel(&positions, RADIUS);
+    let len = positions.len();
+    for column in 0..len {
+        assert!(
+            !matrix[2 * len + column],
+            "row 2 column {column} should be false"
+        );
+        assert!(
+            !matrix[column * len + 2],
+            "row {column} column 2 should be false"
+        );
+    }
+}