Module cholesky 

computes the Cholesky decomposition (either $LL^\top$, $LTL^\top$, or $LBL^\top$) of a given sparse matrix. see crate::linalg::cholesky for more info.

the entry point in this module is SymbolicCholesky and factorize_symbolic_cholesky.

note

the functions in this module accept unsorted input, producing a sorted decomposition factor (simplicial).

example (low level api)

simplicial:

fn simplicial_cholesky() -> Result<(), faer::sparse::FaerError> {
	use faer::Par;
	use faer::dyn_stack::{MemBuffer, MemStack, StackReq};
	use faer::reborrow::*;

	use faer::linalg::cholesky::ldlt::factor::LdltRegularization;
	use faer::linalg::cholesky::llt::factor::LltRegularization;

	use faer::sparse::linalg::amd;
	use faer::sparse::linalg::cholesky::simplicial;
	use faer::sparse::{CreationError, SparseColMat, SymbolicSparseColMat, Triplet};

	use rand::prelude::*;

	// the simplicial cholesky api takes an upper triangular matrix as input, to be
	// interpreted as self-adjoint.
	let dim = 4;
	let A_upper = match SparseColMat::<usize, f64>::try_new_from_triplets(
		dim,
		dim,
		&[
			// diagonal entries
			Triplet::new(0, 0, 10.0),
			Triplet::new(1, 1, 11.0),
			Triplet::new(2, 2, 12.0),
			Triplet::new(3, 3, 13.0),
			// non diagonal entries
			Triplet::new(0, 1, 1.0),
			Triplet::new(0, 3, 1.5),
			Triplet::new(1, 3, -3.2),
		],
	) {
		Ok(A) => Ok(A),
		Err(CreationError::Generic(err)) => Err(err),
		Err(CreationError::OutOfBounds { .. }) => panic!(),
	}?;

	let mut A = faer::sparse::ops::add(A_upper.rb(), A_upper.to_row_major()?.rb().transpose())?;
	for i in 0..dim {
		A[(i, i)] /= 2.0;
	}

	let A_nnz = A_upper.compute_nnz();
	let mut rhs = Vec::new();
	let mut rng = StdRng::seed_from_u64(0);
	rhs.try_reserve_exact(dim)?;
	rhs.resize_with(dim, || rng.gen::<f64>());

	let mut sol = Vec::new();
	sol.try_reserve_exact(dim)?;
	sol.resize(dim, 0.0f64);

	let rhs = faer::MatRef::from_column_major_slice(&rhs, dim, 1);
	let mut sol = faer::MatMut::from_column_major_slice_mut(&mut sol, dim, 1);

	// optional: fill reducing permutation
	let (perm, perm_inv) = {
		let mut perm = Vec::new();
		let mut perm_inv = Vec::new();
		perm.try_reserve_exact(dim)?;
		perm_inv.try_reserve_exact(dim)?;
		perm.resize(dim, 0usize);
		perm_inv.resize(dim, 0usize);

		let mut mem = MemBuffer::try_new(amd::order_scratch::<usize>(dim, A_nnz))?;
		amd::order(
			&mut perm,
			&mut perm_inv,
			A_upper.symbolic(),
			amd::Control::default(),
			MemStack::new(&mut mem),
		)?;

		(perm, perm_inv)
	};

	let perm = unsafe { faer::perm::PermRef::new_unchecked(&perm, &perm_inv, dim) };

	let A_perm_upper = {
		let mut A_perm_col_ptrs = Vec::new();
		let mut A_perm_row_indices = Vec::new();
		let mut A_perm_values = Vec::new();

		A_perm_col_ptrs.try_reserve_exact(dim + 1)?;
		A_perm_col_ptrs.resize(dim + 1, 0usize);
		A_perm_row_indices.try_reserve_exact(A_nnz)?;
		A_perm_row_indices.resize(A_nnz, 0usize);
		A_perm_values.try_reserve_exact(A_nnz)?;
		A_perm_values.resize(A_nnz, 0.0f64);

		let mut mem = MemBuffer::try_new(faer::sparse::utils::permute_self_adjoint_scratch::<
			usize,
		>(dim))?;
		faer::sparse::utils::permute_self_adjoint_to_unsorted(
			&mut A_perm_values,
			&mut A_perm_col_ptrs,
			&mut A_perm_row_indices,
			A_upper.rb(),
			perm,
			faer::Side::Upper,
			faer::Side::Upper,
			MemStack::new(&mut mem),
		);

		SparseColMat::<usize, f64>::new(
			unsafe {
				SymbolicSparseColMat::new_unchecked(
					dim,
					dim,
					A_perm_col_ptrs,
					None,
					A_perm_row_indices,
				)
			},
			A_perm_values,
		)
	};

	// symbolic analysis
	let symbolic = {
		let mut mem = MemBuffer::try_new(StackReq::any_of(&[
			simplicial::prefactorize_symbolic_cholesky_scratch::<usize>(dim, A_nnz),
			simplicial::factorize_simplicial_symbolic_cholesky_scratch::<usize>(dim),
		]))?;
		let stack = MemStack::new(&mut mem);

		let mut etree = Vec::new();
		let mut col_counts = Vec::new();
		etree.try_reserve_exact(dim)?;
		etree.resize(dim, 0isize);
		col_counts.try_reserve_exact(dim)?;
		col_counts.resize(dim, 0usize);

		simplicial::prefactorize_symbolic_cholesky(
			&mut etree,
			&mut col_counts,
			A_perm_upper.symbolic(),
			stack,
		);
		simplicial::factorize_simplicial_symbolic_cholesky(
			A_perm_upper.symbolic(),
			// SAFETY: `etree` was filled correctly by
			// `simplicial::prefactorize_symbolic_cholesky`.
			unsafe { simplicial::EliminationTreeRef::from_inner(&etree) },
			&col_counts,
			stack,
		)?
	};

	// numerical factorization
	let mut mem = MemBuffer::try_new(StackReq::all_of(&[
		simplicial::factorize_simplicial_numeric_llt_scratch::<usize, f64>(dim),
		simplicial::factorize_simplicial_numeric_ldlt_scratch::<usize, f64>(dim),
		faer::perm::permute_rows_in_place_scratch::<usize, f64>(dim, 1),
		symbolic.solve_in_place_scratch::<f64>(dim),
	]))?;
	let stack = MemStack::new(&mut mem);

	// numerical llt factorization
	{
		let mut L_values = Vec::new();
		L_values.try_reserve_exact(symbolic.len_val())?;
		L_values.resize(symbolic.len_val(), 0.0f64);

		match simplicial::factorize_simplicial_numeric_llt::<usize, f64>(
			&mut L_values,
			A_perm_upper.rb(),
			LltRegularization::default(),
			&symbolic,
			stack,
		) {
			Ok(_) => {},
			Err(err) => panic!("matrix is not positive definite: {err}"),
		};

		let llt = simplicial::SimplicialLltRef::<'_, usize, f64>::new(&symbolic, &L_values);

		sol.copy_from(rhs);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm, stack);
		llt.solve_in_place_with_conj(faer::Conj::No, sol.rb_mut(), faer::Par::Seq, stack);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm.inverse(), stack);

		assert!((&A * &sol - &rhs).norm_max() <= 1e-14);
	}

	// numerical ldlt factorization
	{
		let mut L_values = Vec::new();
		L_values.try_reserve_exact(symbolic.len_val())?;
		L_values.resize(symbolic.len_val(), 0.0f64);

		simplicial::factorize_simplicial_numeric_ldlt::<usize, f64>(
			&mut L_values,
			A_perm_upper.rb(),
			LdltRegularization::default(),
			&symbolic,
			stack,
		);

		let ldlt = simplicial::SimplicialLdltRef::<'_, usize, f64>::new(&symbolic, &L_values);

		sol.copy_from(rhs);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm, stack);
		ldlt.solve_in_place_with_conj(faer::Conj::No, sol.rb_mut(), faer::Par::Seq, stack);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm.inverse(), stack);

		assert!((&A * &sol - &rhs).norm_max() <= 1e-14);
	}
	Ok(())
}
simplicial_cholesky().unwrap()

