faer/linalg/

solvers.rs

use crate::internal_prelude::*;
use crate::{assert, get_global_parallelism};
use alloc::vec;
use alloc::vec::Vec;
use dyn_stack::MemBuffer;
use faer_traits::{ComplexConj, math_utils};
use linalg::svd::ComputeSvdVectors;

pub use linalg::cholesky::ldlt::factor::LdltError;
pub use linalg::cholesky::llt::factor::LltError;
pub use linalg::evd::EvdError;
pub use linalg::gevd::{GevdError, SelfAdjointGevdError};
pub use linalg::svd::SvdError;

/// shape info of a linear system solver
pub trait ShapeCore {
	/// returns the number of rows of the matrix
	fn nrows(&self) -> usize;
	/// returns the number of columns of the matrix
	fn ncols(&self) -> usize;
}

/// linear system solver implementation
pub trait SolveCore<T: ComplexField>: ShapeCore {
	/// solves the equation `self × x = rhs`, implicitly conjugating `self` if needed, and stores
	/// the result in `rhs`
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>);
	/// solves the equation `self.transpose() × x = rhs`, implicitly conjugating `self` if needed,
	/// and stores the result in `rhs`
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>);
}
/// least squares linear system solver implementation
pub trait SolveLstsqCore<T: ComplexField>: ShapeCore {
	/// solves the equation `self × x = rhs` in the sense of least squares, implicitly conjugating
	/// `self` if needed, and stores the result in the top rows of `rhs`
	fn solve_lstsq_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>);
}
/// dense linear system solver
pub trait DenseSolveCore<T: ComplexField>: SolveCore<T> {
	/// returns an approximation of the matrix that was used to create the decomposition
	fn reconstruct(&self) -> Mat<T>;
	/// returns an approximation of the inverse of the matrix that was used to create the
	/// decomposition
	fn inverse(&self) -> Mat<T>;
}

impl<S: ?Sized + ShapeCore> ShapeCore for &S {
	#[inline]
	fn nrows(&self) -> usize {
		(**self).nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		(**self).ncols()
	}
}

impl<T: ComplexField, S: ?Sized + SolveCore<T>> SolveCore<T> for &S {
	#[inline]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		(**self).solve_in_place_with_conj(conj, rhs)
	}

	#[inline]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		(**self).solve_transpose_in_place_with_conj(conj, rhs)
	}
}

impl<T: ComplexField, S: ?Sized + SolveLstsqCore<T>> SolveLstsqCore<T> for &S {
	#[inline]
	fn solve_lstsq_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		(**self).solve_lstsq_in_place_with_conj(conj, rhs)
	}
}

impl<T: ComplexField, S: ?Sized + DenseSolveCore<T>> DenseSolveCore<T> for &S {
	#[inline]
	fn reconstruct(&self) -> Mat<T> {
		(**self).reconstruct()
	}

	#[inline]
	fn inverse(&self) -> Mat<T> {
		(**self).inverse()
	}
}

/// [`SolveCore`] extension trait
pub trait Solve<T: ComplexField>: SolveCore<T> {
	#[track_caller]
	#[inline]
	/// solves $A x = b$
	fn solve_in_place(&self, rhs: impl AsMatMut<T = T, Rows = usize>) {
		self.solve_in_place_with_conj(Conj::No, { rhs }.as_mat_mut().as_dyn_cols_mut());
	}
	#[track_caller]
	#[inline]
	/// solves $\bar A x = b$
	fn solve_conjugate_in_place(&self, rhs: impl AsMatMut<T = T, Rows = usize>) {
		self.solve_in_place_with_conj(Conj::Yes, { rhs }.as_mat_mut().as_dyn_cols_mut());
	}

	#[track_caller]
	#[inline]
	/// solves $A^\top x = b$
	fn solve_transpose_in_place(&self, rhs: impl AsMatMut<T = T, Rows = usize>) {
		self.solve_transpose_in_place_with_conj(Conj::No, { rhs }.as_mat_mut().as_dyn_cols_mut());
	}
	#[track_caller]
	#[inline]
	/// solves $A^H x = b$
	fn solve_adjoint_in_place(&self, rhs: impl AsMatMut<T = T, Rows = usize>) {
		self.solve_transpose_in_place_with_conj(Conj::Yes, { rhs }.as_mat_mut().as_dyn_cols_mut());
	}

	#[track_caller]
	#[inline]
	/// solves $x A = b$
	fn rsolve_in_place(&self, lhs: impl AsMatMut<T = T, Cols = usize>) {
		self.solve_transpose_in_place_with_conj(Conj::No, { lhs }.as_mat_mut().as_dyn_rows_mut().transpose_mut());
	}
	#[track_caller]
	#[inline]
	/// solves $x \bar A = b$
	fn rsolve_conjugate_in_place(&self, lhs: impl AsMatMut<T = T, Cols = usize>) {
		self.solve_transpose_in_place_with_conj(Conj::Yes, { lhs }.as_mat_mut().as_dyn_rows_mut().transpose_mut());
	}

	#[track_caller]
	#[inline]
	/// solves $x A^\top = b$
	fn rsolve_transpose_in_place(&self, lhs: impl AsMatMut<T = T, Cols = usize>) {
		self.solve_in_place_with_conj(Conj::No, { lhs }.as_mat_mut().as_dyn_rows_mut().transpose_mut());
	}
	#[track_caller]
	#[inline]
	/// solves $x A^H = b$
	fn rsolve_adjoint_in_place(&self, lhs: impl AsMatMut<T = T, Cols = usize>) {
		self.solve_in_place_with_conj(Conj::Yes, { lhs }.as_mat_mut().as_dyn_rows_mut().transpose_mut());
	}

	#[track_caller]
	#[inline]
	/// solves $A x = b$
	fn solve<Rhs: AsMatRef<T = T, Rows = usize>>(&self, rhs: Rhs) -> Rhs::Owned {
		let rhs = rhs.as_mat_ref();
		let mut out = Rhs::Owned::zeros(rhs.nrows(), rhs.ncols());
		out.as_mat_mut().copy_from(rhs);
		self.solve_in_place(&mut out);
		out
	}
	#[track_caller]
	#[inline]
	/// solves $\bar A x = b$
	fn solve_conjugate<Rhs: AsMatRef<T = T, Rows = usize>>(&self, rhs: Rhs) -> Rhs::Owned {
		let rhs = rhs.as_mat_ref();
		let mut out = Rhs::Owned::zeros(rhs.nrows(), rhs.ncols());
		out.as_mat_mut().copy_from(rhs);
		self.solve_conjugate_in_place(&mut out);
		out
	}

	#[track_caller]
	#[inline]
	/// solves $A^\top x = b$
	fn solve_transpose<Rhs: AsMatRef<T = T, Rows = usize>>(&self, rhs: Rhs) -> Rhs::Owned {
		let rhs = rhs.as_mat_ref();
		let mut out = Rhs::Owned::zeros(rhs.nrows(), rhs.ncols());
		out.as_mat_mut().copy_from(rhs);
		self.solve_transpose_in_place(&mut out);
		out
	}
	#[track_caller]
	#[inline]
	/// solves $A^H x = b$
	fn solve_adjoint<Rhs: AsMatRef<T = T, Rows = usize>>(&self, rhs: Rhs) -> Rhs::Owned {
		let rhs = rhs.as_mat_ref();
		let mut out = Rhs::Owned::zeros(rhs.nrows(), rhs.ncols());
		out.as_mat_mut().copy_from(rhs);
		self.solve_adjoint_in_place(&mut out);
		out
	}

	#[track_caller]
	#[inline]
	/// solves $x A = b$
	fn rsolve<Lhs: AsMatRef<T = T, Cols = usize>>(&self, lhs: Lhs) -> Lhs::Owned {
		let lhs = lhs.as_mat_ref();
		let mut out = Lhs::Owned::zeros(lhs.nrows(), lhs.ncols());
		out.as_mat_mut().copy_from(lhs);
		self.rsolve_in_place(&mut out);
		out
	}
	#[track_caller]
	#[inline]
	/// solves $x \bar A = b$
	fn rsolve_conjugate<Lhs: AsMatRef<T = T, Cols = usize>>(&self, lhs: Lhs) -> Lhs::Owned {
		let lhs = lhs.as_mat_ref();
		let mut out = Lhs::Owned::zeros(lhs.nrows(), lhs.ncols());
		out.as_mat_mut().copy_from(lhs);
		self.rsolve_conjugate_in_place(&mut out);
		out
	}

	#[track_caller]
	#[inline]
	/// solves $x A^\top = b$
	fn rsolve_transpose<Lhs: AsMatRef<T = T, Cols = usize>>(&self, lhs: Lhs) -> Lhs::Owned {
		let lhs = lhs.as_mat_ref();
		let mut out = Lhs::Owned::zeros(lhs.nrows(), lhs.ncols());
		out.as_mat_mut().copy_from(lhs);
		self.rsolve_transpose_in_place(&mut out);
		out
	}
	#[track_caller]
	#[inline]
	/// solves $x A^H = b$
	fn rsolve_adjoint<Lhs: AsMatRef<T = T, Cols = usize>>(&self, lhs: Lhs) -> Lhs::Owned {
		let lhs = lhs.as_mat_ref();
		let mut out = Lhs::Owned::zeros(lhs.nrows(), lhs.ncols());
		out.as_mat_mut().copy_from(lhs);
		self.rsolve_adjoint_in_place(&mut out);
		out
	}
}

