faer 0.24.0

use super::*;
use crate::assert;
use crate::linalg::matmul::matmul;
use linalg::evd::schur;

const MIN_DIM: usize = 32;

/// partial eigendecomposition tuning parameters.
#[derive(Debug, Copy, Clone)]
pub struct PartialEigenParams {
	/// minimum projection subspace dimension.
	pub min_dim: usize,
	/// maximum projection subspace dimension.
	pub max_dim: usize,
	/// maximum number of algorithm restarts.
	pub max_restarts: usize,
	#[doc(hidden)]
	pub non_exhaustive: NonExhaustive,
}

/// partial eigendecomposition tuning parameters.
#[derive(Debug, Copy, Clone)]
pub struct PartialEigenInfo {
	/// number of converged eigenvalues and eigenvectors.
	pub n_converged_eigen: usize,
	#[doc(hidden)]
	pub non_exhaustive: NonExhaustive,
}

impl Default for PartialEigenParams {
	fn default() -> Self {
		Self {
			min_dim: 0,
			max_dim: 0,
			max_restarts: 1000,
			non_exhaustive: NonExhaustive(()),
		}
	}
}

fn iterate_arnoldi<T: ComplexField>(
	A: &dyn LinOp<T>,
	H: MatMut<'_, T>,
	V: MatMut<'_, T>,
	start: usize,
	end: usize,
	par: Par,
	stack: &mut MemStack,
) {
	let mut V = V;
	let mut H = H;
	for j in start..end + 1 {
		let mut H = H.rb_mut().col_mut(j - 1);
		H.fill(zero());
		let (V, Vnext) = V.rb_mut().split_at_col_mut(j);
		let V = V.rb();
		let mut Vnext = Vnext.col_mut(0);
		A.apply(
			Vnext.rb_mut().as_mat_mut(),
			V.col(j - 1).as_mat(),
			par,
			stack,
		);
		let (mut converged, _) = stack.collect(core::iter::repeat_n(false, j));
		let mut h = H.rb_mut().get_mut(..j);
		for i in 0..j {
			let r = V.col(i).adjoint() * Vnext.rb();
			zip!(Vnext.rb_mut(), V.col(i))
				.for_each(|unzip!(y, x): Zip!(&mut T, &T)| *y -= &r * x);
			h[i] = r;
		}
		let ref f =
			from_f64::<T::Real>(Ord::max(j, 8) as f64) * eps::<T::Real>();
		loop {
			let mut all_true = true;
			for i in 0..j {
				if !converged[i] {
					all_true = false;
					let ref r = V.col(i).adjoint() * Vnext.rb();
					zip!(Vnext.rb_mut(), V.col(i))
						.for_each(|unzip!(y, x): Zip!(&mut T, &T)| *y -= r * x);
					h[i] += r;
					converged[i] = r.abs() < f * Vnext.norm_l2();
				}
			}
			if all_true {
				break;
			}
		}
		let norm = Vnext.norm_l2();
		if norm > zero() {
			let norm_inv = norm.recip();
			zip!(&mut Vnext).for_each(|unzip!(v)| *v = v.mul_real(&norm_inv));
		}
		H[j] = norm.to_cplx();
	}
}

fn schur_swap<T: ComplexField>(
	a: MatMut<T>,
	q: Option<MatMut<T>>,
	j0: usize,
	n1: usize,
	n2: usize,
) -> isize {
	if const { T::IS_REAL } {
		unsafe {
			schur::real_schur::schur_swap::<T::Real>(
				core::mem::transmute(a),
				core::mem::transmute(q),
				j0,
				n1,
				n2,
			)
		}
	} else {
		assert!(all(n1 == 1, n2 == 1));
		schur::complex_schur::schur_swap(a, q, j0)
	}
}

fn reorder_schur<T: ComplexField>(
	mut A: MatMut<'_, T>,
	mut Q: Option<MatMut<'_, T>>,
	mut ifst: usize,
	mut ilst: usize,
) {
	let zero = zero::<T>();
	let n = A.nrows();
	let mut nbf = 1;
	if const { T::IS_REAL } {
		if ifst > 0 {
			if A[(ifst, ifst - 1)] != zero {
				ifst -= 1;
			}
		}
		if ifst < n - 1 {
			if A[(ifst + 1, ifst)] != zero {
				nbf = 2;
			}
		}
	}
	let mut nbl = 1;
	if const { T::IS_REAL } {
		if ilst > 0 {
			if A[(ilst, ilst - 1)] != zero {
				ilst = ilst - 1;
			}
		}
		if ilst < n - 1 {
			if A[(ilst + 1, ilst)] != zero {
				nbl = 2;
			}
		}
	}
	if ifst == ilst {
		return;
	}
	if ifst < ilst {
		if nbf == 2 && nbl == 1 {
			ilst -= 1;
		}
		if nbf == 1 && nbl == 2 {
			ilst += 1;
		}
		let mut here = ifst;
		loop {
			if nbf == 1 || nbf == 2 {
				let mut nbnext = 1;
				if const { T::IS_REAL } {
					if here + nbf + 1 <= n - 1 {
						if A[(here + nbf + 1, here + nbf)] != zero {
							nbnext = 2;
						}
					}
				}
				schur_swap(A.rb_mut(), Q.rb_mut(), here, nbf, nbnext);
				here += nbnext;
				if const { T::IS_REAL } {
					if nbf == 2 {
						if A[(here + 1, here)] == zero {
							nbf = 3;
						}
					}
				}
			} else if const { T::IS_REAL } {
				let mut nbnext = 1;
				if here + 3 <= n - 1 {
					if A[(here + 3, here + 2)] != zero {
						nbnext = 2;
					}
				}
				schur_swap(A.rb_mut(), Q.rb_mut(), here + 1, 1, nbnext);
				if nbnext == 1 {
					schur_swap(A.rb_mut(), Q.rb_mut(), here, 1, 1);
					here += 1;
				} else {
					if A[(here + 2, here + 1)] == zero {
						nbnext = 1;
					}
					if nbnext == 2 {
						schur_swap(A.rb_mut(), Q.rb_mut(), here, 1, nbnext);
						here += 2;
					} else {
						schur_swap(A.rb_mut(), Q.rb_mut(), here, 1, 1);
						here += 1;
						schur_swap(A.rb_mut(), Q.rb_mut(), here, 1, 1);
						here += 1;
					}
				}
			}
			if here >= ilst {
				break;
			}
		}
	} else {
		let mut here = ifst;
		loop {
			if nbf == 1 || nbf == 2 {
				let mut nbnext = 1;
				if const { T::IS_REAL } {
					if here >= 2 {
						if A[(here - 1, here - 2)] != zero {
							nbnext = 2;
						}
					}
				}
				schur_swap(A.rb_mut(), Q.rb_mut(), here - nbnext, nbnext, nbf);
				here -= nbnext;
				if const { T::IS_REAL } {
					if nbf == 2 {
						if A[(here + 1, here)] == zero {
							nbf = 3;
						}
					}
				}
			} else if const { T::IS_REAL } {
				let mut nbnext = 1;
				if here >= 2 {
					if A[(here - 1, here - 2)] != zero {
						nbnext = 2;
					}
				}
				schur_swap(A.rb_mut(), Q.rb_mut(), here - nbnext, nbnext, 1);
				if nbnext == 1 {
					schur_swap(A.rb_mut(), Q.rb_mut(), here, nbnext, 1);
					here -= 1;
				} else {
					if A[(here, here - 1)] == zero {
						nbnext = 1;
					}
					if nbnext == 2 {
						schur_swap(A.rb_mut(), Q.rb_mut(), here - 1, 2, 1);
						here -= 2;
					} else {
						schur_swap(A.rb_mut(), Q.rb_mut(), here, 1, 1);
						here -= 1;
						schur_swap(A.rb_mut(), Q.rb_mut(), here, 1, 1);
						here -= 1;
					}
				}
			}
			if here <= ilst {
				break;
			}
		}
	}
}

