faer/linalg/matmul/
triangular.rs

1use super::*;
2use crate::linalg::temp_mat_uninit;
3use crate::linalg::zip::Diag;
4use crate::mat::{AsMatMut, MatMut, MatRef};
5use crate::utils::thread::join_raw;
6use crate::{assert, debug_assert};
7
8#[repr(u8)]
9#[derive(Copy, Clone, Debug)]
10pub(crate) enum DiagonalKind {
11	Zero,
12	Unit,
13	Generic,
14}
15#[inline]
16fn pointer_offset<T>(ptr: *const T) -> usize {
17	if const {
18		core::mem::size_of::<T>().is_power_of_two()
19			&& core::mem::size_of::<T>() <= 64
20	} {
21		ptr.align_offset(64).wrapping_neg() % 16
22	} else {
23		0
24	}
25}
26fn copy_lower<'N, T: ComplexField>(
27	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
28	src: MatRef<'_, T, Dim<'N>, Dim<'N>>,
29	src_diag: DiagonalKind,
30) {
31	let N = dst.nrows();
32	let mut dst = dst;
33	match src_diag {
34		DiagonalKind::Zero => {
35			dst.copy_from_strict_triangular_lower(src);
36			for j in N.indices() {
37				let zero = zero();
38				dst[(j, j)] = zero;
39			}
40		},
41		DiagonalKind::Unit => {
42			dst.copy_from_strict_triangular_lower(src);
43			for j in N.indices() {
44				let one = one();
45				dst[(j, j)] = one;
46			}
47		},
48		DiagonalKind::Generic => dst.copy_from_triangular_lower(src),
49	}
50	zip!(dst)
51		.for_each_triangular_upper(Diag::Skip, |unzip!(dst)| *dst = zero());
52}
53fn accum_lower<'N, T: ComplexField>(
54	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
55	src: MatRef<'_, T, Dim<'N>, Dim<'N>>,
56	skip_diag: bool,
57	beta: Accum,
58) {
59	let N = dst.nrows();
60	debug_assert!(N == dst.nrows());
61	debug_assert!(N == dst.ncols());
62	debug_assert!(N == src.nrows());
63	debug_assert!(N == src.ncols());
64	match beta {
65		Accum::Add => {
66			zip!(dst, src).for_each_triangular_lower(
67				if skip_diag { Diag::Skip } else { Diag::Include },
68				|unzip!(dst, src)| *dst += src,
69			);
70		},
71		Accum::Replace => {
72			zip!(dst, src).for_each_triangular_lower(
73				if skip_diag { Diag::Skip } else { Diag::Include },
74				|unzip!(dst, src)| *dst = src.copy(),
75			);
76		},
77	}
78}
79fn copy_upper<'N, T: ComplexField>(
80	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
81	src: MatRef<'_, T, Dim<'N>, Dim<'N>>,
82	src_diag: DiagonalKind,
83) {
84	copy_lower(dst.transpose_mut(), src.transpose(), src_diag)
85}
86#[repr(align(64))]
87struct Storage<T>([T; 32 * 16]);
88macro_rules! stack_mat_16x16 {
89	($name:ident, $n:expr, $offset:expr, $rs:expr, $cs:expr, $T:ty $(,)?) => {
90		let mut __tmp = core::mem::MaybeUninit::<Storage<$T>>::uninit();
91		let __stack = MemStack::new_any(core::slice::from_mut(&mut __tmp));
92		let mut $name = unsafe { temp_mat_uninit(32, $n, __stack) }.0;
93		let mut $name = $name.as_mat_mut().subrows_mut($offset, $n);
94		if $cs.unsigned_abs() == 1 {
95			$name = $name.transpose_mut();
96			if $cs == 1 {
97				$name = $name.transpose_mut().reverse_cols_mut();
98			}
99		} else if $rs == -1 {
100			$name = $name.reverse_rows_mut();
101		}
102	};
103}
104fn mat_x_lower_impl_unchecked<'M, 'N, T: ComplexField>(
105	dst: MatMut<'_, T, Dim<'M>, Dim<'N>>,
106	beta: Accum,
107	lhs: MatRef<'_, T, Dim<'M>, Dim<'N>>,
108	rhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
109	rhs_diag: DiagonalKind,
110	alpha: &T,
111	conj_lhs: Conj,
112	conj_rhs: Conj,
113	mut par: Par,
114) {
115	let N = rhs.nrows();
116	let M = lhs.nrows();
117	let n = N.unbound();
118	let m = M.unbound();
119	debug_assert!(M == lhs.nrows());
120	debug_assert!(N == lhs.ncols());
121	debug_assert!(N == rhs.nrows());
122	debug_assert!(N == rhs.ncols());
123	debug_assert!(M == dst.nrows());
124	debug_assert!(N == dst.ncols());
125	// One section per level (the join of the two halves); the unit is a multiply-add of the
126	// triangular product, `m*n*n/2`, at the GEMM rate.
127	let join_parallelism = if crate::utils::cost::join_pays::<T>(
128		&par,
129		m.saturating_mul(n).saturating_mul(n) / 2,
130		crate::utils::cost::rate::gemm_mac(),
131		1,
132	) {
133		par.reborrow()
134	} else {
135		Par::Seq
136	};
137	if n <= 16 {
138		let op = {
139			#[inline(never)]
140			|| {
141				stack_mat_16x16!(
142					temp_rhs,
143					N,
144					pointer_offset(rhs.as_ptr()),
145					rhs.row_stride(),
146					rhs.col_stride(),
147					T
148				);
149				copy_lower(temp_rhs.rb_mut(), rhs, rhs_diag);
150				let mut dst = dst;
151				super::matmul_with_conj(
152					dst.rb_mut(),
153					beta,
154					lhs,
155					conj_lhs,
156					temp_rhs.rb(),
157					conj_rhs,
158					alpha.clone(),
159					par,
160				);
161			}
162		};
163		op();
164	} else {
165		make_guard!(HEAD);
166		make_guard!(TAIL);
167		let bs = N.partition(N.checked_idx_inc(N.unbound() / 2), HEAD, TAIL);
168		let (rhs_top_left, _, rhs_bot_left, rhs_bot_right) =
169			rhs.split_with(bs, bs);
170		let (lhs_left, lhs_right) = lhs.split_cols_with(bs);
171		let (mut dst_left, mut dst_right) = dst.split_cols_with_mut(bs);
172		{
173			join_raw(
174				|par| {
175					mat_x_lower_impl_unchecked(
176						dst_left.rb_mut(),
177						beta,
178						lhs_left,
179						rhs_top_left,
180						rhs_diag,
181						alpha,
182						conj_lhs,
183						conj_rhs,
184						par,
185					)
186				},
187				|par| {
188					mat_x_lower_impl_unchecked(
189						dst_right.rb_mut(),
190						beta,
191						lhs_right,
192						rhs_bot_right,
193						rhs_diag,
194						alpha,
195						conj_lhs,
196						conj_rhs,
197						par,
198					)
199				},
200				join_parallelism,
201			)
202		};
203		super::matmul_with_conj(
204			dst_left,
205			Accum::Add,
206			lhs_right,
207			conj_lhs,
208			rhs_bot_left,
209			conj_rhs,
210			alpha.clone(),
211			par,
212		);
213	}
214}
215fn lower_x_lower_into_lower_impl_unchecked<'N, T: ComplexField>(
216	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
217	beta: Accum,
218	skip_diag: bool,
219	lhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
220	lhs_diag: DiagonalKind,
221	rhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
222	rhs_diag: DiagonalKind,
223	alpha: &T,
224	conj_lhs: Conj,
225	conj_rhs: Conj,
226	mut par: Par,
227) {
228	let N = dst.nrows();
229	let n = N.unbound();
230	debug_assert!(N == lhs.nrows());
231	debug_assert!(N == lhs.ncols());
232	debug_assert!(N == rhs.nrows());
233	debug_assert!(N == rhs.ncols());
234	debug_assert!(N == dst.nrows());
235	debug_assert!(N == dst.ncols());
236	if n <= 16 {
237		let op = {
238			#[inline(never)]
239			|| {
240				stack_mat_16x16!(
241					temp_dst,
242					N,
243					pointer_offset(dst.as_ptr()),
244					dst.row_stride(),
245					dst.col_stride(),
246					T
247				);
248				stack_mat_16x16!(
249					temp_lhs,
250					N,
251					pointer_offset(lhs.as_ptr()),
252					lhs.row_stride(),
253					lhs.col_stride(),
254					T
255				);
256				stack_mat_16x16!(
257					temp_rhs,
258					N,
259					pointer_offset(rhs.as_ptr()),
260					rhs.row_stride(),
261					rhs.col_stride(),
262					T
263				);
264				copy_lower(temp_lhs.rb_mut(), lhs, lhs_diag);
265				copy_lower(temp_rhs.rb_mut(), rhs, rhs_diag);
266				super::matmul_with_conj(
267					temp_dst.rb_mut(),
268					Accum::Replace,
269					temp_lhs.rb(),
270					conj_lhs,
271					temp_rhs.rb(),
272					conj_rhs,
273					alpha.clone(),
274					par,
275				);
276				accum_lower(dst, temp_dst.rb(), skip_diag, beta);
277			}
278		};
279		op();
280	} else {
281		make_guard!(HEAD);
282		make_guard!(TAIL);
283		let bs = N.partition(N.checked_idx_inc(N.unbound() / 2), HEAD, TAIL);
284		let (dst_top_left, _, mut dst_bot_left, dst_bot_right) =
285			dst.split_with_mut(bs, bs);
286		let (lhs_top_left, _, lhs_bot_left, lhs_bot_right) =
287			lhs.split_with(bs, bs);
288		let (rhs_top_left, _, rhs_bot_left, rhs_bot_right) =
289			rhs.split_with(bs, bs);
290		lower_x_lower_into_lower_impl_unchecked(
291			dst_top_left,
292			beta,
293			skip_diag,
294			lhs_top_left,
295			lhs_diag,
296			rhs_top_left,
297			rhs_diag,
298			alpha,
299			conj_lhs,
300			conj_rhs,
301			par.reborrow(),
302		);
303		mat_x_lower_impl_unchecked(
304			dst_bot_left.rb_mut(),
305			beta,
306			lhs_bot_left,
307			rhs_top_left,
308			rhs_diag,
309			alpha,
310			conj_lhs,
311			conj_rhs,
312			par.reborrow(),
313		);
314		mat_x_lower_impl_unchecked(
315			dst_bot_left.reverse_rows_and_cols_mut().transpose_mut(),
316			Accum::Add,
317			rhs_bot_left.reverse_rows_and_cols().transpose(),
318			lhs_bot_right.reverse_rows_and_cols().transpose(),
319			lhs_diag,
320			alpha,
321			conj_rhs,
322			conj_lhs,
323			par.reborrow(),
324		);
325		lower_x_lower_into_lower_impl_unchecked(
326			dst_bot_right,
327			beta,
328			skip_diag,
329			lhs_bot_right,
330			lhs_diag,
331			rhs_bot_right,
332			rhs_diag,
333			alpha,
334			conj_lhs,
335			conj_rhs,
336			par,
337		)
338	}
339}
340fn upper_x_lower_impl_unchecked<'N, T: ComplexField>(
341	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
342	beta: Accum,
343	lhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
344	lhs_diag: DiagonalKind,
345	rhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
346	rhs_diag: DiagonalKind,
347	alpha: &T,
348	conj_lhs: Conj,
349	conj_rhs: Conj,
350	par: Par,
351) {
352	let N = dst.nrows();
353	let n = N.unbound();
354	debug_assert!(N == lhs.nrows());
355	debug_assert!(N == lhs.ncols());
356	debug_assert!(N == rhs.nrows());
357	debug_assert!(N == rhs.ncols());
358	debug_assert!(N == dst.nrows());
359	debug_assert!(N == dst.ncols());
360	if n <= 16 {
361		let op = {
362			#[inline(never)]
363			|| {
364				stack_mat_16x16!(
365					temp_lhs,
366					N,
367					pointer_offset(lhs.as_ptr()),
368					lhs.row_stride(),
369					lhs.col_stride(),
370					T
371				);
372				stack_mat_16x16!(
373					temp_rhs,
374					N,
375					pointer_offset(rhs.as_ptr()),
376					rhs.row_stride(),
377					rhs.col_stride(),
378					T
379				);
380				copy_upper(temp_lhs.rb_mut(), lhs, lhs_diag);
381				copy_lower(temp_rhs.rb_mut(), rhs, rhs_diag);
382				super::matmul_with_conj(
383					dst,
384					beta,
385					temp_lhs.rb(),
386					conj_lhs,
387					temp_rhs.rb(),
388					conj_rhs,
389					alpha.clone(),
390					par,
391				);
392			}
393		};
394		op();
395	} else {
396		make_guard!(HEAD);
397		make_guard!(TAIL);
398		let bs = N.partition(N.checked_idx_inc(N.unbound() / 2), HEAD, TAIL);
399		let (mut dst_top_left, dst_top_right, dst_bot_left, dst_bot_right) =
400			dst.split_with_mut(bs, bs);
401		let (lhs_top_left, lhs_top_right, _, lhs_bot_right) =
402			lhs.split_with(bs, bs);
403		let (rhs_top_left, _, rhs_bot_left, rhs_bot_right) =
404			rhs.split_with(bs, bs);
405		join_raw(
406			|mut par| {
407				super::matmul_with_conj(
408					dst_top_left.rb_mut(),
409					beta,
410					lhs_top_right,
411					conj_lhs,
412					rhs_bot_left,
413					conj_rhs,
414					alpha.clone(),
415					par.reborrow(),
416				);
417				upper_x_lower_impl_unchecked(
418					dst_top_left,
419					Accum::Add,
420					lhs_top_left,
421					lhs_diag,
422					rhs_top_left,
423					rhs_diag,
424					alpha,
425					conj_lhs,
426					conj_rhs,
427					par,
428				)
429			},
430			|mut par| {
431				join_raw(
432					|par| {
433						mat_x_lower_impl_unchecked(
434							dst_top_right,
435							beta,
436							lhs_top_right,
437							rhs_bot_right,
438							rhs_diag,
439							alpha,
440							conj_lhs,
441							conj_rhs,
442							par,
443						)
444					},
445					|par| {
446						mat_x_lower_impl_unchecked(
447							dst_bot_left.transpose_mut(),
448							beta,
449							rhs_bot_left.transpose(),
450							lhs_bot_right.transpose(),
451							lhs_diag,
452							alpha,
453							conj_rhs,
454							conj_lhs,
455							par,
456						)
457					},
458					par.reborrow(),
459				);
460				upper_x_lower_impl_unchecked(
461					dst_bot_right,
462					beta,
463					lhs_bot_right,
464					lhs_diag,
465					rhs_bot_right,
466					rhs_diag,
467					alpha,
468					conj_lhs,
469					conj_rhs,
470					par,
471				)
472			},
473			par,
474		);
475	}
476}
477fn upper_x_lower_into_lower_impl_unchecked<'N, T: ComplexField>(
478	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
479	beta: Accum,
480	skip_diag: bool,
481	lhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
482	lhs_diag: DiagonalKind,
483	rhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
484	rhs_diag: DiagonalKind,
485	alpha: &T,
486	conj_lhs: Conj,
487	conj_rhs: Conj,
488	par: Par,
489) {
490	let N = dst.nrows();
491	let n = N.unbound();
492	debug_assert!(N == lhs.nrows());
493	debug_assert!(N == lhs.ncols());
494	debug_assert!(N == rhs.nrows());
495	debug_assert!(N == rhs.ncols());
496	debug_assert!(N == dst.nrows());
497	debug_assert!(N == dst.ncols());
498	if n <= 16 {
499		let op = {
500			#[inline(never)]
501			|| {
502				stack_mat_16x16!(
503					temp_dst,
504					N,
505					pointer_offset(dst.as_ptr()),
506					dst.row_stride(),
507					dst.col_stride(),
508					T
509				);
510				stack_mat_16x16!(
511					temp_lhs,
512					N,
513					pointer_offset(lhs.as_ptr()),
514					lhs.row_stride(),
515					lhs.col_stride(),
516					T
517				);
518				stack_mat_16x16!(
519					temp_rhs,
520					N,
521					pointer_offset(rhs.as_ptr()),
522					rhs.row_stride(),
523					rhs.col_stride(),
524					T
525				);
526				copy_upper(temp_lhs.rb_mut(), lhs, lhs_diag);
527				copy_lower(temp_rhs.rb_mut(), rhs, rhs_diag);
528				super::matmul_with_conj(
529					temp_dst.rb_mut(),
530					Accum::Replace,
531					temp_lhs.rb(),
532					conj_lhs,
533					temp_rhs.rb(),
534					conj_rhs,
535					alpha.clone(),
536					par,
537				);
538				accum_lower(dst, temp_dst.rb(), skip_diag, beta);
539			}
540		};
541		op();
542	} else {
543		make_guard!(HEAD);
544		make_guard!(TAIL);
545		let bs = N.partition(N.checked_idx_inc(N.unbound() / 2), HEAD, TAIL);
546		let (mut dst_top_left, _, dst_bot_left, dst_bot_right) =
547			dst.split_with_mut(bs, bs);
548		let (lhs_top_left, lhs_top_right, _, lhs_bot_right) =
549			lhs.split_with(bs, bs);
550		let (rhs_top_left, _, rhs_bot_left, rhs_bot_right) =
551			rhs.split_with(bs, bs);
552		join_raw(
553			|mut par| {
554				mat_x_mat_into_lower_impl_unchecked(
555					dst_top_left.rb_mut(),
556					beta,
557					skip_diag,
558					lhs_top_right,
559					rhs_bot_left,
560					alpha,
561					conj_lhs,
562					conj_rhs,
563					par.reborrow(),
564				);
565				upper_x_lower_into_lower_impl_unchecked(
566					dst_top_left,
567					Accum::Add,
568					skip_diag,
569					lhs_top_left,
570					lhs_diag,
571					rhs_top_left,
572					rhs_diag,
573					alpha,
574					conj_lhs,
575					conj_rhs,
576					par,
577				)
578			},
579			|mut par| {
580				mat_x_lower_impl_unchecked(
581					dst_bot_left.transpose_mut(),
582					beta,
583					rhs_bot_left.transpose(),
584					lhs_bot_right.transpose(),
585					lhs_diag,
586					alpha,
587					conj_rhs,
588					conj_lhs,
589					par.reborrow(),
590				);
591				upper_x_lower_into_lower_impl_unchecked(
592					dst_bot_right,
593					beta,
594					skip_diag,
595					lhs_bot_right,
596					lhs_diag,
597					rhs_bot_right,
598					rhs_diag,
599					alpha,
600					conj_lhs,
601					conj_rhs,
602					par,
603				)
604			},
605			par,
606		);
607	}
608}
609fn mat_x_mat_into_lower_impl_unchecked<'N, 'K, T: ComplexField>(
610	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
611	beta: Accum,
612	skip_diag: bool,
613	lhs: MatRef<'_, T, Dim<'N>, Dim<'K>>,
614	rhs: MatRef<'_, T, Dim<'K>, Dim<'N>>,
615	alpha: &T,
616	conj_lhs: Conj,
617	conj_rhs: Conj,
618	par: Par,
619) {
620	// Two sections per level (the outer join and the join nested in its second half); the
621	// unit is a multiply-add of the product into one triangle, `n*n*k/2`, at the GEMM rate.
622	let pays = crate::utils::cost::join_pays::<T>(
623		&par,
624		dst.nrows()
625			.unbound()
626			.saturating_mul(dst.nrows().unbound())
627			.saturating_mul(lhs.ncols().unbound())
628			/ 2,
629		crate::utils::cost::rate::gemm_mac(),
630		2,
631	);
632	// On x86_64 the `private_gemm_x86` kernel does the product into the triangle. A `Team` token
633	// runs it on its own members (`Fork` of the team), gated by the kernel's two sections at the
634	// width `width_for` chooses; Rayon and Seq hand it the degree. Without the kernel (other
635	// targets, other dtypes) a team that pays runs the recursion below on the team.
636	#[cfg(all(target_arch = "x86_64", feature = "std"))]
637	if const {
638		T::IS_NATIVE_F64
639			|| T::IS_NATIVE_F32
640			|| T::IS_NATIVE_C64
641			|| T::IS_NATIVE_C32
642	} {
643		use private_gemm_x86::*;
644		let feat = if std::arch::is_x86_feature_detected!("avx512f") {
645			Some(InstrSet::Avx512)
646		} else if std::arch::is_x86_feature_detected!("avx2")
647			&& std::arch::is_x86_feature_detected!("fma")
648		{
649			Some(InstrSet::Avx256)
650		} else {
651			None
652		};
653		if *dst.nrows() > 0 && *dst.ncols() > 0 && *lhs.ncols() > 0 {
654			if let Some(feat) = feat {
655				unsafe {
656					let (dst, lhs) = if skip_diag {
657						(
658							dst.as_dyn_mut().get_mut(1.., ..),
659							lhs.as_dyn().get(1.., ..),
660						)
661					} else {
662						(dst.as_dyn_mut(), lhs.as_dyn())
663					};
664					macro_rules! go {
665						($call:path, $last:expr) => {
666							$call(
667								const {
668									if T::IS_NATIVE_F64 {
669										DType::F64
670									} else if T::IS_NATIVE_F32 {
671										DType::F32
672									} else if T::IS_NATIVE_C64 {
673										DType::C64
674									} else {
675										DType::C32
676									}
677								},
678								const { IType::U32 },
679								feat,
680								dst.nrows(),
681								dst.ncols(),
682								lhs.ncols(),
683								dst.as_ptr_mut() as _,
684								dst.row_stride(),
685								dst.col_stride(),
686								core::ptr::null(),
687								core::ptr::null(),
688								DstKind::Lower,
689								match beta {
690									crate::Accum::Replace => Accum::Replace,
691									crate::Accum::Add => Accum::Add,
692								},
693								lhs.as_ptr() as _,
694								lhs.row_stride(),
695								lhs.col_stride(),
696								conj_lhs == Conj::Yes,
697								core::ptr::null(),
698								0,
699								rhs.as_ptr() as _,
700								rhs.row_stride(),
701								rhs.col_stride(),
702								conj_rhs == Conj::Yes,
703								alpha as *const T as *const (),
704								$last,
705							)
706						};
707					}
708					#[cfg(feature = "threadteam")]
709					if let Par::Team(team) = par {
710						let estimate = crate::utils::cost::serial_time(
711							dst.nrows()
712								.saturating_mul(dst.ncols())
713								.saturating_mul(lhs.ncols())
714								/ 2,
715							crate::utils::cost::rate::gemm_mac().of::<T>(),
716						);
717						match crate::utils::team::kernel::kernel_width(
718							&team, estimate,
719						) {
720							Some(width) => {
721								crate::utils::team::kernel::with_fork(
722									team,
723									width,
724									|fork| {
725										go!(private_gemm_x86::gemm_with, fork)
726									},
727								)
728							},
729							None => go!(private_gemm_x86::gemm, 1),
730						}
731						return;
732					}
733					go!(private_gemm_x86::gemm, par.degree());
734					return;
735				}
736			}
737		}
738	}
739	let N = dst.nrows();
740	let n = N.unbound();
741	debug_assert!(dst.nrows() == dst.ncols());
742	debug_assert!(dst.nrows() == lhs.nrows());
743	debug_assert!(dst.ncols() == rhs.ncols());
744	debug_assert!(lhs.ncols() == rhs.nrows());
745	let par = if pays { par } else { Par::Seq };
746	if n <= 16 {
747		let op = {
748			#[inline(never)]
749			|| {
750				stack_mat_16x16!(
751					temp_dst,
752					N,
753					pointer_offset(dst.as_ptr()),
754					dst.row_stride(),
755					dst.col_stride(),
756					T
757				);
758				super::matmul_with_conj(
759					temp_dst.rb_mut(),
760					Accum::Replace,
761					lhs,
762					conj_lhs,
763					rhs,
764					conj_rhs,
765					alpha.clone(),
766					par,
767				);
768				accum_lower(dst, temp_dst.rb(), skip_diag, beta);
769			}
770		};
771		op();
772	} else {
773		make_guard!(HEAD);
774		make_guard!(TAIL);
775		let bs = N.partition(N.checked_idx_inc(N.unbound() / 2), HEAD, TAIL);
776		let (dst_top_left, _, dst_bot_left, dst_bot_right) =
777			dst.split_with_mut(bs, bs);
778		let (lhs_top, lhs_bot) = lhs.split_rows_with(bs);
779		let (rhs_left, rhs_right) = rhs.split_cols_with(bs);
780		join_raw(
781			|par| {
782				super::matmul_with_conj(
783					dst_bot_left,
784					beta,
785					lhs_bot,
786					conj_lhs,
787					rhs_left,
788					conj_rhs,
789					alpha.clone(),
790					par,
791				)
792			},
793			|par| {
794				join_raw(
795					|par| {
796						mat_x_mat_into_lower_impl_unchecked(
797							dst_top_left,
798							beta,
799							skip_diag,
800							lhs_top,
801							rhs_left,
802							alpha,
803							conj_lhs,
804							conj_rhs,
805							par,
806						)
807					},
808					|par| {
809						mat_x_mat_into_lower_impl_unchecked(
810							dst_bot_right,
811							beta,
812							skip_diag,
813							lhs_bot,
814							rhs_right,
815							alpha,
816							conj_lhs,
817							conj_rhs,
818							par,
819						)
820					},
821					par,
822				)
823			},
824			par,
825		);
826	}
827}
828fn mat_x_lower_into_lower_impl_unchecked<'N, T: ComplexField>(
829	dst: MatMut<'_, T, Dim<'N>, Dim<'N>>,
830	beta: Accum,
831	skip_diag: bool,
832	lhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
833	rhs: MatRef<'_, T, Dim<'N>, Dim<'N>>,
834	rhs_diag: DiagonalKind,
835	alpha: &T,
836	conj_lhs: Conj,
837	conj_rhs: Conj,
838	mut par: Par,
839) {
840	let N = dst.nrows();
841	let n = N.unbound();
842	debug_assert!(N == dst.nrows());
843	debug_assert!(N == dst.ncols());
844	debug_assert!(N == lhs.nrows());
845	debug_assert!(N == lhs.ncols());
846	debug_assert!(N == rhs.nrows());
847	debug_assert!(N == rhs.ncols());
848	if n <= 16 {
849		let op = {
850			#[inline(never)]
851			|| {
852				stack_mat_16x16!(
853					temp_dst,
854					N,
855					pointer_offset(dst.as_ptr()),
856					dst.row_stride(),
857					dst.col_stride(),
858					T
859				);
860				stack_mat_16x16!(
861					temp_rhs,
862					N,
863					pointer_offset(rhs.as_ptr()),
864					rhs.row_stride(),
865					rhs.col_stride(),
866					T
867				);
868				copy_lower(temp_rhs.rb_mut(), rhs, rhs_diag);
869				super::matmul_with_conj(
870					temp_dst.rb_mut(),
871					Accum::Replace,
872					lhs,
873					conj_lhs,
874					temp_rhs.rb(),
875					conj_rhs,
876					alpha.clone(),
877					par,
878				);
879				accum_lower(dst, temp_dst.rb(), skip_diag, beta);
880			}
881		};
882		op();
883	} else {
884		make_guard!(HEAD);
885		make_guard!(TAIL);
886		let bs = N.partition(N.checked_idx_inc(N.unbound() / 2), HEAD, TAIL);
887		let (mut dst_top_left, _, mut dst_bot_left, dst_bot_right) =
888			dst.split_with_mut(bs, bs);
889		let (lhs_top_left, lhs_top_right, lhs_bot_left, lhs_bot_right) =
890			lhs.split_with(bs, bs);
891		let (rhs_top_left, _, rhs_bot_left, rhs_bot_right) =
892			rhs.split_with(bs, bs);
893		super::matmul_with_conj(
894			dst_bot_left.rb_mut(),
895			beta,
896			lhs_bot_right,
897			conj_lhs,
898			rhs_bot_left,
899			conj_rhs,
900			alpha.clone(),
901			par.reborrow(),
902		);
903		mat_x_lower_into_lower_impl_unchecked(
904			dst_bot_right,
905			beta,
906			skip_diag,
907			lhs_bot_right,
908			rhs_bot_right,
909			rhs_diag,
910			alpha,
911			conj_lhs,
912			conj_rhs,
913			par.reborrow(),
914		);
915		mat_x_lower_into_lower_impl_unchecked(
916			dst_top_left.rb_mut(),
917			beta,
918			skip_diag,
919			lhs_top_left,
920			rhs_top_left,
921			rhs_diag,
922			alpha,
923			conj_lhs,
924			conj_rhs,
925			par.reborrow(),
926		);
927		mat_x_mat_into_lower_impl_unchecked(
928			dst_top_left,
929			Accum::Add,
930			skip_diag,
931			lhs_top_right,
932			rhs_bot_left,
933			alpha,
934			conj_lhs,
935			conj_rhs,
936			par.reborrow(),
937		);
938		mat_x_lower_impl_unchecked(
939			dst_bot_left,
940			Accum::Add,
941			lhs_bot_left,
942			rhs_top_left,
943			rhs_diag,
944			alpha,
945			conj_lhs,
946			conj_rhs,
947			par,
948		);
949	}
950}
951/// describes the parts of the matrix that must be accessed.
952#[derive(Debug, Clone, Copy, PartialEq, Eq)]
953pub enum BlockStructure {
954	/// the full matrix is accessed.
955	Rectangular,
956	/// the lower triangular half (including the diagonal) is accessed.
957	TriangularLower,
958	/// the lower triangular half (excluding the diagonal) is accessed.
959	StrictTriangularLower,
960	/// the lower triangular half (excluding the diagonal, which is assumed to
961	/// be equal to `1.0`) is accessed.
962	UnitTriangularLower,
963	/// the upper triangular half (including the diagonal) is accessed.
964	TriangularUpper,
965	/// the upper triangular half (excluding the diagonal) is accessed.
966	StrictTriangularUpper,
967	/// the upper triangular half (excluding the diagonal, which is assumed to
968	/// be equal to `1.0`) is accessed.
969	UnitTriangularUpper,
970}
971impl BlockStructure {
972	/// checks if `self` is full.
973	#[inline]
974	pub fn is_dense(self) -> bool {
975		matches!(self, BlockStructure::Rectangular)
976	}
977
978	/// checks if `self` is triangular lower (either inclusive or exclusive).
979	#[inline]
980	pub fn is_lower(self) -> bool {
981		use BlockStructure::*;
982		matches!(
983			self,
984			TriangularLower | StrictTriangularLower | UnitTriangularLower
985		)
986	}
987
988	/// checks if `self` is triangular upper (either inclusive or exclusive).
989	#[inline]
990	pub fn is_upper(self) -> bool {
991		use BlockStructure::*;
992		matches!(
993			self,
994			TriangularUpper | StrictTriangularUpper | UnitTriangularUpper
995		)
996	}
997
998	/// returns the block structure corresponding to the transposed matrix.
999	#[inline]
1000	pub fn transpose(self) -> Self {
1001		use BlockStructure::*;
1002		match self {
1003			Rectangular => Rectangular,
1004			TriangularLower => TriangularUpper,
1005			StrictTriangularLower => StrictTriangularUpper,
1006			UnitTriangularLower => UnitTriangularUpper,
1007			TriangularUpper => TriangularLower,
1008			StrictTriangularUpper => StrictTriangularLower,
1009			UnitTriangularUpper => UnitTriangularLower,
1010		}
1011	}
1012
1013	#[inline]
1014	pub(crate) fn diag_kind(self) -> DiagonalKind {
1015		use BlockStructure::*;
1016		match self {
1017			Rectangular | TriangularLower | TriangularUpper => {
1018				DiagonalKind::Generic
1019			},
1020			StrictTriangularLower | StrictTriangularUpper => DiagonalKind::Zero,
1021			UnitTriangularLower | UnitTriangularUpper => DiagonalKind::Unit,
1022		}
1023	}
1024}
1025#[track_caller]
1026fn precondition<M: Shape, N: Shape, K: Shape>(
1027	dst_nrows: M,
1028	dst_ncols: N,
1029	dst_structure: BlockStructure,
1030	lhs_nrows: M,
1031	lhs_ncols: K,
1032	lhs_structure: BlockStructure,
1033	rhs_nrows: K,
1034	rhs_ncols: N,
1035	rhs_structure: BlockStructure,
1036) {
1037	assert!(all(
1038		dst_nrows == lhs_nrows,
1039		dst_ncols == rhs_ncols,
1040		lhs_ncols == rhs_nrows,
1041	));
1042	let dst_nrows = dst_nrows.unbound();
1043	let dst_ncols = dst_ncols.unbound();
1044	let lhs_nrows = lhs_nrows.unbound();
1045	let lhs_ncols = lhs_ncols.unbound();
1046	let rhs_nrows = rhs_nrows.unbound();
1047	let rhs_ncols = rhs_ncols.unbound();
1048	if !dst_structure.is_dense() {
1049		assert!(dst_nrows == dst_ncols);
1050	}
1051	if !lhs_structure.is_dense() {
1052		assert!(lhs_nrows == lhs_ncols);
1053	}
1054	if !rhs_structure.is_dense() {
1055		assert!(rhs_nrows == rhs_ncols);
1056	}
1057}
1058/// computes the matrix product `[beta * acc] + alpha * lhs * rhs` (implicitly
1059/// conjugating the operands if needed) and stores the result in `acc`
1060///
1061/// performs the operation:
1062/// - `acc = alpha * lhs * rhs` if `beta` is `accum::replace` (in this case, the
1063///   preexisting
1064/// values in `acc` are not read)
1065/// - `acc = acc + alpha * lhs * rhs` if `beta` is `accum::add`
1066///
1067/// the left hand side and right hand side may be interpreted as triangular
1068/// depending on the given corresponding matrix structure.
1069///
1070/// for the destination matrix, the result is:
1071/// - fully computed if the structure is rectangular,
1072/// - only the triangular half (including the diagonal) is computed if the
1073///   structure is
1074/// triangular
1075/// - only the strict triangular half (excluding the diagonal) is computed if
1076///   the structure is
1077/// strictly triangular or unit triangular
1078///
1079/// # panics
1080///
1081/// panics if the matrix dimensions are not compatible for matrix
1082/// multiplication. i.e.
1083///  - `acc.nrows() == lhs.nrows()`
1084///  - `acc.ncols() == rhs.ncols()`
1085///  - `lhs.ncols() == rhs.nrows()`
1086///
1087/// additionally, matrices that are marked as triangular must be square, i.e.,
1088/// they must have the same number of rows and columns.
1089///
1090/// # example
1091///
1092/// ```
1093/// use faer::linalg::matmul::triangular::{BlockStructure, matmul_with_conj};
1094/// use faer::{Accum, Conj, Mat, Par, mat, unzip, zip};
1095/// let lhs = mat![[0.0, 2.0], [1.0, 3.0]];
1096/// let rhs = mat![[4.0, 6.0], [5.0, 7.0]];
1097/// let mut acc = Mat::<f64>::zeros(2, 2);
1098/// let target = mat![
1099/// 	[
1100/// 		2.5 * (lhs[(0, 0)] * rhs[(0, 0)] + lhs[(0, 1)] * rhs[(1, 0)]),
1101/// 		0.0,
1102/// 	],
1103/// 	[
1104/// 		2.5 * (lhs[(1, 0)] * rhs[(0, 0)] + lhs[(1, 1)] * rhs[(1, 0)]),
1105/// 		2.5 * (lhs[(1, 0)] * rhs[(0, 1)] + lhs[(1, 1)] * rhs[(1, 1)]),
1106/// 	],
1107/// ];
1108/// matmul_with_conj(
1109/// 	&mut acc,
1110/// 	BlockStructure::TriangularLower,
1111/// 	Accum::Replace,
1112/// 	&lhs,
1113/// 	BlockStructure::Rectangular,
1114/// 	Conj::No,
1115/// 	&rhs,
1116/// 	BlockStructure::Rectangular,
1117/// 	Conj::No,
1118/// 	2.5,
1119/// 	Par::Seq,
1120/// );
1121/// zip!(&acc, &target)
1122/// 	.for_each(|unzip!(acc, target)| assert!((acc - target).abs() < 1e-10));
1123/// ```
1124#[track_caller]
1125#[inline]
1126pub fn matmul_with_conj<T: ComplexField, M: Shape, N: Shape, K: Shape>(
1127	dst: impl AsMatMut<T = T, Rows = M, Cols = N>,
1128	dst_structure: BlockStructure,
1129	beta: Accum,
1130	lhs: impl AsMatRef<T = T, Rows = M, Cols = K>,
1131	lhs_structure: BlockStructure,
1132	conj_lhs: Conj,
1133	rhs: impl AsMatRef<T = T, Rows = K, Cols = N>,
1134	rhs_structure: BlockStructure,
1135	conj_rhs: Conj,
1136	alpha: T,
1137	par: Par,
1138) {
1139	let mut dst = dst;
1140	let dst = dst.as_mat_mut();
1141	let lhs = lhs.as_mat_ref();
1142	let rhs = rhs.as_mat_ref();
1143	precondition(
1144		dst.nrows(),
1145		dst.ncols(),
1146		dst_structure,
1147		lhs.nrows(),
1148		lhs.ncols(),
1149		lhs_structure,
1150		rhs.nrows(),
1151		rhs.ncols(),
1152		rhs_structure,
1153	);
1154	make_guard!(M);
1155	make_guard!(N);
1156	make_guard!(K);
1157	let M = dst.nrows().bind(M);
1158	let N = dst.ncols().bind(N);
1159	let K = lhs.ncols().bind(K);
1160	matmul_imp(
1161		dst.as_dyn_stride_mut().as_shape_mut(M, N),
1162		dst_structure,
1163		beta,
1164		lhs.as_dyn_stride().canonical().as_shape(M, K),
1165		lhs_structure,
1166		conj_lhs,
1167		rhs.as_dyn_stride().canonical().as_shape(K, N),
1168		rhs_structure,
1169		conj_rhs,
1170		&alpha,
1171		par,
1172	);
1173}
1174/// computes the matrix product `[beta * acc] + alpha * lhs * rhs` (implicitly
1175/// conjugating the operands if needed) and stores the result in `acc`
1176///
1177/// performs the operation:
1178/// - `acc = alpha * lhs * rhs` if `beta` is `accum::replace` (in this case, the
1179///   preexisting
1180/// values in `acc` are not read)
1181/// - `acc = acc + alpha * lhs * rhs` if `beta` is `accum::add`
1182///
1183/// the left hand side and right hand side may be interpreted as triangular
1184/// depending on the given corresponding matrix structure.
1185///
1186/// for the destination matrix, the result is:
1187/// - fully computed if the structure is rectangular,
1188/// - only the triangular half (including the diagonal) is computed if the
1189///   structure is
1190/// triangular
1191/// - only the strict triangular half (excluding the diagonal) is computed if
1192///   the structure is
1193/// strictly triangular or unit triangular
1194///
1195/// # panics
1196///
1197/// panics if the matrix dimensions are not compatible for matrix
1198/// multiplication. i.e.
1199///  - `acc.nrows() == lhs.nrows()`
1200///  - `acc.ncols() == rhs.ncols()`
1201///  - `lhs.ncols() == rhs.nrows()`
1202///
1203/// additionally, matrices that are marked as triangular must be square, i.e.,
1204/// they must have the same number of rows and columns.
1205///
1206/// # example
1207///
1208/// ```
1209/// use faer::linalg::matmul::triangular::{BlockStructure, matmul};
1210/// use faer::{Accum, Conj, Mat, Par, mat, unzip, zip};
1211/// let lhs = mat![[0.0, 2.0], [1.0, 3.0]];
1212/// let rhs = mat![[4.0, 6.0], [5.0, 7.0]];
1213/// let mut acc = Mat::<f64>::zeros(2, 2);
1214/// let target = mat![
1215/// 	[
1216/// 		2.5 * (lhs[(0, 0)] * rhs[(0, 0)] + lhs[(0, 1)] * rhs[(1, 0)]),
1217/// 		0.0,
1218/// 	],
1219/// 	[
1220/// 		2.5 * (lhs[(1, 0)] * rhs[(0, 0)] + lhs[(1, 1)] * rhs[(1, 0)]),
1221/// 		2.5 * (lhs[(1, 0)] * rhs[(0, 1)] + lhs[(1, 1)] * rhs[(1, 1)]),
1222/// 	],
1223/// ];
1224/// matmul(
1225/// 	&mut acc,
1226/// 	BlockStructure::TriangularLower,
1227/// 	Accum::Replace,
1228/// 	&lhs,
1229/// 	BlockStructure::Rectangular,
1230/// 	&rhs,
1231/// 	BlockStructure::Rectangular,
1232/// 	2.5,
1233/// 	Par::Seq,
1234/// );
1235/// zip!(&acc, &target)
1236/// 	.for_each(|unzip!(acc, target)| assert!((acc - target).abs() < 1e-10));
1237/// ```
1238#[track_caller]
1239#[inline]
1240pub fn matmul<
1241	T: ComplexField,
1242	LhsT: Conjugate<Canonical = T>,
1243	RhsT: Conjugate<Canonical = T>,
1244	M: Shape,
1245	N: Shape,
1246	K: Shape,
1247>(
1248	dst: impl AsMatMut<T = T, Rows = M, Cols = N>,
1249	dst_structure: BlockStructure,
1250	beta: Accum,
1251	lhs: impl AsMatRef<T = LhsT, Rows = M, Cols = K>,
1252	lhs_structure: BlockStructure,
1253	rhs: impl AsMatRef<T = RhsT, Rows = K, Cols = N>,
1254	rhs_structure: BlockStructure,
1255	alpha: T,
1256	par: Par,
1257) {
1258	let mut dst = dst;
1259	let dst = dst.as_mat_mut();
1260	let lhs = lhs.as_mat_ref();
1261	let rhs = rhs.as_mat_ref();
1262	precondition(
1263		dst.nrows(),
1264		dst.ncols(),
1265		dst_structure,
1266		lhs.nrows(),
1267		lhs.ncols(),
1268		lhs_structure,
1269		rhs.nrows(),
1270		rhs.ncols(),
1271		rhs_structure,
1272	);
1273	make_guard!(M);
1274	make_guard!(N);
1275	make_guard!(K);
1276	let M = dst.nrows().bind(M);
1277	let N = dst.ncols().bind(N);
1278	let K = lhs.ncols().bind(K);
1279	matmul_imp(
1280		dst.as_dyn_stride_mut().as_shape_mut(M, N),
1281		dst_structure,
1282		beta,
1283		lhs.as_dyn_stride().canonical().as_shape(M, K),
1284		lhs_structure,
1285		const { Conj::get::<LhsT>() },
1286		rhs.as_dyn_stride().canonical().as_shape(K, N),
1287		rhs_structure,
1288		const { Conj::get::<RhsT>() },
1289		alpha.by_ref(),
1290		par,
1291	);
1292}
1293fn matmul_imp<'M, 'N, 'K, T: ComplexField>(
1294	dst: MatMut<'_, T, Dim<'M>, Dim<'N>>,
1295	dst_structure: BlockStructure,
1296	beta: Accum,
1297	lhs: MatRef<'_, T, Dim<'M>, Dim<'K>>,
1298	lhs_structure: BlockStructure,
1299	conj_lhs: Conj,
1300	rhs: MatRef<'_, T, Dim<'K>, Dim<'N>>,
1301	rhs_structure: BlockStructure,
1302	conj_rhs: Conj,
1303	alpha: &T,
1304	par: Par,
1305) {
1306	let mut acc = dst.as_dyn_mut();
1307	let mut lhs = lhs.as_dyn();
1308	let mut rhs = rhs.as_dyn();
1309	let mut acc_structure = dst_structure;
1310	let mut lhs_structure = lhs_structure;
1311	let mut rhs_structure = rhs_structure;
1312	let mut conj_lhs = conj_lhs;
1313	let mut conj_rhs = conj_rhs;
1314	if rhs_structure.is_lower() {
1315		false
1316	} else if rhs_structure.is_upper() {
1317		acc = acc.reverse_rows_and_cols_mut();
1318		lhs = lhs.reverse_rows_and_cols();
1319		rhs = rhs.reverse_rows_and_cols();
1320		acc_structure = acc_structure.transpose();
1321		lhs_structure = lhs_structure.transpose();
1322		rhs_structure = rhs_structure.transpose();
1323		false
1324	} else if lhs_structure.is_lower() {
1325		acc = acc.reverse_rows_and_cols_mut().transpose_mut();
1326		(lhs, rhs) = (
1327			rhs.reverse_rows_and_cols().transpose(),
1328			lhs.reverse_rows_and_cols().transpose(),
1329		);
1330		(conj_lhs, conj_rhs) = (conj_rhs, conj_lhs);
1331		(lhs_structure, rhs_structure) = (rhs_structure, lhs_structure);
1332		true
1333	} else if lhs_structure.is_upper() {
1334		acc_structure = acc_structure.transpose();
1335		acc = acc.transpose_mut();
1336		(lhs, rhs) = (rhs.transpose(), lhs.transpose());
1337		(conj_lhs, conj_rhs) = (conj_rhs, conj_lhs);
1338		(lhs_structure, rhs_structure) =
1339			(rhs_structure.transpose(), lhs_structure.transpose());
1340		true
1341	} else {
1342		false
1343	};
1344	make_guard!(M);
1345	make_guard!(N);
1346	make_guard!(K);
1347	let M = acc.nrows().bind(M);
1348	let N = acc.ncols().bind(N);
1349	let K = lhs.ncols().bind(K);
1350	let clear_upper = |acc: MatMut<'_, T>, skip_diag: bool| match &beta {
1351		Accum::Add => {},
1352		Accum::Replace => zip!(acc).for_each_triangular_upper(
1353			if skip_diag { Diag::Skip } else { Diag::Include },
1354			|unzip!(acc)| *acc = zero(),
1355		),
1356	};
1357	let skip_diag = matches!(
1358		acc_structure,
1359		BlockStructure::StrictTriangularLower
1360			| BlockStructure::StrictTriangularUpper
1361			| BlockStructure::UnitTriangularLower
1362			| BlockStructure::UnitTriangularUpper
1363	);
1364	let lhs_diag = lhs_structure.diag_kind();
1365	let rhs_diag = rhs_structure.diag_kind();
1366	if acc_structure.is_dense() {
1367		if lhs_structure.is_dense() && rhs_structure.is_dense() {
1368			super::matmul_with_conj(
1369				acc,
1370				beta,
1371				lhs,
1372				conj_lhs,
1373				rhs,
1374				conj_rhs,
1375				alpha.clone(),
1376				par,
1377			);
1378		} else {
1379			debug_assert!(rhs_structure.is_lower());
1380			if lhs_structure.is_dense() {
1381				mat_x_lower_impl_unchecked(
1382					acc.as_shape_mut(M, N),
1383					beta,
1384					lhs.as_shape(M, N),
1385					rhs.as_shape(N, N),
1386					rhs_diag,
1387					alpha,
1388					conj_lhs,
1389					conj_rhs,
1390					par,
1391				)
1392			} else if lhs_structure.is_lower() {
1393				clear_upper(acc.rb_mut(), true);
1394				lower_x_lower_into_lower_impl_unchecked(
1395					acc.as_shape_mut(N, N),
1396					beta,
1397					false,
1398					lhs.as_shape(N, N),
1399					lhs_diag,
1400					rhs.as_shape(N, N),
1401					rhs_diag,
1402					alpha,
1403					conj_lhs,
1404					conj_rhs,
1405					par,
1406				);
1407			} else {
1408				debug_assert!(lhs_structure.is_upper());
1409				upper_x_lower_impl_unchecked(
1410					acc.as_shape_mut(N, N),
1411					beta,
1412					lhs.as_shape(N, N),
1413					lhs_diag,
1414					rhs.as_shape(N, N),
1415					rhs_diag,
1416					alpha,
1417					conj_lhs,
1418					conj_rhs,
1419					par,
1420				)
1421			}
1422		}
1423	} else if acc_structure.is_lower() {
1424		if lhs_structure.is_dense() && rhs_structure.is_dense() {
1425			mat_x_mat_into_lower_impl_unchecked(
1426				acc.as_shape_mut(N, N),
1427				beta,
1428				skip_diag,
1429				lhs.as_shape(N, K),
1430				rhs.as_shape(K, N),
1431				alpha,
1432				conj_lhs,
1433				conj_rhs,
1434				par,
1435			)
1436		} else {
1437			debug_assert!(rhs_structure.is_lower());
1438			if lhs_structure.is_dense() {
1439				mat_x_lower_into_lower_impl_unchecked(
1440					acc.as_shape_mut(N, N),
1441					beta,
1442					skip_diag,
1443					lhs.as_shape(N, N),
1444					rhs.as_shape(N, N),
1445					rhs_diag,
1446					alpha,
1447					conj_lhs,
1448					conj_rhs,
1449					par,
1450				);
1451			} else if lhs_structure.is_lower() {
1452				lower_x_lower_into_lower_impl_unchecked(
1453					acc.as_shape_mut(N, N),
1454					beta,
1455					skip_diag,
1456					lhs.as_shape(N, N),
1457					lhs_diag,
1458					rhs.as_shape(N, N),
1459					rhs_diag,
1460					alpha,
1461					conj_lhs,
1462					conj_rhs,
1463					par,
1464				)
1465			} else {
1466				upper_x_lower_into_lower_impl_unchecked(
1467					acc.as_shape_mut(N, N),
1468					beta,
1469					skip_diag,
1470					lhs.as_shape(N, N),
1471					lhs_diag,
1472					rhs.as_shape(N, N),
1473					rhs_diag,
1474					alpha,
1475					conj_lhs,
1476					conj_rhs,
1477					par,
1478				)
1479			}
1480		}
1481	} else if lhs_structure.is_dense() && rhs_structure.is_dense() {
1482		mat_x_mat_into_lower_impl_unchecked(
1483			acc.as_shape_mut(N, N).transpose_mut(),
1484			beta,
1485			skip_diag,
1486			rhs.transpose().as_shape(N, K),
1487			lhs.transpose().as_shape(K, N),
1488			alpha,
1489			conj_rhs,
1490			conj_lhs,
1491			par,
1492		)
1493	} else {
1494		debug_assert!(rhs_structure.is_lower());
1495		if lhs_structure.is_dense() {
1496			upper_x_lower_into_lower_impl_unchecked(
1497				acc.as_shape_mut(N, N).transpose_mut(),
1498				beta,
1499				skip_diag,
1500				rhs.transpose().as_shape(N, N),
1501				rhs_diag,
1502				lhs.transpose().as_shape(N, N),
1503				lhs_diag,
1504				alpha,
1505				conj_rhs,
1506				conj_lhs,
1507				par,
1508			)
1509		} else if lhs_structure.is_lower() {
1510			if !skip_diag {
1511				match beta {
1512					Accum::Add => {
1513						for j in 0..N.unbound() {
1514							acc[(j, j)] += alpha * &lhs[(j, j)] * &rhs[(j, j)];
1515						}
1516					},
1517					Accum::Replace => {
1518						for j in 0..N.unbound() {
1519							acc[(j, j)] = alpha * &lhs[(j, j)] * &rhs[(j, j)];
1520						}
1521					},
1522				}
1523			}
1524			clear_upper(acc.rb_mut(), true);
1525		} else {
1526			debug_assert!(lhs_structure.is_upper());
1527			upper_x_lower_into_lower_impl_unchecked(
1528				acc.as_shape_mut(N, N).transpose_mut(),
1529				beta,
1530				skip_diag,
1531				rhs.transpose().as_shape(N, N),
1532				rhs_diag,
1533				lhs.transpose().as_shape(N, N),
1534				lhs_diag,
1535				alpha,
1536				conj_rhs,
1537				conj_lhs,
1538				par,
1539			)
1540		}
1541	}
1542}
faer/linalg/matmul/triangular.rs

faer/linalg/matmul/
triangular.rs