impl<C: Conjugate, Inner: for<'short> Reborrow<'short, Target = mat::Ref<'short, C>>> mat::generic::Mat<Inner> {
	#[track_caller]
	/// returns the $LU$ decomposition of `self` with partial (row) pivoting
	pub fn partial_piv_lu(&self) -> PartialPivLu<C::Canonical> {
		PartialPivLu::new(self.rb())
	}

	#[track_caller]
	/// returns the $LU$ decomposition of `self` with full pivoting
	pub fn full_piv_lu(&self) -> FullPivLu<C::Canonical> {
		FullPivLu::new(self.rb())
	}

	#[track_caller]
	/// returns the $QR$ decomposition of `self`
	pub fn qr(&self) -> Qr<C::Canonical> {
		Qr::new(self.rb())
	}

	#[track_caller]
	/// returns the $QR$ decomposition of `self` with column pivoting
	pub fn col_piv_qr(&self) -> ColPivQr<C::Canonical> {
		ColPivQr::new(self.rb())
	}

	#[track_caller]
	/// returns the svd of `self`
	///
	/// singular values are nonnegative and sorted in nonincreasing order
	pub fn svd(&self) -> Result<Svd<C::Canonical>, SvdError> {
		Svd::new(self.rb())
	}

	#[track_caller]
	/// returns the thin svd of `self`
	///
	/// singular values are nonnegative and sorted in nonincreasing order
	pub fn thin_svd(&self) -> Result<Svd<C::Canonical>, SvdError> {
		Svd::new_thin(self.rb())
	}

	#[track_caller]
	/// returns the $L L^\top$ decomposition of `self`
	pub fn llt(&self, side: Side) -> Result<Llt<C::Canonical>, LltError> {
		Llt::new(self.rb(), side)
	}

	#[track_caller]
	/// returns the $L D L^\top$ decomposition of `self`
	pub fn ldlt(&self, side: Side) -> Result<Ldlt<C::Canonical>, LdltError> {
		Ldlt::new(self.rb(), side)
	}

	#[track_caller]
	/// returns the $LBL^\top$ decomposition of `self`
	pub fn lblt(&self, side: Side) -> Lblt<C::Canonical> {
		Lblt::new(self.rb(), side)
	}

	#[track_caller]
	/// returns the eigendecomposition of `self`, assuming it is self-adjoint
	///
	/// eigenvalues sorted in nondecreasing order
	pub fn self_adjoint_eigen(&self, side: Side) -> Result<SelfAdjointEigen<C::Canonical>, EvdError> {
		SelfAdjointEigen::new(self.rb(), side)
	}

	#[track_caller]
	/// returns the eigenvalues of `self`, assuming it is self-adjoint
	///
	/// eigenvalues sorted in nondecreasing order
	pub fn self_adjoint_eigenvalues(&self, side: Side) -> Result<Vec<Real<C>>, EvdError> {
		#[track_caller]
		pub fn imp<T: ComplexField>(mut A: MatRef<'_, T>, side: Side) -> Result<Vec<T::Real>, EvdError> {
			assert!(A.nrows() == A.ncols());
			if side == Side::Upper {
				A = A.transpose();
			}
			let par = get_global_parallelism();
			let n = A.nrows();

			let mut s = Diag::<T>::zeros(n);

			linalg::evd::self_adjoint_evd(
				A,
				s.as_mut(),
				None,
				par,
				MemStack::new(&mut MemBuffer::new(linalg::evd::self_adjoint_evd_scratch::<T>(
					n,
					linalg::evd::ComputeEigenvectors::No,
					par,
					default(),
				))),
				default(),
			)?;

			Ok(s.column_vector().iter().map(|x| real(x)).collect())
		}

		imp(self.rb().canonical(), side)
	}

	#[track_caller]
	/// returns the singular values of `self`
	///
	/// singular values are nonnegative and sorted in nonincreasing order
	pub fn singular_values(&self) -> Result<Vec<Real<C>>, SvdError> {
		pub fn imp<T: ComplexField>(A: MatRef<'_, T>) -> Result<Vec<T::Real>, SvdError> {
			let par = get_global_parallelism();
			let m = A.nrows();
			let n = A.ncols();

			let mut s = Diag::<T>::zeros(Ord::min(m, n));

			linalg::svd::svd(
				A,
				s.as_mut(),
				None,
				None,
				par,
				MemStack::new(&mut MemBuffer::new(linalg::svd::svd_scratch::<T>(
					m,
					n,
					linalg::svd::ComputeSvdVectors::No,
					linalg::svd::ComputeSvdVectors::No,
					par,
					default(),
				))),
				default(),
			)?;

			Ok(s.column_vector().iter().map(|x| real(x)).collect())
		}

		imp(self.rb().canonical())
	}
}