fn partial_schur_real_imp<T: RealField>(
	eigvecs: MatMut<'_, Complex<T>>,
	eigvals: &mut [Complex<T>],
	A: &dyn LinOp<T>,
	v0: ColRef<'_, T>,
	min_dim: usize,
	max_dim: usize,
	n_eigval: usize,
	tol: T,
	restarts: usize,
	par: Par,
	stack: &mut MemStack,
) -> usize {
	let n = A.nrows();
	let (mut H, stack) =
		temp_mat_zeroed::<T, _, _>(max_dim + 1, max_dim, stack);
	let mut H = H.as_mat_mut();
	let (mut V, stack) = temp_mat_zeroed::<T, _, _>(n, max_dim + 1, stack);
	let mut V = V.as_mat_mut();
	let (mut X, stack) = temp_mat_zeroed::<T, _, _>(max_dim, max_dim, stack);
	let mut X = X.as_mat_mut();
	let (mut vecs, stack) = temp_mat_zeroed::<T, _, _>(max_dim, max_dim, stack);
	let mut vecs = vecs.as_mat_mut();
	let (mut tmp, stack) = temp_mat_zeroed::<T, _, _>(n, max_dim, stack);
	let mut tmp = tmp.as_mat_mut();
	let mut active = 0usize;
	if max_dim < n {
		let (mut Q, stack) =
			temp_mat_zeroed::<T, _, _>(max_dim, max_dim, stack);
		let mut Q = Q.as_mat_mut();
		let (mut residual, stack) =
			temp_mat_zeroed::<T, _, _>(max_dim, 1, stack);
		let mut residual = residual.as_mat_mut().col_mut(0);
		let (mut w, stack) = temp_mat_zeroed::<T, _, _>(max_dim, 2, stack);
		let mut w = w.as_mat_mut();
		let block_size =
			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
				max_dim, max_dim,
			);
		let (mut householder, stack) =
			temp_mat_zeroed::<T, _, _>(block_size, max_dim, stack);
		let mut householder = householder.as_mat_mut();
		let f = v0.norm_l2();
		if f > min_positive() {
			let ref f = f.recip();
			zip!(V.rb_mut().col_mut(0), v0)
				.for_each(|unzip!(y, x): Zip!(&mut T, &T)| *y = f * x);
		} else {
			let n0 = n as u32;
			let n1 = (n >> 32) as u32;
			let n = from_f64::<T>(n0 as f64) + from_f64::<T>(n1 as f64);
			let f = n.sqrt().recip();
			zip!(V.rb_mut().col_mut(0)).for_each(|unzip!(y)| *y = f.copy());
		}
		iterate_arnoldi(A, H.as_mut(), V.as_mut(), 1, min_dim, par, stack);
		let mut k = min_dim;
		for iter in 0..restarts {
			_ = iter;
			iterate_arnoldi(
				A,
				H.as_mut(),
				V.as_mut(),
				k + 1,
				max_dim,
				par,
				stack,
			);
			let ref Hmm = H[(max_dim, max_dim - 1)].abs();
			let n = max_dim - active;
			let (mut w_re, mut w_im) =
				w.rb_mut().get_mut(active..max_dim, ..).two_cols_mut(0, 1);
			Q.fill(zero());
			Q.rb_mut().diagonal_mut().fill(one());
			let mut Q_slice =
				Q.rb_mut().get_mut(active..max_dim, active..max_dim);
			let mut H_slice =
				H.rb_mut().get_mut(active..max_dim, active..max_dim);
			{
				let n = max_dim - active;
				let block_size =
					linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
						n, n,
					);
				let mut householder =
					householder.rb_mut().get_mut(..block_size, ..n - 1);
				linalg::evd::hessenberg::hessenberg_in_place(
					H_slice.rb_mut(),
					householder.rb_mut(),
					par,
					stack,
					default(),
				);
				linalg::householder::apply_block_householder_sequence_on_the_right_in_place_with_conj(
					H_slice.rb().submatrix(1, 0, n - 1, n - 1),
					householder.rb(),
					Conj::No,
					Q_slice.rb_mut().submatrix_mut(1, 1, n - 1, n - 1),
					par,
					stack,
				);
				for j in 0..n {
					for i in j + 2..n {
						H_slice[(i, j)] = zero();
					}
				}
			}
			schur::real_schur::multishift_qr(
				true,
				H_slice.rb_mut(),
				Some(Q_slice.rb_mut()),
				w_re.rb_mut(),
				w_im.rb_mut(),
				0,
				n,
				par,
				stack,
				auto!(T),
			);
			let mut j = 0usize;
			while j < n {
				let mut i = j;
				let mut idx = i;
				let mut max = zero::<T>();
				while i < n {
					let cplx = i + 1 < n && H_slice[(i + 1, i)] != zero::<T>();
					let (v, bs) = if cplx {
						(w_re[i].hypot(&w_im[i]), 2)
					} else {
						(w_re[i].abs(), 1)
					};
					if v > max {
						max = v;
						idx = i;
					}
					i += bs;
				}
				let i = idx;
				let cplx = i + 1 < n && H_slice[(i + 1, i)] != zero::<T>();
				let bs = if cplx { 2 } else { 1 };
				if i != j {
					reorder_schur(
						H_slice.rb_mut(),
						Some(Q_slice.rb_mut()),
						i,
						j,
					);
					let x_re = w_re
						.rb_mut()
						.try_as_col_major_mut()
						.unwrap()
						.as_slice_mut();
					let x_im = w_im
						.rb_mut()
						.try_as_col_major_mut()
						.unwrap()
						.as_slice_mut();
					for x in [x_re, x_im] {
						x[j..i + bs].rotate_right(bs)
					}
				}
				j += bs;
			}
			let mut X = X.rb_mut().get_mut(..n, ..n);
			linalg::evd::evd_from_real_schur_imp(
				H_slice.rb(),
				X.as_mut(),
				par,
				auto!(T),
			);
			let mut vecs = vecs.rb_mut().get_mut(..n, ..n);
			matmul(
				vecs.rb_mut(),
				Accum::Replace,
				Q_slice.rb(),
				X.rb(),
				one(),
				par,
			);
			let vecs = vecs.rb();
			let mut H_tmp = tmp.rb_mut().get_mut(..active, ..);
			matmul(
				H_tmp.rb_mut(),
				Accum::Replace,
				H.rb().get(..active, ..),
				Q.rb(),
				one(),
				par,
			);
			H.rb_mut().get_mut(..active, ..).copy_from(&H_tmp);
			let mut j = 0usize;
			while j < n {
				let re = &vecs[(max_dim - active - 1, j)];
				if w_im[j] != zero::<T>() {
					let im = &vecs[(max_dim - active - 1, j + 1)];
					let res = Hmm * re.hypot(im);
					residual[active + j] = res.copy();
					residual[active + j + 1] = res;
					j += 2;
				} else {
					residual[active + j] = Hmm * re.abs();
					j += 1;
				}
			}
			#[derive(Copy, Clone, PartialEq, Eq, Debug)]
			enum Group {
				Lock,
				Retain,
				Purge,
			}
			let mut groups = alloc::vec![Group::Purge; max_dim];
			let mut nev = n_eigval;
			if w_im[nev - active - 1] > zero::<T>() {
				nev += 1;
			}
			let mut nlock = 0usize;
			for j in 0..nev {
				if residual[j] <= tol {
					groups[j] = Group::Lock;
					nlock += 1;
				} else {
					groups[j] = Group::Retain;
				}
			}
			let ideal_size = Ord::min(nlock + min_dim, (min_dim + max_dim) / 2);
			k = nev;
			let mut i = nev;
			while i < max_dim {
				let cplx = w[(i, 1)] != zero::<T>();
				let bs = if cplx { 2 } else { 1 };
				let group;
				if k < ideal_size && residual[i] > tol {
					group = Group::Retain;
					k += bs;
				} else {
					group = Group::Purge;
				}
				for k in 0..bs {
					groups[i + k] = group;
				}
				i += bs;
			}
			let mut purge = 0usize;
			while purge < active && groups[purge] == Group::Lock {
				purge += 1;
			}
			let mut lo = 0usize;
			let mut mi = 0usize;
			let mut hi = 0usize;
			while hi < max_dim {
				let cplx = hi + 1 < max_dim && H[(hi + 1, hi)] != zero::<T>();
				let bs = if cplx { 2 } else { 1 };
				match groups[hi] {
					Group::Lock => {
						reorder_schur(
							H.rb_mut().get_mut(..max_dim, ..max_dim),
							Some(Q.as_mut()),
							hi,
							lo,
						);
						for k in 0..bs {
							residual[lo + k] = residual[hi + k].copy();
						}
						lo += bs;
						mi += bs;
						hi += bs;
					},
					Group::Retain => {
						reorder_schur(
							H.rb_mut().get_mut(..max_dim, ..max_dim),
							Some(Q.as_mut()),
							hi,
							mi,
						);
						mi += bs;
						hi += bs;
					},
					Group::Purge => {
						hi += bs;
					},
				}
			}
			let mut V_tmp = tmp.rb_mut().get_mut(.., purge..k);
			matmul(
				V_tmp.rb_mut(),
				Accum::Replace,
				V.rb().get(.., purge..max_dim),
				Q.rb().get(purge..max_dim, purge..k),
				one(),
				par,
			);
			V.rb_mut().get_mut(.., purge..k).copy_from(&V_tmp);
			let mut b_tmp = tmp.rb_mut().get_mut(0, ..);
			matmul(
				b_tmp.rb_mut(),
				Accum::Replace,
				H.rb().get(max_dim, ..),
				Q.rb(),
				one(),
				par,
			);
			H.rb_mut().get_mut(max_dim, ..).copy_from(b_tmp);
			let (mut x, y) = V.rb_mut().two_cols_mut(k, max_dim);
			x.copy_from(&y);
			let (mut x, mut y) = H.rb_mut().two_rows_mut(k, max_dim);
			x.copy_from(&y);
			y.fill(zero());
			active = nlock;
			if nlock >= n_eigval {
				break;
			}
		}
	} else {
		let mut H = H.rb_mut().get_mut(..n, ..n);
		let mut V = V.rb_mut().get_mut(..n, ..n);
		V.rb_mut().diagonal_mut().fill(one());
		A.apply(H.rb_mut(), V.rb(), par, stack);
		let (mut w, stack) = temp_mat_zeroed::<T, _, _>(n, 2, stack);
		let mut w = w.as_mat_mut();
		let (mut w_re, mut w_im) =
			w.rb_mut().get_mut(active..max_dim, ..).two_cols_mut(0, 1);
		{
			let block_size =
				linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
					n, n,
				);
			let (mut householder, stack) =
				temp_mat_zeroed::<T, _, _>(block_size, n - 1, stack);
			let mut householder = householder.as_mat_mut();
			linalg::evd::hessenberg::hessenberg_in_place(
				H.rb_mut(),
				householder.rb_mut(),
				par,
				stack,
				default(),
			);
			linalg::householder::apply_block_householder_sequence_on_the_right_in_place_with_conj(
				H.rb().submatrix(1, 0, n - 1, n - 1),
				householder.rb(),
				Conj::No,
				V.rb_mut().submatrix_mut(1, 1, n - 1, n - 1),
				par,
				stack,
			);
			for j in 0..n {
				for i in j + 2..n {
					H[(i, j)] = zero();
				}
			}
		}
		schur::real_schur::multishift_qr(
			true,
			H.rb_mut(),
			Some(V.rb_mut()),
			w_re.rb_mut(),
			w_im.rb_mut(),
			0,
			n,
			par,
			stack,
			auto!(T),
		);
		active = n;
	}
	let n = active;
	let H = H.rb().get(..n, ..n);
	let V = V.rb().get(.., ..n);
	let mut X = X.rb_mut().get_mut(..n, ..n);
	linalg::evd::evd_from_real_schur_imp(H, X.as_mut(), par, auto!(T));
	let mut vecs = tmp.rb_mut().get_mut(.., ..n);
	matmul(vecs.rb_mut(), Accum::Replace, V, X.rb(), one(), par);
	let V = vecs.rb();
	let (mut norms, stack) = stack.make_with(n, |_| zero::<T::Real>());
	let (mut perm, stack) = stack.make_with(n, |_| 0usize);
	let _ = stack;
	let perm = &mut *perm;
	let norms = &mut *norms;
	let mut j = 0usize;
	while j < n {
		let cplx = j + 1 < n && H[(j + 1, j)] != zero::<T>();
		let bs = if cplx { 2 } else { 1 };
		let re = &H[(j, j)];
		if cplx {
			let im = H[(j + 1, j)].abs().sqrt() * H[(j, j + 1)].abs().sqrt();
			norms[j] = re.hypot(&im);
			norms[j + 1] = norms[j].copy();
			perm[j] = j;
			perm[j + 1] = j;
		} else {
			norms[j] = re.abs();
			perm[j] = j;
		}
		j += bs;
	}
	perm.sort_unstable_by(|&i, &j| {
		if norms[i] > norms[j] {
			core::cmp::Ordering::Less
		} else if norms[i] < norms[j] {
			core::cmp::Ordering::Greater
		} else {
			core::cmp::Ordering::Equal
		}
	});
	let mut eigvecs = eigvecs;
	let limit = Ord::min(n, n_eigval);
	let mut idx = 0usize;
	while idx < limit {
		let j = perm[idx];
		let cplx = j + 1 < n && H[(j + 1, j)] != zero::<T>();
		let bs = if cplx { 2 } else { 1 };
		let re = &H[(j, j)];
		let v_re = V.col(j);
		if cplx {
			let v_im = V.col(j + 1);
			let im = H[(j + 1, j)].abs().sqrt() * H[(j, j + 1)].abs().sqrt();
			if idx + 1 < limit {
				eigvals[idx + 1] = Complex::new(re.copy(), -&im);
			}
			eigvals[idx] = Complex::new(re.copy(), im);
			if idx + 1 < limit {
				let (ej, ej1) = eigvecs.rb_mut().two_cols_mut(idx, idx + 1);
				zip!(ej, ej1, v_re, v_im).for_each(
					|unzip!(y0, y1, re, im): Zip!(
						&mut Complex<T>,
						&mut Complex<T>,
						&T,
						&T
					)| {
						*y0 = Complex::new(re.copy(), im.copy());
						*y1 = Complex::new(re.copy(), -im);
					},
				);
			} else {
				let ej = eigvecs.rb_mut().col_mut(idx);
				zip!(ej, v_re, v_im).for_each(
					|unzip!(y0, re, im): Zip!(&mut Complex<T>, &T, &T)| {
						*y0 = Complex::new(re.copy(), im.copy());
					},
				);
			}
		} else {
			eigvals[idx] = Complex::new(re.copy(), zero());
			zip!(eigvecs.rb_mut().col_mut(idx), v_re).for_each(
				|unzip!(y, x): Zip!(&mut Complex<T>, &T)| {
					*y = Complex::new(x.copy(), zero());
				},
			);
		}
		idx += bs;
	}
	limit
}

