diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5dfb2f8..9e7904e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -41,6 +41,12 @@ written.
   `GeometricIntegratorsBase` do not join `[deps]`: an integrator is the caller's, and the test
   dependencies supply one.
 
+- **`get_DEIM_indices(Ψ)` is exported.** It returns the `k` rows that the greedy DEIM algorithm
+  selects on the basis `Ψ`, and is the index form of `get_DEIM_interpolation_matrix`:
+  `Π' * A == A[p, :]` for `p = get_DEIM_indices(Ψ)`, so a caller samples the rows without forming
+  the `N × k` matrix `Π`. It was internal under another name; its selection and its
+  `ArgumentError`s are unchanged.
+
 ## [0.2.0] — 2026-10-12
 
 ### New Features
diff --git a/src/ReducedBasisMethods.jl b/src/ReducedBasisMethods.jl
index d3a41ea..d69e17a 100644
--- a/src/ReducedBasisMethods.jl
+++ b/src/ReducedBasisMethods.jl
@@ -29,7 +29,7 @@ export lift, projection_error, reduction_error
 
 include("algorithms/deim.jl")
 
-export get_DEIM_interpolation_matrix
+export get_DEIM_interpolation_matrix, get_DEIM_indices
 
 include("algorithms/symplectic_deim.jl")
 
diff --git a/src/algorithms/deim.jl b/src/algorithms/deim.jl
index f7c58a5..89c6ef1 100644
--- a/src/algorithms/deim.jl
+++ b/src/algorithms/deim.jl
@@ -1,9 +1,13 @@
 """
-    _deim_indices(Ψ)
+    get_DEIM_indices(Ψ)
 
-The rows selected by the greedy DEIM algorithm (Algorithm 1 of
-[doi:10.1137/090766498](https://doi.org/10.1137/090766498)) on the basis `Ψ`, whose `N × k`
-columns are the basis vectors.
+Returns the `k` rows that the greedy DEIM algorithm (Algorithm 1 of
+[doi:10.1137/090766498](https://doi.org/10.1137/090766498)) selects on the basis `Ψ`, whose
+`N × k` columns are the basis vectors.
+
+This is the index form of [`get_DEIM_interpolation_matrix`](@ref): for
+`p = get_DEIM_indices(Ψ)` and `Π = get_DEIM_interpolation_matrix(Ψ)`, `Π' * A == A[p, :]` for
+every finite `N`-row matrix `A`, so `A[p, :]` evaluates `Π' * A` without forming `Π`.
 
 The index of column `i` is the row of the largest magnitude of the residual
 `r = Ψ[:, i] - Ψ[:, 1:i-1] * c`, where `c` solves the interpolation problem of the columns
@@ -15,7 +19,7 @@ is not both positive and finite — `ρ == 0` for a zero column or one exactly d
 ones before it, and `NaN` or `Inf` for a `NaN` or an `Inf` in `Ψ`. Every entry of a column
 enters that column's residual, so the test of each column covers the whole matrix.
 """
-function _deim_indices(Ψ::AbstractMatrix)
+function get_DEIM_indices(Ψ::AbstractMatrix)
     N, k = size(Ψ)
     k ≤ N || throw(ArgumentError(
         "the DEIM basis has $k columns but only $N rows, so k ≤ N is required"))
@@ -57,15 +61,15 @@ Calculates the DEIM interpolation matrix as described in Algorithm 1 in
 
 The input matrix Ψ holds the eigenvectors of the snapshot matrix fᵢ(xⱼ), 1 ≤ i ≤ N, 1 ≤ j ≤ k.
 Returns the `N × k` matrix `Π` whose column `i` is the unit vector of the row that
-[`_deim_indices`](@ref) selects for that column, so that `f(x) ≈ Ψ * inv(Π' * Ψ) * Π' * f(x)`.
+[`get_DEIM_indices`](@ref) selects for that column, so that `f(x) ≈ Ψ * inv(Π' * Ψ) * Π' * f(x)`.
 
 Its element type is `float(eltype(Ψ))`, so a `Float32` basis gives a `Float32` interpolation
 matrix. The selection, and with it the `ArgumentError` for `k > N` and for a column whose
-largest residual magnitude is not finite and positive, is [`_deim_indices`](@ref)'s.
+largest residual magnitude is not finite and positive, is [`get_DEIM_indices`](@ref)'s.
 """
 function get_DEIM_interpolation_matrix(Ψ::AbstractMatrix)
     N, k = size(Ψ)
-    p = _deim_indices(Ψ)
+    p = get_DEIM_indices(Ψ)
     Π = zeros(float(eltype(Ψ)), N, k)
     for i in 1:k
         Π[p[i], i] = one(eltype(Π))
diff --git a/src/algorithms/symplectic_deim.jl b/src/algorithms/symplectic_deim.jl
index 795eb5d..0563b10 100644
--- a/src/algorithms/symplectic_deim.jl
+++ b/src/algorithms/symplectic_deim.jl
@@ -31,7 +31,7 @@ supported.
 Throws an `ArgumentError` where `m > M`, and for a column of the density basis `U` whose largest
 DEIM residual magnitude `ρ` is not both positive and finite: `ρ == 0` where the column is zero
 or exactly dependent on the ones before it, and `ρ` is `NaN` or `Inf` where the column holds a
-`NaN` or an `Inf`. [`_deim_indices`](@ref) is the check.
+`NaN` or an `Inf`. [`get_DEIM_indices`](@ref) is the check.
 """
 struct SymplecticDEIM{T, AT <: AbstractMatrix, ST, HT, GT}
     A::AT