impl<C: Conjugate> MatRef<'_, C> {
	#[track_caller]
	fn eigen_imp(&self) -> Result<Eigen<Real<C>>, EvdError> {
		if const { C::Canonical::IS_REAL } {
			Eigen::new_from_real(unsafe { crate::hacks::coerce(*self) })
		} else if const { C::IS_CANONICAL } {
			Eigen::new(unsafe { crate::hacks::coerce::<_, MatRef<'_, Complex<Real<C>>>>(*self) })
		} else {
			Eigen::new(unsafe { crate::hacks::coerce::<_, MatRef<'_, ComplexConj<Real<C>>>>(*self) })
		}
	}

	#[track_caller]
	fn gen_eigen_imp(&self, B: MatRef<'_, C>) -> Result<GeneralizedEigen<Real<C>>, GevdError> {
		if const { C::Canonical::IS_REAL } {
			GeneralizedEigen::new_from_real(unsafe { crate::hacks::coerce(*self) }, unsafe { crate::hacks::coerce(B) })
		} else if const { C::IS_CANONICAL } {
			GeneralizedEigen::new(unsafe { crate::hacks::coerce::<_, MatRef<'_, Complex<Real<C>>>>(*self) }, unsafe {
				crate::hacks::coerce::<_, MatRef<'_, Complex<Real<C>>>>(B)
			})
		} else {
			GeneralizedEigen::new(unsafe { crate::hacks::coerce::<_, MatRef<'_, ComplexConj<Real<C>>>>(*self) }, unsafe {
				crate::hacks::coerce::<_, MatRef<'_, ComplexConj<Real<C>>>>(B)
			})
		}
	}

	#[track_caller]
	fn eigenvalues_imp(&self) -> Result<Vec<Complex<Real<C>>>, EvdError> {
		let par = get_global_parallelism();

		if const { C::Canonical::IS_REAL } {
			let A = unsafe { crate::hacks::coerce::<_, MatRef<'_, Real<C>>>(*self) };
			assert!(A.nrows() == A.ncols());
			let n = A.nrows();

			let mut s_re = Diag::<Real<C>>::zeros(n);
			let mut s_im = Diag::<Real<C>>::zeros(n);

			linalg::evd::evd_real(
				A,
				s_re.as_mut(),
				s_im.as_mut(),
				None,
				None,
				par,
				MemStack::new(&mut MemBuffer::new(linalg::evd::evd_scratch::<Real<C>>(
					n,
					linalg::evd::ComputeEigenvectors::No,
					linalg::evd::ComputeEigenvectors::No,
					par,
					default(),
				))),
				default(),
			)?;

			Ok(s_re
				.column_vector()
				.iter()
				.zip(s_im.column_vector().iter())
				.map(|(re, im)| Complex::new(re.clone(), im.clone()))
				.collect())
		} else {
			let A = unsafe { crate::hacks::coerce::<_, MatRef<'_, Complex<Real<C>>>>(self.canonical()) };
			assert!(A.nrows() == A.ncols());
			let n = A.nrows();

			let mut s = Diag::<Complex<Real<C>>>::zeros(n);

			linalg::evd::evd_cplx(
				A,
				s.as_mut(),
				None,
				None,
				par,
				MemStack::new(&mut MemBuffer::new(linalg::evd::evd_scratch::<Complex<Real<C>>>(
					n,
					linalg::evd::ComputeEigenvectors::No,
					linalg::evd::ComputeEigenvectors::No,
					par,
					default(),
				))),
				default(),
			)?;

			if const { C::IS_CANONICAL } {
				Ok(s.column_vector().iter().cloned().collect())
			} else {
				Ok(s.column_vector().iter().map(conj).collect())
			}
		}
	}
}

impl<T: Conjugate, Inner: for<'short> Reborrow<'short, Target = mat::Ref<'short, T>>> mat::generic::Mat<Inner> {
	/// returns the generalized_eigendecomposition of `(self, B)`
	#[track_caller]
	pub fn generalized_eigen(&self, B: impl AsMatRef<T = T, Rows = usize, Cols = usize>) -> Result<GeneralizedEigen<Real<T>>, GevdError> {
		self.rb().gen_eigen_imp(B.as_mat_ref())
	}

	/// returns the eigendecomposition of `self`
	#[track_caller]
	pub fn eigen(&self) -> Result<Eigen<Real<T>>, EvdError> {
		self.rb().eigen_imp()
	}

	/// returns the eigenvalues of `self`
	#[track_caller]
	pub fn eigenvalues(&self) -> Result<Vec<Complex<Real<T>>>, EvdError> {
		self.rb().eigenvalues_imp()
	}
}

/// [`SolveLstsqCore`] extension trait
pub trait SolveLstsq<T: ComplexField>: SolveLstsqCore<T> {
	#[track_caller]
	#[inline]
	/// solves $A x = b$ in the sense of least squares.
	fn solve_lstsq_in_place(&self, rhs: impl AsMatMut<T = T, Rows = usize>) {
		self.solve_lstsq_in_place_with_conj(Conj::No, { rhs }.as_mat_mut().as_dyn_cols_mut());
	}

	#[track_caller]
	#[inline]
	/// solves $\bar A x = b$ in the sense of least squares.
	fn solve_conjugate_lstsq_in_place(&self, rhs: impl AsMatMut<T = T, Rows = usize>) {
		self.solve_lstsq_in_place_with_conj(Conj::Yes, { rhs }.as_mat_mut().as_dyn_cols_mut());
	}

	#[track_caller]
	#[inline]
	/// solves $A x = b$ in the sense of least squares.
	fn solve_lstsq<Rhs: AsMatRef<T = T, Rows = usize>>(&self, rhs: Rhs) -> Rhs::Owned {
		let rhs = rhs.as_mat_ref();
		let mut out = Rhs::Owned::zeros(rhs.nrows(), rhs.ncols());
		out.as_mat_mut().copy_from(rhs);
		self.solve_lstsq_in_place(&mut out);
		out.truncate(self.ncols(), rhs.ncols());
		out
	}
	#[track_caller]
	#[inline]
	/// solves $\bar A x = b$ in the sense of least squares.
	fn solve_conjugate_lstsq<Rhs: AsMatRef<T = T, Rows = usize>>(&self, rhs: Rhs) -> Rhs::Owned {
		let rhs = rhs.as_mat_ref();
		let mut out = Rhs::Owned::zeros(rhs.nrows(), rhs.ncols());
		out.as_mat_mut().copy_from(rhs);
		self.solve_conjugate_lstsq_in_place(&mut out);
		out.truncate(self.ncols(), rhs.ncols());
		out
	}
}
/// [`DenseSolveCore`] extension trait
pub trait DenseSolve<T: ComplexField>: DenseSolveCore<T> {}

impl<T: ComplexField, S: ?Sized + SolveCore<T>> Solve<T> for S {}
impl<T: ComplexField, S: ?Sized + SolveLstsqCore<T>> SolveLstsq<T> for S {}
impl<T: ComplexField, S: ?Sized + DenseSolveCore<T>> DenseSolve<T> for S {}

/// $L L^\top$ decomposition
#[derive(Clone, Debug)]
pub struct Llt<T> {
	L: Mat<T>,
}

/// $L D L^\top$ decomposition
#[derive(Clone, Debug)]
pub struct Ldlt<T> {
	L: Mat<T>,
	D: Diag<T>,
}

/// $LBL^\top$ decomposition
#[derive(Clone, Debug)]
pub struct Lblt<T> {
	L: Mat<T>,
	B_diag: Diag<T>,
	B_subdiag: Diag<T>,
	P: Perm<usize>,
}

/// $LU$ decomposition with partial (row) pivoting
#[derive(Clone, Debug)]
pub struct PartialPivLu<T> {
	L: Mat<T>,
	U: Mat<T>,
	P: Perm<usize>,
}

/// $LU$ decomposition with full pivoting
#[derive(Clone, Debug)]
pub struct FullPivLu<T> {
	L: Mat<T>,
	U: Mat<T>,
	P: Perm<usize>,
	Q: Perm<usize>,
}

/// $QR$ decomposition
#[derive(Clone, Debug)]
pub struct Qr<T> {
	Q_basis: Mat<T>,
	Q_coeff: Mat<T>,
	R: Mat<T>,
}

/// $QR$ decomposition with column pivoting
#[derive(Clone, Debug)]
pub struct ColPivQr<T> {
	Q_basis: Mat<T>,
	Q_coeff: Mat<T>,
	R: Mat<T>,
	P: Perm<usize>,
}

/// svd decomposition (either full or thin)
#[derive(Clone, Debug)]
pub struct Svd<T> {
	U: Mat<T>,
	V: Mat<T>,
	S: Diag<T>,
}

/// self-adjoint eigendecomposition
#[derive(Clone, Debug)]
pub struct SelfAdjointEigen<T> {
	U: Mat<T>,
	S: Diag<T>,
}

/// eigendecomposition
#[derive(Clone, Debug)]
pub struct Eigen<T> {
	U: Mat<Complex<T>>,
	S: Diag<Complex<T>>,
}

/// eigendecomposition
#[derive(Clone, Debug)]
pub struct GeneralizedEigen<T> {
	U: Mat<Complex<T>>,
	S_a: Diag<Complex<T>>,
	S_b: Diag<Complex<T>>,
}

impl<T: ComplexField> Llt<T> {
	/// returns the $L L^\top$ decomposition of $A$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>, side: Side) -> Result<Self, LltError> {
		assert!(all(A.nrows() == A.ncols()));
		let n = A.nrows();

		let mut L = Mat::zeros(n, n);
		match side {
			Side::Lower => L.copy_from_triangular_lower(A),
			Side::Upper => L.copy_from_triangular_lower(A.adjoint()),
		}

		Self::new_imp(L)
	}

	#[track_caller]
	fn new_imp(mut L: Mat<T>) -> Result<Self, LltError> {
		let par = get_global_parallelism();

		let n = L.nrows();

		let mut mem = MemBuffer::new(linalg::cholesky::llt::factor::cholesky_in_place_scratch::<T>(n, par, default()));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::llt::factor::cholesky_in_place(L.as_mut(), Default::default(), par, stack, default())?;
		z!(&mut L).for_each_triangular_upper(linalg::zip::Diag::Skip, |uz!(x)| *x = zero());

		Ok(Self { L })
	}

	/// returns the $L$ factor
	pub fn L(&self) -> MatRef<'_, T> {
		self.L.as_ref()
	}
}

impl<T: ComplexField> Ldlt<T> {
	/// returns the $L D L^\top$ decomposition of $A$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>, side: Side) -> Result<Self, LdltError> {
		assert!(all(A.nrows() == A.ncols()));
		let n = A.nrows();

		let mut L = Mat::zeros(n, n);
		match side {
			Side::Lower => L.copy_from_triangular_lower(A),
			Side::Upper => L.copy_from_triangular_lower(A.adjoint()),
		}

		Self::new_imp(L)
	}

	#[track_caller]
	fn new_imp(mut L: Mat<T>) -> Result<Self, LdltError> {
		let par = get_global_parallelism();

		let n = L.nrows();
		let mut D = Diag::zeros(n);

		let mut mem = MemBuffer::new(linalg::cholesky::ldlt::factor::cholesky_in_place_scratch::<T>(n, par, default()));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::ldlt::factor::cholesky_in_place(L.as_mut(), Default::default(), par, stack, default())?;

		D.copy_from(L.diagonal());
		L.diagonal_mut().fill(one());
		z!(&mut L).for_each_triangular_upper(linalg::zip::Diag::Skip, |uz!(x)| *x = zero());

		Ok(Self { L, D })
	}

	/// returns the $L$ factor
	pub fn L(&self) -> MatRef<'_, T> {
		self.L.as_ref()
	}

	/// returns the $D$ factor
	pub fn D(&self) -> DiagRef<'_, T> {
		self.D.as_ref()
	}
}

impl<T: ComplexField> Lblt<T> {
	/// returns the $LBL^\top$ decomposition of $A$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>, side: Side) -> Self {
		assert!(all(A.nrows() == A.ncols()));
		let n = A.nrows();

		let mut L = Mat::zeros(n, n);
		match side {
			Side::Lower => L.copy_from_triangular_lower(A),
			Side::Upper => L.copy_from_triangular_lower(A.adjoint()),
		}
		Self::new_imp(L)
	}

	#[track_caller]
	fn new_imp(mut L: Mat<T>) -> Self {
		let par = get_global_parallelism();

		let n = L.nrows();

		let mut diag = Diag::zeros(n);
		let mut subdiag = Diag::zeros(n);
		let mut perm_fwd = vec![0usize; n];
		let mut perm_bwd = vec![0usize; n];

		let mut mem = MemBuffer::new(linalg::cholesky::lblt::factor::cholesky_in_place_scratch::<usize, T>(n, par, default()));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::lblt::factor::cholesky_in_place(L.as_mut(), subdiag.as_mut(), &mut perm_fwd, &mut perm_bwd, par, stack, default());

		diag.copy_from(L.diagonal());
		L.diagonal_mut().fill(one());
		z!(&mut L).for_each_triangular_upper(linalg::zip::Diag::Skip, |uz!(x)| *x = zero());

		Self {
			L,
			B_diag: diag,
			B_subdiag: subdiag,
			P: unsafe { Perm::new_unchecked(perm_fwd.into_boxed_slice(), perm_bwd.into_boxed_slice()) },
		}
	}

	/// returns the $L$ factor
	pub fn L(&self) -> MatRef<'_, T> {
		self.L.as_ref()
	}

	/// returns the diagonal of the $B$ factor
	pub fn B_diag(&self) -> DiagRef<'_, T> {
		self.B_diag.as_ref()
	}

	/// returns the subdiagonal of the $B$ factor
	pub fn B_subdiag(&self) -> DiagRef<'_, T> {
		self.B_subdiag.as_ref()
	}

	/// returns the pivoting permutation $P$
	pub fn P(&self) -> PermRef<'_, usize> {
		self.P.as_ref()
	}
}

fn split_LU<T: ComplexField>(LU: Mat<T>) -> (Mat<T>, Mat<T>) {
	let (m, n) = LU.shape();
	let size = Ord::min(m, n);

	let (L, U) = if m >= n {
		let mut L = LU;
		let mut U = Mat::zeros(size, size);

		U.copy_from_triangular_upper(L.get(..size, ..size));

		z!(&mut L).for_each_triangular_upper(linalg::zip::Diag::Skip, |uz!(x)| *x = zero());
		L.diagonal_mut().fill(one());

		(L, U)
	} else {
		let mut U = LU;
		let mut L = Mat::zeros(size, size);

		L.copy_from_strict_triangular_lower(U.get(..size, ..size));

		z!(&mut U).for_each_triangular_lower(linalg::zip::Diag::Skip, |uz!(x)| *x = zero());
		L.diagonal_mut().fill(one());

		(L, U)
	};
	(L, U)
}

