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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139

use crate::algebra::linear::{Vector, Matrix};
use super::Substitute;
use crate::algebra::abstr::{Field, Scalar, AbsDiffEq};

impl<T> Substitute<Vector<T>> for Matrix<T> where T: Field + Scalar + AbsDiffEq
{
    fn substitute_forward(&self, a: Vector<T>) -> Result<Vector<T>, ()>
    {
        let mut b: Vector<T> = a;
        for k in 0..self.n
        {

            if self[[k, k]].abs_diff_eq(&T::zero(), T::default_epsilon())
            {
                // free variable
                // if b.get(k).abs_diff_eq(&T::zero(), T::default_epsilon())
                // {
                //     *b.get_mut(k) = T::one();
                // }
                // else
                // {
                //     return Err(())
                // }
            }
            else
            {
                for l in 0..k
                {
                    b[k] = b[k] - self[[k, l]] * b[l];
                }
                b[k] /= self[[k, k]];
            }
        }

        Ok(b)
    }

    fn substitute_backward(&self, c: Vector<T>) -> Result<Vector<T>, ()>
    {
        let mut b: Vector<T> = c;

        for k in (0..self.n).rev()
        {
            if self[[k, k]].abs_diff_eq(&T::zero(), T::default_epsilon())
            {
                // free variable
                // if b.get(k).abs_diff_eq(&T::zero(), T::default_epsilon())
                // {
                //     *b.get_mut(k) = T::one();
                // }
                // else
                // {
                //     return Err(());
                // }
            }
            else
            {
                for l in (k + 1..self.n).rev()
                {
                    b[k] = b[k] - self[[k, l]] * b[l];
                }
                b[k] /= self[[k, k]];
            }
        }

        Ok(b)
    }
}

impl<T> Substitute<Matrix<T>> for Matrix<T> where T: Field + Scalar + AbsDiffEq
{
    fn substitute_forward(&self, a: Matrix<T>) -> Result<Matrix<T>, ()>
    {
        let mut b: Matrix<T> = a;
        let min: usize = std::cmp::min( self.m, self.n);
        for k in 0..min
        {
            if self[[k, k]].abs_diff_eq(&T::zero(), T::default_epsilon())
            {
                //free variable
                // if b.get(k).abs_diff_eq(&T::zero(), T::default_epsilon())
                // {
                //     *b.get_mut(k) = T::one();
                // }
                // else
                // {
                //     return Err(());
                // }
            }
            else
            {
                for l in 0..k
                {
                    b.set_row( & (b.get_row(k) - (b.get_row(l) * self[[k, l]])), k);
                }
                b.set_row( & (b.get_row(k) / self[[k, k]]), k);
            }

        }

        Ok(b)
    }

    fn substitute_backward(&self, a: Matrix<T>) -> Result<Matrix<T>, ()>
    {
        let mut b: Matrix<T> = a;
        let min = std::cmp::min(self.m, self.n);

        for k in (0..min).rev()
        {

            if self[[k, k]].abs_diff_eq(&T::zero(), T::default_epsilon())
            {
                //free variable
                // if b.get(k).abs_diff_eq(&T::zero(), T::default_epsilon())
                // {
                //     for m in 0..min {
                //         *b.get_mut(k, m) = T::one();
                //     }
                // }
                // else
                // {
                //     return Err(())
                // }
            }
            else
            {
                for l in (k + 1..min).rev()
                {
                    b.set_row(&(b.get_row(k) - (b.get_row(l) * self[[k, l]])), k);
                }

                b.set_row(&(b.get_row(k) / self[[k, k]]), k);
            }
        }

        Ok(b)
    }
}