fn partial_schur_cplx_imp<T: ComplexField>(
	eigvecs: MatMut<'_, Complex<T::Real>>,
	eigvals: &mut [Complex<T::Real>],
	A: &dyn LinOp<T>,
	v0: ColRef<'_, T>,
	min_dim: usize,
	max_dim: usize,
	n_eigval: usize,
	tol: T::Real,
	restarts: usize,
	par: Par,
	stack: &mut MemStack,
) -> usize {
	let n = A.nrows();
	let (mut H, stack) =
		temp_mat_zeroed::<T, _, _>(max_dim + 1, max_dim, stack);
	let mut H = H.as_mat_mut();
	let (mut V, stack) = temp_mat_zeroed::<T, _, _>(n, max_dim + 1, stack);
	let mut V = V.as_mat_mut();
	let (mut X, stack) = temp_mat_zeroed::<T, _, _>(max_dim, max_dim, stack);
	let mut X = X.as_mat_mut();
	let (mut vecs, stack) = temp_mat_zeroed::<T, _, _>(max_dim, max_dim, stack);
	let mut vecs = vecs.as_mat_mut();
	let (mut tmp, stack) = temp_mat_zeroed::<T, _, _>(n, max_dim, stack);
	let mut tmp = tmp.as_mat_mut();
	let mut active = 0usize;
	if max_dim < n {
		let (mut w, stack) = temp_mat_zeroed::<T, _, _>(max_dim, 1, stack);
		let mut w = w.as_mat_mut().col_mut(0);
		let (mut residual, stack) =
			temp_mat_zeroed::<T::Real, _, _>(max_dim, 1, stack);
		let mut residual = residual.as_mat_mut().col_mut(0);
		let (mut Q, stack) =
			temp_mat_zeroed::<T, _, _>(max_dim, max_dim, stack);
		let mut Q = Q.as_mat_mut();
		let block_size =
			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
				max_dim, max_dim,
			);
		let (mut householder, stack) =
			temp_mat_zeroed::<T, _, _>(block_size, max_dim, stack);
		let mut householder = householder.as_mat_mut();
		let f = v0.norm_l2();
		if f > min_positive() {
			let f = f.recip();
			zip!(V.rb_mut().col_mut(0), v0)
				.for_each(|unzip!(y, x): Zip!(&mut T, &T)| *y = x.mul_real(&f));
		} else {
			let n0 = n as u32;
			let n1 = (n >> 32) as u32;
			let n = from_f64::<T>(n0 as f64) + from_f64::<T>(n1 as f64);
			let f = n.sqrt().recip();
			zip!(V.rb_mut().col_mut(0)).for_each(|unzip!(y)| *y = f.copy());
		}
		iterate_arnoldi(A, H.as_mut(), V.as_mut(), 1, min_dim, par, stack);
		let mut k = min_dim;
		for iter in 0..restarts {
			_ = iter;
			iterate_arnoldi(
				A,
				H.as_mut(),
				V.as_mut(),
				k + 1,
				max_dim,
				par,
				stack,
			);
			let ref Hmm = H[(max_dim, max_dim - 1)].copy();
			let n = max_dim - active;
			let mut w = w.rb_mut().get_mut(active..max_dim);
			Q.fill(zero());
			Q.rb_mut().diagonal_mut().fill(one());
			let mut Q_slice =
				Q.rb_mut().get_mut(active..max_dim, active..max_dim);
			let mut H_slice =
				H.rb_mut().get_mut(active..max_dim, active..max_dim);
			{
				let n = max_dim - active;
				let block_size =
					linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
						n, n,
					);
				let mut householder =
					householder.rb_mut().get_mut(..block_size, ..n - 1);
				linalg::evd::hessenberg::hessenberg_in_place(
					H_slice.rb_mut(),
					householder.rb_mut(),
					par,
					stack,
					default(),
				);
				linalg::householder::apply_block_householder_sequence_on_the_right_in_place_with_conj(
					H_slice.rb().submatrix(1, 0, n - 1, n - 1),
					householder.rb(),
					Conj::No,
					Q_slice.rb_mut().submatrix_mut(1, 1, n - 1, n - 1),
					par,
					stack,
				);
				for j in 0..n {
					for i in j + 2..n {
						H_slice[(i, j)] = zero();
					}
				}
			}
			schur::complex_schur::multishift_qr(
				true,
				H_slice.rb_mut(),
				Some(Q_slice.rb_mut()),
				w.rb_mut(),
				0,
				n,
				par,
				stack,
				auto!(T),
			);
			for j in 0..n {
				let mut idx = j;
				let mut max = zero::<T::Real>();
				for i in j..n {
					let v = w[i].abs();
					if v > max {
						max = v;
						idx = i;
					}
				}
				let i = idx;
				if i != j {
					reorder_schur(
						H_slice.rb_mut(),
						Some(Q_slice.rb_mut()),
						i,
						j,
					);
					w.rb_mut().try_as_col_major_mut().unwrap().as_slice_mut()
						[j..i + 1]
						.rotate_right(1);
				}
			}
			let mut X = X.rb_mut().get_mut(..n, ..n);
			linalg::evd::evd_from_cplx_schur_imp(
				H_slice.rb(),
				Conj::No,
				X.as_mut(),
				par,
				auto!(T),
			);
			let mut vecs = vecs.rb_mut().get_mut(..n, ..n);
			matmul(
				vecs.rb_mut(),
				Accum::Replace,
				Q_slice.rb(),
				X.rb(),
				one(),
				par,
			);
			let vecs = vecs.rb();
			let mut H_tmp = tmp.rb_mut().get_mut(..active, ..);
			matmul(
				H_tmp.rb_mut(),
				Accum::Replace,
				H.rb().get(..active, ..),
				Q.rb(),
				one(),
				par,
			);
			H.rb_mut().get_mut(..active, ..).copy_from(&H_tmp);
			let ref Hmm_abs = Hmm.abs();
			for j in 0..n {
				residual[active + j] =
					Hmm_abs * vecs[(max_dim - active - 1, j)].abs();
			}
			#[derive(Copy, Clone, PartialEq, Eq, Debug)]
			enum Group {
				Lock,
				Retain,
				Purge,
			}
			let mut groups = alloc::vec![Group::Purge; max_dim];
			let nev = n_eigval;
			let mut nlock = 0usize;
			for j in 0..nev {
				if residual[j] <= tol {
					groups[j] = Group::Lock;
					nlock += 1;
				} else {
					groups[j] = Group::Retain;
				}
			}
			let ideal_size = Ord::min(nlock + min_dim, (min_dim + max_dim) / 2);
			k = nev;
			for i in nev..max_dim {
				let group;
				if k < ideal_size && residual[i] > tol {
					group = Group::Retain;
					k += 1;
				} else {
					group = Group::Purge;
				}
				groups[i] = group;
			}
			let mut purge = 0usize;
			while purge < active && groups[purge] == Group::Lock {
				purge += 1;
			}
			let mut lo = 0usize;
			let mut mi = 0usize;
			let mut hi = 0usize;
			while hi < max_dim {
				match groups[hi] {
					Group::Lock => {
						reorder_schur(
							H.rb_mut().get_mut(..max_dim, ..max_dim),
							Some(Q.as_mut()),
							hi,
							lo,
						);
						residual[lo] = residual[hi].copy();
						lo += 1;
						mi += 1;
						hi += 1;
					},
					Group::Retain => {
						reorder_schur(
							H.rb_mut().get_mut(..max_dim, ..max_dim),
							Some(Q.as_mut()),
							hi,
							mi,
						);
						mi += 1;
						hi += 1;
					},
					Group::Purge => {
						hi += 1;
					},
				}
			}
			let mut V_tmp = tmp.rb_mut().get_mut(.., purge..k);
			matmul(
				V_tmp.rb_mut(),
				Accum::Replace,
				V.rb().get(.., purge..max_dim),
				Q.rb().get(purge..max_dim, purge..k),
				one(),
				par,
			);
			V.rb_mut().get_mut(.., purge..k).copy_from(&V_tmp);
			let mut b_tmp = tmp.rb_mut().get_mut(0, ..);
			matmul(
				b_tmp.rb_mut(),
				Accum::Replace,
				H.rb().get(max_dim, ..),
				Q.rb(),
				one(),
				par,
			);
			H.rb_mut().get_mut(max_dim, ..).copy_from(b_tmp);
			let (mut x, y) = V.rb_mut().two_cols_mut(k, max_dim);
			x.copy_from(&y);
			let (mut x, mut y) = H.rb_mut().two_rows_mut(k, max_dim);
			x.copy_from(&y);
			y.fill(zero());
			active = nlock;
			if nlock >= n_eigval {
				break;
			}
		}
	} else {
		let mut H = H.rb_mut().get_mut(..n, ..n);
		let mut V = V.rb_mut().get_mut(..n, ..n);
		V.rb_mut().diagonal_mut().fill(one());
		A.apply(H.rb_mut(), V.rb(), par, stack);
		let (mut w, stack) = temp_mat_zeroed::<T, _, _>(n, 1, stack);
		let mut w = w.as_mat_mut();
		{
			let block_size =
				linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
					n, n,
				);
			let (mut householder, stack) =
				temp_mat_zeroed::<T, _, _>(block_size, n - 1, stack);
			let mut householder = householder.as_mat_mut();
			linalg::evd::hessenberg::hessenberg_in_place(
				H.rb_mut(),
				householder.rb_mut(),
				par,
				stack,
				default(),
			);
			linalg::householder::apply_block_householder_sequence_on_the_right_in_place_with_conj(
				H.rb().submatrix(1, 0, n - 1, n - 1),
				householder.rb(),
				Conj::No,
				V.rb_mut().submatrix_mut(1, 1, n - 1, n - 1),
				par,
				stack,
			);
			for j in 0..n {
				for i in j + 2..n {
					H[(i, j)] = zero();
				}
			}
		}
		schur::complex_schur::multishift_qr(
			true,
			H.rb_mut(),
			Some(V.rb_mut()),
			w.rb_mut().col_mut(0),
			0,
			n,
			par,
			stack,
			auto!(T),
		);
		active = n;
	}
	let n = active;
	let H = H.rb().get(..n, ..n);
	let V = V.rb().get(.., ..n);
	let mut X = X.rb_mut().get_mut(..n, ..n);
	linalg::evd::evd_from_cplx_schur_imp(
		H,
		Conj::No,
		X.as_mut(),
		par,
		auto!(T),
	);
	let mut vecs = tmp.rb_mut().get_mut(.., ..n);
	matmul(vecs.rb_mut(), Accum::Replace, V, X.rb(), one(), par);
	let V = vecs.rb();
	let (mut norms, stack) = stack.make_with(n, |j| H[(j, j)].abs());
	let (mut perm, stack) = stack.make_with(n, |j| j);
	let _ = stack;
	let perm = &mut *perm;
	let norms = &mut *norms;
	perm.sort_unstable_by(|&i, &j| {
		if norms[i] > norms[j] {
			core::cmp::Ordering::Less
		} else if norms[i] < norms[j] {
			core::cmp::Ordering::Greater
		} else {
			core::cmp::Ordering::Equal
		}
	});
	let mut eigvecs = eigvecs;
	let limit = Ord::min(n, n_eigval);
	for idx in 0..limit {
		let j = perm[idx];
		let w = &H[(j, j)];
		let v = V.col(j);
		eigvals[idx] = Complex::new(w.real(), w.imag());
		zip!(eigvecs.rb_mut().col_mut(idx), v).for_each(
			|unzip!(y, x): Zip!(&mut Complex<T::Real>, &T)| {
				*y = Complex::new(x.real(), x.imag());
			},
		);
	}
	limit
}