impl<T: ComplexField> PartialPivLu<T> {
	/// returns the $LU$ decomposition of $A$ with partial pivoting
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>) -> Self {
		let LU = A.to_owned();
		Self::new_imp(LU)
	}

	#[track_caller]
	fn new_imp(mut LU: Mat<T>) -> Self {
		let par = get_global_parallelism();

		let (m, n) = LU.shape();
		let mut row_perm_fwd = vec![0usize; m];
		let mut row_perm_bwd = vec![0usize; m];

		linalg::lu::partial_pivoting::factor::lu_in_place(
			LU.as_mut(),
			&mut row_perm_fwd,
			&mut row_perm_bwd,
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::lu::partial_pivoting::factor::lu_in_place_scratch::<usize, T>(m, n, par, default()),
			)),
			default(),
		);

		let (L, U) = split_LU(LU);

		Self {
			L,
			U,
			P: unsafe { Perm::new_unchecked(row_perm_fwd.into_boxed_slice(), row_perm_bwd.into_boxed_slice()) },
		}
	}

	/// returns the $L$ factor
	pub fn L(&self) -> MatRef<'_, T> {
		self.L.as_ref()
	}

	/// returns the $U$ factor
	pub fn U(&self) -> MatRef<'_, T> {
		self.U.as_ref()
	}

	/// returns the row pivoting permutation $P$
	pub fn P(&self) -> PermRef<'_, usize> {
		self.P.as_ref()
	}
}

impl<T: ComplexField> FullPivLu<T> {
	/// returns the $LU$ decomposition of $A$ with full pivoting
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>) -> Self {
		let LU = A.to_owned();
		Self::new_imp(LU)
	}

	#[track_caller]
	fn new_imp(mut LU: Mat<T>) -> Self {
		let par = get_global_parallelism();

		let (m, n) = LU.shape();
		let mut row_perm_fwd = vec![0usize; m];
		let mut row_perm_bwd = vec![0usize; m];
		let mut col_perm_fwd = vec![0usize; n];
		let mut col_perm_bwd = vec![0usize; n];

		linalg::lu::full_pivoting::factor::lu_in_place(
			LU.as_mut(),
			&mut row_perm_fwd,
			&mut row_perm_bwd,
			&mut col_perm_fwd,
			&mut col_perm_bwd,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::lu::full_pivoting::factor::lu_in_place_scratch::<usize, T>(
				m,
				n,
				par,
				default(),
			))),
			default(),
		);

		let (L, U) = split_LU(LU);

		Self {
			L,
			U,
			P: unsafe { Perm::new_unchecked(row_perm_fwd.into_boxed_slice(), row_perm_bwd.into_boxed_slice()) },
			Q: unsafe { Perm::new_unchecked(col_perm_fwd.into_boxed_slice(), col_perm_bwd.into_boxed_slice()) },
		}
	}

	/// returns the factor $L$
	pub fn L(&self) -> MatRef<'_, T> {
		self.L.as_ref()
	}

	/// returns the factor $U$
	pub fn U(&self) -> MatRef<'_, T> {
		self.U.as_ref()
	}

	/// returns the row pivoting permutation $P$
	pub fn P(&self) -> PermRef<'_, usize> {
		self.P.as_ref()
	}

	/// returns the column pivoting permutation $P$
	pub fn Q(&self) -> PermRef<'_, usize> {
		self.Q.as_ref()
	}
}

impl<T: ComplexField> Qr<T> {
	/// returns the $QR$ decomposition of $A$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>) -> Self {
		let QR = A.to_owned();
		Self::new_imp(QR)
	}

	#[track_caller]
	fn new_imp(mut QR: Mat<T>) -> Self {
		let par = get_global_parallelism();

		let (m, n) = QR.shape();
		let size = Ord::min(m, n);

		let block_size = linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
		let mut Q_coeff = Mat::zeros(block_size, size);

		linalg::qr::no_pivoting::factor::qr_in_place(
			QR.as_mut(),
			Q_coeff.as_mut(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::no_pivoting::factor::qr_in_place_scratch::<T>(
				m,
				n,
				block_size,
				par,
				default(),
			))),
			default(),
		);

		let (Q_basis, R) = split_LU(QR);

		Self { Q_basis, Q_coeff, R }
	}

	/// returns the householder basis of $Q$
	pub fn Q_basis(&self) -> MatRef<'_, T> {
		self.Q_basis.as_ref()
	}

	/// returns the householder coefficients of $Q$
	pub fn Q_coeff(&self) -> MatRef<'_, T> {
		self.Q_coeff.as_ref()
	}

	/// returns the factor $R$
	pub fn R(&self) -> MatRef<'_, T> {
		self.R.as_ref()
	}

	/// returns the upper trapezoidal part of $R$
	pub fn thin_R(&self) -> MatRef<'_, T> {
		let size = Ord::min(self.nrows(), self.ncols());
		self.R.get(..size, ..)
	}

	/// computes the factor $Q$
	pub fn compute_Q(&self) -> Mat<T> {
		let mut Q = Mat::identity(self.nrows(), self.nrows());
		let par = get_global_parallelism();
		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			Conj::No,
			Q.rb_mut(),
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
					self.nrows(),
					self.Q_coeff.nrows(),
					self.nrows(),
				),
			)),
		);
		Q
	}

	/// computes the first $\min(\text{nrows}, \text{ncols})$ columns of the factor $Q$
	pub fn compute_thin_Q(&self) -> Mat<T> {
		let size = Ord::min(self.nrows(), self.ncols());
		let mut Q = Mat::identity(self.nrows(), size);
		let par = get_global_parallelism();
		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			Conj::No,
			Q.rb_mut(),
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(self.nrows(), self.Q_coeff.nrows(), size),
			)),
		);
		Q
	}
}

impl<T: ComplexField> ColPivQr<T> {
	/// returns the $QR$ decomposition of $A$ with column pivoting
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>) -> Self {
		let QR = A.to_owned();
		Self::new_imp(QR)
	}

	#[track_caller]
	fn new_imp(mut QR: Mat<T>) -> Self {
		let par = get_global_parallelism();

		let (m, n) = QR.shape();
		let size = Ord::min(m, n);

		let mut col_perm_fwd = vec![0usize; n];
		let mut col_perm_bwd = vec![0usize; n];

		let block_size = linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
		let mut Q_coeff = Mat::zeros(block_size, size);

		linalg::qr::col_pivoting::factor::qr_in_place(
			QR.as_mut(),
			Q_coeff.as_mut(),
			&mut col_perm_fwd,
			&mut col_perm_bwd,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::col_pivoting::factor::qr_in_place_scratch::<usize, T>(
				m,
				n,
				block_size,
				par,
				default(),
			))),
			default(),
		);

		let (Q_basis, R) = split_LU(QR);

		Self {
			Q_basis,
			Q_coeff,
			R,
			P: unsafe { Perm::new_unchecked(col_perm_fwd.into_boxed_slice(), col_perm_bwd.into_boxed_slice()) },
		}
	}

	/// returns the householder basis of $Q$
	pub fn Q_basis(&self) -> MatRef<'_, T> {
		self.Q_basis.as_ref()
	}

	/// returns the householder coefficients of $Q$
	pub fn Q_coeff(&self) -> MatRef<'_, T> {
		self.Q_coeff.as_ref()
	}

	/// returns the factor $R$
	pub fn R(&self) -> MatRef<'_, T> {
		self.R.as_ref()
	}

	/// returns the upper trapezoidal part of $R$
	pub fn thin_R(&self) -> MatRef<'_, T> {
		let size = Ord::min(self.nrows(), self.ncols());
		self.R.get(..size, ..)
	}

	/// computes the factor $Q$
	pub fn compute_Q(&self) -> Mat<T> {
		let mut Q = Mat::identity(self.nrows(), self.nrows());
		let par = get_global_parallelism();
		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			Conj::No,
			Q.rb_mut(),
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
					self.nrows(),
					self.Q_coeff.nrows(),
					self.nrows(),
				),
			)),
		);
		Q
	}

	/// computes the first $\min(\text{nrows}, \text{ncols})$ columns of the factor $Q$
	pub fn compute_thin_Q(&self) -> Mat<T> {
		let size = Ord::min(self.nrows(), self.ncols());
		let mut Q = Mat::identity(self.nrows(), size);
		let par = get_global_parallelism();
		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			Conj::No,
			Q.rb_mut(),
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(self.nrows(), self.Q_coeff.nrows(), size),
			)),
		);
		Q
	}

	/// returns the column pivoting permutation $P$
	pub fn P(&self) -> PermRef<'_, usize> {
		self.P.as_ref()
	}
}

impl<T: ComplexField> Svd<T> {
	/// returns the svd of $A$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>) -> Result<Self, SvdError> {
		Self::new_imp(A.canonical(), Conj::get::<C>(), false)
	}

	/// returns the thin svd of $A$
	#[track_caller]
	pub fn new_thin<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>) -> Result<Self, SvdError> {
		Self::new_imp(A.canonical(), Conj::get::<C>(), true)
	}

	#[track_caller]
	fn new_imp(A: MatRef<'_, T>, conj: Conj, thin: bool) -> Result<Self, SvdError> {
		let par = get_global_parallelism();

		let (m, n) = A.shape();
		let size = Ord::min(m, n);

		let mut U = Mat::zeros(m, if thin { size } else { m });
		let mut V = Mat::zeros(n, if thin { size } else { n });
		let mut S = Diag::zeros(size);

		let compute = if thin { ComputeSvdVectors::Thin } else { ComputeSvdVectors::Full };

		linalg::svd::svd(
			A,
			S.as_mut(),
			Some(U.as_mut()),
			Some(V.as_mut()),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::svd::svd_scratch::<T>(m, n, compute, compute, par, default()))),
			default(),
		)?;

		if conj == Conj::Yes {
			for c in U.col_iter_mut() {
				for x in c.iter_mut() {
					*x = math_utils::conj(x);
				}
			}
			for c in V.col_iter_mut() {
				for x in c.iter_mut() {
					*x = math_utils::conj(x);
				}
			}
		}

		Ok(Self { U, V, S })
	}

	/// returns the factor $U$
	pub fn U(&self) -> MatRef<'_, T> {
		self.U.as_ref()
	}

	/// returns the factor $V$
	pub fn V(&self) -> MatRef<'_, T> {
		self.V.as_ref()
	}

	/// returns the factor $S$
	pub fn S(&self) -> DiagRef<'_, T> {
		self.S.as_ref()
	}

	/// returns the pseudoinverse of the original matrix $A$.
	pub fn pseudoinverse(&self) -> Mat<T> {
		let U = self.U();
		let V = self.V();
		let S = self.S();
		let par = get_global_parallelism();
		let stack = &mut MemBuffer::new(linalg::svd::pseudoinverse_from_svd_scratch::<T>(self.nrows(), self.ncols(), par));
		let mut pinv = Mat::zeros(self.ncols(), self.nrows());
		linalg::svd::pseudoinverse_from_svd(pinv.rb_mut(), S, U, V, par, MemStack::new(stack));
		pinv
	}
}

