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

oxiblas_matrix/
triangular.rs
