oxiblas_matrix/
triangular.rs

1//! Triangular matrix types.
2//!
3//! Triangular matrices are square matrices where all elements above (lower)
4//! or below (upper) the main diagonal are zero. They arise frequently in
5//! matrix decompositions (LU, Cholesky, QR) and are important in BLAS.
6//!
7//! # Storage
8//!
9//! This module provides two storage options:
10//! - Full storage: Uses a standard dense matrix, treating the non-stored
11//!   triangle as zeros implicitly.
12//! - Packed storage: Uses `PackedMat` to store only the triangular portion.
13//!
14//! # Unit Triangular
15//!
16//! A unit triangular matrix has ones on the diagonal. This is common in LU
17//! decomposition where L is unit lower triangular.
18
19#[cfg(not(feature = "std"))]
20use alloc::{vec, vec::Vec};
21
22use crate::packed::{PackedMat, PackedMut, PackedRef, TriangularKind};
23use crate::{Mat, MatMut, MatRef};
24use oxiblas_core::scalar::Scalar;
25
26/// Diagonal type for triangular matrices.
27#[derive(Debug, Clone, Copy, PartialEq, Eq)]
28pub enum DiagonalKind {
29    /// Non-unit diagonal (general values on diagonal).
30    NonUnit,
31    /// Unit diagonal (ones on diagonal, not stored).
32    Unit,
33}
34
35/// A triangular matrix view over dense storage.
36///
37/// This is a logical view that treats the non-stored triangle as zeros.
38/// The underlying storage is a full dense matrix, but operations only
39/// access the stored triangle.
40///
41/// # Example
42///
43/// ```
44/// use oxiblas_matrix::{Mat, triangular::{TriangularView, DiagonalKind}};
45/// use oxiblas_matrix::packed::TriangularKind;
46///
47/// let mut m: Mat<f64> = Mat::eye(3);
48/// m[(0, 1)] = 2.0;
49/// m[(0, 2)] = 3.0;
50/// m[(1, 2)] = 4.0;
51///
52/// // Create upper triangular view
53/// let tri = TriangularView::new(m.as_ref(), TriangularKind::Upper, DiagonalKind::NonUnit);
54///
55/// // Access upper triangle
56/// assert_eq!(tri.get(0, 1), Some(&2.0));
57/// // Lower triangle returns zero
58/// assert_eq!(tri.get(1, 0), None);
59/// ```
60#[derive(Clone, Copy)]
61pub struct TriangularView<'a, T: Scalar> {
62    /// Underlying matrix view.
63    inner: MatRef<'a, T>,
64    /// Upper or lower triangular.
65    uplo: TriangularKind,
66    /// Unit or non-unit diagonal.
67    diag: DiagonalKind,
68}
69
70impl<'a, T: Scalar> TriangularView<'a, T> {
71    /// Creates a new triangular view over a matrix.
72    ///
73    /// # Panics
74    /// Panics if the matrix is not square.
75    #[inline]
76    pub fn new(mat: MatRef<'a, T>, uplo: TriangularKind, diag: DiagonalKind) -> Self {
77        assert!(mat.is_square(), "Triangular matrix must be square");
78        TriangularView {
79            inner: mat,
80            uplo,
81            diag,
82        }
83    }
84
85    /// Returns the matrix dimension.
86    #[inline]
87    pub fn dim(&self) -> usize {
88        self.inner.nrows()
89    }
90
91    /// Returns the shape (n, n).
92    #[inline]
93    pub fn shape(&self) -> (usize, usize) {
94        self.inner.shape()
95    }
96
97    /// Returns the triangular kind (upper/lower).
98    #[inline]
99    pub fn uplo(&self) -> TriangularKind {
100        self.uplo
101    }
102
103    /// Returns the diagonal kind (unit/non-unit).
104    #[inline]
105    pub fn diag(&self) -> DiagonalKind {
106        self.diag
107    }
108
109    /// Returns true if element (row, col) is in the stored triangle.
110    #[inline]
111    pub fn in_triangle(&self, row: usize, col: usize) -> bool {
112        match self.uplo {
113            TriangularKind::Upper => row <= col,
114            TriangularKind::Lower => row >= col,
115        }
116    }
117
118    /// Returns a reference to the element at (row, col).
119    ///
120    /// Returns `None` if the element is in the non-stored triangle.
121    /// For unit triangular matrices, diagonal elements return `None`
122    /// (they are implicitly one).
123    #[inline]
124    pub fn get(&self, row: usize, col: usize) -> Option<&T> {
125        if !self.in_triangle(row, col) {
126            return None;
127        }
128
129        if self.diag == DiagonalKind::Unit && row == col {
130            return None; // Diagonal is implicitly one
131        }
132
133        self.inner.get(row, col)
134    }
135
136    /// Returns the underlying matrix reference.
137    #[inline]
138    pub fn as_inner(&self) -> MatRef<'a, T> {
139        self.inner
140    }
141
142    /// Returns a pointer to the matrix data.
143    #[inline]
144    pub fn as_ptr(&self) -> *const T {
145        self.inner.as_ptr()
146    }
147
148    /// Returns the row stride.
149    #[inline]
150    pub fn row_stride(&self) -> usize {
151        self.inner.row_stride()
152    }
153
154    /// Converts to a full dense matrix.
155    pub fn to_dense(&self) -> Mat<T>
156    where
157        T: bytemuck::Zeroable,
158    {
159        let n = self.dim();
160        let mut mat = Mat::zeros(n, n);
161
162        for j in 0..n {
163            for i in 0..n {
164                if self.in_triangle(i, j) {
165                    if self.diag == DiagonalKind::Unit && i == j {
166                        mat[(i, j)] = T::one();
167                    } else {
168                        mat[(i, j)] = self.inner[(i, j)];
169                    }
170                }
171            }
172        }
173
174        mat
175    }
176
177    /// Converts to packed storage.
178    pub fn to_packed(&self) -> PackedMat<T>
179    where
180        T: bytemuck::Zeroable,
181    {
182        PackedMat::from_dense(&self.inner, self.uplo)
183    }
184}
185
186/// A mutable triangular matrix view over dense storage.
187pub struct TriangularViewMut<'a, T: Scalar> {
188    /// Underlying mutable matrix view.
189    inner: MatMut<'a, T>,
190    /// Upper or lower triangular.
191    uplo: TriangularKind,
192    /// Unit or non-unit diagonal.
193    diag: DiagonalKind,
194}
195
196impl<'a, T: Scalar> TriangularViewMut<'a, T> {
197    /// Creates a new mutable triangular view over a matrix.
198    #[inline]
199    pub fn new(mat: MatMut<'a, T>, uplo: TriangularKind, diag: DiagonalKind) -> Self {
200        assert!(mat.is_square(), "Triangular matrix must be square");
201        TriangularViewMut {
202            inner: mat,
203            uplo,
204            diag,
205        }
206    }
207
208    /// Returns the matrix dimension.
209    #[inline]
210    pub fn dim(&self) -> usize {
211        self.inner.nrows()
212    }
213
214    /// Returns the shape (n, n).
215    #[inline]
216    pub fn shape(&self) -> (usize, usize) {
217        self.inner.shape()
218    }
219
220    /// Returns the triangular kind.
221    #[inline]
222    pub fn uplo(&self) -> TriangularKind {
223        self.uplo
224    }
225
226    /// Returns the diagonal kind.
227    #[inline]
228    pub fn diag(&self) -> DiagonalKind {
229        self.diag
230    }
231
232    /// Returns true if element is in stored triangle.
233    #[inline]
234    pub fn in_triangle(&self, row: usize, col: usize) -> bool {
235        match self.uplo {
236            TriangularKind::Upper => row <= col,
237            TriangularKind::Lower => row >= col,
238        }
239    }
240
241    /// Returns a reference to element.
242    #[inline]
243    pub fn get(&self, row: usize, col: usize) -> Option<&T> {
244        if !self.in_triangle(row, col) {
245            return None;
246        }
247        if self.diag == DiagonalKind::Unit && row == col {
248            return None;
249        }
250        self.inner.get(row, col)
251    }
252
253    /// Returns a mutable reference to element.
254    #[inline]
255    pub fn get_mut(&mut self, row: usize, col: usize) -> Option<&mut T> {
256        if !self.in_triangle(row, col) {
257            return None;
258        }
259        if self.diag == DiagonalKind::Unit && row == col {
260            return None;
261        }
262        self.inner.get_mut(row, col)
263    }
264
265    /// Sets an element in the stored triangle.
266    ///
267    /// # Panics
268    /// Panics if the element is outside the stored triangle or on the
269    /// diagonal for unit triangular matrices.
270    #[inline]
271    pub fn set(&mut self, row: usize, col: usize, value: T) {
272        assert!(
273            self.in_triangle(row, col),
274            "Element outside stored triangle"
275        );
276        assert!(
277            !(self.diag == DiagonalKind::Unit && row == col),
278            "Cannot set diagonal of unit triangular matrix"
279        );
280        self.inner.set(row, col, value);
281    }
282
283    /// Creates an immutable reborrow.
284    #[inline]
285    pub fn rb(&self) -> TriangularView<'_, T> {
286        TriangularView {
287            inner: self.inner.rb(),
288            uplo: self.uplo,
289            diag: self.diag,
290        }
291    }
292
293    /// Creates a mutable reborrow.
294    #[inline]
295    pub fn rb_mut(&mut self) -> TriangularViewMut<'_, T> {
296        TriangularViewMut {
297            inner: self.inner.rb_mut(),
298            uplo: self.uplo,
299            diag: self.diag,
300        }
301    }
302
303    /// Returns a mutable pointer to the data.
304    #[inline]
305    pub fn as_mut_ptr(&mut self) -> *mut T {
306        self.inner.as_mut_ptr()
307    }
308
309    /// Fills the stored triangle with a value.
310    ///
311    /// Does not modify the diagonal for unit triangular matrices.
312    pub fn fill(&mut self, value: T) {
313        let n = self.dim();
314        for j in 0..n {
315            for i in 0..n {
316                if self.in_triangle(i, j) {
317                    if self.diag == DiagonalKind::Unit && i == j {
318                        continue;
319                    }
320                    self.inner.set(i, j, value);
321                }
322            }
323        }
324    }
325
326    /// Scales the stored triangle by a scalar.
327    pub fn scale(&mut self, alpha: T) {
328        let n = self.dim();
329        for j in 0..n {
330            for i in 0..n {
331                if self.in_triangle(i, j) {
332                    if self.diag == DiagonalKind::Unit && i == j {
333                        continue;
334                    }
335                    if let Some(val) = self.inner.get(i, j) {
336                        self.inner.set(i, j, *val * alpha);
337                    }
338                }
339            }
340        }
341    }
342
343    /// Clears the non-stored triangle to zero.
344    pub fn zero_non_triangle(&mut self)
345    where
346        T: num_traits::Zero,
347    {
348        let n = self.dim();
349        for j in 0..n {
350            for i in 0..n {
351                if !self.in_triangle(i, j) {
352                    self.inner.set(i, j, T::zero());
353                }
354            }
355        }
356    }
357}
358
359/// Triangular matrix using packed storage.
360///
361/// This is a wrapper around `PackedMat` that provides triangular matrix
362/// semantics with efficient packed storage.
363#[derive(Clone)]
364pub struct TriangularMat<T: Scalar> {
365    /// Packed storage.
366    packed: PackedMat<T>,
367    /// Unit or non-unit diagonal.
368    diag: DiagonalKind,
369}
370
371impl<T: Scalar> TriangularMat<T> {
372    /// Creates a new triangular matrix filled with zeros.
373    pub fn zeros(n: usize, uplo: TriangularKind, diag: DiagonalKind) -> Self
374    where
375        T: bytemuck::Zeroable,
376    {
377        TriangularMat {
378            packed: PackedMat::zeros(n, uplo),
379            diag,
380        }
381    }
382
383    /// Creates a unit triangular matrix (identity-like).
384    ///
385    /// The diagonal is implicitly one and not stored.
386    pub fn unit_zeros(n: usize, uplo: TriangularKind) -> Self
387    where
388        T: bytemuck::Zeroable,
389    {
390        Self::zeros(n, uplo, DiagonalKind::Unit)
391    }
392
393    /// Creates from a packed matrix.
394    #[inline]
395    pub fn from_packed(packed: PackedMat<T>, diag: DiagonalKind) -> Self {
396        TriangularMat { packed, diag }
397    }
398
399    /// Creates from a dense matrix view.
400    pub fn from_dense(mat: &MatRef<'_, T>, uplo: TriangularKind, diag: DiagonalKind) -> Self
401    where
402        T: bytemuck::Zeroable,
403    {
404        TriangularMat {
405            packed: PackedMat::from_dense(mat, uplo),
406            diag,
407        }
408    }
409
410    /// Returns the matrix dimension.
411    #[inline]
412    pub fn dim(&self) -> usize {
413        self.packed.dim()
414    }
415
416    /// Returns the shape (n, n).
417    #[inline]
418    pub fn shape(&self) -> (usize, usize) {
419        let n = self.dim();
420        (n, n)
421    }
422
423    /// Returns the triangular kind.
424    #[inline]
425    pub fn uplo(&self) -> TriangularKind {
426        self.packed.kind()
427    }
428
429    /// Returns the diagonal kind.
430    #[inline]
431    pub fn diag(&self) -> DiagonalKind {
432        self.diag
433    }
434
435    /// Returns the packed storage length.
436    #[inline]
437    pub fn len(&self) -> usize {
438        self.packed.len()
439    }
440
441    /// Returns true if empty.
442    #[inline]
443    pub fn is_empty(&self) -> bool {
444        self.packed.is_empty()
445    }
446
447    /// Returns true if element is in the stored triangle.
448    #[inline]
449    pub fn in_triangle(&self, row: usize, col: usize) -> bool {
450        self.packed.packed_index(row, col).is_some()
451    }
452
453    /// Returns a reference to element.
454    ///
455    /// For unit triangular matrices, returns `None` for diagonal elements.
456    #[inline]
457    pub fn get(&self, row: usize, col: usize) -> Option<&T> {
458        if self.diag == DiagonalKind::Unit && row == col {
459            return None;
460        }
461        self.packed.get(row, col)
462    }
463
464    /// Returns a mutable reference to element.
465    #[inline]
466    pub fn get_mut(&mut self, row: usize, col: usize) -> Option<&mut T> {
467        if self.diag == DiagonalKind::Unit && row == col {
468            return None;
469        }
470        self.packed.get_mut(row, col)
471    }
472
473    /// Sets an element.
474    ///
475    /// # Panics
476    /// Panics if setting diagonal on unit triangular matrix.
477    #[inline]
478    pub fn set(&mut self, row: usize, col: usize, value: T) {
479        assert!(
480            !(self.diag == DiagonalKind::Unit && row == col),
481            "Cannot set diagonal of unit triangular matrix"
482        );
483        self.packed.set(row, col, value);
484    }
485
486    /// Returns a pointer to the packed data.
487    #[inline]
488    pub fn as_ptr(&self) -> *const T {
489        self.packed.as_ptr()
490    }
491
492    /// Returns a mutable pointer to the packed data.
493    #[inline]
494    pub fn as_mut_ptr(&mut self) -> *mut T {
495        self.packed.as_mut_ptr()
496    }
497
498    /// Returns the packed data as a slice.
499    #[inline]
500    pub fn as_slice(&self) -> &[T] {
501        self.packed.as_slice()
502    }
503
504    /// Returns the packed data as a mutable slice.
505    #[inline]
506    pub fn as_slice_mut(&mut self) -> &mut [T] {
507        self.packed.as_slice_mut()
508    }
509
510    /// Returns a reference to the underlying packed matrix.
511    #[inline]
512    pub fn as_packed(&self) -> &PackedMat<T> {
513        &self.packed
514    }
515
516    /// Returns a mutable reference to the underlying packed matrix.
517    #[inline]
518    pub fn as_packed_mut(&mut self) -> &mut PackedMat<T> {
519        &mut self.packed
520    }
521
522    /// Converts to a full dense matrix.
523    pub fn to_dense(&self) -> Mat<T>
524    where
525        T: bytemuck::Zeroable,
526    {
527        let n = self.dim();
528        let mut mat = Mat::zeros(n, n);
529
530        for j in 0..n {
531            for i in 0..n {
532                if let Some(&val) = self.packed.get(i, j) {
533                    if self.diag == DiagonalKind::Unit && i == j {
534                        mat[(i, j)] = T::one();
535                    } else {
536                        mat[(i, j)] = val;
537                    }
538                } else if self.diag == DiagonalKind::Unit && i == j {
539                    mat[(i, j)] = T::one();
540                }
541            }
542        }
543
544        mat
545    }
546
547    /// Returns the diagonal elements.
548    ///
549    /// For unit triangular matrices, returns a vector of ones.
550    pub fn diagonal(&self) -> Vec<T> {
551        let n = self.dim();
552        if self.diag == DiagonalKind::Unit {
553            vec![T::one(); n]
554        } else {
555            self.packed.diagonal()
556        }
557    }
558
559    /// Sets the diagonal elements.
560    ///
561    /// # Panics
562    /// Panics for unit triangular matrices.
563    pub fn set_diagonal(&mut self, diag: &[T]) {
564        assert!(
565            self.diag != DiagonalKind::Unit,
566            "Cannot set diagonal of unit triangular matrix"
567        );
568        self.packed.set_diagonal(diag);
569    }
570
571    /// Fills the stored triangle with a value.
572    pub fn fill(&mut self, value: T) {
573        if self.diag == DiagonalKind::Unit {
574            // Fill off-diagonal only
575            let n = self.dim();
576            for j in 0..n {
577                for i in 0..n {
578                    if self.in_triangle(i, j) && i != j {
579                        self.packed.set(i, j, value);
580                    }
581                }
582            }
583        } else {
584            self.packed.fill(value);
585        }
586    }
587
588    /// Scales the stored triangle by a scalar.
589    pub fn scale(&mut self, alpha: T) {
590        if self.diag == DiagonalKind::Unit {
591            // Scale off-diagonal only
592            let n = self.dim();
593            for j in 0..n {
594                for i in 0..n {
595                    if self.in_triangle(i, j) && i != j {
596                        if let Some(val) = self.packed.get_mut(i, j) {
597                            *val *= alpha;
598                        }
599                    }
600                }
601            }
602        } else {
603            self.packed.scale(alpha);
604        }
605    }
606
607    /// Returns the transpose (flips upper/lower).
608    pub fn transpose(&self) -> Self
609    where
610        T: bytemuck::Zeroable,
611    {
612        TriangularMat {
613            packed: self.packed.transpose(),
614            diag: self.diag,
615        }
616    }
617}
618
619impl<T: Scalar + core::fmt::Debug> core::fmt::Debug for TriangularMat<T> {
620    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
621        let n = self.dim();
622        writeln!(
623            f,
624            "TriangularMat {}×{} ({:?}, {:?}) {{",
625            n,
626            n,
627            self.uplo(),
628            self.diag
629        )?;
630
631        let max_dim = 8.min(n);
632        for i in 0..max_dim {
633            write!(f, "  [")?;
634            for j in 0..max_dim {
635                if j > 0 {
636                    write!(f, ", ")?;
637                }
638                if self.in_triangle(i, j) {
639                    if self.diag == DiagonalKind::Unit && i == j {
640                        write!(f, "{:8.4?}", T::one())?;
641                    } else if let Some(v) = self.packed.get(i, j) {
642                        write!(f, "{:8.4?}", v)?;
643                    } else {
644                        write!(f, "      * ")?;
645                    }
646                } else {
647                    write!(f, "      0 ")?;
648                }
649            }
650            if n > max_dim {
651                write!(f, ", ...")?;
652            }
653            writeln!(f, "]")?;
654        }
655        if n > max_dim {
656            writeln!(f, "  ...")?;
657        }
658        write!(f, "}}")
659    }
660}
661
662/// A reference to triangular packed data.
663#[derive(Clone, Copy)]
664pub struct TriangularRef<'a, T: Scalar> {
665    /// Packed reference.
666    packed: PackedRef<'a, T>,
667    /// Diagonal kind.
668    diag: DiagonalKind,
669}
670
671impl<'a, T: Scalar> TriangularRef<'a, T> {
672    /// Creates a new triangular reference.
673    #[inline]
674    pub fn new(packed: PackedRef<'a, T>, diag: DiagonalKind) -> Self {
675        TriangularRef { packed, diag }
676    }
677
678    /// Creates from a slice.
679    #[inline]
680    pub fn from_slice(data: &'a [T], n: usize, uplo: TriangularKind, diag: DiagonalKind) -> Self {
681        TriangularRef {
682            packed: PackedRef::from_slice(data, n, uplo),
683            diag,
684        }
685    }
686
687    /// Returns the dimension.
688    #[inline]
689    pub fn dim(&self) -> usize {
690        self.packed.dim()
691    }
692
693    /// Returns the triangular kind.
694    #[inline]
695    pub fn uplo(&self) -> TriangularKind {
696        self.packed.kind()
697    }
698
699    /// Returns the diagonal kind.
700    #[inline]
701    pub fn diag(&self) -> DiagonalKind {
702        self.diag
703    }
704
705    /// Returns element at (row, col).
706    #[inline]
707    pub fn get(&self, row: usize, col: usize) -> Option<&T> {
708        if self.diag == DiagonalKind::Unit && row == col {
709            return None;
710        }
711        self.packed.get(row, col)
712    }
713
714    /// Returns the underlying packed reference.
715    #[inline]
716    pub fn as_packed(&self) -> PackedRef<'a, T> {
717        self.packed
718    }
719}
720
721/// A mutable reference to triangular packed data.
722pub struct TriangularMut<'a, T: Scalar> {
723    /// Packed mutable reference.
724    packed: PackedMut<'a, T>,
725    /// Diagonal kind.
726    diag: DiagonalKind,
727}
728
729impl<'a, T: Scalar> TriangularMut<'a, T> {
730    /// Creates a new mutable triangular reference.
731    #[inline]
732    pub fn new(packed: PackedMut<'a, T>, diag: DiagonalKind) -> Self {
733        TriangularMut { packed, diag }
734    }
735
736    /// Creates from a mutable slice.
737    #[inline]
738    pub fn from_slice(
739        data: &'a mut [T],
740        n: usize,
741        uplo: TriangularKind,
742        diag: DiagonalKind,
743    ) -> Self {
744        TriangularMut {
745            packed: PackedMut::from_slice(data, n, uplo),
746            diag,
747        }
748    }
749
750    /// Returns the dimension.
751    #[inline]
752    pub fn dim(&self) -> usize {
753        self.packed.dim()
754    }
755
756    /// Returns the triangular kind.
757    #[inline]
758    pub fn uplo(&self) -> TriangularKind {
759        self.packed.kind()
760    }
761
762    /// Returns the diagonal kind.
763    #[inline]
764    pub fn diag(&self) -> DiagonalKind {
765        self.diag
766    }
767
768    /// Returns element at (row, col).
769    #[inline]
770    pub fn get(&self, row: usize, col: usize) -> Option<&T> {
771        if self.diag == DiagonalKind::Unit && row == col {
772            return None;
773        }
774        self.packed.get(row, col)
775    }
776
777    /// Returns mutable element at (row, col).
778    #[inline]
779    pub fn get_mut(&mut self, row: usize, col: usize) -> Option<&mut T> {
780        if self.diag == DiagonalKind::Unit && row == col {
781            return None;
782        }
783        self.packed.get_mut(row, col)
784    }
785
786    /// Sets element at (row, col).
787    #[inline]
788    pub fn set(&mut self, row: usize, col: usize, value: T) {
789        assert!(
790            !(self.diag == DiagonalKind::Unit && row == col),
791            "Cannot set diagonal of unit triangular matrix"
792        );
793        self.packed.set(row, col, value);
794    }
795
796    /// Creates an immutable reborrow.
797    #[inline]
798    pub fn rb(&self) -> TriangularRef<'_, T> {
799        TriangularRef {
800            packed: self.packed.rb(),
801            diag: self.diag,
802        }
803    }
804
805    /// Creates a mutable reborrow.
806    #[inline]
807    pub fn rb_mut(&mut self) -> TriangularMut<'_, T> {
808        TriangularMut {
809            packed: self.packed.rb_mut(),
810            diag: self.diag,
811        }
812    }
813}
814
815#[cfg(test)]
816mod tests {
817    use super::*;
818
819    #[test]
820    fn test_triangular_view_upper() {
821        let mut m: Mat<f64> = Mat::zeros(3, 3);
822        m[(0, 0)] = 1.0;
823        m[(0, 1)] = 2.0;
824        m[(0, 2)] = 3.0;
825        m[(1, 1)] = 4.0;
826        m[(1, 2)] = 5.0;
827        m[(2, 2)] = 6.0;
828        // Set some values below diagonal (should be ignored)
829        m[(1, 0)] = 99.0;
830        m[(2, 0)] = 99.0;
831        m[(2, 1)] = 99.0;
832
833        let tri = TriangularView::new(m.as_ref(), TriangularKind::Upper, DiagonalKind::NonUnit);
834
835        // Upper triangle
836        assert_eq!(tri.get(0, 0), Some(&1.0));
837        assert_eq!(tri.get(0, 1), Some(&2.0));
838        assert_eq!(tri.get(0, 2), Some(&3.0));
839        assert_eq!(tri.get(1, 1), Some(&4.0));
840        assert_eq!(tri.get(1, 2), Some(&5.0));
841        assert_eq!(tri.get(2, 2), Some(&6.0));
842
843        // Lower triangle returns None
844        assert_eq!(tri.get(1, 0), None);
845        assert_eq!(tri.get(2, 0), None);
846        assert_eq!(tri.get(2, 1), None);
847
848        // Test to_dense
849        let dense = tri.to_dense();
850        assert_eq!(dense[(0, 0)], 1.0);
851        assert_eq!(dense[(0, 1)], 2.0);
852        assert_eq!(dense[(1, 0)], 0.0); // Zero, not 99
853        assert_eq!(dense[(2, 1)], 0.0);
854    }
855
856    #[test]
857    fn test_triangular_view_unit() {
858        let mut m: Mat<f64> = Mat::zeros(3, 3);
859        m[(0, 0)] = 99.0; // Should be treated as 1
860        m[(0, 1)] = 2.0;
861        m[(0, 2)] = 3.0;
862        m[(1, 1)] = 99.0; // Should be treated as 1
863        m[(1, 2)] = 4.0;
864        m[(2, 2)] = 99.0; // Should be treated as 1
865
866        let tri = TriangularView::new(m.as_ref(), TriangularKind::Upper, DiagonalKind::Unit);
867
868        // Diagonal returns None (implicitly one)
869        assert_eq!(tri.get(0, 0), None);
870        assert_eq!(tri.get(1, 1), None);
871        assert_eq!(tri.get(2, 2), None);
872
873        // Off-diagonal upper elements
874        assert_eq!(tri.get(0, 1), Some(&2.0));
875        assert_eq!(tri.get(0, 2), Some(&3.0));
876        assert_eq!(tri.get(1, 2), Some(&4.0));
877
878        // to_dense should have ones on diagonal
879        let dense = tri.to_dense();
880        assert_eq!(dense[(0, 0)], 1.0);
881        assert_eq!(dense[(1, 1)], 1.0);
882        assert_eq!(dense[(2, 2)], 1.0);
883        assert_eq!(dense[(0, 1)], 2.0);
884    }
885
886    #[test]
887    fn test_triangular_view_lower() {
888        let mut m: Mat<f64> = Mat::zeros(3, 3);
889        m[(0, 0)] = 1.0;
890        m[(1, 0)] = 2.0;
891        m[(1, 1)] = 3.0;
892        m[(2, 0)] = 4.0;
893        m[(2, 1)] = 5.0;
894        m[(2, 2)] = 6.0;
895
896        let tri = TriangularView::new(m.as_ref(), TriangularKind::Lower, DiagonalKind::NonUnit);
897
898        // Lower triangle
899        assert_eq!(tri.get(0, 0), Some(&1.0));
900        assert_eq!(tri.get(1, 0), Some(&2.0));
901        assert_eq!(tri.get(1, 1), Some(&3.0));
902        assert_eq!(tri.get(2, 0), Some(&4.0));
903        assert_eq!(tri.get(2, 1), Some(&5.0));
904        assert_eq!(tri.get(2, 2), Some(&6.0));
905
906        // Upper triangle returns None
907        assert_eq!(tri.get(0, 1), None);
908        assert_eq!(tri.get(0, 2), None);
909        assert_eq!(tri.get(1, 2), None);
910    }
911
912    #[test]
913    fn test_triangular_mat_packed() {
914        let mut tri: TriangularMat<f64> =
915            TriangularMat::zeros(3, TriangularKind::Upper, DiagonalKind::NonUnit);
916
917        tri.set(0, 0, 1.0);
918        tri.set(0, 1, 2.0);
919        tri.set(0, 2, 3.0);
920        tri.set(1, 1, 4.0);
921        tri.set(1, 2, 5.0);
922        tri.set(2, 2, 6.0);
923
924        assert_eq!(tri.get(0, 0), Some(&1.0));
925        assert_eq!(tri.get(0, 1), Some(&2.0));
926        assert_eq!(tri.get(1, 2), Some(&5.0));
927
928        // Cannot access lower triangle
929        assert_eq!(tri.get(1, 0), None);
930
931        let diag = tri.diagonal();
932        assert_eq!(diag, vec![1.0, 4.0, 6.0]);
933    }
934
935    #[test]
936    fn test_triangular_mat_unit() {
937        let mut tri: TriangularMat<f64> = TriangularMat::unit_zeros(3, TriangularKind::Lower);
938
939        // Off-diagonal elements
940        tri.set(1, 0, 2.0);
941        tri.set(2, 0, 3.0);
942        tri.set(2, 1, 4.0);
943
944        // Diagonal returns None (implicit one)
945        assert_eq!(tri.get(0, 0), None);
946        assert_eq!(tri.get(1, 1), None);
947        assert_eq!(tri.get(2, 2), None);
948
949        // Off-diagonal elements
950        assert_eq!(tri.get(1, 0), Some(&2.0));
951        assert_eq!(tri.get(2, 0), Some(&3.0));
952        assert_eq!(tri.get(2, 1), Some(&4.0));
953
954        // Diagonal should be ones
955        let diag = tri.diagonal();
956        assert_eq!(diag, vec![1.0, 1.0, 1.0]);
957
958        // to_dense
959        let dense = tri.to_dense();
960        assert_eq!(dense[(0, 0)], 1.0);
961        assert_eq!(dense[(1, 1)], 1.0);
962        assert_eq!(dense[(2, 2)], 1.0);
963        assert_eq!(dense[(1, 0)], 2.0);
964        assert_eq!(dense[(0, 1)], 0.0);
965    }
966
967    #[test]
968    fn test_triangular_mat_transpose() {
969        let mut tri: TriangularMat<f64> =
970            TriangularMat::zeros(3, TriangularKind::Upper, DiagonalKind::NonUnit);
971
972        tri.set(0, 0, 1.0);
973        tri.set(0, 1, 2.0);
974        tri.set(0, 2, 3.0);
975        tri.set(1, 1, 4.0);
976        tri.set(1, 2, 5.0);
977        tri.set(2, 2, 6.0);
978
979        let tri_t = tri.transpose();
980        assert_eq!(tri_t.uplo(), TriangularKind::Lower);
981
982        // Elements should be transposed
983        assert_eq!(tri_t.get(0, 0), Some(&1.0));
984        assert_eq!(tri_t.get(1, 0), Some(&2.0)); // Was (0, 1)
985        assert_eq!(tri_t.get(2, 0), Some(&3.0)); // Was (0, 2)
986        assert_eq!(tri_t.get(1, 1), Some(&4.0));
987        assert_eq!(tri_t.get(2, 1), Some(&5.0)); // Was (1, 2)
988        assert_eq!(tri_t.get(2, 2), Some(&6.0));
989    }
990
991    #[test]
992    fn test_triangular_view_mut() {
993        let mut m: Mat<f64> = Mat::zeros(3, 3);
994        let mut view =
995            TriangularViewMut::new(m.as_mut(), TriangularKind::Upper, DiagonalKind::NonUnit);
996
997        view.set(0, 0, 1.0);
998        view.set(0, 1, 2.0);
999        view.set(1, 1, 3.0);
1000
1001        assert_eq!(view.get(0, 0), Some(&1.0));
1002        assert_eq!(view.get(0, 1), Some(&2.0));
1003        assert_eq!(view.get(1, 1), Some(&3.0));
1004
1005        // zero_non_triangle
1006        view.zero_non_triangle();
1007        // Check that lower triangle is zeroed (it was already zero)
1008
1009        // Fill
1010        view.fill(5.0);
1011        assert_eq!(view.get(0, 0), Some(&5.0));
1012        assert_eq!(view.get(0, 1), Some(&5.0));
1013    }
1014
1015    #[test]
1016    fn test_triangular_from_dense() {
1017        let m = Mat::from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0], &[7.0, 8.0, 9.0]]);
1018
1019        let tri =
1020            TriangularMat::from_dense(&m.as_ref(), TriangularKind::Upper, DiagonalKind::NonUnit);
1021
1022        assert_eq!(tri.get(0, 0), Some(&1.0));
1023        assert_eq!(tri.get(0, 1), Some(&2.0));
1024        assert_eq!(tri.get(0, 2), Some(&3.0));
1025        assert_eq!(tri.get(1, 1), Some(&5.0));
1026        assert_eq!(tri.get(1, 2), Some(&6.0));
1027        assert_eq!(tri.get(2, 2), Some(&9.0));
1028
1029        // Lower part is not stored
1030        assert_eq!(tri.get(1, 0), None);
1031        assert_eq!(tri.get(2, 0), None);
1032        assert_eq!(tri.get(2, 1), None);
1033    }
1034
1035    #[test]
1036    fn test_triangular_ref_mut() {
1037        let mut data = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
1038        let mut tmut =
1039            TriangularMut::from_slice(&mut data, 3, TriangularKind::Upper, DiagonalKind::NonUnit);
1040
1041        assert_eq!(tmut.get(0, 0), Some(&1.0));
1042        assert_eq!(tmut.get(0, 1), Some(&2.0));
1043
1044        tmut.set(0, 0, 10.0);
1045        assert_eq!(tmut.get(0, 0), Some(&10.0));
1046        assert_eq!(data[0], 10.0);
1047    }
1048}
oxiblas_matrix/triangular.rs

oxiblas_matrix/
triangular.rs