@@ -49,7 +49,7 @@ function SymplecticDEIM(
 
     # The weights express the sum of the DEIM-interpolated densities using sampled values:
     # 𝟙ᵀU(ΠᵀU)⁻¹Πᵀh = wᵀh[p].
-    p = _deim_indices(U)
+    p = get_DEIM_indices(U)
     w = convert(Vector{T}, (U[p, :])' \ (U' * ones(T, size(U, 1))))
 
     # `L̃ = Aᵀ·sym(L)·A`, symmetrised once more after the product, so that the Jacobian of the
@@ -123,7 +123,7 @@ the write reads 0 B, and anything a callback allocates is the caller's.
 The write gives no guarantee for an aliased array: `g` and `z` must be distinct, because the
 quadratic part is written with `mul!`, which requires its output and its input to differ.
 
-The constructor validates the density basis through [`_deim_indices`](@ref): `m > M`, and a
+The constructor validates the density basis through [`get_DEIM_indices`](@ref): `m > M`, and a
 column whose largest DEIM residual magnitude is zero or non-finite, raise `ArgumentError`.
 This write propagates errors from array operations and callbacks. A `NaN` or an `Inf` in `z`,
 or from `∇h`, enters `g` with no check.
diff --git a/test/algorithms/deim.jl b/test/algorithms/deim.jl
index 3b8a356..34b9221 100644
--- a/test/algorithms/deim.jl
+++ b/test/algorithms/deim.jl
@@ -2,7 +2,6 @@ using LinearAlgebra
 using Random
 using ReducedBasisMethods
 using ReducedComplexityModeling: Parameter, ParameterSpace