supernodal:

fn supernodal_cholesky() -> Result<(), faer::sparse::FaerError> {
	use faer::Par;
	use faer::dyn_stack::{MemBuffer, MemStack, StackReq};
	use faer::reborrow::*;

	use faer::linalg::cholesky::ldlt::factor::LdltRegularization;

	use faer::sparse::linalg::amd;
	use faer::sparse::linalg::cholesky::{simplicial, supernodal};
	use faer::sparse::{CreationError, SparseColMat, SymbolicSparseColMat, Triplet};
	use rand::prelude::*;

	// the supernodal cholesky api takes a lower triangular matrix as input, to be
	// interpreted as self-adjoint.
	let dim = 8;
	let A_lower = match SparseColMat::<usize, f64>::try_new_from_triplets(
		dim,
		dim,
		&[
			// diagonal entries
			Triplet::new(0, 0, 11.0),
			Triplet::new(1, 1, 11.0),
			Triplet::new(2, 2, 12.0),
			Triplet::new(3, 3, 13.0),
			Triplet::new(4, 4, 14.0),
			Triplet::new(5, 5, 16.0),
			Triplet::new(6, 6, 16.0),
			Triplet::new(7, 7, 16.0),
			// non diagonal entries
			Triplet::new(1, 0, 10.0),
			Triplet::new(3, 0, 10.5),
			Triplet::new(4, 0, 10.0),
			Triplet::new(7, 0, 10.5),
			Triplet::new(3, 1, 10.5),
			Triplet::new(4, 1, 10.0),
			Triplet::new(7, 1, 10.5),
			Triplet::new(3, 2, 10.5),
			Triplet::new(4, 2, 10.0),
			Triplet::new(7, 2, 10.0),
		],
	) {
		Ok(A) => Ok(A),
		Err(CreationError::Generic(err)) => Err(err),
		Err(CreationError::OutOfBounds { .. }) => panic!(),
	}?;

	let mut A = faer::sparse::ops::add(A_lower.rb(), A_lower.to_row_major()?.rb().transpose())?;
	for i in 0..dim {
		A[(i, i)] /= 2.0;
	}

	let A_nnz = A_lower.compute_nnz();
	let mut rhs = Vec::new();
	let mut rng = StdRng::seed_from_u64(0);
	rhs.try_reserve_exact(dim)?;
	rhs.resize_with(dim, || rng.gen::<f64>());

	let mut sol = Vec::new();
	sol.try_reserve_exact(dim)?;
	sol.resize(dim, 0.0f64);

	let rhs = faer::MatRef::from_column_major_slice(&rhs, dim, 1);
	let mut sol = faer::MatMut::from_column_major_slice_mut(&mut sol, dim, 1);

	// optional: fill reducing permutation
	let (perm, perm_inv) = {
		let mut perm = Vec::new();
		let mut perm_inv = Vec::new();
		perm.try_reserve_exact(dim)?;
		perm_inv.try_reserve_exact(dim)?;
		perm.resize(dim, 0usize);
		perm_inv.resize(dim, 0usize);

		let mut mem = MemBuffer::try_new(amd::order_scratch::<usize>(dim, A_nnz))?;
		amd::order(
			&mut perm,
			&mut perm_inv,
			A_lower.symbolic(),
			amd::Control::default(),
			MemStack::new(&mut mem),
		)?;

		(perm, perm_inv)
	};

	let perm = unsafe { faer::perm::PermRef::new_unchecked(&perm, &perm_inv, dim) };

	let A_perm_lower = {
		let mut A_perm_col_ptrs = Vec::new();
		let mut A_perm_row_indices = Vec::new();
		let mut A_perm_values = Vec::new();

		A_perm_col_ptrs.try_reserve_exact(dim + 1)?;
		A_perm_col_ptrs.resize(dim + 1, 0usize);
		A_perm_row_indices.try_reserve_exact(A_nnz)?;
		A_perm_row_indices.resize(A_nnz, 0usize);
		A_perm_values.try_reserve_exact(A_nnz)?;
		A_perm_values.resize(A_nnz, 0.0f64);

		let mut mem = MemBuffer::try_new(faer::sparse::utils::permute_self_adjoint_scratch::<
			usize,
		>(dim))?;
		faer::sparse::utils::permute_self_adjoint_to_unsorted(
			&mut A_perm_values,
			&mut A_perm_col_ptrs,
			&mut A_perm_row_indices,
			A_lower.rb(),
			perm,
			faer::Side::Lower,
			faer::Side::Lower,
			MemStack::new(&mut mem),
		);

		SparseColMat::<usize, f64>::new(
			unsafe {
				SymbolicSparseColMat::new_unchecked(
					dim,
					dim,
					A_perm_col_ptrs,
					None,
					A_perm_row_indices,
				)
			},
			A_perm_values,
		)
	};

	let A_perm_upper = A_perm_lower.rb().transpose().symbolic().to_col_major()?;

	// symbolic analysis
	let symbolic = {
		let mut mem = MemBuffer::try_new(StackReq::any_of(&[
			simplicial::prefactorize_symbolic_cholesky_scratch::<usize>(dim, A_nnz),
			supernodal::factorize_supernodal_symbolic_cholesky_scratch::<usize>(dim),
		]))?;
		let stack = MemStack::new(&mut mem);

		let mut etree = Vec::new();
		let mut col_counts = Vec::new();
		etree.try_reserve_exact(dim)?;
		etree.resize(dim, 0isize);
		col_counts.try_reserve_exact(dim)?;
		col_counts.resize(dim, 0usize);

		simplicial::prefactorize_symbolic_cholesky(
			&mut etree,
			&mut col_counts,
			A_perm_upper.rb(),
			stack,
		);
		supernodal::factorize_supernodal_symbolic_cholesky(
			A_perm_upper.rb(),
			// SAFETY: `etree` was filled correctly by
			// `simplicial::prefactorize_symbolic_cholesky`.
			unsafe { simplicial::EliminationTreeRef::from_inner(&etree) },
			&col_counts,
			stack,
			faer::sparse::linalg::SymbolicSupernodalParams {
				relax: Some(&[(usize::MAX, 1.0)]),
			},
		)?
	};

	// numerical factorization
	let mut mem = MemBuffer::try_new(StackReq::any_of(&[
		supernodal::factorize_supernodal_numeric_llt_scratch::<usize, f64>(
			&symbolic,
			faer::Par::Seq,
			Default::default(),
		),
		supernodal::factorize_supernodal_numeric_ldlt_scratch::<usize, f64>(
			&symbolic,
			faer::Par::Seq,
			Default::default(),
		),
		supernodal::factorize_supernodal_numeric_intranode_lblt_scratch::<usize, f64>(
			&symbolic,
			faer::Par::Seq,
			Default::default(),
		),
		faer::perm::permute_rows_in_place_scratch::<usize, f64>(dim, 1),
		symbolic.solve_in_place_scratch::<f64>(dim, Par::Seq),
	]))?;
	let stack = MemStack::new(&mut mem);

	// llt skipped since a is not positive-definite

	// numerical ldlt factorization
	{
		let mut L_values = Vec::new();
		L_values.try_reserve_exact(symbolic.len_val())?;
		L_values.resize(symbolic.len_val(), 0.0f64);

		supernodal::factorize_supernodal_numeric_ldlt::<usize, f64>(
			&mut L_values,
			A_perm_lower.rb(),
			LdltRegularization::default(),
			&symbolic,
			faer::Par::Seq,
			stack,
			Default::default(),
		);

		let ldlt = supernodal::SupernodalLdltRef::<'_, usize, f64>::new(&symbolic, &L_values);

		sol.copy_from(rhs);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm, stack);
		ldlt.solve_in_place_with_conj(faer::Conj::No, sol.rb_mut(), faer::Par::Seq, stack);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm.inverse(), stack);

		assert!((&A * &sol - &rhs).norm_max() <= 1e-14);
	}

	// numerical intranodal LBLT factorization
	{
		let mut L_values = Vec::new();
		let mut subdiag = Vec::new();
		let mut pivot_perm = Vec::new();
		let mut pivot_perm_inv = Vec::new();

		L_values.try_reserve_exact(symbolic.len_val())?;
		L_values.resize(symbolic.len_val(), 0.0f64);
		subdiag.try_reserve_exact(dim)?;
		subdiag.resize(dim, 0.0f64);
		pivot_perm.try_reserve(dim)?;
		pivot_perm.resize(dim, 0usize);
		pivot_perm_inv.try_reserve(dim)?;
		pivot_perm_inv.resize(dim, 0usize);

		supernodal::factorize_supernodal_numeric_intranode_lblt::<usize, f64>(
			&mut L_values,
			&mut subdiag,
			&mut pivot_perm,
			&mut pivot_perm_inv,
			A_perm_lower.rb(),
			&symbolic,
			faer::Par::Seq,
			stack,
			Default::default(),
		);

		let piv_perm =
			unsafe { faer::perm::PermRef::new_unchecked(&pivot_perm, &pivot_perm_inv, dim) };
		let lblt = supernodal::SupernodalIntranodeLbltRef::<'_, usize, f64>::new(
			&symbolic, &L_values, &subdiag, piv_perm,
		);

		sol.copy_from(rhs);
		// we can merge these two permutations if we want to be optimal
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm, stack);
		faer::perm::permute_rows_in_place(sol.rb_mut(), piv_perm, stack);

		lblt.solve_in_place_no_numeric_permute_with_conj(
			faer::Conj::No,
			sol.rb_mut(),
			faer::Par::Seq,
			stack,
		);

		// we can also merge these two permutations if we want to be optimal
		faer::perm::permute_rows_in_place(sol.rb_mut(), piv_perm.inverse(), stack);
		faer::perm::permute_rows_in_place(sol.rb_mut(), perm.inverse(), stack);

		assert!((&A * &sol - &rhs).norm_max() <= 1e-14);
	}
	Ok(())
}

supernodal_cholesky().unwrap()

Modules

simplicial

simplicial factorization module

supernodal

supernodal factorization module

Structs

CholeskySymbolicParams

tuning parameters for the symbolic cholesky factorization

IntranodeLbltRef

sparse intranodal $LBL^\top$ factorization wrapper

LdltRef

sparse $LDL^H$ factorization wrapper

LltRef

sparse $LL^H$ factorization wrapper

SymbolicCholesky

the symbolic structure of a sparse cholesky decomposition

Enums

SymbolicCholeskyRaw

the inner factorization used for the symbolic cholesky, either simplicial or symbolic

SymmetricOrdering

fill reducing ordering to use for the cholesky factorization

Functions

factorize_symbolic_cholesky

computes the symbolic cholesky factorization of the matrix $A$, or returns an error if the operation could not be completed