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

faer/linalg/matmul/
triangular.rs