impl<T: ComplexField> SelfAdjointEigen<T> {
	/// returns the eigendecomposition of $A$, assuming it is self-adjoint
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = T>>(A: MatRef<'_, C>, side: Side) -> Result<Self, EvdError> {
		assert!(A.nrows() == A.ncols());

		match side {
			Side::Lower => Self::new_imp(A.canonical(), Conj::get::<C>()),
			Side::Upper => Self::new_imp(A.adjoint().canonical(), Conj::get::<C::Conj>()),
		}
	}

	#[track_caller]
	fn new_imp(A: MatRef<'_, T>, conj: Conj) -> Result<Self, EvdError> {
		let par = get_global_parallelism();

		let n = A.nrows();

		let mut U = Mat::zeros(n, n);
		let mut S = Diag::zeros(n);

		linalg::evd::self_adjoint_evd(
			A,
			S.as_mut(),
			Some(U.as_mut()),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::evd::self_adjoint_evd_scratch::<T>(
				n,
				linalg::evd::ComputeEigenvectors::Yes,
				par,
				default(),
			))),
			default(),
		)?;

		if conj == Conj::Yes {
			for c in U.col_iter_mut() {
				for x in c.iter_mut() {
					*x = math_utils::conj(x);
				}
			}
		}

		Ok(Self { U, S })
	}

	/// returns the factor $U$
	pub fn U(&self) -> MatRef<'_, T> {
		self.U.as_ref()
	}

	/// returns the factor $S$
	pub fn S(&self) -> DiagRef<'_, T> {
		self.S.as_ref()
	}

	/// returns the pseudoinverse of the original matrix $A$.
	pub fn pseudoinverse(&self) -> Mat<T> {
		let U = self.U();
		let S = self.S();
		let par = get_global_parallelism();
		let stack = &mut MemBuffer::new(linalg::evd::pseudoinverse_from_self_adjoint_evd_scratch::<T>(self.nrows(), par));
		let mut pinv = Mat::zeros(self.ncols(), self.nrows());
		linalg::evd::pseudoinverse_from_self_adjoint_evd(pinv.rb_mut(), S, U, par, MemStack::new(stack));
		pinv
	}
}

fn real_to_cplx<T: RealField>(
	mut U: MatMut<'_, Complex<T>>,
	mut S: DiagMut<'_, Complex<T>>,
	U_real: MatRef<'_, T>,
	S_re: DiagRef<'_, T>,
	S_im: DiagRef<'_, T>,
) {
	let n = U.ncols();

	let mut j = 0;
	while j < n {
		if S_im[j] == zero() {
			S[j] = Complex::new(S_re[j].clone(), zero());

			for i in 0..n {
				U[(i, j)] = Complex::new(U_real[(i, j)].clone(), zero());
			}

			j += 1;
		} else {
			S[j] = Complex::new(S_re[j].clone(), S_im[j].clone());
			S[j + 1] = Complex::new(S_re[j].clone(), neg(&S_im[j]));

			for i in 0..n {
				U[(i, j)] = Complex::new(U_real[(i, j)].clone(), U_real[(i, j + 1)].clone());
				U[(i, j + 1)] = Complex::new(U_real[(i, j)].clone(), neg(&U_real[(i, j + 1)]));
			}

			j += 2;
		}
	}
}

impl<T: RealField> Eigen<T> {
	/// returns the eigendecomposition of $A$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = Complex<T>>>(A: MatRef<'_, C>) -> Result<Self, EvdError> {
		assert!(A.nrows() == A.ncols());
		Self::new_imp(A.canonical(), Conj::get::<C>())
	}

	/// returns the eigendecomposition of $A$
	#[track_caller]
	pub fn new_from_real(A: MatRef<'_, T>) -> Result<Self, EvdError> {
		assert!(A.nrows() == A.ncols());

		let par = get_global_parallelism();

		let n = A.nrows();

		let mut U_real = Mat::zeros(n, n);
		let mut S_re = Diag::zeros(n);
		let mut S_im = Diag::zeros(n);

		linalg::evd::evd_real(
			A,
			S_re.as_mut(),
			S_im.as_mut(),
			None,
			Some(U_real.as_mut()),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::evd::evd_scratch::<T>(
				n,
				linalg::evd::ComputeEigenvectors::No,
				linalg::evd::ComputeEigenvectors::Yes,
				par,
				default(),
			))),
			default(),
		)?;

		let mut U = Mat::zeros(n, n);
		let mut S = Diag::zeros(n);

		real_to_cplx(U.as_mut(), S.as_mut(), U_real.as_ref(), S_re.as_ref(), S_im.as_ref());

		Ok(Self { U, S })
	}

	fn new_imp(A: MatRef<'_, Complex<T>>, conj: Conj) -> Result<Self, EvdError> {
		let par = get_global_parallelism();

		let n = A.nrows();

		let mut U = Mat::zeros(n, n);
		let mut S = Diag::zeros(n);

		linalg::evd::evd_cplx(
			A,
			S.as_mut(),
			None,
			Some(U.as_mut()),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::evd::evd_scratch::<Complex<T>>(
				n,
				linalg::evd::ComputeEigenvectors::No,
				linalg::evd::ComputeEigenvectors::Yes,
				par,
				default(),
			))),
			default(),
		)?;

		if conj == Conj::Yes {
			zip!(&mut U).for_each(|unzip!(c)| *c = math_utils::conj(c));
			zip!(&mut S).for_each(|unzip!(c)| *c = math_utils::conj(c));
		}

		Ok(Self { U, S })
	}

	/// returns the factor $U$
	pub fn U(&self) -> MatRef<'_, Complex<T>> {
		self.U.as_ref()
	}

	/// returns the factor $S$
	pub fn S(&self) -> DiagRef<'_, Complex<T>> {
		self.S.as_ref()
	}
}

impl<T: RealField> GeneralizedEigen<T> {
	/// returns the generalized eigendecomposition of $(A, B)$
	#[track_caller]
	pub fn new<C: Conjugate<Canonical = Complex<T>>>(A: MatRef<'_, C>, B: MatRef<'_, C>) -> Result<Self, GevdError> {
		let n = A.nrows();
		assert!(all(A.nrows() == n, A.ncols() == n, B.nrows() == n, B.ncols() == n));
		Self::new_imp(A.canonical(), B.canonical(), Conj::get::<C>())
	}

	/// returns the generalized eigendecomposition of $(A, B)$
	#[track_caller]
	pub fn new_from_real(A: MatRef<'_, T>, B: MatRef<'_, T>) -> Result<Self, GevdError> {
		let n = A.nrows();
		assert!(all(A.nrows() == n, A.ncols() == n, B.nrows() == n, B.ncols() == n));

		let par = get_global_parallelism();

		let mut U_real = Mat::zeros(n, n);
		let mut S_re = Diag::zeros(n);
		let mut S_im = Diag::zeros(n);
		let mut S_b = Diag::zeros(n);
		let A = &mut A.cloned();
		let B = &mut B.cloned();

		linalg::gevd::gevd_real(
			A.as_mut(),
			B.as_mut(),
			S_re.as_mut(),
			S_im.as_mut(),
			S_b.as_mut(),
			None,
			Some(U_real.as_mut()),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::gevd::gevd_scratch::<T>(
				n,
				linalg::evd::ComputeEigenvectors::No,
				linalg::evd::ComputeEigenvectors::Yes,
				par,
				default(),
			))),
			default(),
		)?;

		let mut U = Mat::zeros(n, n);
		let mut S_a = Diag::zeros(n);
		let S_b = zip!(&S_b).map(|unzip!(x)| Complex::new(x.clone(), zero()));

		real_to_cplx(U.as_mut(), S_a.as_mut(), U_real.as_ref(), S_re.as_ref(), S_im.as_ref());

		Ok(Self { U, S_a, S_b })
	}

	fn new_imp(A: MatRef<'_, Complex<T>>, B: MatRef<'_, Complex<T>>, conj: Conj) -> Result<Self, GevdError> {
		let par = get_global_parallelism();

		let n = A.nrows();

		let mut U = Mat::zeros(n, n);
		let mut S_a = Diag::zeros(n);
		let mut S_b = Diag::zeros(n);
		let A = &mut A.cloned();
		let B = &mut B.cloned();

		linalg::gevd::gevd_cplx(
			A.as_mut(),
			B.as_mut(),
			S_a.as_mut(),
			S_b.as_mut(),
			None,
			Some(U.as_mut()),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::gevd::gevd_scratch::<Complex<T>>(
				n,
				linalg::evd::ComputeEigenvectors::No,
				linalg::evd::ComputeEigenvectors::Yes,
				par,
				default(),
			))),
			default(),
		)?;

		if conj == Conj::Yes {
			zip!(&mut U).for_each(|unzip!(c)| *c = math_utils::conj(c));
			zip!(&mut S_a).for_each(|unzip!(c)| *c = math_utils::conj(c));
			zip!(&mut S_b).for_each(|unzip!(c)| *c = math_utils::conj(c));
		}

		Ok(Self { U, S_a, S_b })
	}

	/// returns the factor $U$
	pub fn U(&self) -> MatRef<'_, Complex<T>> {
		self.U.as_ref()
	}

	/// returns the factor $S_a$
	pub fn S_a(&self) -> DiagRef<'_, Complex<T>> {
		self.S_a.as_ref()
	}

	/// returns the factor $S_b$
	pub fn S_b(&self) -> DiagRef<'_, Complex<T>> {
		self.S_b.as_ref()
	}
}