/// computes the layout of required workspace for computing the `n_eigval`
/// eigenvalues (and corresponding eigenvectors) of $A$ with the largest
/// magnitude.
pub fn partial_eigen_scratch<T: ComplexField>(
	A: &dyn LinOp<T>,
	n_eigval: usize,
	par: Par,
	params: PartialEigenParams,
) -> StackReq {
	let n = A.nrows();
	assert!(A.ncols() == n);
	if n == 0 {
		return StackReq::EMPTY;
	}
	let n_eigval = Ord::min(n_eigval, n);
	let max_dim = Ord::min(
		Ord::max(params.max_dim, Ord::max(2 * MIN_DIM, 2 * n_eigval)),
		n,
	);
	let w = temp_mat_scratch::<T>(max_dim, if T::IS_REAL { 2 } else { 1 });
	let residual = temp_mat_scratch::<T::Real>(max_dim, 1);
	let H = temp_mat_scratch::<T>(max_dim + 1, max_dim);
	let V = temp_mat_scratch::<T>(n, max_dim + 1);
	let Q = temp_mat_scratch::<T>(max_dim, max_dim);
	let X = Q;
	let vecs = Q;
	let tmp = temp_mat_scratch::<T>(n, max_dim);
	let block_size = linalg::qr::no_pivoting::factor::recommended_block_size::<T>(
		max_dim, max_dim,
	);
	let householder = temp_mat_scratch::<T>(block_size, max_dim);
	let arnoldi = A.apply_scratch(1, par).or(StackReq::new::<bool>(max_dim));
	let hess = linalg::evd::hessenberg::hessenberg_in_place_scratch::<T>(
		max_dim,
		block_size,
		par,
		default(),
	);
	let apply_house = linalg::householder::apply_block_householder_sequence_on_the_right_in_place_scratch::<T>(max_dim - 1, block_size, max_dim - 1);
	let schur = schur::multishift_qr_scratch::<T>(
		max_dim,
		max_dim,
		true,
		true,
		par,
		auto!(T),
	);
	StackReq::all_of(&[
		w,
		residual,
		H,
		V,
		Q,
		X,
		vecs,
		tmp,
		householder,
		StackReq::any_of(&[hess, apply_house, schur, arnoldi]),
	])
}

