faer/sparse/linalg/

amd.rs

//! approximate minimum degree ordering.
use crate::internal_prelude_sp::*;
use crate::{assert, debug_assert};
#[inline]
fn post_tree<'n, I: Index>(
	root: Idx<'n, usize>,
	k: usize,
	child: &mut Array<'n, MaybeIdx<'n, I>>,
	sibling: &Array<'n, MaybeIdx<'n, I>>,
	order: &mut Array<'n, I::Signed>,
	stack: &mut Array<'n, I::Signed>,
) -> usize {
	let N = child.len();
	let I = I::Signed::truncate;
	let mut top = 1usize;
	stack[N.check(0)] = I(*root);
	let mut k = k;
	while top > 0 {
		let i = N.check(stack[N.check(top - 1)].zx());
		if let Some(child_) = child[i].idx() {
			let mut f = child_.zx();
			loop {
				top += 1;
				if let Some(new_f) = sibling[f].idx() {
					f = new_f.zx();
				} else {
					break;
				}
			}
			let mut t = top;
			let mut f = child_.zx();
			loop {
				t -= 1;
				stack[N.check(t)] = I(*f);
				if let Some(new_f) = sibling[f].idx() {
					f = new_f.zx();
				} else {
					break;
				}
			}
			child[i] = MaybeIdx::none();
		} else {
			top -= 1;
			order[i] = I(k);
			k += 1;
		}
	}
	k
}
#[inline]
fn postorder<'n, 'out, I: Index>(
	order: &'out mut Array<'n, I::Signed>,
	etree: &Array<'n, MaybeIdx<'n, I>>,
	nv: &Array<'n, I::Signed>,
	f_size: &Array<'n, I::Signed>,
	stack: &mut MemStack,
) {
	let N = order.len();
	let n = *N;
	if n == 0 {
		return;
	}
	let I = I::Signed::truncate;
	let zero = I(0);
	let none = I(NONE);
	let (ref mut child, stack) =
		stack.collect(repeat_n!(MaybeIdx::<'_, I>::none(), n));
	let (ref mut sibling, stack) =
		stack.collect(repeat_n!(MaybeIdx::<'_, I>::none(), n));
	let (ref mut stack, _) = stack.collect(repeat_n!(I(0), n));
	let child = Array::from_mut(child, N);
	let sibling = Array::from_mut(sibling, N);
	let stack = Array::from_mut(stack, N);
	for j in N.indices().rev() {
		if nv[j] > zero {
			if let Some(parent) = etree[j].idx() {
				let parent = parent.zx();
				sibling[j] = child[parent];
				child[parent] = MaybeIdx::from_index(j.truncate());
			}
		}
	}
	for i in N.indices() {
		if nv[i] > zero {
			if let Some(child_) = child[i].idx() {
				let child_ = child_.zx();
				let mut fprev = MaybeIdx::<'n>::none();
				let mut bigfprev = MaybeIdx::<'n>::none();
				let mut bigf = MaybeIdx::<'n>::none();
				let mut maxfrsize = none;
				let mut f = MaybeIdx::from_index(child_);
				while let Some(f_) = f.idx() {
					let frsize = f_size[f_];
					if frsize >= maxfrsize {
						maxfrsize = frsize;
						bigfprev = fprev;
						bigf = f;
					}
					fprev = f;
					f = sibling[f_].sx();
				}
				let bigf = bigf.idx().unwrap();
				let fnext = sibling[bigf];
				if let Some(fnext) = fnext.idx() {
					if let Some(bigfprev) = bigfprev.idx() {
						sibling[bigfprev] = MaybeIdx::from_index(fnext);
					} else {
						child[i] = MaybeIdx::from_index(fnext);
					}
					let fprev = fprev.idx().unwrap();
					sibling[bigf] = MaybeIdx::none();
					sibling[fprev] = MaybeIdx::from_index(bigf.truncate());
				}
			}
		}
	}
	order.as_mut().fill(none);
	let mut k = 0usize;
	for i in N.indices() {
		if etree[i].idx().is_none() && nv[i] > zero {
			k = post_tree(i, k, child, sibling, order, stack);
		}
	}
}
#[inline(always)]
fn flip<I: SignedIndex>(i: I) -> I {
	-I::truncate(2) - i
}
#[inline]
fn clear_flag<I: SignedIndex>(wflg: I, wbig: I, w: &mut [I]) -> I {
	let I = I::truncate;
	let zero = I(0);
	let one = I(1);
	if wflg < I(2) || wflg >= wbig {
		for x in w {
			if *x != zero {
				*x = one;
			}
		}
		return I(2);
	}
	wflg
}
#[allow(clippy::comparison_chain)]
fn amd_2<I: Index>(
	pe: &mut [I::Signed],
	iw: &mut [I::Signed],
	len: &mut [I::Signed],
	pfree: usize,
	next: &mut [I::Signed],
	last: &mut [I::Signed],
	control: Control,
	stack: &mut MemStack,
) -> FlopCount {
	let n = pe.len();
	assert!(n < I::Signed::MAX.zx());
	let mut pfree = pfree;
	let iwlen = iw.len();
	let I = I::Signed::truncate;
	let none = I(NONE);
	let zero = I(0);
	let one = I(1);
	let alpha = control.dense;
	let aggressive = control.aggressive;
	let mut mindeg = 0usize;
	let mut ncmpa = 0usize;
	let mut lemax = 0usize;
	let mut ndiv = 0.0;
	let mut nms_lu = 0.0;
	let mut nms_ldl = 0.0;
	let mut nel = 0usize;
	let mut me = none;
	let dense = if alpha < 0.0 {
		n - 2
	} else {
		(alpha * ComplexFieldExt::sqrt(&(n as f64))) as usize
	};
	let dense = Ord::max(dense, 16);
	let dense = Ord::min(dense, n);
	let (mut w, stack) = stack.collect(repeat_n!(one, n));
	let (mut nv, stack) = stack.collect(repeat_n!(one, n));
	let (mut elen, stack) = stack.collect(repeat_n!(zero, n));
	let nv = &mut *nv;
	let elen = &mut *elen;
	let w = &mut *w;
	let wbig = I::Signed::MAX - I(n);
	let mut wflg = I(0);
	let mut ndense = zero;
	{
		let (mut head, stack) = stack.collect(repeat_n!(none, n));
		let (mut degree, _) = stack.collect(repeat_n!(none, n));
		let head = &mut *head;
		let degree = &mut *degree;
		last.fill(none);
		head.fill(none);
		next.fill(none);
		degree.copy_from_slice(len);
		for i in 0..n {
			let deg = degree[i].zx();
			if deg == 0 {
				elen[i] = flip(one);
				nel += 1;
				pe[i] = none;
				w[i] = zero;
			} else if deg > dense {
				ndense += one;
				nv[i] = zero;
				elen[i] = none;
				pe[i] = none;
				nel += 1;
			} else {
				let inext = head[deg];
				if inext != none {
					last[inext.zx()] = I(i);
				}
				next[i] = inext;
				head[deg] = I(i);
			}
		}
		while nel < n {
			let mut deg = mindeg;
			while deg < n {
				me = head[deg];
				if me != none {
					break;
				}
				deg += 1;
			}
			mindeg = deg;
			let me = me.zx();
			let inext = next[me];
			if inext != none {
				last[inext.zx()] = none;
			}
			head[deg] = inext;
			let elenme = elen[me];
			let mut nvpiv = nv[me];
			nel += nvpiv.zx();
			nv[me] = -nvpiv;
			let mut degme = 0usize;
			let mut pme1;
			let mut pme2;
			if elenme == zero {
				pme1 = pe[me];
				pme2 = pme1 - one;
				for p in pme1.zx()..(pme1 + len[me]).zx() {
					let i = iw[p].zx();
					let nvi = nv[i];
					if nvi > zero {
						degme += nvi.zx();
						nv[i] = -nvi;
						pme2 += one;
						iw[pme2.zx()] = I(i);
						let ilast = last[i];
						let inext = next[i];
						if inext != none {
							last[inext.zx()] = ilast;
						}
						if ilast != none {
							next[ilast.zx()] = inext;
						} else {
							head[degree[i].zx()] = inext;
						}
					}
				}
			} else {
				let mut p = pe[me].zx();
				pme1 = I(pfree);
				let slenme = (len[me] - elenme).zx();
				for knt1 in 1..elenme.zx() + 2 {
					let e;
					let mut pj;
					let ln;
					if I(knt1) > elenme {
						e = me;
						pj = I(p);
						ln = slenme;
					} else {
						e = iw[p].zx();
						p += 1;
						pj = pe[e];
						ln = len[e].zx();
					}
					for knt2 in 1..ln + 1 {
						let i = iw[pj.zx()].zx();
						pj += one;
						let nvi = nv[i];
						if nvi > zero {
							if pfree >= iwlen {
								pe[me] = I(p);
								len[me] -= I(knt1);
								if len[me] == zero {
									pe[me] = none;
								}
								pe[e] = pj;
								len[e] = I(ln - knt2);
								if len[e] == zero {
									pe[e] = none;
								}
								ncmpa += 1;
								for (j, pe) in pe.iter_mut().enumerate() {
									let pn = *pe;
									if pn >= zero {
										let pn = pn.zx();
										*pe = iw[pn];
										iw[pn] = flip(I(j));
									}
								}
								let mut psrc = 0usize;
								let mut pdst = 0usize;
								let pend = pme1.zx();
								while psrc < pend {
									let j = flip(iw[psrc]);
									psrc += 1;
									if j >= zero {
										let j = j.zx();
										iw[pdst] = pe[j];
										pe[j] = I(pdst);
										pdst += 1;
										let lenj = len[j].zx();
										if lenj > 0 {
											iw.copy_within(
												psrc..psrc + lenj - 1,
												pdst,
											);
											psrc += lenj - 1;
											pdst += lenj - 1;
										}
									}
								}
								let p1 = pdst;
								iw.copy_within(pme1.zx()..pfree, pdst);
								pdst += pfree - pme1.zx();
								pme1 = I(p1);
								pfree = pdst;
								pj = pe[e];
								p = pe[me].zx();
							}
							degme += nvi.zx();
							nv[i] = -nvi;
							iw[pfree] = I(i);
							pfree += 1;
							let ilast = last[i];
							let inext = next[i];
							if inext != none {
								last[inext.zx()] = ilast;
							}
							if ilast != none {
								next[ilast.zx()] = inext;
							} else {
								head[degree[i].zx()] = inext;
							}
						}
					}
					if e != me {
						pe[e] = flip(I(me));
						w[e] = zero;
					}
				}
				pme2 = I(pfree - 1);
			}
			degree[me] = I(degme);
			pe[me] = pme1;
			len[me] = pme2 - pme1 + one;
			elen[me] = flip(nvpiv + I(degme));
			wflg = clear_flag(wflg, wbig, w);
			assert!(pme1 >= zero);
			for pme in pme1.zx()..(pme2 + one).zx() {
				let i = iw[pme].zx();
				let eln = elen[i];
				if eln > zero {
					let nvi = -nv[i];
					let wnvi = wflg - nvi;
					for iw in iw[pe[i].zx()..][..eln.zx()].iter() {
						let e = iw.zx();
						let mut we = w[e];
						if we >= wflg {
							we -= nvi;
						} else if we != zero {
							we = degree[e] + wnvi;
						}
						w[e] = we;
					}
				}
			}
			for pme in pme1.zx()..(pme2 + one).zx() {
				let i = iw[pme].zx();
				let p1 = pe[i].zx();
				let p2 = p1 + elen[i].zx();
				let mut pn = p1;
				let mut hash = 0usize;
				deg = 0usize;
				if aggressive {
					for p in p1..p2 {
						let e = iw[p].zx();
						let we = w[e];
						if we != zero {
							let dext = we - wflg;
							if dext > zero {
								deg += dext.zx();
								iw[pn] = I(e);
								pn += 1;
								hash = hash.wrapping_add(e);
							} else {
								pe[e] = flip(I(me));
								w[e] = zero;
							}
						}
					}
				} else {
					for p in p1..p2 {
						let e = iw[p].zx();
						let we = w[e];
						if we != zero {
							let dext = (we - wflg).zx();
							deg += dext;
							iw[pn] = I(e);
							pn += 1;
							hash = hash.wrapping_add(e);
						}
					}
				}
				elen[i] = I(pn - p1 + 1);
				let p3 = pn;
				let p4 = p1 + len[i].zx();
				for p in p2..p4 {
					let j = iw[p].zx();
					let nvj = nv[j];
					if nvj > zero {
						deg += nvj.zx();
						iw[pn] = I(j);
						pn += 1;
						hash = hash.wrapping_add(j);
					}
				}
				if elen[i] == one && p3 == pn {
					pe[i] = flip(I(me));
					let nvi = -nv[i];
					assert!(nvi >= zero);
					degme -= nvi.zx();
					nvpiv += nvi;
					nel += nvi.zx();
					nv[i] = zero;
					elen[i] = none;
				} else {
					degree[i] = Ord::min(degree[i], I(deg));
					iw[pn] = iw[p3];
					iw[p3] = iw[p1];
					iw[p1] = I(me);
					len[i] = I(pn - p1 + 1);
					hash %= n;
					let j = head[hash];
					if j <= none {
						next[i] = flip(j);
						head[hash] = flip(I(i));
					} else {
						next[i] = last[j.zx()];
						last[j.zx()] = I(i);
					}
					last[i] = I(hash);
				}
			}
			degree[me] = I(degme);
			lemax = Ord::max(lemax, degme);
			wflg += I(lemax);
			wflg = clear_flag(wflg, wbig, w);
			for pme in pme1.zx()..(pme2 + one).zx() {
				let mut i = iw[pme].zx();
				if nv[i] < zero {
					let hash = last[i].zx();
					let j = head[hash];
					if j == none {
						i = NONE;
					} else if j < none {
						i = flip(j).zx();
						head[hash] = none;
					} else {
						i = last[j.zx()].sx();
						last[j.zx()] = none;
					}
					while i != NONE && next[i] != none {
						let ln = len[i];
						let eln = elen[i];
						for p in (pe[i] + one).zx()..(pe[i] + ln).zx() {
							w[iw[p].zx()] = wflg;
						}
						let mut jlast = i;
						let mut j = next[i].sx();
						while j != NONE {
							let mut ok = len[j] == ln && elen[j] == eln;
							if ok {
								for p in (pe[j] + one).zx()..(pe[j] + ln).zx() {
									if w[iw[p].zx()] != wflg {
										ok = false;
										break;
									}
								}
							}
							if ok {
								pe[j] = flip(I(i));
								nv[i] += nv[j];
								nv[j] = zero;
								elen[j] = none;
								j = next[j].sx();
								next[jlast] = I(j);
							} else {
								jlast = j;
								j = next[j].sx();
							}
						}
						wflg += one;
						i = next[i].sx();
					}
				}
			}
			let mut p = pme1.zx();
			let nleft = n - nel;
			for pme in pme1.zx()..(pme2 + one).zx() {
				let i = iw[pme].zx();
				let nvi = -nv[i];
				if nvi > zero {
					nv[i] = nvi;
					deg = degree[i].zx() + degme - nvi.zx();
					deg = Ord::min(deg, nleft - nvi.zx());
					let inext = head[deg];
					if inext != none {
						last[inext.zx()] = I(i);
					}
					next[i] = inext;
					last[i] = none;
					head[deg] = I(i);
					mindeg = Ord::min(mindeg, deg);
					degree[i] = I(deg);
					iw[p] = I(i);
					p += 1;
				}
			}
			nv[me] = nvpiv;
			len[me] = I(p) - pme1;
			if len[me] == zero {
				pe[me] = none;
				w[me] = zero;
			}
			if elenme != zero {
				pfree = p;
			}
			{
				let f = nvpiv.sx() as isize as f64;
				let r = degme as f64 + ndense.sx() as isize as f64;
				let lnzme = f * r + (f - 1.0) * f / 2.0;
				ndiv += lnzme;
				let s = f * r * r
					+ r * (f - 1.0) * f
					+ (f - 1.0) * f * (2.0 * f - 1.0) / 6.0;
				nms_lu += s;
				nms_ldl += (s + lnzme) / 2.0;
			}
		}
		{
			let f = ndense.sx() as isize as f64;
			let lnzme = (f - 1.0) * f / 2.0;
			ndiv += lnzme;
			let s = (f - 1.0) * f * (2.0 * f - 1.0) / 6.0;
			nms_lu += s;
			nms_ldl += (s + lnzme) / 2.0;
		}
		for pe in pe.iter_mut() {
			*pe = flip(*pe);
		}
		for elen in elen.iter_mut() {
			*elen = flip(*elen);
		}
		for i in 0..n {
			if nv[i] == zero {
				let mut j = pe[i].sx();
				if j == NONE {
					continue;
				}
				while nv[j] == zero {
					j = pe[j].zx();
				}
				let e = I(j);
				let mut j = i;
				while nv[j] == zero {
					let jnext = pe[j];
					pe[j] = e;
					j = jnext.zx();
				}
			}
		}
	}
	with_dim!(N, n);
	postorder(
		Array::from_mut(w, N),
		Array::from_ref(MaybeIdx::<'_, I>::from_slice_ref_checked(pe, N), N),
		Array::from_ref(nv, N),
		Array::from_ref(elen, N),
		stack,
	);
	let (ref mut head, _) = stack.collect(repeat_n!(none, n));
	next.fill(none);
	for (e, &k) in w.iter().enumerate() {
		if k != none {
			head[k.zx()] = I(e);
		}
	}
	nel = 0;
	for &e in head.iter() {
		if e == none {
			break;
		}
		let e = e.zx();
		next[e] = I(nel);
		nel += nv[e].zx();
	}
	assert!(nel == n - ndense.zx());
	for i in 0..n {
		if nv[i] == zero {
			let e = pe[i];
			if e != none {
				let e = e.zx();
				next[i] = next[e];
				next[e] += one;
			} else {
				next[i] = I(nel);
				nel += 1;
			}
		}
	}
	assert!(nel == n);
	for i in 0..n {
		last[next[i].zx()] = I(i);
	}
	let _ = ncmpa;
	FlopCount {
		n_div: ndiv,
		n_mult_subs_ldl: nms_ldl,
		n_mult_subs_lu: nms_lu,
	}
}
#[allow(clippy::comparison_chain)]
fn amd_1<I: Index>(
	perm: &mut [I::Signed],
	perm_inv: &mut [I::Signed],
	A: SymbolicSparseColMatRef<'_, I>,
	len: &mut [I::Signed],
	iwlen: usize,
	control: Control,
	stack: &mut MemStack,
) -> FlopCount {
	let n = perm.len();
	let I = I::Signed::truncate;
	let zero = I(0);
	let one = I(1);
	let (p_e, stack) = unsafe { stack.make_raw::<I::Signed>(n) };
	let (s_p, stack) = unsafe { stack.make_raw::<I::Signed>(n) };
	let (i_w, stack) = unsafe { stack.make_raw::<I::Signed>(iwlen) };
	let mut pfree = zero;
	for j in 0..n {
		p_e[j] = pfree;
		s_p[j] = pfree;
		pfree += len[j];
	}
	let pfree = pfree.zx();
	assert!(iwlen >= pfree + n);
	{
		with_dim!(N, n);
		let (t_p, _) = unsafe { stack.make_raw::<I::Signed>(n) };
		let A = A.as_shape(N, N);
		let s_p = Array::from_mut(s_p, N);
		let t_p = Array::from_mut(t_p, N);
		for k in N.indices() {
			let mut seen = zero;
			let mut j_prev = usize::MAX;
			for j in A.row_idx_of_col(k) {
				if j_prev == *j {
					j_prev = *j;
					continue;
				}
				j_prev = *j;
				if j < k {
					i_w[s_p[j].zx()] = I(*k);
					s_p[j] += one;
					i_w[s_p[k].zx()] = I(*j);
					s_p[k] += one;
					seen += one;
				} else {
					if j == k {
						seen += one;
					}
					break;
				}
				let mut seen_j = zero;
				let mut i_prev = usize::MAX;
				for i in &A.row_idx_of_col_raw(j)[t_p[j].zx()..] {
					if i_prev == *i.zx() {
						i_prev = *i.zx();
						continue;
					}
					i_prev = *i.zx();
					let i = i.zx();
					if i < k {
						i_w[s_p[i].zx()] = I(*j);
						s_p[i] += one;
						i_w[s_p[j].zx()] = I(*i);
						s_p[j] += one;
						seen_j += one;
					} else {
						if i == k {
							seen_j += one;
						}
						break;
					}
				}
				t_p[j] += seen_j;
			}
			t_p[k] = seen;
		}
		for j in N.indices() {
			let mut i_prev = usize::MAX;
			for i in &A.row_idx_of_col_raw(j)[t_p[j].zx()..] {
				if i_prev == *i.zx() {
					i_prev = *i.zx();
					continue;
				}
				i_prev = *i.zx();
				let i = i.zx();
				i_w[s_p[i].zx()] = I(*j);
				s_p[i] += one;
				i_w[s_p[j].zx()] = I(*i);
				s_p[j] += one;
			}
		}
	}
	debug_assert!(s_p[n - 1] == I(pfree));
	amd_2::<I>(p_e, i_w, len, pfree, perm_inv, perm, control, stack)
}
fn preprocess<'out, I: Index>(
	new_col_ptr: &'out mut [I],
	new_row_idx: &'out mut [I],
	A: SymbolicSparseColMatRef<'_, I>,
	stack: &mut MemStack,
) -> SymbolicSparseColMatRef<'out, I> {
	let n = A.nrows();
	with_dim!(N, n);
	let I = I::Signed::truncate;
	let zero = I(0);
	let one = I(1);
	let (ref mut w, stack) = stack.collect(repeat_n!(I(0), n));
	let (ref mut flag, _) = stack.collect(repeat_n!(I(NONE), n));
	let w = Array::from_mut(w, N);
	let flag = Array::from_mut(flag, N);
	let A = A.as_shape(N, N);
	for j in N.indices() {
		let j_ = I(*j);
		for i in A.row_idx_of_col(j) {
			if flag[i] != j_ {
				w[i] += one;
				flag[i] = j_;
			}
		}
	}
	new_col_ptr[0] = I::from_signed(zero);
	for (i, [r, r_next]) in iter::zip(
		N.indices(),
		windows2(Cell::as_slice_of_cells(Cell::from_mut(new_col_ptr))),
	) {
		r_next.set(r.get() + I::from_signed(w[i]));
	}
	w.as_mut()
		.copy_from_slice(bytemuck::cast_slice(&new_col_ptr[..n]));
	flag.as_mut().fill(I(NONE));
	for j in N.indices() {
		let j_ = I(*j);
		for i in A.row_idx_of_col(j) {
			if flag[i] != j_ {
				new_row_idx[w[i].zx()] = I::from_signed(j_);
				w[i] += one;
				flag[i] = j_;
			}
		}
	}
	unsafe {
		SymbolicSparseColMatRef::new_unchecked(
			n,
			n,
			&*new_col_ptr,
			None,
			&new_row_idx[..new_col_ptr[n].zx()],
		)
	}
}
#[allow(clippy::comparison_chain)]
fn aat<I: Index>(
	len: &mut [I::Signed],
	A: SymbolicSparseColMatRef<'_, I>,
	stack: &mut MemStack,
) -> Result<usize, FaerError> {
	{
		with_dim!(N, A.nrows());
		let I = I::Signed::truncate;
		let zero = I(0);
		let one = I(1);
		let A = A.as_shape(N, N);
		let n = *N;
		let t_p = &mut *unsafe { stack.make_raw::<I::Signed>(n).0 };
		let len = Array::from_mut(len, N);
		let t_p = Array::from_mut(t_p, N);
		for k in N.indices() {
			let mut seen = zero;
			let mut j_prev = usize::MAX;
			for j in A.row_idx_of_col(k) {
				if j_prev == *j {
					j_prev = *j;
					continue;
				}
				j_prev = *j;
				if j < k {
					seen += one;
					len[j] += one;
					len[k] += one;
				} else {
					if j == k {
						seen += one;
					}
					break;
				}
				let mut seen_j = zero;
				let mut i_prev = usize::MAX;
				for i in &A.row_idx_of_col_raw(j)[t_p[j].zx()..] {
					if i_prev == *i.zx() {
						i_prev = *i.zx();
						continue;
					}
					i_prev = *i.zx();
					let i = i.zx();
					if i < k {
						len[i] += one;
						len[j] += one;
						seen_j += one;
					} else {
						if i == k {
							seen_j += one;
						}
						break;
					}
				}
				t_p[j] += seen_j;
			}
			t_p[k] = seen;
		}
		for j in N.indices() {
			let mut i_prev = usize::MAX;
			for i in &A.row_idx_of_col_raw(j)[t_p[j].zx()..] {
				if i_prev == *i.zx() {
					i_prev = *i.zx();
					continue;
				}
				i_prev = *i.zx();
				len[i.zx()] += one;
				len[j] += one;
			}
		}
	}
	let nzaat = I::Signed::sum_nonnegative(len).map(I::from_signed);
	nzaat.ok_or(FaerError::IndexOverflow).map(I::zx)
}
/// computes the layout of required workspace for computing the amd ordering of
/// a sorted matrix
pub fn order_scratch<I: Index>(n: usize, nnz_upper: usize) -> StackReq {
	let n_scratch = StackReq::new::<I>(n);
	let nzaat = nnz_upper.checked_mul(2).unwrap_or(usize::MAX);
	StackReq::all_of(&[
		n_scratch,
		StackReq::new::<I>(nzaat.checked_add(nzaat / 5).unwrap_or(usize::MAX)),
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
		n_scratch,
	])
}
/// computes the layout of required workspace for computing the amd ordering of
/// an unsorted matrix
pub fn order_maybe_unsorted_scratch<I: Index>(
	n: usize,
	nnz_upper: usize,
) -> StackReq {
	StackReq::all_of(&[
		order_scratch::<I>(n, nnz_upper),
		StackReq::new::<I>(n + 1),
		StackReq::new::<I>(nnz_upper),
	])
}
/// computes the approximate minimum degree ordering for reducing the fill-in
/// during the sparse cholesky factorization of a matrix with the sparsity
/// pattern of $A + A^\top$
pub fn order<I: Index>(
	perm: &mut [I],
	perm_inv: &mut [I],
	A: SymbolicSparseColMatRef<'_, I>,
	control: Control,
	stack: &mut MemStack,
) -> Result<FlopCount, FaerError> {
	let n = perm.len();
	if n == 0 {
		return Ok(FlopCount {
			n_div: 0.0,
			n_mult_subs_ldl: 0.0,
			n_mult_subs_lu: 0.0,
		});
	}
	let (ref mut len, stack) =
		stack.collect(repeat_n!(I::Signed::truncate(0), n));
	let nzaat = aat(len, A, stack)?;
	let iwlen = nzaat
		.checked_add(nzaat / 5)
		.and_then(|x| x.checked_add(n))
		.ok_or(FaerError::IndexOverflow)?;
	Ok(amd_1::<I>(
		bytemuck::cast_slice_mut(perm),
		bytemuck::cast_slice_mut(perm_inv),
		A,
		len,
		iwlen,
		control,
		stack,
	))
}
/// computes the approximate minimum degree ordering for reducing the fill-in
/// during the sparse cholesky factorization of a matrix with the sparsity
/// pattern of $A + A^\top$
///
/// # note
/// allows unsorted matrices
pub fn order_maybe_unsorted<I: Index>(
	perm: &mut [I],
	perm_inv: &mut [I],
	A: SymbolicSparseColMatRef<'_, I>,
	control: Control,
	stack: &mut MemStack,
) -> Result<FlopCount, FaerError> {
	let n = perm.len();
	if n == 0 {
		return Ok(FlopCount {
			n_div: 0.0,
			n_mult_subs_ldl: 0.0,
			n_mult_subs_lu: 0.0,
		});
	}
	let nnz = A.compute_nnz();
	let (new_col_ptr, stack) = unsafe { stack.make_raw::<I>(n + 1) };
	let (new_row_idx, stack) = unsafe { stack.make_raw::<I>(nnz) };
	let A = preprocess(new_col_ptr, new_row_idx, A, stack);
	order(perm, perm_inv, A, control, stack)
}
/// tuning parameters for the amd implementation
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Control {
	/// "dense" if degree > dense * sqrt(n)
	pub dense: f64,
	/// do aggressive absorption
	pub aggressive: bool,
}
impl Default for Control {
	#[inline]
	fn default() -> Self {
		Self {
			dense: 10.0,
			aggressive: true,
		}
	}
}
/// flop count of the ldlt and lu factorizations if the provided ordering is
/// used
#[derive(Default, Debug, Copy, Clone, PartialEq)]
pub struct FlopCount {
	/// number of division
	pub n_div: f64,
	/// number of multiplications and subtractions for the ldlt factorization
	pub n_mult_subs_ldl: f64,
	/// number of multiplications and subtractions for the lu factorization
	pub n_mult_subs_lu: f64,
}

#[cfg(test)]
mod tests {
	use super::*;
	use alloc::alloc::Layout;
	use core::slice;
	use std::alloc::System;

	struct GlobalAlloc;
	#[global_allocator]
	static GLOBAL_ALLOC: GlobalAlloc = GlobalAlloc;
	unsafe impl core::alloc::GlobalAlloc for GlobalAlloc {
		unsafe fn alloc(&self, layout: Layout) -> *mut u8 {
			let result = unsafe { System.alloc(layout) };
			// Fill the memory with garbage.
			unsafe {
				slice::from_raw_parts_mut(result, layout.size()).fill(0x7F)
			};
			result
		}

		unsafe fn realloc(
			&self,
			ptr: *mut u8,
			layout: Layout,
			new_size: usize,
		) -> *mut u8 {
			unsafe { System.realloc(ptr, layout, new_size) }
		}

		unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) {
			unsafe { System.dealloc(ptr, layout) }
		}
	}

	#[test]
	fn test_gh_258() {
		use faer::Side;
		use faer::sparse::linalg::cholesky::{
			SymmetricOrdering, factorize_symbolic_cholesky,
		};
		use faer::sparse::{SparseColMat, Triplet};
		use std::fs::read_to_string;

		let data =
			read_to_string("test_data/sparse_cholesky/gh_258.txt").unwrap();
		let mut triplets = Vec::<Triplet<u32, u32, f64>>::new();
		for line in data.lines() {
			let [i, j, a] = Vec::from_iter(line.split(" ")).try_into().unwrap();
			triplets.push(Triplet::new(
				i.parse().unwrap(),
				j.parse().unwrap(),
				a.parse().unwrap(),
			));
		}
		let mat =
			SparseColMat::try_new_from_triplets(52, 52, &triplets).unwrap();
		_ = factorize_symbolic_cholesky(
			mat.symbolic(),
			Side::Lower,
			SymmetricOrdering::Amd,
			<_>::default(),
		);
	}
}