impl<T: ComplexField> ShapeCore for Llt<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.L().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.L().ncols()
	}
}
impl<T: ComplexField> ShapeCore for Ldlt<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.L().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.L().ncols()
	}
}
impl<T: ComplexField> ShapeCore for Lblt<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.L().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.L().ncols()
	}
}
impl<T: ComplexField> ShapeCore for PartialPivLu<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.L().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.U().ncols()
	}
}
impl<T: ComplexField> ShapeCore for FullPivLu<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.L().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.U().ncols()
	}
}
impl<T: ComplexField> ShapeCore for Qr<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.Q_basis().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.R().ncols()
	}
}
impl<T: ComplexField> ShapeCore for ColPivQr<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.Q_basis().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.R().ncols()
	}
}
impl<T: ComplexField> ShapeCore for Svd<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.U().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.V().nrows()
	}
}
impl<T: ComplexField> ShapeCore for SelfAdjointEigen<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.U().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.U().nrows()
	}
}
impl<T: RealField> ShapeCore for Eigen<T> {
	#[inline]
	fn nrows(&self) -> usize {
		self.U().nrows()
	}

	#[inline]
	fn ncols(&self) -> usize {
		self.U().nrows()
	}
}

impl<T: ComplexField> SolveCore<T> for Llt<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		let mut mem = MemBuffer::new(linalg::cholesky::llt::solve::solve_in_place_scratch::<T>(
			self.L.nrows(),
			rhs.ncols(),
			par,
		));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::llt::solve::solve_in_place_with_conj(self.L.as_ref(), conj, rhs, par, stack);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		let mut mem = MemBuffer::new(linalg::cholesky::llt::solve::solve_in_place_scratch::<T>(
			self.L.nrows(),
			rhs.ncols(),
			par,
		));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::llt::solve::solve_in_place_with_conj(self.L.as_ref(), conj.compose(Conj::Yes), rhs, par, stack);
	}
}

#[math]
fn make_self_adjoint<T: ComplexField>(mut A: MatMut<'_, T>) {
	assert!(A.nrows() == A.ncols());
	let n = A.nrows();
	for j in 0..n {
		A[(j, j)] = from_real(real(A[(j, j)]));
		for i in 0..j {
			A[(i, j)] = conj(A[(j, i)]);
		}
	}
}

impl<T: ComplexField> DenseSolveCore<T> for Llt<T> {
	#[track_caller]
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();

		let n = self.L.nrows();
		let mut out = Mat::zeros(n, n);

		let mut mem = MemBuffer::new(linalg::cholesky::llt::reconstruct::reconstruct_scratch::<T>(n, par));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::llt::reconstruct::reconstruct(out.as_mut(), self.L(), par, stack);

		make_self_adjoint(out.as_mut());
		out
	}

	#[track_caller]
	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		let n = self.L.nrows();
		let mut out = Mat::zeros(n, n);

		let mut mem = MemBuffer::new(linalg::cholesky::llt::inverse::inverse_scratch::<T>(n, par));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::llt::inverse::inverse(out.as_mut(), self.L(), par, stack);

		make_self_adjoint(out.as_mut());
		out
	}
}

impl<T: ComplexField> SolveCore<T> for Ldlt<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		let mut mem = MemBuffer::new(linalg::cholesky::ldlt::solve::solve_in_place_scratch::<T>(
			self.L.nrows(),
			rhs.ncols(),
			par,
		));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::ldlt::solve::solve_in_place_with_conj(self.L.as_ref(), self.D.as_ref(), conj, rhs, par, stack);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		let mut mem = MemBuffer::new(linalg::cholesky::ldlt::solve::solve_in_place_scratch::<T>(
			self.L.nrows(),
			rhs.ncols(),
			par,
		));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::ldlt::solve::solve_in_place_with_conj(self.L(), self.D(), conj.compose(Conj::Yes), rhs, par, stack);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for Ldlt<T> {
	#[track_caller]
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();

		let n = self.L.nrows();
		let mut out = Mat::zeros(n, n);

		let mut mem = MemBuffer::new(linalg::cholesky::ldlt::reconstruct::reconstruct_scratch::<T>(n, par));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::ldlt::reconstruct::reconstruct(out.as_mut(), self.L(), self.D(), par, stack);

		make_self_adjoint(out.as_mut());
		out
	}

	#[track_caller]
	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		let n = self.L.nrows();
		let mut out = Mat::zeros(n, n);

		let mut mem = MemBuffer::new(linalg::cholesky::ldlt::inverse::inverse_scratch::<T>(n, par));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::ldlt::inverse::inverse(out.as_mut(), self.L(), self.D(), par, stack);

		make_self_adjoint(out.as_mut());
		out
	}
}

impl<T: ComplexField> SolveCore<T> for Lblt<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		let mut mem = MemBuffer::new(linalg::cholesky::lblt::solve::solve_in_place_scratch::<usize, T>(
			self.L.nrows(),
			rhs.ncols(),
			par,
		));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::lblt::solve::solve_in_place_with_conj(self.L.as_ref(), self.B_diag(), self.B_subdiag(), conj, self.P(), rhs, par, stack);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		let mut mem = MemBuffer::new(linalg::cholesky::lblt::solve::solve_in_place_scratch::<usize, T>(
			self.L.nrows(),
			rhs.ncols(),
			par,
		));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::lblt::solve::solve_in_place_with_conj(
			self.L(),
			self.B_diag(),
			self.B_subdiag(),
			conj.compose(Conj::Yes),
			self.P(),
			rhs,
			par,
			stack,
		);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for Lblt<T> {
	#[track_caller]
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();

		let n = self.L.nrows();
		let mut out = Mat::zeros(n, n);

		let mut mem = MemBuffer::new(linalg::cholesky::lblt::reconstruct::reconstruct_scratch::<usize, T>(n, par));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::lblt::reconstruct::reconstruct(out.as_mut(), self.L(), self.B_diag(), self.B_subdiag(), self.P(), par, stack);

		make_self_adjoint(out.as_mut());
		out
	}

	#[track_caller]
	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		let n = self.L.nrows();
		let mut out = Mat::zeros(n, n);

		let mut mem = MemBuffer::new(linalg::cholesky::lblt::inverse::inverse_scratch::<usize, T>(n, par));
		let stack = MemStack::new(&mut mem);

		linalg::cholesky::lblt::inverse::inverse(out.as_mut(), self.L(), self.B_diag(), self.B_subdiag(), self.P(), par, stack);

		make_self_adjoint(out.as_mut());
		out
	}
}

impl<T: ComplexField> SolveCore<T> for PartialPivLu<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.nrows() == rhs.nrows(),));

		let k = rhs.ncols();

		linalg::lu::partial_pivoting::solve::solve_in_place_with_conj(
			self.L(),
			self.U(),
			self.P(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::lu::partial_pivoting::solve::solve_in_place_scratch::<usize, T>(self.nrows(), k, par),
			)),
		);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.ncols() == rhs.nrows(),));

		let k = rhs.ncols();

		linalg::lu::partial_pivoting::solve::solve_transpose_in_place_with_conj(
			self.L(),
			self.U(),
			self.P(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::lu::partial_pivoting::solve::solve_transpose_in_place_scratch::<usize, T>(self.nrows(), k, par),
			)),
		);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for PartialPivLu<T> {
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();
		let m = self.nrows();
		let n = self.ncols();

		let mut out = Mat::zeros(m, n);

		linalg::lu::partial_pivoting::reconstruct::reconstruct(
			out.as_mut(),
			self.L(),
			self.U(),
			self.P(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::lu::partial_pivoting::reconstruct::reconstruct_scratch::<
				usize,
				T,
			>(m, n, par))),
		);

		out
	}

	#[track_caller]
	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		assert!(self.nrows() == self.ncols());

		let n = self.ncols();

		let mut out = Mat::zeros(n, n);

		linalg::lu::partial_pivoting::inverse::inverse(
			out.as_mut(),
			self.L(),
			self.U(),
			self.P(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::lu::partial_pivoting::inverse::inverse_scratch::<usize, T>(
				n, par,
			))),
		);

		out
	}
}

impl<T: ComplexField> SolveCore<T> for FullPivLu<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.nrows() == rhs.nrows(),));

		let k = rhs.ncols();

		linalg::lu::full_pivoting::solve::solve_in_place_with_conj(
			self.L(),
			self.U(),
			self.P(),
			self.Q(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::lu::full_pivoting::solve::solve_in_place_scratch::<usize, T>(
				self.nrows(),
				k,
				par,
			))),
		);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.ncols() == rhs.nrows(),));

		let k = rhs.ncols();

		linalg::lu::full_pivoting::solve::solve_transpose_in_place_with_conj(
			self.L(),
			self.U(),
			self.P(),
			self.Q(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::lu::full_pivoting::solve::solve_transpose_in_place_scratch::<
				usize,
				T,
			>(self.nrows(), k, par))),
		);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for FullPivLu<T> {
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();
		let m = self.nrows();
		let n = self.ncols();

		let mut out = Mat::zeros(m, n);

		linalg::lu::full_pivoting::reconstruct::reconstruct(
			out.as_mut(),
			self.L(),
			self.U(),
			self.P(),
			self.Q(),
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::lu::full_pivoting::reconstruct::reconstruct_scratch::<usize, T>(m, n, par),
			)),
		);

		out
	}

	#[track_caller]
	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		assert!(self.nrows() == self.ncols());

		let n = self.ncols();

		let mut out = Mat::zeros(n, n);

		linalg::lu::full_pivoting::inverse::inverse(
			out.as_mut(),
			self.L(),
			self.U(),
			self.P(),
			self.Q(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::lu::full_pivoting::inverse::inverse_scratch::<usize, T>(
				n, par,
			))),
		);

		out
	}
}