/// computes the layout of required workspace for computing the `n_eigval`
/// eigenvalues (and corresponding eigenvectors) of $A$ with the largest
/// magnitude.
pub fn partial_self_adjoint_eigen_scratch<T: ComplexField>(
	A: &dyn LinOp<T>,
	n_eigval: usize,
	par: Par,
	params: PartialEigenParams,
) -> StackReq {
	let n = A.nrows();
	assert!(A.ncols() == n);
	if n == 0 {
		return StackReq::EMPTY;
	}

	let n_eigval = Ord::min(n_eigval, n);
	let max_dim = Ord::min(
		Ord::max(params.max_dim, Ord::max(2 * MIN_DIM, 2 * n_eigval)),
		n,
	);
	let residual = temp_mat_scratch::<T::Real>(max_dim, 1);
	let H = temp_mat_scratch::<T>(max_dim + 1, max_dim);
	let V = temp_mat_scratch::<T>(n, max_dim + 1);
	let tmp = temp_mat_scratch::<T>(n, max_dim);
	let Q = temp_mat_scratch::<T>(max_dim, max_dim);

	let X = Q;
	let s = temp_mat_scratch::<T>(max_dim, 1);
	let evd = linalg::evd::evd_scratch::<T>(
		max_dim,
		linalg::evd::ComputeEigenvectors::Yes,
		linalg::evd::ComputeEigenvectors::Yes,
		par,
		default(),
	);
	let norms = StackReq::new::<T::Real>(max_dim);
	let perm = StackReq::new::<usize>(max_dim);

	let arnoldi = StackReq::any_of(&[
		A.apply_scratch(1, par),
		StackReq::new::<bool>(max_dim),
	]);

	StackReq::all_of(&[
		H,
		V,
		Q,
		tmp,
		X,
		residual,
		StackReq::any_of(&[
			arnoldi,
			StackReq::all_of(&[s, evd]),
			StackReq::all_of(&[norms, perm]),
		]),
	])
}

