@@ -9,6 +9,17 @@
//! unstructured problems for parallel processing*, Computing Systems in Engineering 2:135–
//! 1991, with vertex separators). Ties are broken by node index, so the ordering is
//! deterministic.
//!
//! A plane cuts a thin curved wall (a shell meshed with solids, a bellows, a revolved part) many
//! times, so its separators are large there. [`graph_dissection`] takes its separators from the
//! graph instead: a middle level of the breadth-first level structure rooted at a
//! pseudo-peripheral node (A. George, J. W. H. Liu, *An automatic nested dissection algorithm for
//! irregular finite element problems*, SIAM J. Numer. Anal. 15(5):1053–
//! pseudo-peripheral node of N. E. Gibbs, W. G. Poole, P. K. Stockmeyer, *An algorithm for
//! reducing the bandwidth and profile of a sparse matrix*, SIAM J. Numer. Anal. 13(2):236–
//! 1976), which cuts such a wall across its surface. [`fill_reducing_order`] computes both and
//! keeps the one whose Cholesky factor has fewer entries; the symbolic count costs O(|L|) time
//! and O(n + |A|) memory, small next to the numeric factorization.
/// Parts this small are ordered as they come.
const LEAF : usize = 16 ;
@@ -86,6 +97,179 @@ pub fn nested_dissection(xyz: &[[f64; 3]], adjacency: &[Vec<usize>]) -> Vec<usiz
order
}
/// The nested dissection order of the nodes that fills the factor less: by coordinate planes or
/// by the graph's level structures. Ties keep the coordinate order.
pub fn fill_reducing_order ( xyz : & [ [ f64 ; 3 ] ] , adjacency : & [ Vec < usize > ] ) -> Vec < usize > {
let by_planes = nested_dissection ( xyz , adjacency ) ;
let by_graph = graph_dissection ( adjacency ) ;
if fill ( adjacency , & by_graph ) < fill ( adjacency , & by_planes ) { by_graph } else { by_planes }
}
/// Searches for a pseudo-peripheral node per part.
const PERIPHERAL_SEARCHES : usize = 8 ;
/// The breadth-first level structure of the nodes with `owner == token` reached from `root`.
/// Marks each reached node's level in `level` (which must be `usize::MAX` there on entry).
fn level_structure ( adjacency : & [ Vec < usize > ] , root : usize , token : usize , owner : & [ usize ] , level : & mut [ usize ] ) -> Vec < Vec < usize > > {
let mut out = vec! [ vec! [ root ] ] ;
if let Some ( l ) = level . get_mut ( root ) {
* l = 0 ;
}
loop {
let depth = out . len ( ) ;
let mut next = Vec ::new ( ) ;
for & v in out . last ( ) . into_iter ( ) . flatten ( ) {
for & w in adjacency . get ( v ) . into_iter ( ) . flatten ( ) {
if owner . get ( w ) = = Some ( & token )
& & let Some ( l ) = level . get_mut ( w )
& & * l = = usize ::MAX
{
* l = depth ;
next . push ( w ) ;
}
}
}
if next . is_empty ( ) {
return out ;
}
out . push ( next ) ;
}
}
/// Clears the levels of `structure`'s nodes.
fn clear ( level : & mut [ usize ] , structure : & [ Vec < usize > ] ) {
for & v in structure . iter ( ) . flatten ( ) {
if let Some ( l ) = level . get_mut ( v ) {
* l = usize ::MAX ;
}
}
}
/// A permutation of the nodes by nested dissection on the graph alone (George & Liu 1978, see
/// the module notes): each part's separator is the nodes of its middle level that touch the
/// next level. Disconnected parts are ordered one after the other. Ties go by node index, so
/// the ordering is deterministic.
pub fn graph_dissection ( adjacency : & [ Vec < usize > ] ) -> Vec < usize > {
let n = adjacency . len ( ) ;
let mut order = Vec ::with_capacity ( n ) ;
// The part each node is in (a token per part), and its level in the current search.
let mut owner = vec! [ 0 usize ; n ] ;
let mut level = vec! [ usize ::MAX ; n ] ;
let mut token = 0 usize ;
enum Task {
Split ( Vec < usize > ) ,
Emit ( Vec < usize > ) ,
}
let mut stack = vec! [ Task ::Split ( ( 0 .. n ) . collect ( ) ) ] ;
while let Some ( task ) = stack . pop ( ) {
let part = match task {
Task ::Emit ( sep ) = > {
order . extend ( sep ) ;
continue ;
}
Task ::Split ( p ) = > p ,
} ;
let Some ( & first ) = part . first ( ) else { continue } ;
if part . len ( ) < = LEAF {
order . extend ( part ) ;
continue ;
}
token + = 1 ;
for & v in & part {
if let Some ( o ) = owner . get_mut ( v ) {
* o = token ;
}
}
let degree = | v : usize , owner : & [ usize ] | adjacency . get ( v ) . map_or ( 0 , | row | row . iter ( ) . filter ( | w | owner . get ( * * w ) = = Some ( & token ) ) . count ( ) ) ;
// A pseudo-peripheral root: from the first node, move to the least connected node of the
// last level while that deepens the structure.
let mut structure = level_structure ( adjacency , first , token , & owner , & mut level ) ;
for _ in 0 .. PERIPHERAL_SEARCHES {
let Some ( candidate ) = structure . last ( ) . and_then ( | last | last . iter ( ) . copied ( ) . min_by_key ( | & v | ( degree ( v , & owner ) , v ) ) ) else { break } ;
clear ( & mut level , & structure ) ;
let deeper = level_structure ( adjacency , candidate , token , & owner , & mut level ) ;
if deeper . len ( ) < = structure . len ( ) {
clear ( & mut level , & deeper ) ;
structure =
level_structure ( adjacency , structure . first ( ) . and_then ( | l | l . first ( ) ) . copied ( ) . unwrap_or ( first ) , token , & owner , & mut level ) ;
break ;
}
structure = deeper ;
}
let reached : usize = structure . iter ( ) . map ( Vec ::len ) . sum ( ) ;
if reached < part . len ( ) {
// Disconnected: the reached component, then the rest.
let rest : Vec < usize > = part . iter ( ) . copied ( ) . filter ( | & v | level . get ( v ) = = Some ( & usize ::MAX ) ) . collect ( ) ;
clear ( & mut level , & structure ) ;
let mut component : Vec < usize > = structure . into_iter ( ) . flatten ( ) . collect ( ) ;
component . sort_unstable ( ) ;
stack . push ( Task ::Split ( rest ) ) ;
stack . push ( Task ::Split ( component ) ) ;
continue ;
}
if structure . len ( ) < 3 {
clear ( & mut level , & structure ) ;
order . extend ( part ) ;
continue ;
}
// The middle level: the first at which the levels up to it hold half the part, kept
// off the last level so both sides are non-empty.
let mut count = 0 ;
let mut middle = 1 ;
for ( i , l ) in structure . iter ( ) . enumerate ( ) . skip ( 1 ) {
count + = l . len ( ) ;
if 2 * count > = part . len ( ) {
middle = i . min ( structure . len ( ) - 2 ) ;
break ;
}
}
let touches_next =
| v : usize | adjacency . get ( v ) . is_some_and ( | row | row . iter ( ) . any ( | & w | owner . get ( w ) = = Some ( & token ) & & level . get ( w ) = = Some ( & ( middle + 1 ) ) ) ) ;
let mut separator = Vec ::new ( ) ;
let mut before : Vec < usize > = structure . iter ( ) . take ( middle ) . flatten ( ) . copied ( ) . collect ( ) ;
for & v in structure . get ( middle ) . into_iter ( ) . flatten ( ) {
if touches_next ( v ) { separator . push ( v ) } else { before . push ( v ) }
}
let mut after : Vec < usize > = structure . iter ( ) . skip ( middle + 1 ) . flatten ( ) . copied ( ) . collect ( ) ;
clear ( & mut level , & structure ) ;
before . sort_unstable ( ) ;
separator . sort_unstable ( ) ;
after . sort_unstable ( ) ;
// Order: before, then after, then the separator (the stack runs in reverse).
stack . push ( Task ::Emit ( separator ) ) ;
stack . push ( Task ::Split ( after ) ) ;
stack . push ( Task ::Split ( before ) ) ;
}
order
}
/// The entries of the Cholesky factor of the node graph `adjacency` eliminated in `order` (one
/// per node pair): the fill an ordering gives. `order` must be a permutation of the nodes.
pub fn fill ( adjacency : & [ Vec < usize > ] , order : & [ usize ] ) -> usize {
let n = adjacency . len ( ) ;
let mut position = vec! [ usize ::MAX ; n ] ;
for ( k , & v ) in order . iter ( ) . enumerate ( ) {
if let Some ( p ) = position . get_mut ( v ) {
* p = k ;
}
}
let mut ptr = Vec ::with_capacity ( n + 1 ) ;
ptr . push ( 0 ) ;
let mut idx = Vec ::new ( ) ;
let mut rows = Vec ::new ( ) ;
for ( k , & v ) in order . iter ( ) . enumerate ( ) {
rows . clear ( ) ;
rows . extend ( adjacency . get ( v ) . into_iter ( ) . flatten ( ) . filter_map ( | w | position . get ( * w ) . copied ( ) ) . filter ( | & r | r < k ) ) ;
rows . push ( k ) ;
rows . sort_unstable ( ) ;
rows . dedup ( ) ;
idx . extend_from_slice ( & rows ) ;
ptr . push ( idx . len ( ) ) ;
}
let val = vec! [ 0.0 ; idx . len ( ) ] ;
crate ::cholesky ::factor_entries ( & crate ::cholesky ::UpperColumns { n : order . len ( ) , ptr , idx , val } )
}
#[ cfg(test) ]
mod tests {
use super ::* ;
@@ -116,4 +300,170 @@ mod tests {
assert! ( last . iter ( ) . all ( | v | v / m = = last [ 0 ] / m ) | | last . iter ( ) . all ( | v | v % m = = last [ 0 ] % m ) , " {last:?} " ) ;
assert_eq! ( order , nested_dissection ( & xyz , & adj ) , " deterministic " ) ;
}
/// The node graph of a structured hex8 grid of `dims` cells (27-point stencil), with the
/// node coordinates from `at(i, j, k)`; `wrap` closes the first index around.
fn hex_grid ( dims : [ usize ; 3 ] , wrap : bool , at : impl Fn ( usize , usize , usize ) -> [ f64 ; 3 ] ) -> ( Vec < [ f64 ; 3 ] > , Vec < Vec < usize > > ) {
let ni = if wrap { dims [ 0 ] } else { dims [ 0 ] + 1 } ;
let ( nj , nk ) = ( dims [ 1 ] + 1 , dims [ 2 ] + 1 ) ;
let id = | i : usize , j : usize , k : usize | ( i % ni ) + ni * ( j + nj * k ) ;
let mut xyz = vec! [ [ 0.0 ; 3 ] ; ni * nj * nk ] ;
let mut adj = vec! [ Vec ::new ( ) ; ni * nj * nk ] ;
for k in 0 .. nk {
for j in 0 .. nj {
for i in 0 .. ni {
xyz [ id ( i , j , k ) ] = at ( i , j , k ) ;
for t in 0 .. 27 {
let ( di , dj , dk ) = ( t % 3 , ( t / 3 ) % 3 , t / 9 ) ;
let ( a , b , c ) = ( i + ni + di - 1 , j + dj , k + dk ) ;
let a = if wrap {
a % ni
} else if a < ni | | a > = 2 * ni {
continue ;
} else {
a - ni
} ;
if ( 1 ..= nj ) . contains ( & b ) & & ( 1 ..= nk ) . contains ( & c ) {
adj [ id ( i , j , k ) ] . push ( id ( a , b - 1 , c - 1 ) ) ;
}
}
}
}
}
for row in & mut adj {
row . sort_unstable ( ) ;
row . dedup ( ) ;
}
( xyz , adj )
}
fn is_permutation ( order : & [ usize ] , n : usize ) -> bool {
let mut sorted = order . to_vec ( ) ;
sorted . sort_unstable ( ) ;
sorted = = ( 0 .. n ) . collect ::< Vec < _ > > ( )
}
/// A thin cylinder wall, two solid elements thick, 160 around and 40 high.
fn cylinder_wall ( ) -> ( Vec < [ f64 ; 3 ] > , Vec < Vec < usize > > ) {
hex_grid ( [ 160 , 2 , 40 ] , true , | i , j , k | {
let t = i as f64 / 160.0 * std ::f64 ::consts ::TAU ;
let r = 50.0 + j as f64 * 0.5 ;
[ r * t . cos ( ) , r * t . sin ( ) , k as f64 * 2.0 ]
} )
}
#[ test ]
fn graph_dissection_is_a_deterministic_permutation ( ) {
let ( _ , adj ) = cylinder_wall ( ) ;
let order = graph_dissection ( & adj ) ;
assert! ( is_permutation ( & order , adj . len ( ) ) ) ;
assert_eq! ( order , graph_dissection ( & adj ) ) ;
// Two disconnected blocks and an isolated node.
let ( _ , mut two ) = hex_grid ( [ 6 , 6 , 6 ] , false , | i , j , k | [ i as f64 , j as f64 , k as f64 ] ) ;
let m = two . len ( ) ;
two . extend ( two . clone ( ) . into_iter ( ) . map ( | row | row . into_iter ( ) . map ( | v | v + m ) . collect ( ) ) ) ;
two . push ( Vec ::new ( ) ) ;
assert! ( is_permutation ( & graph_dissection ( & two ) , 2 * m + 1 ) ) ;
assert! ( graph_dissection ( & [ ] ) . is_empty ( ) ) ;
}
#[ test ]
fn the_graph_cuts_a_thin_curved_wall_with_less_fill ( ) {
let ( xyz , adj ) = cylinder_wall ( ) ;
let planes = fill ( & adj , & nested_dissection ( & xyz , & adj ) ) ;
let graph = fill ( & adj , & graph_dissection ( & adj ) ) ;
assert! ( graph < planes , " graph {graph} planes {planes} " ) ;
assert_eq! ( fill ( & adj , & fill_reducing_order ( & xyz , & adj ) ) , graph ) ;
}
#[ test ]
fn a_block_keeps_the_better_ordering ( ) {
let ( xyz , adj ) = hex_grid ( [ 16 , 16 , 16 ] , false , | i , j , k | [ i as f64 , j as f64 , k as f64 ] ) ;
let planes = fill ( & adj , & nested_dissection ( & xyz , & adj ) ) ;
let chosen = fill ( & adj , & fill_reducing_order ( & xyz , & adj ) ) ;
assert! ( chosen < = planes , " chosen {chosen} planes {planes} " ) ;
}
/// Factor entries and factor time (s) of a diagonally dominant 3-unknowns-per-node matrix on
/// the graph `adj`, eliminated in `order`.
fn factor_3dof ( adj : & [ Vec < usize > ] , order : & [ usize ] ) -> ( usize , f64 ) {
let mut place = vec! [ 0 ; adj . len ( ) ] ;
for ( p , & v ) in order . iter ( ) . enumerate ( ) {
place [ v ] = p ;
}
let neq = 3 * adj . len ( ) ;
let mut cols : Vec < Vec < ( usize , f64 ) > > = vec! [ Vec ::new ( ) ; neq ] ;
for ( v , row ) in adj . iter ( ) . enumerate ( ) {
for & w in row {
for a in 0 .. 3 {
for b in 0 .. 3 {
let ( r , c ) = ( 3 * place [ w ] + b , 3 * place [ v ] + a ) ;
if r > c {
continue ;
}
let value = if v = = w {
if a = = b { 30.0 } else { 0.5 }
} else if a = = b {
- 1.0
} else {
0.01 * ( ( v + w ) % 3 ) as f64
} ;
cols [ c ] . push ( ( r , value ) ) ;
}
}
}
}
let mut ptr = vec! [ 0 ] ;
let ( mut idx , mut val ) = ( Vec ::new ( ) , Vec ::new ( ) ) ;
for mut col in cols {
col . sort_by_key ( | e | e . 0 ) ;
idx . extend ( col . iter ( ) . map ( | e | e . 0 ) ) ;
val . extend ( col . iter ( ) . map ( | e | e . 1 ) ) ;
ptr . push ( idx . len ( ) ) ;
}
let upper = crate ::cholesky ::UpperColumns { n : neq , ptr , idx , val } ;
let t = std ::time ::Instant ::now ( ) ;
let factor = crate ::cholesky ::Factor ::new ( & upper , usize ::MAX ) . unwrap ( ) ;
( factor . entries ( ) , t . elapsed ( ) . as_secs_f64 ( ) )
}
/// The measurement table of the ordering choice (fill and time per ordering).
#[ test ]
#[ ignore = " measurement; run with --ignored --nocapture " ]
fn measure_orderings ( ) {
let cases = [
( " cylinder wall 160x2x40 " , cylinder_wall ( ) ) ,
(
" bellows 120x2x120 " ,
hex_grid ( [ 120 , 2 , 120 ] , true , | i , j , k | {
let t = i as f64 / 120.0 * std ::f64 ::consts ::TAU ;
let z = k as f64 ;
// Convolutions of the wall in z, the wall thickness along the normal.
let r = 40.0 + 12.0 * ( z / 120.0 * 8.0 * std ::f64 ::consts ::TAU ) . sin ( ) + j as f64 * 0.5 ;
[ r * t . cos ( ) , r * t . sin ( ) , z * 1.0 ]
} ) ,
) ,
( " block 24^3 " , hex_grid ( [ 24 , 24 , 24 ] , false , | i , j , k | [ i as f64 , j as f64 , k as f64 ] ) ) ,
] ;
for ( name , ( xyz , adj ) ) in cases {
let t = std ::time ::Instant ::now ( ) ;
let planes = nested_dissection ( & xyz , & adj ) ;
let tp = t . elapsed ( ) ;
let t = std ::time ::Instant ::now ( ) ;
let graph = graph_dissection ( & adj ) ;
let tg = t . elapsed ( ) ;
let t = std ::time ::Instant ::now ( ) ;
let ( fp , fg ) = ( fill ( & adj , & planes ) , fill ( & adj , & graph ) ) ;
let tf = t . elapsed ( ) ;
let ( ep , sp ) = factor_3dof ( & adj , & planes ) ;
let ( eg , sg ) = factor_3dof ( & adj , & graph ) ;
println! (
" | {name} | {} | {fp} / {fg} | {ep} / {eg} | {sp:.3} / {sg:.3} | {:.1} / {:.1} / {:.1} | " ,
3 * adj . len ( ) ,
tp . as_secs_f64 ( ) * 1e3 ,
tg . as_secs_f64 ( ) * 1e3 ,
tf . as_secs_f64 ( ) * 1e3
) ;
}
}
}
