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