impl<T: ComplexField> SolveCore<T> for Qr<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.nrows() == rhs.nrows(),));

		let n = self.nrows();
		let block_size = self.Q_coeff().nrows();
		let k = rhs.ncols();

		linalg::qr::no_pivoting::solve::solve_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::no_pivoting::solve::solve_in_place_scratch::<T>(
				n, block_size, k, par,
			))),
		);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.ncols() == rhs.nrows(),));

		let n = self.nrows();
		let block_size = self.Q_coeff().nrows();
		let k = rhs.ncols();

		linalg::qr::no_pivoting::solve::solve_transpose_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::qr::no_pivoting::solve::solve_transpose_in_place_scratch::<T>(n, block_size, k, par),
			)),
		);
	}
}

impl<T: ComplexField> SolveLstsqCore<T> for Qr<T> {
	#[track_caller]
	fn solve_lstsq_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == rhs.nrows(), self.nrows() >= self.ncols(),));

		let m = self.nrows();
		let n = self.ncols();
		let block_size = self.Q_coeff().nrows();
		let k = rhs.ncols();

		linalg::qr::no_pivoting::solve::solve_lstsq_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::no_pivoting::solve::solve_lstsq_in_place_scratch::<T>(
				m, n, block_size, k, par,
			))),
		);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for Qr<T> {
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();
		let m = self.nrows();
		let n = self.ncols();
		let block_size = self.Q_coeff().nrows();

		let mut out = Mat::zeros(m, n);

		linalg::qr::no_pivoting::reconstruct::reconstruct(
			out.as_mut(),
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::no_pivoting::reconstruct::reconstruct_scratch::<T>(
				m, n, block_size, par,
			))),
		);

		out
	}

	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();
		assert!(self.nrows() == self.ncols());

		let n = self.ncols();
		let block_size = self.Q_coeff().nrows();

		let mut out = Mat::zeros(n, n);

		linalg::qr::no_pivoting::inverse::inverse(
			out.as_mut(),
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::no_pivoting::inverse::inverse_scratch::<T>(
				n, block_size, par,
			))),
		);

		out
	}
}

impl<T: ComplexField> SolveCore<T> for ColPivQr<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.nrows() == rhs.nrows(),));

		let n = self.nrows();
		let block_size = self.Q_coeff().nrows();
		let k = rhs.ncols();

		linalg::qr::col_pivoting::solve::solve_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			self.P(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::col_pivoting::solve::solve_in_place_scratch::<usize, T>(
				n, block_size, k, par,
			))),
		);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.ncols() == rhs.nrows(),));

		let n = self.nrows();
		let block_size = self.Q_coeff().nrows();
		let k = rhs.ncols();

		linalg::qr::col_pivoting::solve::solve_transpose_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			self.P(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::col_pivoting::solve::solve_transpose_in_place_scratch::<
				usize,
				T,
			>(n, block_size, k, par))),
		);
	}
}

impl<T: ComplexField> SolveLstsqCore<T> for ColPivQr<T> {
	#[track_caller]
	fn solve_lstsq_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == rhs.nrows(), self.nrows() >= self.ncols(),));

		let m = self.nrows();
		let n = self.ncols();
		let block_size = self.Q_coeff().nrows();
		let k = rhs.ncols();

		linalg::qr::col_pivoting::solve::solve_lstsq_in_place_with_conj(
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			self.P(),
			conj,
			rhs,
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::col_pivoting::solve::solve_lstsq_in_place_scratch::<
				usize,
				T,
			>(m, n, block_size, k, par))),
		);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for ColPivQr<T> {
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();
		let m = self.nrows();
		let n = self.ncols();
		let block_size = self.Q_coeff().nrows();

		let mut out = Mat::zeros(m, n);

		linalg::qr::col_pivoting::reconstruct::reconstruct(
			out.as_mut(),
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			self.P(),
			par,
			MemStack::new(&mut MemBuffer::new(
				linalg::qr::col_pivoting::reconstruct::reconstruct_scratch::<usize, T>(m, n, block_size, par),
			)),
		);

		out
	}

	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();
		assert!(self.nrows() == self.ncols());

		let n = self.ncols();
		let block_size = self.Q_coeff().nrows();

		let mut out = Mat::zeros(n, n);

		linalg::qr::col_pivoting::inverse::inverse(
			out.as_mut(),
			self.Q_basis(),
			self.Q_coeff(),
			self.R(),
			self.P(),
			par,
			MemStack::new(&mut MemBuffer::new(linalg::qr::col_pivoting::inverse::inverse_scratch::<usize, T>(
				n, block_size, par,
			))),
		);

		out
	}
}

impl<T: ComplexField> SolveCore<T> for Svd<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.nrows() == rhs.nrows(),));

		let mut rhs = rhs;
		let n = self.nrows();
		let k = rhs.ncols();
		let mut tmp = Mat::zeros(n, k);

		linalg::matmul::matmul_with_conj(
			tmp.as_mut(),
			Accum::Replace,
			self.U().transpose(),
			conj.compose(Conj::Yes),
			rhs.as_ref(),
			Conj::No,
			one(),
			par,
		);

		for j in 0..k {
			for i in 0..n {
				let s = recip(&real(&self.S()[i]));
				tmp[(i, j)] = mul_real(&tmp[(i, j)], &s);
			}
		}

		linalg::matmul::matmul_with_conj(rhs.as_mut(), Accum::Replace, self.V(), conj, tmp.as_ref(), Conj::No, one(), par);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.ncols() == rhs.nrows(),));

		let mut rhs = rhs;
		let n = self.nrows();
		let k = rhs.ncols();
		let mut tmp = Mat::zeros(n, k);

		linalg::matmul::matmul_with_conj(
			tmp.as_mut(),
			Accum::Replace,
			self.V().transpose(),
			conj,
			rhs.as_ref(),
			Conj::No,
			one(),
			par,
		);

		for j in 0..k {
			for i in 0..n {
				let s = recip(&real(&self.S()[i]));
				tmp[(i, j)] = mul_real(&tmp[(i, j)], &s);
			}
		}

		linalg::matmul::matmul_with_conj(
			rhs.as_mut(),
			Accum::Replace,
			self.U(),
			conj.compose(Conj::Yes),
			tmp.as_ref(),
			Conj::No,
			one(),
			par,
		);
	}
}

impl<T: ComplexField> SolveLstsqCore<T> for Svd<T> {
	#[track_caller]
	fn solve_lstsq_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == rhs.nrows(), self.nrows() >= self.ncols(),));

		let m = self.nrows();
		let n = self.ncols();

		let size = Ord::min(m, n);

		let U = self.U().get(.., ..size);
		let V = self.V().get(.., ..size);

		let k = rhs.ncols();

		let mut tmp = Mat::zeros(size, k);

		linalg::matmul::matmul_with_conj(
			tmp.as_mut(),
			Accum::Replace,
			U.transpose(),
			conj.compose(Conj::Yes),
			rhs.as_ref(),
			Conj::No,
			one(),
			par,
		);

		for j in 0..k {
			for i in 0..size {
				let s = recip(&real(&self.S()[i]));
				tmp[(i, j)] = mul_real(&tmp[(i, j)], &s);
			}
		}

		linalg::matmul::matmul_with_conj(rhs.get_mut(..size, ..), Accum::Replace, V, conj, tmp.as_ref(), Conj::No, one(), par);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for Svd<T> {
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();
		let m = self.nrows();
		let n = self.ncols();

		let size = Ord::min(m, n);

		let U = self.U().get(.., ..size);
		let V = self.V().get(.., ..size);
		let S = self.S();

		let mut UxS = Mat::zeros(m, size);
		for j in 0..size {
			let s = real(&S[j]);
			for i in 0..m {
				UxS[(i, j)] = mul_real(&U[(i, j)], &s);
			}
		}

		let mut out = Mat::zeros(m, n);

		linalg::matmul::matmul(out.as_mut(), Accum::Replace, UxS.as_ref(), V.adjoint(), one(), par);

		out
	}

	#[track_caller]
	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		assert!(self.nrows() == self.ncols());
		let n = self.nrows();

		let U = self.U();
		let V = self.V();
		let S = self.S();

		let mut VxS = Mat::zeros(n, n);
		for j in 0..n {
			let s = recip(&real(&S[j]));

			for i in 0..n {
				VxS[(i, j)] = mul_real(&V[(i, j)], &s);
			}
		}

		let mut out = Mat::zeros(n, n);

		linalg::matmul::matmul(out.as_mut(), Accum::Replace, VxS.as_ref(), U.adjoint(), one(), par);

		out
	}
}