/// computes the layout of required workspace for computing the `n_sval`
/// singular values (and corresponding singular vectors) of $A$ with the largest
/// magnitude.
pub fn partial_svd_scratch<T: ComplexField>(
	A: &dyn BiLinOp<T>,
	n_sval: usize,
	par: Par,
	params: PartialEigenParams,
) -> StackReq {
	let m = A.nrows();
	let n = A.ncols();

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

	if k == 0 {
		return StackReq::EMPTY;
	}

	let n_sval = Ord::min(n_sval, k);
	let max_dim = Ord::min(
		Ord::max(params.max_dim, Ord::max(2 * MIN_DIM, 2 * n_sval)),
		k,
	);

	let H = temp_mat_scratch::<T>(max_dim, max_dim + 1);
	let P = temp_mat_scratch::<T>(m, max_dim);
	let Q = temp_mat_scratch::<T>(n, max_dim + 1);
	let tmp = temp_mat_scratch::<T>(Ord::max(m, n), max_dim);

	let X = temp_mat_scratch::<T>(max_dim, max_dim);
	let Y = X;
	let residual = temp_mat_scratch::<T>(max_dim, 1);

	let s = temp_mat_scratch::<T>(max_dim, 1);
	let svd = linalg::svd::svd_scratch::<T>(
		max_dim,
		max_dim,
		linalg::svd::ComputeSvdVectors::Full,
		linalg::svd::ComputeSvdVectors::Full,
		par,
		default(),
	);
	let norms = StackReq::new::<T::Real>(max_dim);
	let perm = StackReq::new::<usize>(max_dim);

	let arnoldi = StackReq::any_of(&[
		A.apply_scratch(1, par),
		A.transpose_apply_scratch(1, par),
		StackReq::new::<bool>(max_dim),
	]);

	StackReq::all_of(&[
		H,
		P,
		Q,
		tmp,
		X,
		Y,
		residual,
		StackReq::any_of(&[
			arnoldi,
			StackReq::all_of(&[s, svd]),
			StackReq::all_of(&[norms, perm]),
		]),
	])
}

