494 lines
15 KiB
Go
494 lines
15 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package linalg
|
||
|
||
import (
|
||
"cmp"
|
||
"slices"
|
||
"sourcedock.dev/petrbalvin/tensor/internal/base"
|
||
)
|
||
|
||
// Fill-reducing orderings for the sparse factorisations. The
|
||
// symmetrised adjacency builder and the orderings live here; the
|
||
// factorisation consumes nothing but a permutation, so a new ordering
|
||
// slots in without touching the numeric code.
|
||
|
||
// symmetrisedAdjacency returns the adjacency lists of the pattern of
|
||
// A + Aᵀ with the diagonal dropped: for every stored entry (i, j),
|
||
// i ≠ j, i sits on j's list and j on i's. Every list is sorted by
|
||
// ascending degree and then index, which fixes the visiting order the
|
||
// breadth first search uses. The lists slice out of one flat arena
|
||
// sized by a counting pass, so the construction allocates twice
|
||
// whatever the vertex count; the degree snapshot counts the entries a
|
||
// list received, duplicates included, exactly as the incremental build
|
||
// measured them, and the degree-then-index sort is a total order, so
|
||
// the fill order inside a list cannot reach the result.
|
||
func symmetrisedAdjacency(c *SparseCSC) ([][]int, error) {
|
||
adj, counts, err := symmetrisedAdjacencyBuild(c)
|
||
if err != nil {
|
||
return nil, err
|
||
}
|
||
for i := range adj {
|
||
slices.SortStableFunc(adj[i], func(x, y int) int {
|
||
if d := cmp.Compare(counts[x], counts[y]); d != 0 {
|
||
return d
|
||
}
|
||
return cmp.Compare(x, y)
|
||
})
|
||
}
|
||
return adj, nil
|
||
}
|
||
|
||
// adjacencyIndexOrdered returns the same adjacency lists in plain
|
||
// ascending index order, the order the minimum-degree absorption
|
||
// merges: it skips the degree-then-index sort the breadth first
|
||
// search needs and no index-ordered consumer does.
|
||
func adjacencyIndexOrdered(c *SparseCSC) ([][]int, error) {
|
||
adj, _, err := symmetrisedAdjacencyBuild(c)
|
||
return adj, err
|
||
}
|
||
|
||
// symmetrisedAdjacencyBuild builds the adjacency lists in ascending
|
||
// index order: one flat arena sized by a counting pass, every list
|
||
// sorted and compacted. The counts it also returns are the raw
|
||
// duplicate-inclusive neighbour counts the build measured, the values
|
||
// the degree-then-index ordering sorts by.
|
||
func symmetrisedAdjacencyBuild(c *SparseCSC) ([][]int, []int, error) {
|
||
counts := make([]int, c.Cols)
|
||
for j := range c.Cols {
|
||
for p := c.ColStart[j]; p < c.ColStart[j+1]; p++ {
|
||
i := c.RowIdx[p]
|
||
if i == j {
|
||
continue
|
||
}
|
||
if i < 0 || i >= c.Cols {
|
||
return nil, nil, base.Errf("adjacency: row index %d out of range for %d columns", i, c.Cols)
|
||
}
|
||
counts[j]++
|
||
counts[i]++
|
||
}
|
||
}
|
||
total := 0
|
||
for i := range c.Cols {
|
||
total += counts[i]
|
||
}
|
||
arena := make([]int, total)
|
||
adj := make([][]int, c.Cols)
|
||
next := make([]int, c.Cols)
|
||
off := 0
|
||
for i := range c.Cols {
|
||
adj[i] = arena[off : off+counts[i]]
|
||
next[i] = off
|
||
off += counts[i]
|
||
}
|
||
for j := range c.Cols {
|
||
for p := c.ColStart[j]; p < c.ColStart[j+1]; p++ {
|
||
i := c.RowIdx[p]
|
||
if i == j {
|
||
continue
|
||
}
|
||
arena[next[j]] = i
|
||
next[j]++
|
||
arena[next[i]] = j
|
||
next[i]++
|
||
}
|
||
}
|
||
for i := range c.Cols {
|
||
slices.Sort(adj[i])
|
||
adj[i] = slices.Compact(adj[i])
|
||
}
|
||
return adj, counts, nil
|
||
}
|
||
|
||
// breadthFirstSearch walks the component of start and returns the
|
||
// visit order, labelling level for every vertex it reaches. It
|
||
// initialises level to the unvisited state itself, so the caller can
|
||
// pass a scratch slice.
|
||
func breadthFirstSearch(adj [][]int, level []int, start int) []int {
|
||
for i := range level {
|
||
level[i] = -1
|
||
}
|
||
order := []int{start}
|
||
level[start] = 0
|
||
for head := 0; head < len(order); head++ {
|
||
v := order[head]
|
||
for _, u := range adj[v] {
|
||
if level[u] < 0 {
|
||
level[u] = level[v] + 1
|
||
order = append(order, u)
|
||
}
|
||
}
|
||
}
|
||
return order
|
||
}
|
||
|
||
// deepestLevel returns the vertices of the largest level in visit
|
||
// order.
|
||
func deepestLevel(order []int, level []int) []int {
|
||
max := 0
|
||
for _, v := range order {
|
||
if level[v] > max {
|
||
max = level[v]
|
||
}
|
||
}
|
||
last := make([]int, 0, len(order))
|
||
for _, v := range order {
|
||
if level[v] == max {
|
||
last = append(last, v)
|
||
}
|
||
}
|
||
return last
|
||
}
|
||
|
||
// pseudoPeripheralStart finds a start vertex whose eccentricity is
|
||
// close to the component's diameter: two sweeps of taking the
|
||
// shallowest-degree vertex of the deepest level, the standard
|
||
// pseudo-peripheral heuristic. Every choice is deterministic, so the
|
||
// ordering is a pure function of the pattern.
|
||
func pseudoPeripheralStart(adj [][]int, level []int, guess int) int {
|
||
start := guess
|
||
for range 2 {
|
||
order := breadthFirstSearch(adj, level, start)
|
||
last := deepestLevel(order, level)
|
||
best := last[0]
|
||
for _, v := range last[1:] {
|
||
if len(adj[v]) < len(adj[best]) {
|
||
best = v
|
||
}
|
||
}
|
||
start = best
|
||
}
|
||
return start
|
||
}
|
||
|
||
// reverseCuthillMcKee returns the elimination order of the symmetric
|
||
// pattern: the vertices of every component in reverse breadth first
|
||
// order from a pseudo-peripheral start. The reverse of the search
|
||
// order is what shrinks the bandwidth, and with it the fill a
|
||
// triangular factorisation produces. Components are taken in
|
||
// ascending order of their smallest unreached vertex, so the result
|
||
// is a pure function of the pattern: position k holds the original
|
||
// index of the row eliminated k-th.
|
||
func reverseCuthillMcKee(c *SparseCSC) ([]int, error) {
|
||
adj, err := symmetrisedAdjacency(c)
|
||
if err != nil {
|
||
return nil, err
|
||
}
|
||
level := make([]int, c.Cols)
|
||
visited := make([]bool, c.Cols)
|
||
order := make([]int, 0, c.Cols)
|
||
for s := range c.Cols {
|
||
if visited[s] {
|
||
continue
|
||
}
|
||
start := pseudoPeripheralStart(adj, level, s)
|
||
component := breadthFirstSearch(adj, level, start)
|
||
for _, v := range component {
|
||
visited[v] = true
|
||
}
|
||
for _, c := range slices.Backward(component) {
|
||
order = append(order, c)
|
||
}
|
||
}
|
||
return order, nil
|
||
}
|
||
|
||
// minimumDegree returns the elimination order of the symmetric
|
||
// pattern by minimum degree: at every step the uneliminated vertex
|
||
// with the fewest remaining neighbours is eliminated, its neighbours
|
||
// absorb its pattern (the element the elimination creates), and
|
||
// degree ties break to the smallest index, so the order is a pure
|
||
// function of the pattern. On irregular patterns it shrinks the fill
|
||
// well below what a bandwidth ordering reaches; on structured meshes
|
||
// the reverse Cuthill-McKee order is its match.
|
||
//
|
||
// The selection runs on a lazy frontier instead of a scan over the
|
||
// surviving vertices: a binary heap of (degree, vertex) pairs where a
|
||
// vertex enters when its degree is first known and its pair is rewritten
|
||
// in place whenever the degree moves, so the heap holds every surviving
|
||
// vertex exactly once and hands out the smallest (degree, index) among
|
||
// them, which is the first vertex of the smallest degree the scan would
|
||
// pick, at a logarithmic price per entry instead of a full sweep per
|
||
// elimination.
|
||
func minimumDegree(c *SparseCSC) ([]int, error) {
|
||
adj, err := adjacencyIndexOrdered(c)
|
||
if err != nil {
|
||
return nil, err
|
||
}
|
||
n := c.Cols
|
||
eliminated := make([]bool, n)
|
||
degree := make([]int, n)
|
||
frontier := °reeFrontier{entries: make([]degreeEntry, 0, n)}
|
||
arena := &intArena{}
|
||
for i := range n {
|
||
degree[i] = len(adj[i])
|
||
frontier.push(degreeEntry{degree[i], i})
|
||
}
|
||
order := make([]int, 0, n)
|
||
for range n {
|
||
if frontier.len() == 0 {
|
||
return nil, base.Errf("minimumDegree: no vertex left to eliminate")
|
||
}
|
||
p := frontier.pop().vertex
|
||
order = append(order, p)
|
||
eliminated[p] = true
|
||
absorbElement(adj, degree, eliminated, p, frontier, arena)
|
||
}
|
||
return order, nil
|
||
}
|
||
|
||
// degreeEntry pairs a vertex with the degree it carried when it
|
||
// entered the frontier.
|
||
type degreeEntry struct {
|
||
degree int
|
||
vertex int
|
||
}
|
||
|
||
// before orders the pairs by degree, then index: the heap's top is
|
||
// the vertex the selection would take.
|
||
func (e degreeEntry) before(other degreeEntry) bool {
|
||
if e.degree != other.degree {
|
||
return e.degree < other.degree
|
||
}
|
||
return e.vertex < other.vertex
|
||
}
|
||
|
||
// degreeFrontier is a binary min-heap of degreeEntry pairs holding
|
||
// every vertex at most once: push rewrites the vertex's pair in place
|
||
// when one is already queued, so the heap never carries a stale pair
|
||
// and its top is always the vertex the selection takes. The absence of
|
||
// stale pairs keeps the heap at the vertex count instead of one entry
|
||
// per degree movement, which is what a sift-down per pop would
|
||
// otherwise charge.
|
||
type degreeFrontier struct {
|
||
entries []degreeEntry
|
||
pos []int
|
||
}
|
||
|
||
// len reports how many pairs the frontier holds.
|
||
func (f *degreeFrontier) len() int { return len(f.entries) }
|
||
|
||
// push records the vertex's current degree: the pair of a vertex
|
||
// already queued moves to its new place, a vertex without a pair
|
||
// enters at the end. pos grows to cover the vertex and fills its new
|
||
// slots with -1, so an absent vertex reads -1 and a queued one reads
|
||
// its heap index.
|
||
func (f *degreeFrontier) push(e degreeEntry) {
|
||
for len(f.pos) <= e.vertex {
|
||
f.pos = append(f.pos, -1)
|
||
}
|
||
if at := f.pos[e.vertex]; at != -1 {
|
||
f.entries[at] = e
|
||
f.fix(at)
|
||
return
|
||
}
|
||
f.pos[e.vertex] = len(f.entries)
|
||
f.entries = append(f.entries, e)
|
||
f.up(len(f.entries) - 1)
|
||
}
|
||
|
||
// pop removes and returns the smallest pair. It is the caller's job to
|
||
// check len first; the returned pair is never stale.
|
||
func (f *degreeFrontier) pop() degreeEntry {
|
||
entries := f.entries
|
||
top := entries[0]
|
||
f.pos[top.vertex] = -1
|
||
last := len(entries) - 1
|
||
entries[0] = entries[last]
|
||
f.entries = entries[:last]
|
||
if last > 0 {
|
||
f.pos[entries[0].vertex] = 0
|
||
f.down(0)
|
||
}
|
||
return top
|
||
}
|
||
|
||
// fix restores the heap order around i after the pair at i changed.
|
||
func (f *degreeFrontier) fix(i int) {
|
||
if !f.up(i) {
|
||
f.down(i)
|
||
}
|
||
}
|
||
|
||
// up lifts the pair at i until its parent stops outranking it, and
|
||
// reports whether it moved at all.
|
||
func (f *degreeFrontier) up(i int) bool {
|
||
entries := f.entries
|
||
moved := false
|
||
for i > 0 {
|
||
parent := (i - 1) / 2
|
||
if !entries[i].before(entries[parent]) {
|
||
break
|
||
}
|
||
entries[i], entries[parent] = entries[parent], entries[i]
|
||
f.pos[entries[i].vertex] = i
|
||
i = parent
|
||
moved = true
|
||
}
|
||
f.pos[entries[i].vertex] = i
|
||
return moved
|
||
}
|
||
|
||
// down drops the pair at i until both children stop outranking it.
|
||
func (f *degreeFrontier) down(i int) {
|
||
entries := f.entries
|
||
for {
|
||
child := 2*i + 1
|
||
if child >= len(entries) {
|
||
break
|
||
}
|
||
if right := child + 1; right < len(entries) && entries[right].before(entries[child]) {
|
||
child = right
|
||
}
|
||
if !entries[child].before(entries[i]) {
|
||
break
|
||
}
|
||
entries[i], entries[child] = entries[child], entries[i]
|
||
f.pos[entries[i].vertex] = i
|
||
i = child
|
||
}
|
||
f.pos[entries[i].vertex] = i
|
||
}
|
||
|
||
// intArena hands out int slices from one bump-allocated backing block,
|
||
// so the minimum-degree absorption merges cost an allocation per block
|
||
// instead of one per merge. A returned slice owns its range: the merge
|
||
// appends into a zero-length view whose capacity is its exact upper
|
||
// bound, and the arena never hands the same range out twice. The
|
||
// capacity a merge leaves over returns to the bump space through
|
||
// release, which must reach the arena before its next alloc, while the
|
||
// released slice is still the most recent one. Backing blocks are kept
|
||
// alive by the adjacency lists that point into them, not by the arena.
|
||
type intArena struct {
|
||
buf []int
|
||
// bump is the offset of the next free entry in buf; the block is
|
||
// replaced once the space behind bump no longer fits a request.
|
||
bump int
|
||
}
|
||
|
||
// arenaBlockInts is the smallest backing block: big enough that the
|
||
// small early merges share one allocation, small enough that a block
|
||
// full of dead lists retires without pinning much memory.
|
||
const arenaBlockInts = 16384
|
||
|
||
// alloc returns a zero-length slice with capacity n. The arena works
|
||
// from an empty buffer as well as a partially consumed one: the bump
|
||
// offset starts at zero and the block check reads it against the block
|
||
// length, so a zero-value arena allocates rather than panics.
|
||
func (a *intArena) alloc(n int) []int {
|
||
if len(a.buf)-a.bump < n {
|
||
a.buf = make([]int, max(arenaBlockInts, n))
|
||
a.bump = 0
|
||
}
|
||
out := a.buf[a.bump : a.bump : a.bump+n]
|
||
a.bump += n
|
||
return out
|
||
}
|
||
|
||
// release folds the capacity a merge left over back into the bump
|
||
// space, so the next alloc reuses it instead of moving past it. The
|
||
// slice must be the arena's most recent allocation and must not be
|
||
// released twice.
|
||
func (a *intArena) release(out []int) {
|
||
a.bump -= cap(out) - len(out)
|
||
}
|
||
|
||
// absorbElement merges the element the eliminated vertex p leaves
|
||
// behind into every surviving neighbour's adjacency list, the fill
|
||
// the elimination creates: adj(j) becomes adj(j) ∪ adj(p) minus j and
|
||
// minus every eliminated vertex. The set difference matters: p sits on
|
||
// j's list and j on p's, so the plain union hands j its own index
|
||
// back, and a self entry both inflates the degree the selection reads
|
||
// and spreads to other lists through later absorptions. A neighbour's
|
||
// degree is the count of surviving entries in the merged list, the
|
||
// number a scan would recount at the same step, kept here where the
|
||
// list changes; a neighbour whose degree moved re-enters the frontier.
|
||
func absorbElement(adj [][]int, degree []int, eliminated []bool, p int, frontier *degreeFrontier, arena *intArena) {
|
||
// The element list is compacted once, in place, before any merge:
|
||
// p is eliminated and its list is dead, no vertex becomes
|
||
// eliminated while the absorb runs, and every neighbour's merge
|
||
// walks the same list, so one filter pass here replaces one per
|
||
// merged element. The compaction is stable, so the merge order is
|
||
// the list's own.
|
||
b := adj[p]
|
||
w := 0
|
||
for _, u := range b {
|
||
if !eliminated[u] {
|
||
b[w] = u
|
||
w++
|
||
}
|
||
}
|
||
b = b[:w]
|
||
for _, j := range b {
|
||
before := degree[j]
|
||
adj[j] = mergedElement(adj[j], b, j, eliminated, degree, arena)
|
||
if degree[j] != before {
|
||
frontier.push(degreeEntry{degree[j], j})
|
||
}
|
||
}
|
||
}
|
||
|
||
// mergedElement merges two sorted unique slices into one sorted unique
|
||
// slice without the self entry and without the eliminated entries, and
|
||
// stores the surviving count into degree[self]. The b side arrives
|
||
// already free of eliminated entries, so only the entries that come
|
||
// from a (and the merged equal pairs) pay the filter; b still carries
|
||
// self exactly once, which its own check drops. The degree the
|
||
// selection reads is the count of surviving entries either way, so the
|
||
// elimination order the frontier produces cannot move. The result lands
|
||
// in an arena slice capped at len(a)+len(b), which the merge cannot
|
||
// exceed: every entry of both lists is emitted at most once and the
|
||
// self entry never is.
|
||
func mergedElement(a, b []int, self int, eliminated []bool, degree []int, arena *intArena) []int {
|
||
n := len(a) + len(b)
|
||
out := arena.alloc(n)
|
||
// The arena reserved exactly n entries, so the merge writes through
|
||
// the stretched view and cuts the result back at the end.
|
||
out = out[:n:n]
|
||
w := 0
|
||
i, j := 0, 0
|
||
for i < len(a) && j < len(b) {
|
||
u := a[i]
|
||
switch v := b[j]; {
|
||
case u < v:
|
||
i++
|
||
if u != self && !eliminated[u] {
|
||
out[w] = u
|
||
w++
|
||
}
|
||
case v < u:
|
||
j++
|
||
if v != self {
|
||
out[w] = v
|
||
w++
|
||
}
|
||
default:
|
||
i++
|
||
j++
|
||
if u != self && !eliminated[u] {
|
||
out[w] = u
|
||
w++
|
||
}
|
||
}
|
||
}
|
||
for ; i < len(a); i++ {
|
||
if u := a[i]; u != self && !eliminated[u] {
|
||
out[w] = u
|
||
w++
|
||
}
|
||
}
|
||
for ; j < len(b); j++ {
|
||
if v := b[j]; v != self {
|
||
out[w] = v
|
||
w++
|
||
}
|
||
}
|
||
out = out[:w]
|
||
degree[self] = w
|
||
// The merge filled w of the capacity it reserved: the tail goes
|
||
// back to the arena before the next merge carves past it.
|
||
arena.release(out)
|
||
return out
|
||
}
|