faer/linalg/lu/full_pivoting/
factor.rs

1use crate::internal_prelude::*;
2use crate::utils::cost::{ParThreshold, rate};
3use crate::perm::{swap_cols_idx, swap_rows_idx};
4use crate::utils::thread::par_split_indices;
5use crate::utils::thread::with_width;
6use faer_traits::{Real, RealReg};
7use linalg::matmul::matmul;
8use pulp::Simd;
9
10#[inline(always)]
11fn best_value<T: ComplexField, S: Simd>(
12	simd: &SimdCtx<T, S>,
13	best_value: RealReg<T::SimdVec<S>>,
14	best_indices: T::SimdIndex<S>,
15	value: T::SimdVec<S>,
16	indices: T::SimdIndex<S>,
17) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
18	let value = simd.abs1(value);
19	let is_better = (**simd).gt(value, best_value);
20	(
21		RealReg(simd.select(is_better, value.0, best_value.0)),
22		simd.iselect(is_better, indices, best_indices),
23	)
24}
25#[inline(always)]
26fn best_score<T: ComplexField, S: Simd>(
27	simd: &SimdCtx<T, S>,
28	best_score: RealReg<T::SimdVec<S>>,
29	best_indices: T::SimdIndex<S>,
30	score: RealReg<T::SimdVec<S>>,
31	indices: T::SimdIndex<S>,
32) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
33	let is_better = (**simd).gt(score, best_score);
34	(
35		RealReg(simd.select(is_better, score.0, best_score.0)),
36		simd.iselect(is_better, indices, best_indices),
37	)
38}
39#[inline(always)]
40fn best_score_2d<T: ComplexField, S: Simd>(
41	simd: &SimdCtx<T, S>,
42	best_score: RealReg<T::SimdVec<S>>,
43	best_row: T::SimdIndex<S>,
44	best_col: T::SimdIndex<S>,
45	score: RealReg<T::SimdVec<S>>,
46	row: T::SimdIndex<S>,
47	col: T::SimdIndex<S>,
48) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>, T::SimdIndex<S>) {
49	let is_better = (**simd).gt(score, best_score);
50	(
51		RealReg(simd.select(is_better, score.0, best_score.0)),
52		simd.iselect(is_better, row, best_row),
53		simd.iselect(is_better, col, best_col),
54	)
55}
56#[inline(always)]
57fn reduce_2d<T: ComplexField, S: Simd>(
58	simd: &SimdCtx<T, S>,
59	best_values: RealReg<T::SimdVec<S>>,
60	best_row: T::SimdIndex<S>,
61	best_col: T::SimdIndex<S>,
62) -> (usize, usize, Real<T>) {
63	let best_val = simd.reduce_max_real(best_values);
64	let best_val_splat = simd.splat_real(&best_val);
65	let is_best = (**simd).ge(best_values, best_val_splat);
66	// The lanes that hold the maximum tie with each other (a matrix of small integers ties a
67	// lot), and which lane comes first in the vector depends on the lane count: with 8 lanes an
68	// AVX-512 build picks another row than an AVX2 build for the same matrix, and the elimination
69	// then follows another path. Break the tie the way the scalar search does - the smallest
70	// column, then the smallest row - so the pivots do not depend on the vector width.
71	let sentinel = simd.isplat(T::Index::truncate(usize::MAX));
72	let rows_tied = simd.iselect(is_best, best_row, sentinel);
73	let cols_tied = simd.iselect(is_best, best_col, sentinel);
74	let row_lanes = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
75		core::slice::from_ref(&rows_tied),
76	);
77	let col_lanes = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
78		core::slice::from_ref(&cols_tied),
79	);
80	let mut best: Option<(usize, usize)> = None;
81	for lane in 0..row_lanes.len() {
82		let (row, col) = (row_lanes[lane].zx(), col_lanes[lane].zx());
83		if row == usize::MAX || col == usize::MAX {
84			continue;
85		}
86		if best.is_none_or(|(brow, bcol)| (col, row) < (bcol, brow)) {
87			best = Some((row, col));
88		}
89	}
90	let (best_row, best_col) = best.expect("the maximum is held by some lane");
91	(best_row, best_col, best_val)
92}
93#[inline(always)]
94fn best_in_col_simd<'M, T: ComplexField, S: Simd>(
95	simd: SimdCtx<'M, T, S>,
96	data: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
97) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
98	let iota = T::simd_iota(&simd.0);
99	let lane_count =
100		core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
101	let inc = simd.isplat(T::Index::truncate(lane_count));
102
103	let indices = simd.batch_indices::<4>();
104	let mut best_val = [RealReg(simd.zero()); 4];
105	let mut best_idx = [simd.izero(); 4];
106
107	let mut idx = simd.iadd(
108		iota,
109		simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
110	);
111	simd_iter!(for (IDX, i) in [indices; 4] {
112		(best_val[IDX], best_idx[IDX]) = best_value(
113			&simd,
114			best_val[IDX],
115			best_idx[IDX],
116			simd.read(data, i),
117			idx,
118		);
119		idx = simd.iadd(idx, inc);
120	});
121	let _ = idx;
122
123	let (best_val0, best_idx0) =
124		best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
125	let (best_val2, best_idx2) =
126		best_score(&simd, best_val[2], best_idx[2], best_val[3], best_idx[3]);
127	best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
128}
129#[inline(always)]
130fn update_and_best_in_col_simd<'M, T: ComplexField, S: Simd>(
131	simd: SimdCtx<'M, T, S>,
132	data: ColMut<'_, T, Dim<'M>, ContiguousFwd>,
133	lhs: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
134	rhs: T,
135) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
136	let mut data = data;
137	let iota = T::simd_iota(&simd.0);
138	let lane_count =
139		core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
140	let inc = simd.isplat(T::Index::truncate(lane_count));
141	let mut best_val = [RealReg(simd.zero()); 3];
142	let mut best_idx = [simd.izero(); 3];
143
144	let mut idx = simd.iadd(
145		iota,
146		simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
147	);
148
149	let indices = simd.batch_indices::<3>();
150	let rhs = simd.splat(&-rhs);
151	simd_iter!(for (IDX, i) in [indices; 3] {
152		let mut x0 = simd.read(data.rb(), i);
153		let l0 = simd.read(lhs, i);
154		x0 = simd.mul_add(l0, rhs, x0);
155		(best_val[IDX], best_idx[IDX]) =
156			best_value(&simd, best_val[IDX], best_idx[IDX], x0, idx);
157		simd.write(data.rb_mut(), i, x0);
158		idx = simd.iadd(idx, inc);
159	});
160	let _ = idx;
161
162	let (best_val0, best_idx0) =
163		best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
164	let (best_val2, best_idx2) = (best_val[2], best_idx[2]);
165	best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
166}
167#[inline(always)]
168fn best_in_mat_simd<T: ComplexField>(
169	data: MatRef<'_, T, usize, usize, ContiguousFwd>,
170) -> (usize, usize, Real<T>) {
171	struct Impl<'a, 'M, 'N, T: ComplexField> {
172		data: MatRef<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
173	}
174	impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
175		type Output = (usize, usize, Real<T>);
176
177		#[inline(always)]
178		fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
179			let Self { data } = self;
180			let M = data.nrows();
181			let N = data.ncols();
182			let simd = SimdCtx::<'_, T, S>::new(T::simd_ctx(simd), M);
183			let mut best_row = simd.isplat(T::Index::truncate(0));
184			let mut best_col = simd.isplat(T::Index::truncate(0));
185			let mut best_val = simd.splat_real(&zero());
186			for j in N.indices() {
187				let col = data.col(j);
188				let (best_val_j, best_row_j) = best_in_col_simd(simd, col);
189				(best_val, best_row, best_col) = best_score_2d(
190					&simd,
191					best_val,
192					best_row,
193					best_col,
194					best_val_j,
195					best_row_j,
196					simd.isplat(T::Index::truncate(*j)),
197				);
198			}
199			reduce_2d(&simd, best_val, best_row, best_col)
200		}
201	}
202	with_dim!(M, data.nrows());
203	with_dim!(N, data.ncols());
204	dispatch!(
205		Impl {
206			data: data.as_shape(M, N)
207		},
208		Impl,
209		T
210	)
211}
212#[inline(always)]
213fn update_and_best_in_mat_simd<T: ComplexField>(
214	data: MatMut<'_, T, usize, usize, ContiguousFwd>,
215	lhs: ColRef<'_, T, usize, ContiguousFwd>,
216	rhs: RowRef<'_, T, usize>,
217	align: usize,
218) -> (usize, usize, Real<T>) {
219	struct Impl<'a, 'M, 'N, T: ComplexField> {
220		data: MatMut<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
221		lhs: ColRef<'a, T, Dim<'M>, ContiguousFwd>,
222		rhs: RowRef<'a, T, Dim<'N>>,
223		align: usize,
224	}
225	impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
226		type Output = (usize, usize, Real<T>);
227
228		#[inline(always)]
229		fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
230			let Self {
231				data,
232				lhs,
233				rhs,
234				align,
235			} = self;
236			let M = data.nrows();
237			let N = data.ncols();
238			let simd =
239				SimdCtx::<'_, T, S>::new_align(T::simd_ctx(simd), M, align);
240			let mut best_row = simd.isplat(T::Index::truncate(0));
241			let mut best_col = simd.isplat(T::Index::truncate(0));
242			let mut best_val = simd.splat_real(&zero());
243			let mut data = data;
244			for j in N.indices() {
245				let data = data.rb_mut().col_mut(j);
246				let rhs = rhs[j].copy();
247				let (best_val_j, best_row_j) =
248					update_and_best_in_col_simd(simd, data, lhs, rhs);
249				(best_val, best_row, best_col) = best_score_2d(
250					&simd,
251					best_val,
252					best_row,
253					best_col,
254					best_val_j,
255					best_row_j,
256					simd.isplat(T::Index::truncate(*j)),
257				);
258			}
259			reduce_2d(&simd, best_val, best_row, best_col)
260		}
261	}
262	with_dim!(M, data.nrows());
263	with_dim!(N, data.ncols());
264	dispatch!(
265		Impl {
266			data: data.as_shape_mut(M, N),
267			lhs: lhs.as_row_shape(M),
268			rhs: rhs.as_col_shape(N),
269			align,
270		},
271		Impl,
272		T
273	)
274}
275fn best_in_matrix_fallback<T: ComplexField>(
276	data: MatRef<'_, T>,
277) -> (usize, usize, Real<T>) {
278	let mut max = zero();
279	let mut row = 0;
280	let mut col = 0;
281	let (m, n) = data.shape();
282	for j in 0..n {
283		for i in 0..m {
284			let abs = data[(i, j)].abs1();
285			if abs > max {
286				row = i;
287				col = j;
288				max = abs;
289			}
290		}
291	}
292	(row, col, max)
293}
294fn best_in_matrix<T: ComplexField>(
295	data: MatRef<'_, T>,
296) -> (usize, usize, Real<T>) {
297	if const { T::SIMD_CAPABILITIES.is_simd() } {
298		if let Some(dst) = data.try_as_col_major() {
299			best_in_mat_simd(dst)
300		} else {
301			best_in_matrix_fallback(data)
302		}
303	} else {
304		best_in_matrix_fallback(data)
305	}
306}
307fn rank_one_update_and_best_in_matrix<T: ComplexField>(
308	mut dst: MatMut<'_, T>,
309	lhs: ColRef<'_, T>,
310	rhs: RowRef<'_, T>,
311	align: usize,
312) -> (usize, usize, Real<T>) {
313	if const { T::SIMD_CAPABILITIES.is_simd() } {
314		if let (Some(dst), Some(lhs)) =
315			(dst.rb_mut().try_as_col_major_mut(), lhs.try_as_col_major())
316		{
317			update_and_best_in_mat_simd(dst, lhs, rhs, align)
318		} else {
319			matmul(
320				dst.rb_mut(),
321				Accum::Add,
322				lhs.as_mat(),
323				rhs.as_mat(),
324				-one::<T>(),
325				Par::Seq,
326			);
327			best_in_matrix(dst.rb())
328		}
329	} else {
330		matmul(
331			dst.rb_mut(),
332			Accum::Add,
333			lhs.as_mat(),
334			rhs.as_mat(),
335			-one::<T>(),
336			Par::Seq,
337		);
338		best_in_matrix(dst.rb())
339	}
340}
341fn lu_in_place_unblocked<T: ComplexField>(
342	A: MatMut<'_, T>,
343	row_trans: &mut [usize],
344	col_trans: &mut [usize],
345	best: &mut [(usize, usize, T::Real)],
346	par: Par,
347	transpose: bool,
348	params: Spec<FullPivLuParams, T>,
349) -> usize {
350	// the parallel arm is the only reader (absent without `rayon`)
351	let _ = &best;
352	let params = params.config;
353	let mut n_trans = 0;
354	let (m, n) = A.shape();
355	if m == 0 || n == 0 {
356		return 0;
357	}
358	let mut par = par;
359	// `Never` keeps the whole call sequential (the steps decide their own width otherwise)
360	if matches!(params.par_threshold, ParThreshold::Never) {
361		par = Par::Seq;
362	}
363	let mut A = A;
364	let (mut max_row, mut max_col, mut max_score) = best_in_matrix(A.rb());
365	for k in 0..Ord::min(m, n) {
366		if max_score < min_positive() {
367			for (i, (row, col)) in
368				core::iter::zip(&mut row_trans[k..], &mut col_trans[k..])
369					.enumerate()
370			{
371				*row = i + k;
372				*col = i + k;
373			}
374			break;
375		}
376		row_trans[k] = max_row;
377		col_trans[k] = max_col;
378		if max_row != k {
379			swap_rows_idx(A.rb_mut(), k, max_row);
380			n_trans += 1;
381		}
382		if max_col != k {
383			swap_cols_idx(A.rb_mut(), k, max_col);
384			n_trans += 1;
385		}
386		let ref inv = A[(k, k)].recip();
387		if transpose {
388			for j in k + 1..n {
389				A[(k, j)] *= inv;
390			}
391		} else {
392			for i in k + 1..m {
393				A[(i, k)] *= inv;
394			}
395		}
396		if k + 1 == Ord::min(m, n) {
397			break;
398		}
399		// One section per step (the fused rank-one update and pivot search over column
400		// blocks); the unit is an element of the trailing block.
401		let width = params.par_threshold.region_width::<T>(
402			&par,
403			(m - k - 1) * (n - k - 1),
404			rate::lu_full_pivot(),
405			1,
406		);
407		let (_, A01, A10, mut A11) = A.rb_mut().split_at_mut(k + 1, k + 1);
408		let lhs = A10.col(k);
409		let rhs = A01.row(k);
410		with_width(par.reborrow(), width, |par| {
411			match par {
412				Par::Seq => {
413					(max_row, max_col, max_score) =
414						rank_one_update_and_best_in_matrix(
415							A11.rb_mut(),
416							lhs,
417							rhs,
418							simd_align(k + 1),
419						);
420				},
421		#[cfg(feature = "rayon")]
422				_ => {
423				let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap();
424					use rayon::prelude::*;
425					let nthreads = nthreads.get();
426					// one entry per block, from the call's scratch (no allocation per step, #111); every
427					// entry is written by its block
428					let best = &mut best[..nthreads];
429					let full_cols = A11.ncols();
430					crate::utils::thread::for_each(par, 
431						nthreads,
432						best.par_iter_mut()
433							.zip_eq(A11.rb_mut().par_col_partition_mut(nthreads))
434							.zip_eq(rhs.par_partition(nthreads))
435							.enumerate(),
436						|(idx, (((max_row, max_col, max_score), A11), rhs))| {
437							(*max_row, *max_col, *max_score) = {
438								let (a, mut b, c) =
439									rank_one_update_and_best_in_matrix(
440										A11,
441										lhs,
442										rhs,
443										simd_align(k + 1),
444									);
445								b += par_split_indices(full_cols, idx, nthreads).0;
446								(a, b, c)
447							};
448						},
449					);
450					max_row = 0;
451					max_col = 0;
452					max_score = zero();
453					for (row, col, val) in best.iter() {
454						if *val > max_score {
455							max_row = *row;
456							max_col = *col;
457							max_score = val.copy();
458						}
459					}
460				},
461			}
462		});
463		max_row += k + 1;
464		max_col += k + 1;
465	}
466	n_trans
467}
468/// $LU$ factorization tuning parameters
469#[derive(Copy, Clone, Debug)]
470pub struct FullPivLuParams {
471	/// floor of the parallel regions: [`ParThreshold::Derived`] (the default) runs a step in
472	/// parallel when its serial time pays for one section; [`ParThreshold::Fixed`] counts
473	/// elements of the trailing block `(m-k-1)*(n-k-1)`
474	pub par_threshold: ParThreshold,
475	#[doc(hidden)]
476	pub non_exhaustive: NonExhaustive,
477}
478impl<T: ComplexField> Auto<T> for FullPivLuParams {
479	#[inline]
480	fn auto() -> Self {
481		Self {
482			par_threshold: ParThreshold::Derived,
483			non_exhaustive: NonExhaustive(()),
484		}
485	}
486}
487#[inline]
488pub fn lu_in_place_scratch<I: Index, T: ComplexField>(
489	nrows: usize,
490	ncols: usize,
491	par: ParShape,
492	params: Spec<FullPivLuParams, T>,
493) -> StackReq {
494	_ = params;
495	let size = Ord::min(nrows, ncols);
496	// the transpositions, and one pivot candidate per block of a parallel step
497	StackReq::new::<usize>(size)
498		.array(2)
499		.and(StackReq::new::<(usize, usize, T::Real)>(par.degree()))
500}
501#[derive(Copy, Clone, Debug)]
502pub struct FullPivLuInfo {
503	pub transposition_count: usize,
504}
505pub fn lu_in_place<'out, I: Index, T: ComplexField>(
506	mat: MatMut<'_, T>,
507	row_perm: &'out mut [I],
508	row_perm_inv: &'out mut [I],
509	col_perm: &'out mut [I],
510	col_perm_inv: &'out mut [I],
511	mut par: Par,
512	stack: &mut MemStack,
513	params: Spec<FullPivLuParams, T>,
514) -> (FullPivLuInfo, PermRef<'out, I>, PermRef<'out, I>) {
515	#[cfg(feature = "perf-warn")]
516	if (mat.col_stride().unsigned_abs() == 1
517		|| mat.row_stride().unsigned_abs() != 1)
518		&& crate::__perf_warn!(LU_WARN)
519	{
520		log::warn!(
521			target : "faer_perf",
522			"LU with full pivoting prefers column-major or row-major matrix. Found matrix with generic strides."
523		);
524	}
525	// Same row-major adaptation as the partial-pivoting factorisation above.
526	let (M, N) = mat.shape();
527	if mat.row_stride().unsigned_abs() != 1
528		&& M >= N
529		&& N > 1
530		&& linalg::row_major_copy_pays(M, N)
531	{
532		let mut work = linalg::column_major_copy(mat.rb(), par.reborrow());
533		let mut work = MatMut::from_column_major_slice_mut(&mut work, M, N);
534		let ret = lu_in_place(
535			work.rb_mut(),
536			row_perm,
537			row_perm_inv,
538			col_perm,
539			col_perm_inv,
540			par.reborrow(),
541			stack,
542			params,
543		);
544		linalg::copy_across_layouts(mat, work.rb(), par);
545		return ret;
546	}
547	let size = Ord::min(M, N);
548	let (mut row_transpositions, stack) = stack.make_with(size, |_| 0);
549	let row_transpositions = row_transpositions.as_mut();
550	let (mut col_transpositions, stack) = stack.make_with(size, |_| 0);
551	let col_transpositions = col_transpositions.as_mut();
552	let (mut best, _) = stack.make_with(par.degree(), |_| (0usize, 0usize, zero::<T::Real>()));
553	let best = best.as_mut();
554	let n_transpositions = if mat.row_stride().abs() < mat.col_stride().abs() {
555		lu_in_place_unblocked(
556			mat,
557			row_transpositions,
558			col_transpositions,
559			best,
560			par,
561			false,
562			params,
563		)
564	} else {
565		lu_in_place_unblocked(
566			mat.transpose_mut(),
567			col_transpositions,
568			row_transpositions,
569			best,
570			par,
571			true,
572			params,
573		)
574	};
575	for i in 0..M {
576		row_perm[i] = I::truncate(i);
577	}
578	for (i, t) in row_transpositions.iter().copied().enumerate() {
579		row_perm.as_mut().swap(i, t);
580	}
581	for i in 0..M {
582		row_perm_inv[row_perm[i].zx()] = I::truncate(i);
583	}
584	for j in 0..N {
585		col_perm[j] = I::truncate(j);
586	}
587	for (i, t) in col_transpositions.iter().copied().enumerate() {
588		col_perm.as_mut().swap(i, t);
589	}
590	for j in 0..N {
591		col_perm_inv[col_perm[j].zx()] = I::truncate(j);
592	}
593	unsafe {
594		(
595			FullPivLuInfo {
596				transposition_count: n_transpositions,
597			},
598			PermRef::new_unchecked(row_perm, row_perm_inv, M),
599			PermRef::new_unchecked(col_perm, col_perm_inv, N),
600		)
601	}
602}
603#[cfg(test)]
604mod tests {
605	use super::*;
606	use crate::stats::prelude::*;
607	use crate::utils::approx::*;
608	use crate::{Mat, assert, c64};
609	use dyn_stack::MemBuffer;
610	#[test]
611	fn test_flu() {
612		let rng = &mut StdRng::seed_from_u64(0);
613		let mut run = |mut par: Par<'_>| {
614			let par_shape = par.shape();
615			for m in [8, 16, 24, 32, 128, 255, 256, 257] {
616				let n = 8;
617				let approx_eq = CwiseMat(ApproxEq {
618					abs_tol: 1e-10,
619					rel_tol: 1e-10,
620				});
621				let A = CwiseMatDistribution {
622					nrows: m,
623					ncols: n,
624					dist: ComplexDistribution::new(
625						StandardNormal,
626						StandardNormal,
627					),
628				}
629				.rand::<Mat<c64>>(rng);
630				let A = A.as_ref();
631				let mut LU = A.cloned();
632				let row_perm = &mut *vec![0usize; m];
633				let row_perm_inv = &mut *vec![0usize; m];
634				let col_perm = &mut *vec![0usize; n];
635				let col_perm_inv = &mut *vec![0usize; n];
636				let (_, p, q) = lu_in_place(
637					LU.as_mut(),
638					row_perm,
639					row_perm_inv,
640					col_perm,
641					col_perm_inv,
642					par.reborrow(),
643					MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
644						usize,
645						c64,
646					>(
647						n, n, par_shape, default()
648					))),
649					default(),
650				);
651				let mut L = LU.as_ref().cloned();
652				let mut U = LU.as_ref().cloned();
653				for j in 0..n {
654					for i in 0..j {
655						L[(i, j)] = c64::ZERO;
656					}
657					L[(j, j)] = c64::ONE;
658				}
659				for j in 0..n {
660					for i in j + 1..m {
661						U[(i, j)] = c64::ZERO;
662					}
663				}
664				let L = L.as_ref();
665				let U = U.as_ref();
666				let U = U.subrows(0, n);
667				assert!(p.inverse() * L * U * q ~ A);
668			}
669			for n in [16, 24, 32, 128, 255, 256, 257] {
670				let approx_eq = CwiseMat(ApproxEq {
671					abs_tol: 1e-10,
672					rel_tol: 1e-10,
673				});
674				let A = CwiseMatDistribution {
675					nrows: n,
676					ncols: n,
677					dist: StandardNormal,
678				}
679				.rand::<Mat<f64>>(rng);
680				let A = A.as_ref();
681				let mut LU = A.cloned();
682				let row_perm = &mut *vec![0usize; n];
683				let row_perm_inv = &mut *vec![0usize; n];
684				let col_perm = &mut *vec![0usize; n];
685				let col_perm_inv = &mut *vec![0usize; n];
686				let (_, p, q) = lu_in_place(
687					LU.as_mut(),
688					row_perm,
689					row_perm_inv,
690					col_perm,
691					col_perm_inv,
692					par.reborrow(),
693					MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
694						usize,
695						f64,
696					>(
697						n, n, par_shape, default()
698					))),
699					default(),
700				);
701				let mut L = LU.as_ref().cloned();
702				let mut U = LU.as_ref().cloned();
703				for j in 0..n {
704					for i in 0..j {
705						L[(i, j)] = 0.0;
706					}
707					L[(j, j)] = 1.0;
708				}
709				for j in 0..n {
710					for i in j + 1..n {
711						U[(i, j)] = 0.0;
712					}
713				}
714				let L = L.as_ref();
715				let U = U.as_ref();
716				assert!(p.inverse() * L * U * q ~ A);
717			}
718		};
719		run(Par::Seq);
720		crate::team::with_team(8, |par| run(par));
721	}
722	/// The pivot search breaks a tie the way the scalar search does: the smallest column, then the
723	/// smallest row. Without that, which lane wins depends on the lane count, so an AVX2 build and
724	/// an AVX-512 build eliminate in a different order for the same matrix (and on a singular
725	/// fixture one of the two ends on an exact zero pivot, which the solve turns into `NaN`).
726	#[test]
727	fn pivot_search_breaks_ties_by_column_then_row() {
728		// column 2 holds the maximum twice (rows 1 and 5): both routes must take row 1, and a
729		// maximum in column 2 must win over an equal one in a later column
730		let a = crate::Mat::<f64>::from_fn(8, 6, |i, j| match (j, i) {
731			(2, 1) | (2, 5) => 3.0,
732			(4, 0) | (4, 7) => 3.0,
733			(1, _) => 0.5,
734			_ => 0.0,
735		});
736		assert_eq!(best_in_matrix(a.as_ref()), (1, 2, 3.0));
737		assert_eq!(best_in_matrix_fallback(a.as_ref()), (1, 2, 3.0));
738
739		// a matrix of small integers, where every column ties: the whole search must agree with
740		// the scalar one, whatever the lane count is
741		let mut rng = core::num::Wrapping(7u64);
742		for (m, n) in [(8usize, 8usize), (27, 13), (16, 40), (1, 9), (9, 1)] {
743			let mut a = crate::Mat::<f64>::zeros(m, n);
744			for j in 0..n {
745				for i in 0..m {
746					rng += core::num::Wrapping(6364136223846793005u64);
747					a[(i, j)] = ((rng.0 >> 33) % 3) as f64;
748				}
749			}
750			assert_eq!(
751				best_in_matrix(a.as_ref()),
752				best_in_matrix_fallback(a.as_ref()),
753				"{m}x{n}"
754			);
755		}
756	}
757
758	#[test]
759	fn test_gh238() {
760		#[rustfmt::skip]
761        let A = [
762            [
763                -2.0,
764                0.0,
765                1.0,
766                0.0,
767                0.0,
768                0.0,
769                0.0,
770                0.0,
771                0.0,
772                0.0,
773                0.0,
774                0.0,
775                0.0,
776                0.0,
777                0.0,
778                0.0,
779                0.0,
780                0.0,
781                0.0,
782                0.0,
783                0.0,
784                0.0,
785                0.0,
786                0.0,
787                0.0,
788                0.0,
789                0.0,
790                0.0,
791            ],
792            [
793                0.0,
794                -2.0,
795                0.0,
796                1.0,
797                0.0,
798                0.0,
799                0.0,
800                0.0,
801                0.0,
802                0.0,
803                0.0,
804                0.0,
805                0.0,
806                0.0,
807                0.0,
808                0.0,
809                0.0,
810                0.0,
811                0.0,
812                0.0,
813                0.0,
814                0.0,
815                0.0,
816                0.0,
817                0.0,
818                0.0,
819                0.0,
820                0.0,
821            ],
822            [
823                0.0,
824                0.0,
825                0.0,
826                0.0,
827                -2.0,
828                0.0,
829                1.0,
830                0.0,
831                0.0,
832                0.0,
833                0.0,
834                0.0,
835                0.0,
836                0.0,
837                0.0,
838                0.0,
839                0.0,
840                0.0,
841                0.0,
842                0.0,
843                0.0,
844                0.0,
845                0.0,
846                0.0,
847                0.0,
848                0.0,
849                0.0,
850                0.0,
851            ],
852            [
853                0.0,
854                0.0,
855                0.0,
856                0.0,
857                0.0,
858                -2.0,
859                0.0,
860                1.0,
861                0.0,
862                0.0,
863                0.0,
864                0.0,
865                0.0,
866                0.0,
867                0.0,
868                0.0,
869                0.0,
870                0.0,
871                0.0,
872                0.0,
873                0.0,
874                0.0,
875                0.0,
876                0.0,
877                0.0,
878                0.0,
879                0.0,
880                0.0,
881            ],
882            [
883                0.0,
884                0.0,
885                0.0,
886                0.0,
887                0.0,
888                0.0,
889                0.0,
890                0.0,
891                -2.0,
892                0.0,
893                1.0,
894                0.0,
895                0.0,
896                0.0,
897                0.0,
898                0.0,
899                0.0,
900                0.0,
901                0.0,
902                0.0,
903                0.0,
904                0.0,
905                0.0,
906                0.0,
907                0.0,
908                0.0,
909                0.0,
910                0.0,
911            ],
912            [
913                0.0,
914                0.0,
915                0.0,
916                0.0,
917                0.0,
918                0.0,
919                0.0,
920                0.0,
921                0.0,
922                -2.0,
923                0.0,
924                1.0,
925                0.0,
926                0.0,
927                0.0,
928                0.0,
929                0.0,
930                0.0,
931                0.0,
932                0.0,
933                0.0,
934                0.0,
935                0.0,
936                0.0,
937                0.0,
938                0.0,
939                0.0,
940                0.0,
941            ],
942            [
943                0.0,
944                0.0,
945                0.0,
946                0.0,
947                0.0,
948                0.0,
949                0.0,
950                0.0,
951                0.0,
952                0.0,
953                0.0,
954                0.0,
955                -2.0,
956                0.0,
957                1.0,
958                0.0,
959                0.0,
960                0.0,
961                0.0,
962                0.0,
963                0.0,
964                0.0,
965                0.0,
966                0.0,
967                0.0,
968                0.0,
969                0.0,
970                0.0,
971            ],
972            [
973                0.0,
974                0.0,
975                0.0,
976                0.0,
977                0.0,
978                0.0,
979                0.0,
980                0.0,
981                0.0,
982                0.0,
983                0.0,
984                0.0,
985                0.0,
986                -2.0,
987                0.0,
988                1.0,
989                0.0,
990                0.0,
991                0.0,
992                0.0,
993                0.0,
994                0.0,
995                0.0,
996                0.0,
997                0.0,
998                0.0,
999                0.0,
1000                0.0,
1001            ],
1002            [
1003                0.0,
1004                0.0,
1005                0.0,
1006                0.0,
1007                0.0,
1008                0.0,
1009                0.0,
1010                0.0,
1011                0.0,
1012                0.0,
1013                0.0,
1014                0.0,
1015                0.0,
1016                0.0,
1017                0.0,
1018                0.0,
1019                -2.0,
1020                0.0,
1021                1.0,
1022                0.0,
1023                0.0,
1024                0.0,
1025                0.0,
1026                0.0,
1027                0.0,
1028                0.0,
1029                0.0,
1030                0.0,
1031            ],
1032            [
1033                0.0,
1034                0.0,
1035                0.0,
1036                0.0,
1037                0.0,
1038                0.0,
1039                0.0,
1040                0.0,
1041                0.0,
1042                0.0,
1043                0.0,
1044                0.0,
1045                0.0,
1046                0.0,
1047                0.0,
1048                0.0,
1049                0.0,
1050                -2.0,
1051                0.0,
1052                1.0,
1053                0.0,
1054                0.0,
1055                0.0,
1056                0.0,
1057                0.0,
1058                0.0,
1059                0.0,
1060                0.0,
1061            ],
1062            [
1063                0.0,
1064                0.0,
1065                0.0,
1066                0.0,
1067                0.0,
1068                0.0,
1069                0.0,
1070                0.0,
1071                0.0,
1072                0.0,
1073                0.0,
1074                0.0,
1075                0.0,
1076                0.0,
1077                0.0,
1078                0.0,
1079                0.0,
1080                0.0,
1081                0.0,
1082                0.0,
1083                -2.0,
1084                0.0,
1085                1.0,
1086                0.0,
1087                0.0,
1088                0.0,
1089                0.0,
1090                0.0,
1091            ],
1092            [
1093                0.0,
1094                0.0,
1095                0.0,
1096                0.0,
1097                0.0,
1098                0.0,
1099                0.0,
1100                0.0,
1101                0.0,
1102                0.0,
1103                0.0,
1104                0.0,
1105                0.0,
1106                0.0,
1107                0.0,
1108                0.0,
1109                0.0,
1110                0.0,
1111                0.0,
1112                0.0,
1113                0.0,
1114                -2.0,
1115                0.0,
1116                1.0,
1117                0.0,
1118                0.0,
1119                0.0,
1120                0.0,
1121            ],
1122            [
1123                0.0,
1124                0.0,
1125                0.0,
1126                0.0,
1127                0.0,
1128                0.0,
1129                0.0,
1130                0.0,
1131                0.0,
1132                0.0,
1133                0.0,
1134                0.0,
1135                0.0,
1136                0.0,
1137                0.0,
1138                0.0,
1139                0.0,
1140                0.0,
1141                0.0,
1142                0.0,
1143                0.0,
1144                0.0,
1145                0.0,
1146                0.0,
1147                -2.0,
1148                0.0,
1149                1.0,
1150                0.0,
1151            ],
1152            [
1153                0.0,
1154                0.0,
1155                0.0,
1156                0.0,
1157                0.0,
1158                0.0,
1159                0.0,
1160                0.0,
1161                0.0,
1162                0.0,
1163                0.0,
1164                0.0,
1165                0.0,
1166                0.0,
1167                0.0,
1168                0.0,
1169                0.0,
1170                0.0,
1171                0.0,
1172                0.0,
1173                0.0,
1174                0.0,
1175                0.0,
1176                0.0,
1177                0.0,
1178                -2.0,
1179                0.0,
1180                1.0,
1181            ],
1182            [
1183                0.0,
1184                0.0,
1185                0.0,
1186                0.0,
1187                0.0,
1188                0.0,
1189                0.0,
1190                0.0,
1191                0.0,
1192                0.0,
1193                0.0,
1194                0.0,
1195                0.0,
1196                0.0,
1197                0.0,
1198                0.0,
1199                0.0,
1200                0.0,
1201                0.0,
1202                0.0,
1203                0.0,
1204                0.0,
1205                0.0,
1206                0.0,
1207                0.0,
1208                0.0,
1209                0.0,
1210                1.0,
1211            ],
1212            [
1213                0.0,
1214                0.0,
1215                0.0,
1216                1.0,
1217                0.0,
1218                0.0,
1219                1.0,
1220                0.0,
1221                0.0,
1222                0.0,
1223                0.0,
1224                0.0,
1225                0.0,
1226                0.0,
1227                0.0,
1228                0.0,
1229                0.0,
1230                0.0,
1231                0.0,
1232                0.0,
1233                0.0,
1234                0.0,
1235                0.0,
1236                0.0,
1237                0.0,
1238                0.0,
1239                0.0,
1240                0.0,
1241            ],
1242            [
1243                0.0,
1244                -1.0,
1245                0.0,
1246                0.0,
1247                1.0,
1248                0.0,
1249                0.0,
1250                0.0,
1251                0.0,
1252                0.0,
1253                0.0,
1254                0.0,
1255                0.0,
1256                0.0,
1257                0.0,
1258                0.0,
1259                0.0,
1260                0.0,
1261                0.0,
1262                0.0,
1263                0.0,
1264                0.0,
1265                0.0,
1266                0.0,
1267                0.0,
1268                0.0,
1269                0.0,
1270                0.0,
1271            ],
1272            [
1273                0.0,
1274                0.0,
1275                0.0,
1276                0.0,
1277                0.0,
1278                0.0,
1279                0.0,
1280                0.0,
1281                0.0,
1282                0.0,
1283                0.0,
1284                0.0,
1285                0.0,
1286                0.0,
1287                0.0,
1288                1.0,
1289                0.0,
1290                0.0,
1291                1.0,
1292                0.0,
1293                0.0,
1294                0.0,
1295                0.0,
1296                0.0,
1297                0.0,
1298                0.0,
1299                0.0,
1300                0.0,
1301            ],
1302            [
1303                0.0,
1304                0.0,
1305                0.0,
1306                0.0,
1307                0.0,
1308                0.0,
1309                0.0,
1310                0.0,
1311                0.0,
1312                0.0,
1313                0.0,
1314                0.0,
1315                0.0,
1316                -1.0,
1317                0.0,
1318                0.0,
1319                1.0,
1320                0.0,
1321                0.0,
1322                0.0,
1323                0.0,
1324                0.0,
1325                0.0,
1326                0.0,
1327                0.0,
1328                0.0,
1329                0.0,
1330                0.0,
1331            ],
1332            [
1333                0.0,
1334                0.0,
1335                0.0,
1336                0.0,
1337                0.0,
1338                0.0,
1339                0.0,
1340                1.0,
1341                0.0,
1342                0.0,
1343                1.0,
1344                0.0,
1345                0.0,
1346                0.0,
1347                0.0,
1348                0.0,
1349                0.0,
1350                0.0,
1351                0.0,
1352                0.0,
1353                0.0,
1354                0.0,
1355                0.0,
1356                0.0,
1357                0.0,
1358                0.0,
1359                0.0,
1360                0.0,
1361            ],
1362            [
1363                0.0,
1364                0.0,
1365                0.0,
1366                0.0,
1367                0.0,
1368                1.0,
1369                0.0,
1370                0.0,
1371                -1.0,
1372                0.0,
1373                0.0,
1374                0.0,
1375                0.0,
1376                0.0,
1377                0.0,
1378                0.0,
1379                0.0,
1380                0.0,
1381                0.0,
1382                0.0,
1383                0.0,
1384                0.0,
1385                0.0,
1386                0.0,
1387                0.0,
1388                0.0,
1389                0.0,
1390                0.0,
1391            ],
1392            [
1393                0.0,
1394                0.0,
1395                0.0,
1396                0.0,
1397                0.0,
1398                0.0,
1399                0.0,
1400                0.0,
1401                0.0,
1402                0.0,
1403                0.0,
1404                1.0,
1405                0.0,
1406                0.0,
1407                1.0,
1408                0.0,
1409                0.0,
1410                0.0,
1411                0.0,
1412                0.0,
1413                0.0,
1414                0.0,
1415                0.0,
1416                0.0,
1417                0.0,
1418                0.0,
1419                0.0,
1420                0.0,
1421            ],
1422            [
1423                0.0,
1424                0.0,
1425                0.0,
1426                0.0,
1427                0.0,
1428                0.0,
1429                0.0,
1430                0.0,
1431                0.0,
1432                1.0,
1433                0.0,
1434                0.0,
1435                -1.0,
1436                0.0,
1437                0.0,
1438                0.0,
1439                0.0,
1440                0.0,
1441                0.0,
1442                0.0,
1443                0.0,
1444                0.0,
1445                0.0,
1446                0.0,
1447                0.0,
1448                0.0,
1449                0.0,
1450                0.0,
1451            ],
1452            [
1453                0.0,
1454                0.0,
1455                0.0,
1456                0.0,
1457                0.0,
1458                0.0,
1459                0.0,
1460                0.0,
1461                0.0,
1462                0.0,
1463                0.0,
1464                0.0,
1465                0.0,
1466                0.0,
1467                0.0,
1468                0.0,
1469                0.0,
1470                0.0,
1471                0.0,
1472                0.0,
1473                0.0,
1474                0.0,
1475                0.0,
1476                1.0,
1477                0.0,
1478                0.0,
1479                1.0,
1480                0.0,
1481            ],
1482            [
1483                0.0,
1484                0.0,
1485                0.0,
1486                0.0,
1487                0.0,
1488                0.0,
1489                0.0,
1490                0.0,
1491                0.0,
1492                0.0,
1493                0.0,
1494                0.0,
1495                0.0,
1496                0.0,
1497                0.0,
1498                0.0,
1499                0.0,
1500                0.0,
1501                0.0,
1502                0.0,
1503                0.0,
1504                -1.0,
1505                0.0,
1506                0.0,
1507                1.0,
1508                0.0,
1509                0.0,
1510                0.0,
1511            ],
1512            [
1513                0.0,
1514                0.0,
1515                1.0,
1516                0.0,
1517                0.0,
1518                0.0,
1519                0.0,
1520                0.0,
1521                0.0,
1522                0.0,
1523                0.0,
1524                0.0,
1525                0.0,
1526                0.0,
1527                0.0,
1528                0.0,
1529                0.0,
1530                0.0,
1531                0.0,
1532                0.0,
1533                0.0,
1534                0.0,
1535                0.0,
1536                0.0,
1537                0.0,
1538                0.0,
1539                0.0,
1540                0.0,
1541            ],
1542            [
1543                0.0,
1544                0.0,
1545                0.0,
1546                0.0,
1547                0.0,
1548                0.0,
1549                0.0,
1550                0.0,
1551                0.0,
1552                0.0,
1553                0.0,
1554                0.0,
1555                0.0,
1556                0.0,
1557                0.0,
1558                0.0,
1559                0.0,
1560                0.0,
1561                0.0,
1562                1.0,
1563                0.0,
1564                0.0,
1565                1.0,
1566                0.0,
1567                0.0,
1568                0.0,
1569                0.0,
1570                0.0,
1571            ],
1572            [
1573                0.0,
1574                0.0,
1575                0.0,
1576                0.0,
1577                0.0,
1578                0.0,
1579                0.0,
1580                0.0,
1581                0.0,
1582                0.0,
1583                0.0,
1584                0.0,
1585                0.0,
1586                0.0,
1587                0.0,
1588                0.0,
1589                0.0,
1590                1.0,
1591                0.0,
1592                0.0,
1593                -1.0,
1594                0.0,
1595                0.0,
1596                0.0,
1597                0.0,
1598                0.0,
1599                0.0,
1600                0.0,
1601            ],
1602        ];
1603		let A = MatRef::from_row_major_array(&A).to_owned();
1604		let b = [
1605			-0.0, -0.0, -2.0, -2.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0,
1606			-0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1607			0.0, 0.0, 0.0, 0.0,
1608		];
1609		let b = ColRef::from_slice(&b);
1610		let lu = A.full_piv_lu(Par::Seq);
1611		assert!(lu.solve(b, Par::Seq).is_all_finite());
1612	}
1613}
faer/linalg/lu/full_pivoting/factor.rs

faer/linalg/lu/full_pivoting/
factor.rs