-using ReducedBasisMethods: _deim_indices
 using Test
 
 """
@@ -93,21 +92,21 @@ end
         # selects row 3, while the largest magnitude selects row 2.
         # The first pick is row 1, so neither selection repeats it.
         Ψ₁ = T[1 0; 0 -2; 0 1] ./ T[1 sqrt(5)]
-        @test _deim_indices(Ψ₁) == [1, 2]
+        @test get_DEIM_indices(Ψ₁) == [1, 2]
         @test get_DEIM_interpolation_matrix(Ψ₁) == T[1 0; 0 1; 0 0]
 
         # The first pick: the largest magnitude is in row 2, the largest value in row 3.
         Ψ₂ = T[0.1 0; -1 -2; 0.2 1] ./ T[norm([0.1, -1, 0.2]) sqrt(5)]
-        @test _deim_indices(Ψ₂) == [2, 3]
+        @test get_DEIM_indices(Ψ₂) == [2, 3]
         # At k = 1 the first pick is the whole selection.
-        @test _deim_indices(Ψ₂[:, 1:1]) == [2]
+        @test get_DEIM_indices(Ψ₂[:, 1:1]) == [2]
 
         # The residual is zero at the selected row and has equal largest magnitudes in
         # rows 2 and 3, so the tie goes to row 2.
         Ψ₃ = T[1 0; 0 -1; 0 -1] ./ T[1 sqrt(2)]
-        @test _deim_indices(Ψ₃) == [1, 2]
+        @test get_DEIM_indices(Ψ₃) == [1, 2]
 
-        p₃ = _deim_indices(Ψ₃)
+        p₃ = get_DEIM_indices(Ψ₃)
         Π₃ = get_DEIM_interpolation_matrix(Ψ₃)
         @test allunique(p₃)
         @test rank(Π₃' * Ψ₃) == 2
@@ -118,15 +117,15 @@ end
     # Π follows the element type of Ψ.
     Ψ₃₂ = Float32[1 0; 0 -1; 0 -1] ./ Float32[1 sqrt(2)]
     @test eltype(get_DEIM_interpolation_matrix(Ψ₃₂)) === Float32
-    @test _deim_indices(Ψ₃₂) == [1, 2]
+    @test get_DEIM_indices(Ψ₃₂) == [1, 2]
 
     # An integer Ψ gives a Float64 Π, never an integer one.
     Ψᵢ = Int32[1 0; 0 -2; 0 1]
-    @test _deim_indices(Ψᵢ) == [1, 2]
+    @test get_DEIM_indices(Ψᵢ) == [1, 2]
     @test eltype(get_DEIM_interpolation_matrix(Ψᵢ)) === Float64
 
     # k = 0: the empty selection, and the N × 0 matrix built from it.
-    @test _deim_indices(zeros(3, 0)) == Int[]
+    @test get_DEIM_indices(zeros(3, 0)) == Int[]
     @test size(get_DEIM_interpolation_matrix(zeros(3, 0))) == (3, 0)
 
     # k > N leaves too few distinct rows for the interpolation.
@@ -142,13 +141,13 @@ end
     for T in (Float32, Float64), bad in (zero(T), T(NaN), T(Inf))
 
         @test_throws ArgumentError get_DEIM_interpolation_matrix(reshape(T[bad, 0, 0], 3, 1))
-        @test_throws ArgumentError _deim_indices(reshape(T[bad, 0, 0], 3, 1))
+        @test_throws ArgumentError get_DEIM_indices(reshape(T[bad, 0, 0], 3, 1))
 
         # A bad first column of a two-column basis is reported as column 1, not as column 2, and
         # a zero first column is rejected before it reaches the solve, which as a singular
         # `PᵀΨ` would raise `SingularException` where the `ArgumentError` is due.
         err = try
-            _deim_indices(T[bad 1; 0 0; 0 1])
+            get_DEIM_indices(T[bad 1; 0 0; 0 1])
             nothing
         catch e
             e
@@ -156,19 +155,33 @@ end
         @test err isa ArgumentError
         @test occursin("column 1", sprint(showerror, err))
     end
-    @test_throws ArgumentError _deim_indices([0.0 1; 0 0; 0 0])
+    @test_throws ArgumentError get_DEIM_indices([0.0 1; 0 0; 0 0])
 
     # The selection is scale-invariant: the residual is compared with zero, never with a
     # cutoff. At 3.7e-8 in Float32 the largest residual is below `eps(Float32)`.
     for T in (Float32, Float64)
         Ψ, _ = pod_basis(example1_training(T); k = 20)
-        p = _deim_indices(Ψ)
+        p = get_DEIM_indices(Ψ)
 
-        @test _deim_indices(T(3.7e8) .* Ψ) == p
-        @test _deim_indices(T(3.7e-8) .* Ψ) == p
+        @test get_DEIM_indices(T(3.7e8) .* Ψ) == p
+        @test get_DEIM_indices(T(3.7e-8) .* Ψ) == p
 
         # A subnormal first column is positive, so it is selected rather than refused: column 1's
         # test is the same `0 < ρ`, with no threshold.
-        @test _deim_indices(reshape(T[nextfloat(zero(T)), zero(T)], 2, 1)) == [1]
+        @test get_DEIM_indices(reshape(T[nextfloat(zero(T)), zero(T)], 2, 1)) == [1]
+    end
+end
+
+@testset "indices are the index form of the interpolation matrix" begin
+    # Π' * A samples the rows of A that get_DEIM_indices selects, for any basis and any finite A.
+    for T in (Float32, Float64)
+        rng = Xoshiro(7)
+        Ψ = Matrix(qr(randn(rng, T, 50, 8)).Q)
+        A = randn(rng, T, 50, 5)
+        p = get_DEIM_indices(Ψ)
+
+        @test length(p) == 8
+        @test allunique(p)
+        @test get_DEIM_interpolation_matrix(Ψ)' * A == A[p, :]
     end
 end
diff --git a/test/algorithms/symplectic_deim.jl b/test/algorithms/symplectic_deim.jl
index da22012..ccc78a6 100644
--- a/test/algorithms/symplectic_deim.jl
+++ b/test/algorithms/symplectic_deim.jl
@@ -3,7 +3,6 @@ using GeometricIntegrators
 using LinearAlgebra
 using Random
 using ReducedBasisMethods
-using ReducedBasisMethods: _deim_indices
 using ReducedComplexityModeling: Parameter, ParameterSpace
 using Test
 
@@ -178,7 +177,7 @@ end
         # At m = M the interpolation is exact for the densities: w = 𝟙 and H̃(z) = H(Az).
         sd = SymplecticDEIM(b.A, b.L, U, stencil_test, h_test, ∇h_test)
 
-        @test _deim_indices(U) != 1:NQ
+        @test get_DEIM_indices(U) != 1:NQ
         @test norm(sd.w .- one(T), Inf) ≤ sqrt(eps(T))
 
         for z in eachslice(randn(Xoshiro(4), T, 2KT, 20); dims = 2)
@@ -213,7 +212,7 @@ end
         ∇hᵢ(xs, i) = [(1 + i / NQ) * sin(xs[1])]
         sdᵢ = SymplecticDEIM(b.A, b.L, U, stencil_test, hᵢ, ∇hᵢ)
 
-        @test _deim_indices(U) != 1:NQ
+        @test get_DEIM_indices(U) != 1:NQ
 
         for z in eachslice(randn(Xoshiro(5), T, 2KT, 20); dims = 2)
             Az = b.A * z