impl<T: ComplexField> SolveCore<T> for SelfAdjointEigen<T> {
	#[track_caller]
	fn solve_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.nrows() == rhs.nrows(),));

		let mut rhs = rhs;
		let n = self.nrows();
		let k = rhs.ncols();
		let mut tmp = Mat::zeros(n, k);

		linalg::matmul::matmul_with_conj(
			tmp.as_mut(),
			Accum::Replace,
			self.U().transpose(),
			conj.compose(Conj::Yes),
			rhs.as_ref(),
			Conj::No,
			one(),
			par,
		);

		for j in 0..k {
			for i in 0..n {
				let s = recip(&real(&self.S()[i]));
				tmp[(i, j)] = mul_real(&tmp[(i, j)], &s);
			}
		}

		linalg::matmul::matmul_with_conj(rhs.as_mut(), Accum::Replace, self.U(), conj, tmp.as_ref(), Conj::No, one(), par);
	}

	#[track_caller]
	fn solve_transpose_in_place_with_conj(&self, conj: Conj, rhs: MatMut<'_, T>) {
		let par = get_global_parallelism();

		assert!(all(self.nrows() == self.ncols(), self.ncols() == rhs.nrows(),));

		let mut rhs = rhs;
		let n = self.nrows();
		let k = rhs.ncols();
		let mut tmp = Mat::zeros(n, k);

		linalg::matmul::matmul_with_conj(
			tmp.as_mut(),
			Accum::Replace,
			self.U().transpose(),
			conj,
			rhs.as_ref(),
			Conj::No,
			one(),
			par,
		);

		for j in 0..k {
			for i in 0..n {
				let s = recip(&real(&self.S()[i]));
				tmp[(i, j)] = mul_real(&tmp[(i, j)], &s);
			}
		}

		linalg::matmul::matmul_with_conj(
			rhs.as_mut(),
			Accum::Replace,
			self.U(),
			conj.compose(Conj::Yes),
			tmp.as_ref(),
			Conj::No,
			one(),
			par,
		);
	}
}

impl<T: ComplexField> DenseSolveCore<T> for SelfAdjointEigen<T> {
	fn reconstruct(&self) -> Mat<T> {
		let par = get_global_parallelism();
		let m = self.nrows();
		let n = self.ncols();

		let size = Ord::min(m, n);

		let U = self.U().get(.., ..size);
		let V = self.U().get(.., ..size);
		let S = self.S();

		let mut UxS = Mat::zeros(m, size);
		for j in 0..size {
			let s = real(&S[j]);
			for i in 0..m {
				UxS[(i, j)] = mul_real(&U[(i, j)], &s);
			}
		}

		let mut out = Mat::zeros(m, n);

		linalg::matmul::matmul(out.as_mut(), Accum::Replace, UxS.as_ref(), V.adjoint(), one(), par);

		out
	}

	fn inverse(&self) -> Mat<T> {
		let par = get_global_parallelism();

		assert!(self.nrows() == self.ncols());
		let n = self.nrows();

		let U = self.U();
		let V = self.U();
		let S = self.S();

		let mut VxS = Mat::zeros(n, n);
		for j in 0..n {
			let s = recip(&real(&S[j]));

			for i in 0..n {
				VxS[(i, j)] = mul_real(&V[(i, j)], &s);
			}
		}

		let mut out = Mat::zeros(n, n);

		linalg::matmul::matmul(out.as_mut(), Accum::Replace, VxS.as_ref(), U.adjoint(), one(), par);

		out
	}
}

#[cfg(test)]
mod tests {
	use super::*;
	use crate::assert;
	use crate::stats::prelude::*;
	use crate::utils::approx::*;

	#[track_caller]
	fn test_solver(A: MatRef<'_, c64>, A_dec: impl SolveCore<c64>) {
		#[track_caller]
		fn test_solver_imp(A: MatRef<'_, c64>, A_dec: &dyn SolveCore<c64>) {
			let rng = &mut StdRng::seed_from_u64(0xC0FFEE);

			let n = A.nrows();
			let approx_eq = CwiseMat(ApproxEq::eps() * 128.0 * (n as f64));

			let k = 3;

			let ref R = CwiseMatDistribution {
				nrows: n,
				ncols: k,
				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
			}
			.rand::<Mat<c64>>(rng);

			let ref L = CwiseMatDistribution {
				nrows: k,
				ncols: n,
				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
			}
			.rand::<Mat<c64>>(rng);

			assert!(A * A_dec.solve(R) ~ R);
			assert!(A.conjugate() * A_dec.solve_conjugate(R) ~ R);
			assert!(A.transpose() * A_dec.solve_transpose(R) ~ R);
			assert!(A.adjoint() * A_dec.solve_adjoint(R) ~ R);

			assert!(A_dec.rsolve(L) * A ~ L);
			assert!(A_dec.rsolve_conjugate(L) * A.conjugate() ~ L);
			assert!(A_dec.rsolve_transpose(L) * A.transpose() ~ L);
			assert!(A_dec.rsolve_adjoint(L) * A.adjoint() ~ L);
		}

		test_solver_imp(A, &A_dec)
	}

	#[test]
	fn test_all_solvers() {
		let rng = &mut StdRng::seed_from_u64(0);
		let n = 50;

		let ref A = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		}
		.rand::<Mat<c64>>(rng);
		let A = A.rb();

		test_solver(A, A.partial_piv_lu());
		test_solver(A, A.full_piv_lu());
		test_solver(A, A.qr());
		test_solver(A, A.col_piv_qr());
		test_solver(A, A.svd().unwrap());

		{
			let ref A = A * A.adjoint();
			let A = A.rb();
			test_solver(A, A.llt(Side::Lower).unwrap());
			test_solver(A, A.ldlt(Side::Lower).unwrap());
		}

		{
			let ref A = A + A.adjoint();
			let A = A.rb();
			test_solver(A, A.lblt(Side::Lower));
			test_solver(A, A.self_adjoint_eigen(Side::Lower).unwrap());
		}
	}

	#[test]
	fn test_eigen_cplx() {
		let rng = &mut StdRng::seed_from_u64(0);
		let n = 50;

		let A = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		}
		.rand::<Mat<c64>>(rng);

		let n = A.nrows();
		let approx_eq = CwiseMat(ApproxEq::eps() * 128.0 * (n as f64));

		{
			let evd = A.eigen().unwrap();
			let e = A.eigenvalues().unwrap();
			assert!(&A * evd.U() ~ evd.U() * evd.S());
			assert!(evd.S().column_vector() ~ ColRef::from_slice(&e));
		}
		{
			let evd = A.conjugate().eigen().unwrap();
			let e = A.conjugate().eigenvalues().unwrap();
			assert!(A.conjugate() * evd.U() ~ evd.U() * evd.S());
			assert!(evd.S().column_vector() ~ ColRef::from_slice(&e));
		}
	}

	#[test]
	fn test_geigen_cplx() {
		let rng = &mut StdRng::seed_from_u64(0);
		let n = 50;

		let A = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		}
		.rand::<Mat<c64>>(rng);

		let B = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		}
		.rand::<Mat<c64>>(rng);

		let n = A.nrows();
		let approx_eq = CwiseMat(ApproxEq::eps() * 128.0 * (n as f64));

		{
			let evd = A.generalized_eigen(&B).unwrap();
			let e = zip!(evd.S_a(), evd.S_b()).map(|unzip!(a, b)| a / b);
			assert!(&A * evd.U() ~ &B * evd.U() * e);
		}

		{
			let evd = A.conjugate().generalized_eigen(B.conjugate()).unwrap();
			let e = zip!(evd.S_a(), evd.S_b()).map(|unzip!(a, b)| a / b);
			assert!(A.conjugate() * evd.U() ~ B.conjugate() * evd.U() * e);
		}
	}

	#[test]
	fn test_eigen_real() {
		let rng = &mut StdRng::seed_from_u64(0);
		let n = 50;

		let A = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: StandardNormal,
		}
		.rand::<Mat<f64>>(rng);

		let n = A.nrows();
		let approx_eq = CwiseMat(ApproxEq::eps() * 128.0 * (n as f64));

		let evd = A.eigen().unwrap();
		let e = A.eigenvalues().unwrap();

		let A = Mat::from_fn(A.nrows(), A.ncols(), |i, j| c64::from(A[(i, j)]));

		assert!(&A * evd.U() ~ evd.U() * evd.S());
		assert!(evd.S().column_vector() ~ ColRef::from_slice(&e));
	}

	#[test]
	fn test_geigen_real() {
		let rng = &mut StdRng::seed_from_u64(0);
		let n = 50;

		let A = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: StandardNormal,
		}
		.rand::<Mat<f64>>(rng);

		let B = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: StandardNormal,
		}
		.rand::<Mat<f64>>(rng);

		let n = A.nrows();
		let approx_eq = CwiseMat(ApproxEq::eps() * 128.0 * (n as f64));

		let Ac = zip!(&A).map(|unzip!(x)| c64::new(*x, 0.0));
		let Bc = zip!(&B).map(|unzip!(x)| c64::new(*x, 0.0));

		{
			let evd = A.generalized_eigen(&B).unwrap();
			let e = zip!(evd.S_a(), evd.S_b()).map(|unzip!(a, b)| a / b);
			assert!(&Ac * evd.U() ~ &Bc * evd.U() * e);
		}
	}

	#[test]
	fn test_svd_solver_for_rectangular_matrix() {
		#[rustfmt::skip]
    	let A = crate::mat![
    	    [4.,   5.,   7.],
    	    [8.,   8.,   2.],
    	    [4.,   0.,   9.],
    	    [2.,   6.,   2.],
    	    [0.,   6.,   0.],
    	];
		#[rustfmt::skip]
    	let B = crate::mat![
        	[105.,    49.],
        	[ 98.,    54.],
        	[113.,    35.],
        	[ 46.,    34.],
        	[ 12.,    24.],
     	];

		#[rustfmt::skip]
	    let X_true= crate::mat![
	      [8.,   2.],
	      [2.,   4.],
	      [9.,   3.],
	    ];

		let approx_eq = CwiseMat(ApproxEq::eps() * 128.0 * (A.nrows() as f64));
		let svd = A.svd().unwrap();
		let mut X = B.cloned();
		svd.solve_lstsq_in_place_with_conj(crate::Conj::No, X.as_mat_mut());
		assert!(X.get(..X_true.nrows(),..) ~ X_true);
	}
}