/// computes an estimate of the eigenvalues (and corresponding eigenvectors) of
/// $A$ with the largest magnitude until the provided outputs are full or the
/// maximum number of algorithm restarts is reached.
pub fn partial_eigen<T: ComplexField>(
	eigvecs: MatMut<'_, Complex<T::Real>>,
	eigvals: &mut [Complex<T::Real>],
	A: &dyn LinOp<T>,
	v0: ColRef<'_, T>,
	tolerance: T::Real,
	par: Par,
	stack: &mut MemStack,
	params: PartialEigenParams,
) -> PartialEigenInfo {
	let n = v0.nrows();
	assert!(all(
		eigvals.len() == eigvecs.ncols(),
		A.nrows() == n,
		A.ncols() == n,
		eigvecs.nrows() == n,
	));
	let n_eigval = eigvals.len();
	let n_eigval = Ord::min(n_eigval, n);
	if n == 0 {
		return PartialEigenInfo {
			n_converged_eigen: 0,
			non_exhaustive: NonExhaustive(()),
		};
	}
	let min_dim =
		Ord::min(Ord::max(params.min_dim, Ord::max(MIN_DIM, n_eigval)), n);
	let max_dim = Ord::min(
		Ord::max(params.max_dim, Ord::max(2 * MIN_DIM, 2 * n_eigval)),
		n,
	);
	let n_eigval = if const { T::IS_REAL } {
		partial_schur_real_imp::<T::Real>(
			eigvecs,
			eigvals,
			unsafe { core::mem::transmute(A) },
			unsafe { core::mem::transmute(v0) },
			min_dim,
			max_dim,
			n_eigval,
			tolerance,
			params.max_restarts,
			par,
			stack,
		)
	} else {
		partial_schur_cplx_imp(
			eigvecs,
			eigvals,
			A,
			v0,
			min_dim,
			max_dim,
			n_eigval,
			tolerance,
			params.max_restarts,
			par,
			stack,
		)
	};
	PartialEigenInfo {
		n_converged_eigen: n_eigval,
		non_exhaustive: NonExhaustive(()),
	}
}

/// computes an estimate of the eigenvalues (and corresponding eigenvectors) of
/// $A$ with the largest magnitude, assuming $A$ is self-adjoint, until the
/// provided outputs are full or the maximum number of algorithm restarts is
/// reached.
pub fn partial_self_adjoint_eigen<T: ComplexField>(
	eigvecs: MatMut<'_, T>,
	eigvals: &mut [T],
	A: &dyn LinOp<T>,
	v0: ColRef<'_, T>,
	tolerance: T::Real,
	par: Par,
	stack: &mut MemStack,
	params: PartialEigenParams,
) -> PartialEigenInfo {
	let n = v0.nrows();
	assert!(all(
		eigvals.len() == eigvecs.ncols(),
		A.nrows() == n,
		A.ncols() == n,
		eigvecs.nrows() == n,
	));
	let n_eigval = eigvals.len();
	let n_eigval = Ord::min(n_eigval, n);
	if n == 0 {
		return PartialEigenInfo {
			n_converged_eigen: 0,
			non_exhaustive: NonExhaustive(()),
		};
	}
	let min_dim =
		Ord::min(Ord::max(params.min_dim, Ord::max(MIN_DIM, n_eigval)), n);
	let max_dim = Ord::min(
		Ord::max(params.max_dim, Ord::max(2 * MIN_DIM, 2 * n_eigval)),
		n,
	);
	let n_eigval = {
		super::self_adjoint_eigen::partial_self_adjoint_eigen_imp(
			eigvecs,
			eigvals,
			A,
			v0,
			min_dim,
			max_dim,
			n_eigval,
			tolerance,
			params.max_restarts,
			par,
			stack,
		)
	};
	PartialEigenInfo {
		n_converged_eigen: n_eigval,
		non_exhaustive: NonExhaustive(()),
	}
}

/// computes an estimate of the singular values (and corresponding singular
/// vectors) of $A$ with the largest magnitude, until the provided outputs are
/// full or the maximum number of algorithm restarts is reached.
pub fn partial_svd<T: ComplexField>(
	left_singular_vecs: MatMut<'_, T>,
	right_singular_vecs: MatMut<'_, T>,
	singular_vals: &mut [T],
	A: &dyn BiLinOp<T>,
	v0: ColRef<'_, T>,
	tolerance: T::Real,
	par: Par,
	stack: &mut MemStack,
	params: PartialEigenParams,
) -> PartialEigenInfo {
	let m = A.nrows();
	let n = A.ncols();
	let k = Ord::min(m, n);

	assert!(all(
		singular_vals.len() == left_singular_vecs.ncols(),
		singular_vals.len() == right_singular_vecs.ncols(),
		v0.nrows() == n,
		left_singular_vecs.nrows() == m,
		right_singular_vecs.nrows() == n,
	));
	let n_sval = singular_vals.len();
	let n_sval = Ord::min(n_sval, k);

	if k == 0 {
		return PartialEigenInfo {
			n_converged_eigen: 0,
			non_exhaustive: NonExhaustive(()),
		};
	}

	let min_dim =
		Ord::min(Ord::max(params.min_dim, Ord::max(MIN_DIM, n_sval)), k);
	let max_dim = Ord::min(
		Ord::max(params.max_dim, Ord::max(2 * MIN_DIM, 2 * n_sval)),
		k,
	);
	let n_sval = super::svd::partial_svd_imp(
		left_singular_vecs,
		right_singular_vecs,
		singular_vals,
		A,
		v0,
		min_dim,
		max_dim,
		n_sval,
		tolerance,
		params.max_restarts,
		par,
		stack,
	);

	PartialEigenInfo {
		n_converged_eigen: n_sval,
		non_exhaustive: NonExhaustive(()),
	}
}

#[cfg(test)]
mod tests {
	use super::*;
	use crate::stats::prelude::*;
	use crate::{Scale, assert};
	use rand::prelude::*;

