maryada 0.2.0

No-std binary64 real interval arithmetic conforming to IEEE 1788.1
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

#![allow(clippy::arithmetic_side_effects, clippy::indexing_slicing)]

use nalgebra::{
    Const, DefaultAllocator, Dim, DimMin, Scalar, Storage,
    allocator::Allocator,
    constraint::{AreMultipliable, ShapeConstraint},
};

use crate::{EpsilonInflationSolver, IntervalOps};

use super::{IntervalMatrix, OIntervalMatrix, OIntervalVector, Solver};

/// Hansen–Bliek–Rohn verified solver for interval systems with an H-matrix.
pub struct HansenBliekRohnSolver;
impl<T, D> Solver<T, D> for HansenBliekRohnSolver
where
    T: IntervalOps + Scalar,
    D: Dim + DimMin<D, Output = D>,
    DefaultAllocator: Allocator<D, D> + Allocator<D>,
    ShapeConstraint: AreMultipliable<D, D, D, D> + AreMultipliable<D, D, D, Const<1>>,
{
    fn solve<SA, SB>(
        &self,
        lhs: &IntervalMatrix<T, D, D, SA>,
        rhs: &IntervalMatrix<T, D, Const<1>, SB>,
    ) -> Option<OIntervalVector<T, D>>
    where
        SA: Storage<T, D, D>,
        SB: Storage<T, D, Const<1>>,
    {
        if !lhs.is_square()
            || lhs.is_empty()
            || lhs.nrows() != rhs.nrows()
            || lhs.has_empty_entries()
            || lhs.has_nai_entries()
            || rhs.has_empty_entries()
            || rhs.has_nai_entries()
            || !lhs.is_h_matrix()
        {
            return None;
        }

        let a_comp = lhs.comparison_matrix()?;
        let (dim, _) = rhs.shape_generic();
        let a_comp_interval = OIntervalMatrix::from_singletons(&a_comp);
        let identity = OIntervalMatrix::<T, D, D>::identity_generic(dim);
        let a_comp_inv =
            EpsilonInflationSolver::default().solve_matrix(&a_comp_interval, &identity)?;

        let b_mag = rhs.mag();
        let b_mag = OIntervalVector::from_singletons(&b_mag);
        let u = &a_comp_inv * &b_mag;
        let d = a_comp_inv.diagonal();
        if d.any(|entry| entry.inf() <= 0.0) {
            return None;
        }

        let a_diag = a_comp_interval.diagonal();
        let d_recip = d.map(IntervalOps::recip);
        let alpha_interval = &a_diag - &d_recip;
        let alpha = alpha_interval.mag();
        let scaled_u = u.component_div(&d);
        let beta_interval = &scaled_u - &b_mag;
        let beta = beta_interval.mag();
        if alpha
            .iter()
            .chain(beta.iter())
            .any(|entry| !entry.is_finite())
        {
            return None;
        }

        for i in 0..lhs.nrows() {
            let denominator = lhs[(i, i)] + T::new(-alpha[i], alpha[i]);
            if denominator.mig() == 0.0 {
                return None;
            }
        }

        let x = OIntervalVector::from_fn_generic(dim, Const::<1>, |i, _| {
            (rhs[i] + T::new(-beta[i], beta[i])) / (lhs[(i, i)] + T::new(-alpha[i], alpha[i]))
        });
        Some(x)
    }
}