Skip to main content

faer/
lib.rs

1//! `faer` is a general-purpose linear algebra library for rust, with a focus on
2//! high performance for algebraic operations on medium/large matrices, as well
3//! as matrix decompositions
4//!
5//! most of the high-level functionality in this library is provided through
6//! associated functions in its vocabulary types: [`Mat`]/[`MatRef`]/[`MatMut`]
7//!
8//! `faer` is recommended for applications that handle medium to large dense
9//! matrices, and its design is not well suited for applications that operate
10//! mostly on low dimensional vectors and matrices such as computer graphics or
11//! game development. for such applications, `nalgebra` and `cgmath` may be
12//! better suited
13//!
14//! # basic usage
15//!
16//! [`Mat`] is a resizable matrix type with dynamic capacity, which can be
17//! created using [`Mat::new`] to produce an empty $0\times 0$ matrix,
18//! [`Mat::zeros`] to create a rectangular matrix filled with zeros,
19//! [`Mat::identity`] to create an identity matrix, or [`Mat::from_fn`]
20//! for the more general case
21//!
22//! Given a `&Mat<T>` (resp. `&mut Mat<T>`), a [`MatRef<'_, T>`](MatRef) (resp.
23//! [`MatMut<'_, T>`](MatMut)) can be created by calling [`Mat::as_ref`] (resp.
24//! [`Mat::as_mut`]), which allow for more flexibility than `Mat` in that they
25//! allow slicing ([`MatRef::get`]) and splitting ([`MatRef::split_at`])
26//!
27//! `MatRef` and `MatMut` are lightweight view objects. the former can be copied
28//! freely while the latter has move and reborrow semantics, as described in its
29//! documentation
30//!
31//! most of the matrix operations can be used through the corresponding math
32//! operators: `+` for matrix addition, `-` for subtraction, `*` for either
33//! scalar or matrix multiplication depending on the types of the operands.
34//!
35//! ## example
36//! ```
37//! use faer::{Mat, Scale, mat};
38//! let a = mat![
39//! 	[1.0, 5.0, 9.0], //
40//! 	[2.0, 6.0, 10.0],
41//! 	[3.0, 7.0, 11.0],
42//! 	[4.0, 8.0, 12.0f64],
43//! ];
44//! let b = Mat::from_fn(4, 3, |i, j| (i + j) as f64);
45//! let add = &a + &b;
46//! let sub = &a - &b;
47//! let scale = Scale(3.0) * &a;
48//! let mul = &a * b.transpose();
49//! let a00 = a[(0, 0)];
50//! ```
51//!
52//! # matrix decompositions
53//! `faer` provides a variety of matrix factorizations, each with its own
54//! advantages and drawbacks:
55//!
56//! ## $LL^\top$ decomposition
57//! [`Mat::llt`] decomposes a self-adjoint positive definite matrix $A$ such
58//! that $$A = LL^H,$$
59//! where $L$ is a lower triangular matrix. this decomposition is highly
60//! efficient and has good stability properties
61//!
62//! [an implementation for sparse matrices is also
63//! available](sparse::linalg::solvers::Llt)
64//!
65//! ## $LBL^\top$ decomposition
66//! [`Mat::lblt`] decomposes a self-adjoint (possibly indefinite) matrix $A$
67//! such that $$P A P^\top = LBL^H,$$
68//! where $P$ is a permutation matrix, $L$ is a lower triangular matrix, and $B$
69//! is a block diagonal matrix, with $1 \times 1$ or $2 \times 2$ diagonal
70//! blocks. this decomposition is efficient and has good stability properties
71//!
72//! ## $LU$ decomposition with partial pivoting
73//! [`Mat::partial_piv_lu`] decomposes a square invertible matrix $A$ into a
74//! lower triangular matrix $L$, a unit upper triangular matrix $U$, and a
75//! permutation matrix $P$, such that $$PA = LU$$
76//! it is used by default for computing the determinant, and is generally the
77//! recommended method for solving a square linear system or computing the
78//! inverse of a matrix (although we generally recommend using a
79//! [`faer::linalg::solvers::Solve`](crate::linalg::solvers::Solve) instead of
80//! computing the inverse explicitly)
81//!
82//! [an implementation for sparse matrices is also
83//! available](sparse::linalg::solvers::Lu)
84//!
85//! ## $LU$ decomposition with full pivoting
86//! [`Mat::full_piv_lu`] decomposes a generic rectangular matrix $A$ into a
87//! lower triangular matrix $L$, a unit upper triangular matrix $U$, and
88//! permutation matrices $P$ and $Q$, such that $$PAQ^\top = LU$$
89//! it can be more stable than the LU decomposition with partial pivoting, in
90//! exchange for being more computationally expensive
91//!
92//! ## $QR$ decomposition
93//! [`Mat::qr`] decomposes a matrix $A$ into the product $$A = QR,$$
94//! where $Q$ is a unitary matrix, and $R$ is an upper trapezoidal matrix. it is
95//! often used for solving least squares problems
96//!
97//! [an implementation for sparse matrices is also
98//! available](sparse::linalg::solvers::Qr)
99//!
100//! ## $QR$ decomposition with column pivoting
101//! ([`Mat::col_piv_qr`]) decomposes a matrix $A$ into the product $$AP^\top =
102//! QR,$$ where $P$ is a permutation matrix, $Q$ is a unitary matrix, and $R$ is
103//! an upper trapezoidal matrix
104//!
105//! it is slower than the version with no pivoting, in exchange for being more
106//! numerically stable for rank-deficient matrices
107//!
108//! ## singular value decomposition
109//! the SVD of a matrix $A$ of shape $(m, n)$ is a decomposition into three
110//! components $U$, $S$, and $V$, such that:
111//!
112//! - $U$ has shape $(m, m)$ and is a unitary matrix,
113//! - $V$ has shape $(n, n)$ and is a unitary matrix,
114//! - $S$ has shape $(m, n)$ and is zero everywhere except the main diagonal,
115//!   with nonnegative
116//! diagonal values in nonincreasing order,
117//! - and finally:
118//!
119//! $$A = U S V^H$$
120//!
121//! the SVD is provided in two forms: either the full matrices $U$ and $V$ are
122//! computed, using [`Mat::svd`], or only their first $\min(m, n)$ columns are
123//! computed, using [`Mat::thin_svd`]
124//!
125//! if only the singular values (elements of $S$) are desired, they can be
126//! obtained in nonincreasing order using [`Mat::singular_values`]
127//!
128//! ## eigendecomposition
129//! **note**: the order of the eigenvalues is currently unspecified and may be
130//! changed in a future release
131//!
132//! the eigenvalue decomposition of a square matrix $A$ of shape $(n, n)$ is a
133//! decomposition into two components $U$, $S$:
134//!
135//! - $U$ has shape $(n, n)$ and is invertible,
136//! - $S$ has shape $(n, n)$ and is a diagonal matrix,
137//! - and finally:
138//!
139//! $$A = U S U^{-1}$$
140//!
141//! if $A$ is self-adjoint, then $U$ can be made unitary ($U^{-1} = U^H$), and
142//! $S$ is real valued. additionally, the eigenvalues are sorted in
143//! nondecreasing order
144//!
145//! Depending on the domain of the input matrix and whether it is self-adjoint,
146//! multiple methods are provided to compute the eigendecomposition:
147//! * [`Mat::self_adjoint_eigen`] can be used with either real or complex
148//!   matrices,
149//! producing an eigendecomposition of the same type,
150//! * [`Mat::eigen`] can be used with real or complex matrices, but always
151//!   produces complex values.
152//!
153//! if only the eigenvalues (elements of $S$) are desired, they can be obtained
154//! using [`Mat::self_adjoint_eigenvalues`] (nondecreasing order),
155//! [`Mat::eigenvalues`] with the same conditions described above.
156//!
157//! # parallelism
158//! there is no process wide parallelism setting. Every function that can run in parallel takes a
159//! [`Par`] argument, the high level ones included ([`Mat::svd`], [`Mat::llt`],
160//! [`Mat::singular_values`], the `solve` family of [`linalg::solvers::Solve`], ...). Pass
161//! [`Par::Seq`] for sequential execution, or [`Par::rayon`] for the global rayon threadpool. The
162//! operator overloads (`&a * &b`, `&a + &b`, ...) have no room for the argument and always run
163//! sequentially; use [`linalg::matmul::matmul`] for a parallel product.
164//!
165//! # crate features
166//!
167//! - `std`: enabled by default. links with the standard library to enable
168//!   additional features such as cpu feature detection at runtime
169//! - `rayon`: enabled by default. enables the `rayon` parallel backend
170//! - `npy`: enables conversions to/from numpy's matrix file format
171//! - `perf-warn`: produces performance warnings when matrix operations are
172//!   called with suboptimal
173//! data layout
174//! - `nightly`: requires the nightly compiler. enables experimental simd
175//!   features such as avx512
176#![cfg_attr(not(feature = "std"), no_std)]
177#![allow(non_snake_case)]
178#![warn(rustdoc::broken_intra_doc_links)]
179extern crate alloc;
180#[cfg(feature = "std")]
181extern crate std;
182/// see: [`generativity::make_guard`]
183#[macro_export]
184macro_rules! make_guard {
185    ($($name:ident),* $(,)?) => {
186        $(#[allow(unused_unsafe)] let $name = unsafe { $crate::generativity::Id::new() };
187        #[allow(unused, unused_unsafe)] let lifetime_brand = unsafe {
188        $crate::generativity::LifetimeBrand::new(&$name) }; #[allow(unused_unsafe)] let
189        $name = unsafe { $crate::generativity::Guard::new($name) };)*
190    };
191}
192macro_rules! repeat_n {
193	($e:expr, $n:expr) => {
194		iter::repeat_n($e, $n)
195	};
196}
197use core::num::NonZeroUsize;
198use equator::{assert, debug_assert};
199use faer_traits::*;
200/// shorthand for `<_ as Auto::<T>>::auto()`
201#[macro_export]
202macro_rules! auto {
203	($ty:ty $(,)?) => {
204		$crate::Auto::<$ty>::auto()
205	};
206}
207macro_rules! dispatch {
208	($imp:expr, $ty:ident, $T:ty $(,)?) => {
209		if const { <$T>::IS_NATIVE_C32 } {
210			unsafe {
211				transmute(
212					<ComplexImpl<f32> as ComplexField>::Arch::default()
213						.dispatch(transmute::<_, $ty<ComplexImpl<f32>>>($imp)),
214				)
215			}
216		} else if const { <$T>::IS_NATIVE_C64 } {
217			unsafe {
218				transmute(
219					<ComplexImpl<f64> as ComplexField>::Arch::default()
220						.dispatch(transmute::<_, $ty<ComplexImpl<f64>>>($imp)),
221				)
222			}
223		} else {
224			<$T>::Arch::default().dispatch($imp)
225		}
226	};
227}
228macro_rules! stack_mat {
229	($name:ident, $m:expr, $n:expr, $A:expr, $N:expr, $T:ty $(,)?) => {
230		let mut __tmp = {
231			#[repr(align(64))]
232			struct __Col<T, const A: usize>([T; A]);
233			struct __Mat<T, const A: usize, const N: usize>([__Col<T, A>; N]);
234			core::mem::MaybeUninit::<__Mat<$T, $A, $N>>::uninit()
235		};
236		let __stack = MemStack::new_any(core::slice::from_mut(&mut __tmp));
237		let mut $name =
238			$crate::linalg::temp_mat_zeroed::<$T, _, _>($m, $n, __stack).0;
239		let mut $name = $crate::mat::AsMatMut::as_mat_mut(&mut $name);
240	};
241	($name:ident, $m:expr, $n:expr, $T:ty $(,)?) => {
242		stack_mat!($name, $m, $n, $m, $n, $T)
243	};
244}
245#[macro_export]
246#[doc(hidden)]
247macro_rules! __dbg {
248    () => {
249        std::eprintln!("[{}:{}:{}]", std::file!(), std::line!(), std::column!())
250    };
251    ($val:expr $(,)?) => {
252        match $val { tmp => { std::eprintln!("[{}:{}:{}] {} = {:20?}", std::file!(),
253        std::line!(), std::column!(), std::stringify!($val), & tmp); tmp } }
254    };
255    ($($val:expr),+ $(,)?) => {
256        ($($crate::__dbg!($val)),+,)
257    };
258}
259#[cfg(feature = "perf-warn")]
260#[macro_export]
261#[doc(hidden)]
262macro_rules! __perf_warn {
263	($name:ident) => {{
264		#[inline(always)]
265		#[allow(non_snake_case)]
266		fn $name() -> &'static ::core::sync::atomic::AtomicBool {
267			static $name: ::core::sync::atomic::AtomicBool =
268				::core::sync::atomic::AtomicBool::new(false);
269			&$name
270		}
271		::core::matches!(
272			$name().compare_exchange(
273				false,
274				true,
275				::core::sync::atomic::Ordering::Relaxed,
276				::core::sync::atomic::Ordering::Relaxed,
277			),
278			Ok(_)
279		)
280	}};
281}
282#[doc(hidden)]
283#[macro_export]
284macro_rules! with_dim {
285    ($name:ident, $value:expr $(,)?) => {
286        let __val__ = $value; $crate::make_guard!($name); let $name =
287        $crate::utils::bound::Dim::new(__val__, $name);
288    };
289    ({ $(let $name:ident = $value:expr;)* }) => {
290        $(let __val__ = $value; $crate::make_guard!($name); let $name =
291        $crate::utils::bound::Dim::new(__val__, $name);)*
292    };
293}
294/// zips together matrix of the same size, so that coefficient-wise operations
295/// can be performed on their elements.
296///
297/// # note
298/// the order in which the matrix elements are traversed is unspecified.
299///
300/// # example
301/// ```
302/// use faer::{Mat, mat, unzip, zip};
303/// let nrows = 2;
304/// let ncols = 3;
305/// let a = mat![[1.0, 3.0, 5.0], [2.0, 4.0, 6.0]];
306/// let b = mat![[7.0, 9.0, 11.0], [8.0, 10.0, 12.0]];
307/// let mut sum = Mat::<f64>::zeros(nrows, ncols);
308/// zip!(&mut sum, &a, &b).for_each(|unzip!(sum, a, b)| {
309/// 	*sum = a + b;
310/// });
311/// for i in 0..nrows {
312/// 	for j in 0..ncols {
313/// 		assert_eq!(sum[(i, j)], a[(i, j)] + b[(i, j)]);
314/// 	}
315/// }
316/// ```
317#[macro_export]
318macro_rules! zip {
319    ($head:expr $(,)?) => {
320        $crate::linalg::zip::LastEq($crate::linalg::zip::IntoView::into_view($head),
321        ::core::marker::PhantomData)
322    };
323    ($head:expr, $($tail:expr),* $(,)?) => {
324        $crate::linalg::zip::ZipEq::new($crate::linalg::zip::IntoView::into_view($head),
325        $crate::zip!($($tail,)*))
326    };
327}
328/// expands to the type of zipped items.
329///
330/// # example
331/// ```
332/// use faer::{Mat, Zip, mat, unzip, zip};
333/// let nrows = 2;
334/// let ncols = 3;
335/// let a = mat![[1.0, 3.0, 5.0], [2.0, 4.0, 6.0]];
336/// let b = mat![[7.0, 9.0, 11.0], [8.0, 10.0, 12.0]];
337/// let mut sum = Mat::<f64>::zeros(nrows, ncols);
338/// zip!(&mut sum, &a, &b).for_each(
339/// 	|unzip!(sum, a, b): Zip!(&mut f64, &f64, &f64)| {
340/// 		*sum = a + b;
341/// 	},
342/// );
343/// for i in 0..nrows {
344/// 	for j in 0..ncols {
345/// 		assert_eq!(sum[(i, j)], a[(i, j)] + b[(i, j)]);
346/// 	}
347/// }
348/// ```
349#[macro_export]
350macro_rules! Zip {
351    (..$(,)?) => {
352        _
353    };
354    ($head:ty $(,)?) => {
355        $crate::linalg::zip::Last::<$head >
356    };
357    ($head:ty, $($tail:tt)*) => {
358        $crate::linalg::zip::Zip::<$head, $crate::Zip!($($tail)*) >
359    };
360}
361/// used to undo the zipping by the [`zip!`] macro.
362///
363/// # example
364/// ```
365/// use faer::{Mat, mat, unzip, zip};
366/// let nrows = 2;
367/// let ncols = 3;
368/// let a = mat![[1.0, 3.0, 5.0], [2.0, 4.0, 6.0]];
369/// let b = mat![[7.0, 9.0, 11.0], [8.0, 10.0, 12.0]];
370/// let mut sum = Mat::<f64>::zeros(nrows, ncols);
371/// zip!(&mut sum, &a, &b).for_each(|unzip!(sum, a, b)| {
372/// 	*sum = a + b;
373/// });
374/// for i in 0..nrows {
375/// 	for j in 0..ncols {
376/// 		assert_eq!(sum[(i, j)], a[(i, j)] + b[(i, j)]);
377/// 	}
378/// }
379/// ```
380#[macro_export]
381macro_rules! unzip {
382    (..$(,)?) => {
383        _
384    };
385    ($head:pat $(,)?) => {
386        $crate::linalg::zip::Last($head)
387    };
388    ($head:pat, $($tail:tt)*) => {
389        $crate::linalg::zip::Zip($head, $crate::unzip!($($tail)*))
390    };
391}
392#[macro_export]
393#[doc(hidden)]
394macro_rules! __transpose_impl {
395    ([$([$($col:expr),*])*] $($v:expr;)*) => {
396        [$([$($col,)*],)* [$($v,)*]]
397    };
398    ([$([$($col:expr),*])*] $($v0:expr, $($v:expr),*;)*) => {
399        $crate::__transpose_impl!([$([$($col),*])* [$($v0),*]] $($($v),*;)*)
400    };
401}
402/// creates a [`Mat`] containing the arguments.
403///
404/// ```
405/// use faer::mat;
406/// let matrix = mat![
407/// 	[1.0, 5.0, 9.0], //
408/// 	[2.0, 6.0, 10.0],
409/// 	[3.0, 7.0, 11.0],
410/// 	[4.0, 8.0, 12.0f64],
411/// ];
412/// assert_eq!(matrix[(0, 0)], 1.0);
413/// assert_eq!(matrix[(1, 0)], 2.0);
414/// assert_eq!(matrix[(2, 0)], 3.0);
415/// assert_eq!(matrix[(3, 0)], 4.0);
416/// assert_eq!(matrix[(0, 1)], 5.0);
417/// assert_eq!(matrix[(1, 1)], 6.0);
418/// assert_eq!(matrix[(2, 1)], 7.0);
419/// assert_eq!(matrix[(3, 1)], 8.0);
420/// assert_eq!(matrix[(0, 2)], 9.0);
421/// assert_eq!(matrix[(1, 2)], 10.0);
422/// assert_eq!(matrix[(2, 2)], 11.0);
423/// assert_eq!(matrix[(3, 2)], 12.0);
424/// ```
425#[macro_export]
426macro_rules! mat {
427    () => {
428        { compile_error!("number of columns in the matrix is ambiguous"); }
429    };
430    ($([$($v:expr),* $(,)?]),* $(,)?) => {
431        { let __data = ::core::mem::ManuallyDrop::new($crate::__transpose_impl!([]
432        $($($v),*;)*)); let __data = &* __data; let __ncols = __data.len(); let __nrows =
433        (* __data.get(0).unwrap()).len(); #[allow(unused_unsafe)] unsafe {
434        $crate::mat::Mat::from_fn(__nrows, __ncols, | i, j | ::core::ptr::from_ref(&
435        __data[j] [i]).read()) } }
436    };
437    (
438        @ alloc $($unsafe:ident)? ($stack:ident) ($var:pat) (uninit::<$ty:ty >, $($arg:expr),+
439        $(,)?)
440    ) => {
441        let (mut __mat__, $stack) = $($unsafe)? { $crate::linalg::temp_mat_uninit::<$ty, _, _
442        > ($($arg,)* $stack,) }; let $var = $crate::mat::AsMatMut::as_mat_mut(&mut __mat__);
443    };
444    (
445        @ alloc $($unsafe:ident)? ($stack:ident) ($var:pat) (zero::<$ty:ty >, $($arg:expr),+
446        $(,)?)
447    ) => {
448        let (mut __mat__, $stack) = $($unsafe)? { $crate::linalg::temp_mat_zeroed::<$ty, _, _
449        > ($($arg,)* $stack,) }; let $var = $crate::mat::AsMatMut::as_mat_mut(&mut __mat__);
450    };
451}
452/// creates a [`col::Col`] containing the arguments
453///
454/// ```
455/// use faer::col;
456/// let col_vec = col![3.0, 5.0, 7.0, 9.0];
457/// assert_eq!(col_vec[0], 3.0);
458/// assert_eq!(col_vec[1], 5.0);
459/// assert_eq!(col_vec[2], 7.0);
460/// assert_eq!(col_vec[3], 9.0);
461/// ```
462#[macro_export]
463macro_rules! col {
464    ($($v:expr),* $(,)?) => {
465        { let __data = ::core::mem::ManuallyDrop::new([$($v,)*]); let __data = &* __data;
466        let __len = __data.len(); #[allow(unused_unsafe)] unsafe {
467        $crate::col::Col::from_fn(__len, | i | ::core::ptr::from_ref(& __data[i]).read())
468        } }
469    };
470    (
471        @ alloc $($unsafe:ident)? ($stack:ident) ($var:pat) (uninit::<$ty:ty >, $($arg:expr),+
472        $(,)?)
473    ) => {
474        let (mut __mat__, $stack) = $($unsafe)? { $crate::linalg::temp_mat_uninit::<$ty, _, _
475        > ($($arg,)* 1, $stack,) }; let $var = $crate::mat::AsMatMut::as_mat_mut(&mut __mat__).col_mut(0);
476    };
477
478    (
479        @ alloc $($unsafe:ident)? ($stack:ident) ($var:pat) (zero::<$ty:ty >, $($arg:expr),+
480        $(,)?)
481    ) => {
482        let (mut __mat__, $stack) = $($unsafe)? { $crate::linalg::temp_mat_zeroed::<$ty, _, _
483        > ($($arg,)* 1, $stack,) }; let $var = $crate::mat::AsMatMut::as_mat_mut(&mut __mat__).col_mut(0);
484    };
485}
486/// creates a [`row::Row`] containing the arguments
487///
488/// ```
489/// use faer::row;
490/// let row_vec = row![3.0, 5.0, 7.0, 9.0];
491/// assert_eq!(row_vec[0], 3.0);
492/// assert_eq!(row_vec[1], 5.0);
493/// assert_eq!(row_vec[2], 7.0);
494/// assert_eq!(row_vec[3], 9.0);
495/// ```
496#[macro_export]
497macro_rules! row {
498    ($($v:expr),* $(,)?) => {
499        { let __data = ::core::mem::ManuallyDrop::new([$($v,)*]); let __data = &* __data;
500        let __len = __data.len(); #[allow(unused_unsafe)] unsafe {
501        $crate::row::Row::from_fn(__len, | i | ::core::ptr::from_ref(& __data[i]).read())
502        } }
503    };
504    (
505        @ alloc $($unsafe:ident)? ($stack:ident) ($var:pat) (uninit::<$ty:ty >, $($arg:expr),+
506        $(,)?)
507    ) => {
508        let (mut __mat__, $stack) = $($unsafe)? { $crate::linalg::temp_mat_uninit::<$ty, _, _
509        > ($($arg,)* 1, $stack,) }; let $var = $crate::mat::AsMatMut::as_mat_mut(&mut __mat__).col_mut(0).transpose_mut();
510    };
511
512    (
513        @ alloc $($unsafe:ident)? ($stack:ident) ($var:pat) (zero::<$ty:ty >, $($arg:expr),+
514        $(,)?)
515    ) => {
516        let (mut __mat__, $stack) = $($unsafe)? { $crate::linalg::temp_mat_zeroed::<$ty, _, _
517        > ($($arg,)* 1, $stack,) }; let $var = $crate::mat::AsMatMut::as_mat_mut(&mut __mat__).col_mut(0).transpose_mut();
518    };
519}
520/// convenience function to concatenate a nested list of matrices into a single
521/// big ['Mat']. concatonation pattern follows the numpy.block convention that
522/// each sub-list must have an equal number of columns (net) but the boundaries
523/// do not need to align. in other words, this sort of thing:
524/// ```notcode
525///   AAAbb
526///   AAAbb
527///   cDDDD
528/// ```
529/// is perfectly acceptable
530#[doc(hidden)]
531#[track_caller]
532pub fn concat_impl<T: ComplexField>(
533	blocks: &[&[(mat::MatRef<'_, T>, Conj)]],
534) -> mat::Mat<T> {
535	#[inline(always)]
536	fn count_total_columns<T: ComplexField>(
537		block_row: &[(mat::MatRef<'_, T>, Conj)],
538	) -> usize {
539		let mut out: usize = 0;
540		for (elem, _) in block_row.iter() {
541			out += elem.ncols();
542		}
543		out
544	}
545	#[inline(always)]
546	#[track_caller]
547	fn count_rows<T: ComplexField>(
548		block_row: &[(mat::MatRef<'_, T>, Conj)],
549	) -> usize {
550		let mut out: usize = 0;
551		for (i, (e, _)) in block_row.iter().enumerate() {
552			if i == 0 {
553				out = e.nrows();
554			} else {
555				assert!(e.nrows() == out);
556			}
557		}
558		out
559	}
560	let mut n: usize = 0;
561	let mut m: usize = 0;
562	for row in blocks.iter() {
563		n += count_rows(row);
564	}
565	for (i, row) in blocks.iter().enumerate() {
566		let cols = count_total_columns(row);
567		if i == 0 {
568			m = cols;
569		} else {
570			assert!(cols == m);
571		}
572	}
573	let mut mat = mat::Mat::<T>::zeros(n, m);
574	let mut ni: usize = 0;
575	let mut mj: usize;
576	for row in blocks.iter() {
577		mj = 0;
578		for (elem, conj) in row.iter() {
579			let mut dst =
580				mat.as_mut()
581					.submatrix_mut(ni, mj, elem.nrows(), elem.ncols());
582			if *conj == Conj::No {
583				dst.copy_from(elem);
584			} else {
585				dst.copy_from(elem.conjugate());
586			}
587			mj += elem.ncols();
588		}
589		ni += row[0].0.nrows();
590	}
591	mat
592}
593/// concatenates the matrices in each row horizontally,
594/// then concatenates the results vertically
595///
596/// `concat![[a0, a1, a2], [b1, b2]]` results in the matrix
597/// ```notcode
598/// [a0 | a1 | a2][b0 | b1]
599/// ```
600#[macro_export]
601macro_rules! concat {
602    () => {
603        { compile_error!("number of columns in the matrix is ambiguous"); }
604    };
605    ($([$($v:expr),* $(,)?]),* $(,)?) => {
606        { $crate::concat_impl(& [$(& [$(($v).as_ref().__canonicalize(),)*],)*]) }
607    };
608}
609/// column vector
610pub mod col;
611/// diagonal matrix
612pub mod diag;
613/// de-serialization from common matrix file formats
614#[cfg(feature = "std")]
615pub mod io;
616pub mod linalg;
617/// rectangular matrix
618pub mod mat;
619#[path = "./operator/mod.rs"]
620pub mod matrix_free;
621/// permutation matrix
622pub mod perm;
623/// row vector
624pub mod row;
625pub mod sparse;
626/// helper utilities
627pub mod utils;
628/// native unsigned integer type
629pub trait Index: traits::IndexCore + traits::Index + seal::Seal {}
630impl<T: faer_traits::Index<Signed: seal::Seal> + seal::Seal> Index for T {}
631mod seal {
632	pub trait Seal {}
633	impl<T: faer_traits::Seal> Seal for T {}
634	impl Seal for crate::utils::bound::Dim<'_> {}
635	impl<I: crate::Index> Seal for crate::utils::bound::Idx<'_, I> {}
636	impl<I: crate::Index> Seal for crate::utils::bound::IdxInc<'_, I> {}
637	impl<I: crate::Index> Seal for crate::utils::bound::MaybeIdx<'_, I> {}
638	impl<I: crate::Index> Seal for crate::utils::bound::IdxIncOne<I> {}
639	impl<I: crate::Index> Seal for crate::utils::bound::MaybeIdxOne<I> {}
640	impl Seal for crate::utils::bound::One {}
641	impl Seal for crate::utils::bound::Zero {}
642	impl Seal for crate::ContiguousFwd {}
643	impl Seal for crate::ContiguousBwd {}
644}
645/// sealed trait for types that can be created from "unbound" values, as long as
646/// their struct preconditions are upheld
647pub trait Unbind<I = usize>:
648	Send + Sync + Copy + core::fmt::Debug + seal::Seal
649{
650	/// creates new value
651	/// # safety
652	/// safety invariants must be upheld
653	unsafe fn new_unbound(idx: I) -> Self;
654	/// returns the unbound value, unconstrained by safety invariants
655	fn unbound(self) -> I;
656}
657/// type that can be used to index into a range
658pub type Idx<Dim, I = usize> = <Dim as ShapeIdx>::Idx<I>;
659/// type that can be used to partition a range
660pub type IdxInc<Dim, I = usize> = <Dim as ShapeIdx>::IdxInc<I>;
661/// either an index or a negative value
662pub type MaybeIdx<Dim, I = usize> = <Dim as ShapeIdx>::MaybeIdx<I>;
663/// base trait for [`Shape`]
664pub trait ShapeIdx {
665	/// type that can be used to index into a range
666	type Idx<I: Index>: Unbind<I> + Ord + Eq;
667	/// type that can be used to partition a range
668	type IdxInc<I: Index>: Unbind<I> + Ord + Eq + From<Idx<Self, I>>;
669	/// either an index or a negative value
670	type MaybeIdx<I: Index>: Unbind<I::Signed> + Ord + Eq;
671}
672/// matrix dimension
673pub trait Shape:
674	Unbind
675	+ Ord
676	+ ShapeIdx<
677		Idx<usize>: Ord + Eq + PartialOrd<Self>,
678		IdxInc<usize>: Ord + Eq + PartialOrd<Self>,
679	>
680{
681	/// whether the types involved have any safety invariants
682	const IS_BOUND: bool = true;
683	/// bind the current value using a invariant lifetime guard
684	#[inline]
685	fn bind<'n>(self, guard: generativity::Guard<'n>) -> utils::bound::Dim<'n> {
686		utils::bound::Dim::new(self.unbound(), guard)
687	}
688	/// cast a slice of bound values to unbound values
689	#[inline]
690	fn cast_idx_slice<I: Index>(slice: &[Idx<Self, I>]) -> &[I] {
691		unsafe { core::slice::from_raw_parts(slice.as_ptr() as _, slice.len()) }
692	}
693	/// cast a slice of bound values to unbound values
694	#[inline]
695	fn cast_idx_inc_slice<I: Index>(slice: &[IdxInc<Self, I>]) -> &[I] {
696		unsafe { core::slice::from_raw_parts(slice.as_ptr() as _, slice.len()) }
697	}
698	/// returns the index `0`, which is always valid
699	#[inline(always)]
700	fn start() -> IdxInc<Self> {
701		unsafe { IdxInc::<Self>::new_unbound(0) }
702	}
703	/// returns the incremented value, as an inclusive index
704	#[inline(always)]
705	fn next(idx: Idx<Self>) -> IdxInc<Self> {
706		unsafe { IdxInc::<Self>::new_unbound(idx.unbound() + 1) }
707	}
708	/// returns the last value, equal to the dimension
709	#[inline(always)]
710	fn end(self) -> IdxInc<Self> {
711		unsafe { IdxInc::<Self>::new_unbound(self.unbound()) }
712	}
713	/// checks if the index is valid, returning `Some(_)` in that case
714	#[inline(always)]
715	fn idx(self, idx: usize) -> Option<Idx<Self>> {
716		if idx < self.unbound() {
717			Some(unsafe { Idx::<Self>::new_unbound(idx) })
718		} else {
719			None
720		}
721	}
722	/// checks if the index is valid, returning `Some(_)` in that case
723	#[inline(always)]
724	fn idx_inc(self, idx: usize) -> Option<IdxInc<Self>> {
725		if idx <= self.unbound() {
726			Some(unsafe { IdxInc::<Self>::new_unbound(idx) })
727		} else {
728			None
729		}
730	}
731	/// checks if the index is valid, and panics otherwise
732	#[inline(always)]
733	fn checked_idx(self, idx: usize) -> Idx<Self> {
734		equator::assert!(idx < self.unbound());
735		unsafe { Idx::<Self>::new_unbound(idx) }
736	}
737	/// checks if the index is valid, and panics otherwise
738	#[inline(always)]
739	fn checked_idx_inc(self, idx: usize) -> IdxInc<Self> {
740		equator::assert!(idx <= self.unbound());
741		unsafe { IdxInc::<Self>::new_unbound(idx) }
742	}
743	/// assumes the index is valid
744	/// # safety
745	/// the index must be valid
746	#[inline(always)]
747	unsafe fn unchecked_idx(self, idx: usize) -> Idx<Self> {
748		equator::debug_assert!(idx < self.unbound());
749		unsafe { Idx::<Self>::new_unbound(idx) }
750	}
751	/// assumes the index is valid
752	/// # safety
753	/// the index must be valid
754	#[inline(always)]
755	unsafe fn unchecked_idx_inc(self, idx: usize) -> IdxInc<Self> {
756		equator::debug_assert!(idx <= self.unbound());
757		unsafe { IdxInc::<Self>::new_unbound(idx) }
758	}
759	/// returns an iterator over the indices between `from` and `to`
760	#[inline(always)]
761	fn indices(
762		from: IdxInc<Self>,
763		to: IdxInc<Self>,
764	) -> impl Clone + ExactSizeIterator + DoubleEndedIterator<Item = Idx<Self>>
765	{
766		(from.unbound()..to.unbound()).map(
767			#[inline(always)]
768			|i| unsafe { Idx::<Self>::new_unbound(i) },
769		)
770	}
771}
772impl<T: Send + Sync + Copy + core::fmt::Debug + faer_traits::Seal> Unbind<T>
773	for T
774{
775	#[inline(always)]
776	unsafe fn new_unbound(idx: T) -> Self {
777		idx
778	}
779
780	#[inline(always)]
781	fn unbound(self) -> T {
782		self
783	}
784}
785impl ShapeIdx for usize {
786	type Idx<I: Index> = I;
787	type IdxInc<I: Index> = I;
788	type MaybeIdx<I: Index> = I::Signed;
789}
790impl Shape for usize {
791	const IS_BOUND: bool = false;
792}
793/// stride distance between two consecutive elements along a given dimension
794pub trait Stride:
795	seal::Seal + core::fmt::Debug + Copy + Send + Sync + 'static
796{
797	/// the reversed stride type
798	type Rev: Stride<Rev = Self>;
799	/// returns the reversed stride
800	fn rev(self) -> Self::Rev;
801	/// returns the stride in elements
802	fn element_stride(self) -> isize;
803}
804impl Stride for isize {
805	type Rev = Self;
806
807	#[inline(always)]
808	fn rev(self) -> Self::Rev {
809		-self
810	}
811
812	#[inline(always)]
813	fn element_stride(self) -> isize {
814		self
815	}
816}
817/// contiguous stride equal to `+1`
818#[derive(Copy, Clone, Debug)]
819pub struct ContiguousFwd;
820/// contiguous stride equal to `-1`
821#[derive(Copy, Clone, Debug)]
822pub struct ContiguousBwd;
823impl Stride for ContiguousFwd {
824	type Rev = ContiguousBwd;
825
826	#[inline(always)]
827	fn rev(self) -> Self::Rev {
828		ContiguousBwd
829	}
830
831	#[inline(always)]
832	fn element_stride(self) -> isize {
833		1
834	}
835}
836impl Stride for ContiguousBwd {
837	type Rev = ContiguousFwd;
838
839	#[inline(always)]
840	fn rev(self) -> Self::Rev {
841		ContiguousFwd
842	}
843
844	#[inline(always)]
845	fn element_stride(self) -> isize {
846		-1
847	}
848}
849/// memory allocation error
850#[derive(Copy, Clone, Debug, PartialEq, Eq)]
851pub enum TryReserveError {
852	///rRequired allocation does not fit within `isize` bytes
853	CapacityOverflow,
854	/// allocator could not provide an allocation with the requested layout
855	AllocError {
856		/// requested layout
857		layout: core::alloc::Layout,
858	},
859}
860/// determines whether the input should be implicitly conjugated or not
861#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
862pub enum Conj {
863	/// no implicit conjugation
864	No,
865	/// implicit conjugation
866	Yes,
867}
868#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
869pub(crate) enum DiagStatus {
870	Unit,
871	Generic,
872}
873/// determines whether to replace or add to the result of a matmul operatio
874#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
875pub enum Accum {
876	/// overwrites the output buffer
877	Replace,
878	/// adds the result to the output buffer
879	Add,
880}
881/// determines which side of a self-adjoint matrix should be accessed
882#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
883pub enum Side {
884	/// lower triangular half
885	Lower,
886	/// upper triangular half
887	Upper,
888}
889impl Conj {
890	/// returns `self == Conj::Yes`
891	#[inline]
892	pub const fn is_conj(self) -> bool {
893		matches!(self, Conj::Yes)
894	}
895
896	/// returns the composition of `self` and `other`
897	#[inline]
898	pub const fn compose(self, other: Self) -> Self {
899		match (self, other) {
900			(Conj::No, Conj::No) => Conj::No,
901			(Conj::Yes, Conj::Yes) => Conj::No,
902			(Conj::No, Conj::Yes) => Conj::Yes,
903			(Conj::Yes, Conj::No) => Conj::Yes,
904		}
905	}
906
907	/// returns `Conj::No` if `T` is the canonical representation, otherwise
908	/// `Conj::Yes`
909	#[inline]
910	pub const fn get<T: Conjugate>() -> Self {
911		if T::IS_CANONICAL { Self::No } else { Self::Yes }
912	}
913
914	#[inline]
915	pub(crate) fn apply<T: Conjugate>(value: &T) -> T::Canonical {
916		let value = unsafe { &*(value as *const T as *const T::Canonical) };
917		if const { matches!(Self::get::<T>(), Conj::Yes) } {
918			T::Canonical::conj_impl(value)
919		} else {
920			T::Canonical::copy_impl(value)
921		}
922	}
923
924	#[inline]
925	pub(crate) fn apply_rt<T: ComplexField>(self, value: &T) -> T {
926		if self.is_conj() {
927			T::conj_impl(value)
928		} else {
929			T::copy_impl(value)
930		}
931	}
932}
933/// determines the parallelization configuration, and carries the right to execute with it
934///
935/// `Par` is neither `Copy` nor `Clone`: a `Team` owns exclusive authority over reserved
936/// threads, so a duplicate could hand the same threads to two concurrent children. Every
937/// function that runs work takes `Par` by value; a caller that needs another phase lends the
938/// parallelism with [`Par::reborrow`], and concurrent children are created by splitting
939/// (see `utils::thread::split`), never by copying. Queries that only need the size (workspace
940/// requirements, thresholds) use [`ParShape`] or [`Par::degree`].
941///
942/// ```
943/// use faer::{Par, ParShape};
944/// // sizes are plain values
945/// fn is_copy<T: Copy>() {}
946/// is_copy::<ParShape>();
947/// // phases lend, they do not copy
948/// fn phase(_: Par<'_>) {}
949/// let mut par = Par::Seq;
950/// phase(par.reborrow());
951/// phase(par);
952/// ```
953///
954/// Not `Copy`:
955///
956/// ```compile_fail
957/// fn is_copy<T: Copy>() {}
958/// is_copy::<faer::Par<'static>>();
959/// ```
960///
961/// Not `Clone`:
962///
963/// ```compile_fail
964/// fn is_clone<T: Clone>() {}
965/// is_clone::<faer::Par<'static>>();
966/// ```
967///
968/// A phase that takes the parallelism consumes it:
969///
970/// ```compile_fail
971/// fn phase(_: faer::Par<'_>) {}
972/// let par = faer::Par::Seq;
973/// phase(par);
974/// phase(par);
975/// ```
976///
977/// Concurrent children cannot both use the parent:
978///
979/// ```compile_fail
980/// fn phase(_: faer::Par<'_>) {}
981/// let par = faer::Par::Seq;
982/// faer::utils::thread::join_raw(|_| phase(par), |_| phase(par), faer::Par::Seq);
983/// ```
984#[derive(Debug)]
985pub enum Par<'team> {
986	/// sequential, non portable across different platforms
987	Seq,
988	/// parallelized using the global rayon threadpool, non portable across
989	/// different platforms
990	#[cfg(feature = "rayon")]
991	Rayon(NonZeroUsize),
992	/// experimental: parallelized on the reserved members of a threadteam team. The value is
993	/// the exclusive authority over them for as long as it lives; the number of participants is
994	/// the size of the team.
995	#[cfg(feature = "threadteam")]
996	Team(threadteam::TeamExec<'team>),
997	#[cfg(not(feature = "threadteam"))]
998	#[doc(hidden)]
999	__Lifetime(core::convert::Infallible, core::marker::PhantomData<&'team ()>),
1000}
1001/// the shape of a [`Par`]: which kind of parallelism and how many threads, without the
1002/// right to run anything. `Copy`, so workspace queries and thresholds can use it freely; it
1003/// cannot be passed to any function that executes work.
1004#[derive(Copy, Clone, Debug, PartialEq, Eq)]
1005pub enum ParShape {
1006	/// sequential
1007	Seq,
1008	/// the global rayon threadpool with this many threads
1009	#[cfg(feature = "rayon")]
1010	Rayon(NonZeroUsize),
1011	/// a threadteam team of this many participants, with the dispatch costs of the pool it
1012	/// belongs to (see [`ParShape::dispatch_cost`])
1013	#[cfg(feature = "threadteam")]
1014	Team(NonZeroUsize, threadteam::DispatchCosts),
1015}
1016impl ParShape {
1017	/// sequential
1018	#[inline]
1019	pub const fn seq() -> Self {
1020		Self::Seq
1021	}
1022
1023	/// `Rayon(nthreads)`, or `Rayon(rayon::current_num_threads())` when `nthreads` is zero
1024	#[inline]
1025	#[cfg(feature = "rayon")]
1026	pub fn rayon(nthreads: usize) -> Self {
1027		match Par::rayon(nthreads) {
1028			Par::Rayon(n) => Self::Rayon(n),
1029			_ => unreachable!(),
1030		}
1031	}
1032
1033	/// `Team(nthreads)` with the host default dispatch costs, zero treated as one
1034	#[inline]
1035	#[cfg(feature = "threadteam")]
1036	pub fn team(nthreads: usize) -> Self {
1037		Self::Team(
1038			NonZeroUsize::new(nthreads.max(1)).unwrap(),
1039			threadteam::DispatchCosts::host_default(),
1040		)
1041	}
1042
1043	/// the number of threads that should ideally execute an operation of this shape
1044	#[inline]
1045	pub fn degree(self) -> usize {
1046		match self {
1047			Self::Seq => 1,
1048			#[cfg(feature = "rayon")]
1049			Self::Rayon(n) => n.get(),
1050			#[cfg(feature = "threadteam")]
1051			Self::Team(n, _) => n.get(),
1052		}
1053	}
1054}
1055impl Par<'_> {
1056	/// lends the parallelism for one phase. The parent is exclusively borrowed and unusable
1057	/// while the lent value lives, so a phase sequence never holds the authority twice.
1058	#[inline]
1059	pub fn reborrow(&mut self) -> Par<'_> {
1060		match self {
1061			Par::Seq => Par::Seq,
1062			#[cfg(feature = "rayon")]
1063			Par::Rayon(n) => Par::Rayon(*n),
1064			#[cfg(feature = "threadteam")]
1065			Par::Team(team) => Par::Team(team.reborrow()),
1066			#[cfg(not(feature = "threadteam"))]
1067			Par::__Lifetime(never, _) => match *never {},
1068		}
1069	}
1070
1071	/// the shape of this parallelism, without the right to execute with it
1072	#[inline]
1073	pub fn shape(&self) -> ParShape {
1074		match self {
1075			Par::Seq => ParShape::Seq,
1076			#[cfg(feature = "rayon")]
1077			Par::Rayon(n) => ParShape::Rayon(*n),
1078			#[cfg(feature = "threadteam")]
1079			Par::Team(team) => ParShape::Team(team.size(), team.dispatch_costs()),
1080			#[cfg(not(feature = "threadteam"))]
1081			Par::__Lifetime(never, _) => match *never {},
1082		}
1083	}
1084
1085	/// returns `Par::Rayon(nthreads)` if `nthreads` is non-zero, or
1086	/// `Par::Rayon(rayon::current_num_threads())` otherwise
1087	#[inline]
1088	#[cfg(feature = "rayon")]
1089	pub fn rayon(nthreads: usize) -> Self {
1090		if nthreads == 0 {
1091			Self::Rayon(
1092				NonZeroUsize::new(rayon::current_num_threads()).unwrap(),
1093			)
1094		} else {
1095			Self::Rayon(NonZeroUsize::new(nthreads).unwrap())
1096		}
1097	}
1098
1099	/// the number of threads that should ideally execute an operation with the
1100	/// given parallelism
1101	#[inline]
1102	pub fn degree(&self) -> usize {
1103		utils::thread::parallelism_degree(self)
1104	}
1105}
1106#[allow(non_camel_case_types)]
1107/// `Complex<f32>`
1108pub type c32 = traits::c32;
1109#[allow(non_camel_case_types)]
1110/// `Complex<f64>`
1111pub type c64 = traits::c64;
1112#[allow(non_camel_case_types)]
1113/// `Complex<f64>`
1114pub type cx128 = traits::cx128;
1115#[allow(non_camel_case_types)]
1116/// `Complex<f64>`
1117pub type fx128 = traits::fx128;
1118pub use col::{Col, ColMut, ColRef};
1119pub use mat::{Mat, MatMut, MatRef};
1120pub use row::{Row, RowMut, RowRef};
1121#[allow(unused_imports, dead_code)]
1122mod internal_prelude {
1123	pub use crate::simd_iter;
1124
1125	pub trait DivCeil: Sized {
1126		fn msrv_div_ceil(self, rhs: Self) -> Self;
1127		fn msrv_next_multiple_of(self, rhs: Self) -> Self;
1128		fn msrv_checked_next_multiple_of(self, rhs: Self) -> Option<Self>;
1129	}
1130	impl DivCeil for usize {
1131		#[inline]
1132		fn msrv_div_ceil(self, rhs: Self) -> Self {
1133			let d = self / rhs;
1134			let r = self % rhs;
1135			if r > 0 { d + 1 } else { d }
1136		}
1137
1138		#[inline]
1139		fn msrv_next_multiple_of(self, rhs: Self) -> Self {
1140			match self % rhs {
1141				0 => self,
1142				r => self + (rhs - r),
1143			}
1144		}
1145
1146		#[inline]
1147		fn msrv_checked_next_multiple_of(self, rhs: Self) -> Option<Self> {
1148			{
1149				match self.checked_rem(rhs)? {
1150					0 => Some(self),
1151					r => self.checked_add(rhs - r),
1152				}
1153			}
1154		}
1155	}
1156	pub(crate) use crate::col::{Col, ColMut, ColRef};
1157	pub(crate) use crate::diag::{Diag, DiagMut, DiagRef};
1158	pub(crate) use crate::hacks::transmute;
1159	pub(crate) use crate::linalg::{
1160		self, temp_mat_scratch, temp_mat_uninit, temp_mat_zeroed,
1161	};
1162	pub(crate) use crate::mat::{
1163		AsMat, AsMatMut, AsMatRef, Mat, MatMut, MatRef,
1164	};
1165	pub(crate) use crate::perm::{Perm, PermRef};
1166	pub(crate) use crate::prelude::*;
1167	pub(crate) use crate::row::{AsRowMut, AsRowRef, Row, RowMut, RowRef};
1168	pub(crate) use crate::utils::bound::{Array, Dim, Idx, IdxInc, MaybeIdx};
1169	pub(crate) use crate::utils::simd::SimdCtx;
1170	pub(crate) use crate::{Auto, NonExhaustive, Side, Spec};
1171	pub use faer_traits::ext::*;
1172	pub use faer_traits::math_utils::*;
1173	pub use faer_traits::{
1174		ComplexField, ComplexImpl, ComplexImplConj, Conjugate, Index,
1175		IndexCore, Real, RealField, SignedIndex, SimdArch, Symbolic,
1176	};
1177	pub use num_complex::Complex;
1178	#[cfg(test)]
1179	pub(crate) use std::dbg;
1180	#[cfg(test)]
1181	pub(crate) use {alloc::boxed::Box, alloc::vec, alloc::vec::Vec};
1182	#[inline]
1183	pub fn simd_align(i: usize) -> usize {
1184		i.wrapping_neg()
1185	}
1186	pub(crate) use self::{unzip as uz, zip as z};
1187	pub(crate) use crate::{
1188		Accum, Conj, ContiguousBwd, ContiguousFwd, DiagStatus, Par, ParShape,
1189		Shape, Stride, Unbind, make_guard, unzip, zip,
1190	};
1191	pub use dyn_stack::{MemStack, StackReq, alloc as alloca};
1192	pub use equator::{
1193		assert, assert as Assert, debug_assert, debug_assert as DebugAssert,
1194	};
1195	pub use reborrow::*;
1196}
1197#[allow(unused_imports)]
1198pub(crate) mod internal_prelude_sp {
1199	pub(crate) use crate::internal_prelude::*;
1200	pub(crate) use crate::sparse::{
1201		FaerError, NONE, Pair, SparseColMat, SparseColMatMut, SparseColMatRef,
1202		SparseRowMat, SparseRowMatMut, SparseRowMatRef, SymbolicSparseColMat,
1203		SymbolicSparseColMatRef, SymbolicSparseRowMat, SymbolicSparseRowMatRef,
1204		Triplet, csc_numeric, csc_symbolic, csr_numeric, csr_symbolic,
1205		linalg as linalg_sp, try_collect, try_zeroed, windows2,
1206	};
1207	pub(crate) use core::cell::Cell;
1208	pub(crate) use core::iter;
1209	pub(crate) use dyn_stack::MemBuffer;
1210}
1211/// useful imports for general usage of the library
1212pub mod prelude {
1213	#[cfg(feature = "linalg")]
1214	pub use super::linalg::solvers::{DenseSolve, Solve, SolveLstsq};
1215	#[cfg(feature = "sparse")]
1216	pub use super::prelude_sp::*;
1217	pub use super::{Par, ParShape, Scale, c32, c64, col, mat, row, unzip, zip};
1218	pub use col::{Col, ColMut, ColRef};
1219	pub use mat::{Mat, MatMut, MatRef};
1220	pub use reborrow::{IntoConst, Reborrow, ReborrowMut};
1221	pub use row::{Row, RowMut, RowRef};
1222	/// see [`Default`]
1223	#[inline]
1224	pub fn default<T: Default>() -> T {
1225		Default::default()
1226	}
1227}
1228#[cfg(feature = "sparse")]
1229mod prelude_sp {
1230	use crate::sparse;
1231	pub use sparse::{
1232		SparseColMat, SparseColMatMut, SparseColMatRef, SparseRowMat,
1233		SparseRowMatMut, SparseRowMatRef,
1234	};
1235}
1236/// scaling factor for multiplying matrices.
1237#[derive(Copy, Clone, Debug)]
1238#[repr(transparent)]
1239pub struct Scale<T>(pub T);
1240impl<T> Scale<T> {
1241	/// create a reference to a scaling factor from a reference to a value.
1242	#[inline(always)]
1243	pub fn from_ref(value: &T) -> &Self {
1244		unsafe { &*(value as *const T as *const Self) }
1245	}
1246
1247	/// create a mutable reference to a scaling factor from a mutable reference
1248	/// to a value.
1249	#[inline(always)]
1250	pub fn from_mut(value: &mut T) -> &mut Self {
1251		unsafe { &mut *(value as *mut T as *mut Self) }
1252	}
1253}
1254#[doc(hidden)]
1255pub mod hacks;
1256/// statistics and randomness functionality
1257pub mod stats;
1258mod non_exhaustive {
1259	#[doc(hidden)]
1260	#[derive(Debug, Copy, Clone, PartialEq, Eq)]
1261	#[repr(transparent)]
1262	pub struct NonExhaustive(pub(crate) ());
1263}
1264pub(crate) use non_exhaustive::NonExhaustive;
1265/// like `Default`, but with an extra type parameter so that algorithm
1266/// hyperparameters can be tuned per scalar type.
1267pub trait Auto<T> {
1268	/// returns the default value for the type `T`
1269	fn auto() -> Self;
1270}
1271/// implements [`Default`] based on `Config`'s [`Auto`] implementation for the
1272/// type `T`.
1273pub struct Spec<Config, T> {
1274	/// wrapped config value
1275	pub config: Config,
1276	__marker: core::marker::PhantomData<fn() -> T>,
1277}
1278impl<Config, T> core::ops::Deref for Spec<Config, T> {
1279	type Target = Config;
1280
1281	#[inline]
1282	fn deref(&self) -> &Self::Target {
1283		&self.config
1284	}
1285}
1286impl<Config, T> core::ops::DerefMut for Spec<Config, T> {
1287	#[inline]
1288	fn deref_mut(&mut self) -> &mut Self::Target {
1289		&mut self.config
1290	}
1291}
1292impl<Config: Copy, T> Copy for Spec<Config, T> {}
1293impl<Config: Clone, T> Clone for Spec<Config, T> {
1294	#[inline]
1295	fn clone(&self) -> Self {
1296		Self::new(self.config.clone())
1297	}
1298}
1299impl<Config: core::fmt::Debug, T> core::fmt::Debug for Spec<Config, T> {
1300	#[inline]
1301	fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
1302		self.config.fmt(f)
1303	}
1304}
1305impl<Config, T> Spec<Config, T> {
1306	/// wraps the given config value
1307	#[inline]
1308	pub fn new(config: Config) -> Self {
1309		Spec {
1310			config,
1311			__marker: core::marker::PhantomData,
1312		}
1313	}
1314}
1315impl<T, Config> From<Config> for Spec<Config, T> {
1316	#[inline]
1317	fn from(config: Config) -> Self {
1318		Spec {
1319			config,
1320			__marker: core::marker::PhantomData,
1321		}
1322	}
1323}
1324impl<T, Config: Auto<T>> Default for Spec<Config, T> {
1325	#[inline]
1326	fn default() -> Self {
1327		Spec {
1328			config: Auto::<T>::auto(),
1329			__marker: core::marker::PhantomData,
1330		}
1331	}
1332}
1333mod into_range {
1334	use super::*;
1335	pub trait IntoRange<I> {
1336		type Len<N: Shape>: Shape;
1337		fn into_range(self, min: I, max: I) -> core::ops::Range<I>;
1338	}
1339	impl<I> IntoRange<I> for core::ops::Range<I> {
1340		type Len<N: Shape> = usize;
1341
1342		#[inline]
1343		fn into_range(self, _: I, _: I) -> core::ops::Range<I> {
1344			self
1345		}
1346	}
1347	impl<I> IntoRange<I> for core::ops::RangeFrom<I> {
1348		type Len<N: Shape> = usize;
1349
1350		#[inline]
1351		fn into_range(self, _: I, max: I) -> core::ops::Range<I> {
1352			self.start..max
1353		}
1354	}
1355	impl<I> IntoRange<I> for core::ops::RangeTo<I> {
1356		type Len<N: Shape> = usize;
1357
1358		#[inline]
1359		fn into_range(self, min: I, _: I) -> core::ops::Range<I> {
1360			min..self.end
1361		}
1362	}
1363	impl<I> IntoRange<I> for core::ops::RangeFull {
1364		type Len<N: Shape> = N;
1365
1366		#[inline]
1367		fn into_range(self, min: I, max: I) -> core::ops::Range<I> {
1368			min..max
1369		}
1370	}
1371}
1372mod sort;
1373pub extern crate dyn_stack;
1374pub extern crate faer_traits as traits;
1375#[doc(hidden)]
1376pub extern crate generativity;
1377pub extern crate num_complex as complex;
1378#[cfg(feature = "rand")]
1379#[cfg_attr(docs_rs, doc(cfg(feature = "rand")))]
1380pub extern crate rand;
1381pub extern crate reborrow;
1382
1383extern crate self as faer;
1384
1385#[macro_export]
1386#[doc(hidden)]
1387macro_rules! simd_iter {
1388	(for ($batch_id:tt, $i:pat $(,)?) in [$indices:expr; $batch_size:expr] $b:block $(,)?) => {#[allow(non_upper_case_globals, unconditional_panic)] 'out:{
1389		let (__head__, __body_batch__, __body_item__, __tail__) = &$indices;
1390		let __body_batch__ = __body_batch__.clone();
1391		let mut __body_item__ = __body_item__.clone();
1392		const {
1393			::core::assert!($batch_size > 0);
1394			::core::assert!($batch_size <= 8);
1395		}
1396		{
1397		if const { $batch_size == 1 } {if $batch_size == 1 { const $batch_id: usize = 0; if let &Some($i) = __head__ $b; }}
1398		if const { $batch_size == 2 } {if $batch_size == 2 { const $batch_id: usize = 1; if let &Some($i) = __head__ $b; }}
1399		if const { $batch_size == 3 } {if $batch_size == 3 { const $batch_id: usize = 2; if let &Some($i) = __head__ $b; }}
1400		if const { $batch_size == 4 } {if $batch_size == 4 { const $batch_id: usize = 3; if let &Some($i) = __head__ $b; }}
1401		if const { $batch_size == 5 } {if $batch_size == 5 { const $batch_id: usize = 4; if let &Some($i) = __head__ $b; }}
1402		if const { $batch_size == 6 } {if $batch_size == 6 { const $batch_id: usize = 5; if let &Some($i) = __head__ $b; }}
1403		if const { $batch_size == 7 } {if $batch_size == 7 { const $batch_id: usize = 6; if let &Some($i) = __head__ $b; }}
1404		if const { $batch_size == 8 } {if $batch_size == 8 { const $batch_id: usize = 7; if let &Some($i) = __head__ $b; }}
1405		}
1406		for __i__ in __body_batch__ {
1407			let __i__: [_; $batch_size] = __i__;
1408			let __i__ = __i__.as_slice();
1409			::core::assert!(__i__.len() == $batch_size);
1410
1411			{
1412			if const { $batch_size > 0 } {if __i__.len() > 0 { const $batch_id: usize = 0; let $i = __i__[0]; $b }}
1413			if const { $batch_size > 1 } {if __i__.len() > 1 { const $batch_id: usize = 1; let $i = __i__[1]; $b }}
1414			if const { $batch_size > 2 } {if __i__.len() > 2 { const $batch_id: usize = 2; let $i = __i__[2]; $b }}
1415			if const { $batch_size > 3 } {if __i__.len() > 3 { const $batch_id: usize = 3; let $i = __i__[3]; $b }}
1416			if const { $batch_size > 4 } {if __i__.len() > 4 { const $batch_id: usize = 4; let $i = __i__[4]; $b }}
1417			if const { $batch_size > 5 } {if __i__.len() > 5 { const $batch_id: usize = 5; let $i = __i__[5]; $b }}
1418			if const { $batch_size > 6 } {if __i__.len() > 6 { const $batch_id: usize = 6; let $i = __i__[6]; $b }}
1419			if const { $batch_size > 7 } {if __i__.len() > 7 { const $batch_id: usize = 7; let $i = __i__[7]; $b }}
1420			}
1421		};
1422
1423		{
1424		if const { $batch_size > 0 } { if $batch_size > 0 { const $batch_id: usize = 0; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1425		if const { $batch_size > 1 } { if $batch_size > 1 { const $batch_id: usize = 1; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1426		if const { $batch_size > 2 } { if $batch_size > 2 { const $batch_id: usize = 2; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1427		if const { $batch_size > 3 } { if $batch_size > 3 { const $batch_id: usize = 3; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1428		if const { $batch_size > 4 } { if $batch_size > 4 { const $batch_id: usize = 4; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1429		if const { $batch_size > 5 } { if $batch_size > 5 { const $batch_id: usize = 5; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1430		if const { $batch_size > 6 } { if $batch_size > 6 { const $batch_id: usize = 6; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1431		if const { $batch_size > 7 } { if $batch_size > 7 { const $batch_id: usize = 7; if let Some($i) = __body_item__.next() $b else { if let &Some($i) = __tail__ $b; break 'out; } } }
1432		}
1433	}};
1434	(for $i:pat in [$indices:expr] $b:block $(,)?) => {#[allow(non_upper_case_globals)]{
1435		let (__head__, __body__, __tail__) = &$indices;
1436		if let &Some($i) = __head__ $b;
1437		for $i in __body__.clone() $b;
1438		if let &Some($i) = __tail__ $b;
1439	}};
1440}