	#[test]
	fn test_arnoldi_real() {
		let rng = &mut StdRng::seed_from_u64(1);
		let n = 100;
		let n_eigval = 20;
		let min_dim = 30;
		let max_dim = 60;
		let restarts = 1000;
		let mat = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: StandardNormal,
		};
		let col = CwiseColDistribution {
			nrows: n,
			dist: StandardNormal,
		};
		let A: Mat<f64> = mat.sample(rng);
		let mut v0: Col<f64> = col.sample(rng);
		v0 /= v0.norm_l2();
		let A = A.as_ref();
		let v0 = v0.as_ref();
		let par = Par::Seq;
		let mem = &mut MemBuffer::new(partial_eigen_scratch(
			&A,
			n_eigval,
			par,
			PartialEigenParams {
				min_dim,
				max_dim,
				max_restarts: restarts,
				..Default::default()
			},
		));
		let mut V = Mat::zeros(n, n_eigval);
		let mut w = vec![c64::ZERO; n_eigval];
		let info = partial_eigen(
			V.rb_mut(),
			&mut w,
			&A,
			v0,
			f64::EPSILON * 128.0,
			par,
			MemStack::new(mem),
			PartialEigenParams {
				min_dim,
				max_dim,
				max_restarts: restarts,
				..Default::default()
			},
		);
		assert!(w.iter().map(|x| x.norm()).is_sorted_by(|x, y| x >= y));
		let A = &zip!(A).map(|unzip!(x)| Complex::from(*x));
		for j in 0..info.n_converged_eigen {
			assert!((A * V.col(j) - Scale(w[j]) * V.col(j)).norm_l2() < 1e-10);
		}
	}
	#[test]
	fn test_arnoldi_cplx() {
		let rng = &mut StdRng::seed_from_u64(1);
		let n = 100;
		let n_eigval = 20;
		let min_dim = 30;
		let max_dim = 60;
		let restarts = 1000;
		let mat = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		};
		let col = CwiseColDistribution {
			nrows: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		};
		let A: Mat<c64> = mat.sample(rng);
		let mut v0: Col<c64> = col.sample(rng);
		v0 /= v0.norm_l2();
		let A = A.as_ref();
		let v0 = v0.as_ref();
		let par = Par::Seq;
		let mem = &mut MemBuffer::new(partial_eigen_scratch(
			&A,
			n_eigval,
			par,
			PartialEigenParams {
				min_dim,
				max_dim,
				max_restarts: restarts,
				..Default::default()
			},
		));
		let mut V = Mat::zeros(n, n_eigval);
		let mut w = vec![c64::ZERO; n_eigval];
		let info = partial_eigen(
			V.rb_mut(),
			&mut w,
			&A,
			v0,
			f64::EPSILON * 128.0,
			par,
			MemStack::new(mem),
			PartialEigenParams {
				min_dim,
				max_dim,
				max_restarts: restarts,
				..Default::default()
			},
		);
		assert!(w.iter().map(|x| x.norm()).is_sorted_by(|x, y| x >= y));
		for j in 0..info.n_converged_eigen {
			assert!((A * V.col(j) - Scale(w[j]) * V.col(j)).norm_l2() < 1e-10);
		}
	}
	#[test]
	fn test_toeplitz() {
		let n = 10;
		let mut mat = Mat::<f64>::zeros(n, n);
		for i in 0..n {
			mat[(i, i)] = 6.0;
			if i + 1 < n {
				mat[(i, i + 1)] = 3.0;
			}
			if i > 0 {
				mat[(i, i - 1)] = 3.0;
			}
		}
		let rng = &mut StdRng::seed_from_u64(1);
		let col = CwiseColDistribution {
			nrows: n,
			dist: StandardNormal,
		};
		let mut v0: Col<f64> = col.sample(rng);
		v0 /= v0.norm_l2();
		let A = mat.as_ref();
		let v0 = v0.as_ref();
		let par = Par::Seq;
		let n_eigval = 5;
		let min_dim = 7;
		let max_dim = 9;
		let mem = &mut MemBuffer::new(partial_eigen_scratch(
			&A,
			n_eigval,
			par,
			PartialEigenParams {
				min_dim,
				max_dim,
				max_restarts: 20,
				..Default::default()
			},
		));
		let mut V = Mat::zeros(n, n_eigval);
		let mut w = vec![c64::ZERO; n_eigval];
		let info = partial_eigen(
			V.rb_mut(),
			&mut w,
			&A,
			v0,
			0.000001,
			par,
			MemStack::new(mem),
			PartialEigenParams {
				min_dim,
				max_dim,
				max_restarts: 20,
				..Default::default()
			},
		);
		assert!(w.iter().map(|x| x.norm()).is_sorted_by(|x, y| x >= y));
		let A = &zip!(A).map(|unzip!(x)| Complex::from(*x));
		for j in 0..info.n_converged_eigen {
			assert!((A * V.col(j) - Scale(w[j]) * V.col(j)).norm_l2() < 1e-10);
		}
	}
	#[test]
	fn test_small_real() {
		let rng = &mut StdRng::seed_from_u64(1);
		let n = 5;
		let n_eigval = 3;
		let mat = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: StandardNormal,
		};
		let col = CwiseColDistribution {
			nrows: n,
			dist: StandardNormal,
		};
		let A: Mat<f64> = mat.sample(rng);
		let mut v0: Col<f64> = col.sample(rng);
		v0 /= v0.norm_l2();
		let A = A.as_ref();
		let v0 = v0.as_ref();
		let par = Par::Seq;
		let mem = &mut MemBuffer::new(partial_eigen_scratch(
			&A,
			n_eigval,
			par,
			default(),
		));
		let mut V = Mat::zeros(n, n_eigval);
		let mut w = vec![c64::ZERO; n_eigval];
		let info = partial_eigen(
			V.rb_mut(),
			&mut w,
			&A,
			v0,
			f64::EPSILON * 128.0,
			par,
			MemStack::new(mem),
			default(),
		);
		assert!(w.iter().map(|x| x.norm()).is_sorted_by(|x, y| x >= y));
		let A = &zip!(A).map(|unzip!(x)| Complex::from(*x));
		for j in 0..info.n_converged_eigen {
			assert!((A * V.col(j) - Scale(w[j]) * V.col(j)).norm_l2() < 1e-10);
		}
	}
	#[test]
	fn test_small_cplx() {
		let rng = &mut StdRng::seed_from_u64(1);
		let n = 5;
		let n_eigval = 3;
		let mat = CwiseMatDistribution {
			nrows: n,
			ncols: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		};
		let col = CwiseColDistribution {
			nrows: n,
			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
		};
		let A: Mat<c64> = mat.sample(rng);
		let mut v0: Col<c64> = col.sample(rng);
		v0 /= v0.norm_l2();
		let A = A.as_ref();
		let v0 = v0.as_ref();
		let par = Par::Seq;
		let mem = &mut MemBuffer::new(partial_eigen_scratch(
			&A,
			n_eigval,
			par,
			default(),
		));
		let mut V = Mat::zeros(n, n_eigval);
		let mut w = vec![c64::ZERO; n_eigval];
		let info = partial_eigen(
			V.rb_mut(),
			&mut w,
			&A,
			v0,
			f64::EPSILON * 128.0,
			par,
			MemStack::new(mem),
			default(),
		);
		assert!(w.iter().map(|x| x.norm()).is_sorted_by(|x, y| x >= y));
		let A = &zip!(A).map(|unzip!(x)| Complex::from(*x));
		for j in 0..info.n_converged_eigen {
			assert!((A * V.col(j) - Scale(w[j]) * V.col(j)).norm_l2() < 1e-10);
		}
	}
}