rlst 0.8.0

A Rust native linear algebra library.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

//! Implement the triangular solver trait.

use crate::{
    UnsafeRandom1DAccessByValue,
    base_types::{RlstResult, UpLo},
    dense::{
        array::{Array, DynArray},
        linalg::lapack::interface::trsm::{TrsmDiag, TrsmSide, TrsmTransA, TrsmUplo},
    },
    traits::{
        base_operations::Shape,
        linalg::{decompositions::SolveTriangular, lapack::Lapack},
    },
};

impl<Item, ArrayImpl, RhsArrayImpl, const NDIM: usize> SolveTriangular<Array<RhsArrayImpl, NDIM>>
    for Array<ArrayImpl, 2>
where
    Item: Lapack,
    ArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<2>,
    RhsArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<NDIM>,
{
    type Output = DynArray<Item, NDIM>;

    fn solve_triangular(
        &self,
        uplo: UpLo,
        rhs: &Array<RhsArrayImpl, NDIM>,
    ) -> RlstResult<Self::Output> {
        let [m, n] = self.shape();

        assert_eq!(
            m, n,
            "Triangular solver requires a square matrix, but got shape: {m}x{n}"
        );

        assert_eq!(
            m,
            rhs.shape()[0],
            "Left-hand side and right-hand side must have the same number of rows. {} != {}",
            m,
            rhs.shape()[0]
        );

        let nrhs = if NDIM == 1 {
            1
        } else if NDIM == 2 {
            rhs.shape()[1]
        } else {
            panic!("The right-hand side must be one or two dimensional: NDIM = {NDIM}");
        };

        let mut a = DynArray::new_from(self);
        let mut b = DynArray::new_from(rhs);

        let uplo = match uplo {
            UpLo::Upper => TrsmUplo::Upper,
            UpLo::Lower => TrsmUplo::Lower,
        };

        Item::trsm(
            TrsmSide::Left,
            uplo,
            TrsmTransA::NoTrans,
            TrsmDiag::NonUnit,
            m,
            nrhs,
            Item::one(),
            a.data_mut().unwrap(),
            m,
            b.data_mut().unwrap(),
            m,
        )?;

        Ok(b)
    }
}

#[cfg(test)]
mod test {

    use super::*;
    use crate::base_types::{c32, c64};
    use crate::dense::array::DynArray;
    use crate::dot_blas;
    use num::Zero;
    use paste::paste;

    macro_rules! implement_triangular_solve_test {
        ($scalar:ty, $tol:expr) => {
            paste! {

            #[test]
            fn [<test_solve_triangular_$scalar>]() {
                let n = 10;
                let mut a = DynArray::<$scalar, 2>::from_shape([n, n]);
                a.fill_from_seed_equally_distributed(0);

                for row in 0..n {
                    for col in 1 + row..n {
                        a[[row, col]] = <$scalar>::zero(); // Make it lower triangular
                    }
                }

                let mut x_actual = DynArray::<$scalar, 2>::from_shape([n, 1]);
                x_actual.fill_from_seed_equally_distributed(1);

                let b = dot_blas!(a.r(), x_actual.r());

                let x = a.solve_triangular(UpLo::Lower, &b).unwrap();

                crate::assert_array_relative_eq!(x_actual, x, $tol);
            }


                    }
        };
    }

    implement_triangular_solve_test!(f32, 5E-3);
    implement_triangular_solve_test!(f64, 1E-10);
    implement_triangular_solve_test!(c32, 5E-3);
    implement_triangular_solve_test!(c64, 1E-10);
}