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
140
141
142
143
144
145
146
147
148
149
150

use matrix::{BaseMatrix, BaseMatrixMut};
use libnum::{Zero, Num};
use utils::in_place_vec_bin_op;

pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M)
    where T: Zero, M: BaseMatrixMut<T> {
    for (i, mut row) in matrix.row_iter_mut().enumerate() {
        for element in row.raw_slice_mut().iter_mut().take(i) {
            *element = T::zero();
        }
    }
}

pub fn nullify_upper_triangular_part<T, M>(matrix: &mut M)
    where T: Zero, M: BaseMatrixMut<T> {
    for (i, mut row) in matrix.row_iter_mut().enumerate() {
        for element in row.raw_slice_mut().iter_mut().skip(i + 1) {
            *element = T::zero();
        }
    }
}

/// Given a vector `x` and a `m x n` matrix `A`, compute
/// `y = A^T x`.
///
/// This is a stopgap solution until we have a more proper
/// BLIS/BLAS-like API.
pub fn transpose_gemv<T, M>(a: &M, x: &[T], y: &mut [T])
    where M: BaseMatrix<T>, T: Num + Copy
{
    let m = a.rows();
    let n = a.cols();

    assert!(x.len() == m, "A and x must be dimensionally compatible.");
    assert!(y.len() == n, "A and y must be dimensionally compatible.");

    for element in y.iter_mut() {
        *element = T::zero();
    }

    for j in 0 .. m {
        let a_j = a.row(j).raw_slice();
        axpy(x[j], a_j, y);
    }
}

// Given scalar `a` and vectors `x` and `y` of same length, computes the
// scalar-vector product `y <- y + a * x`.
//
/// This function is a stopgap solution until we have a more proper
/// BLIS/BLAS-like API.
pub fn axpy<T>(a: T, x: &[T], y: &mut [T])
    where T: Num + Copy
{
    assert!(x.len() == y.len());
    in_place_vec_bin_op(y, x, |y, x| {
        *y = y.clone() + a * x.clone();
    });
}

/// Given a `m x n` matrix `A` and vectors `x` and `y` and
/// a scalar `alpha`, perform the rank-1 update
/// `A <- A + alpha * x y^T`.
///
/// This function is a stopgap solution until we have a more proper
/// BLIS/BLAS-like API.
pub fn ger<T, M>(a: &mut M, alpha: T, x: &[T], y: &[T])
    where M: BaseMatrixMut<T>, T: Num + Copy
{
    let m = a.rows();
    let n = a.cols();

    assert!(x.len() == m);
    assert!(y.len() == n);

    for i in 0 .. m {
        let mut a_i = a.row_mut(i).raw_slice_mut();
        // Let a_i be the ith row. Then
        // a_i <- a_i + alpha * x[i] * y
        // is just an axpy operation
        axpy(alpha * x[i], &y, &mut a_i);
    }
}

#[cfg(test)]
mod tests {
    use vector::Vector;

    use super::nullify_lower_triangular_part;
    use super::nullify_upper_triangular_part;
    use super::transpose_gemv;
    use super::ger;

    #[test]
    fn nullify_lower_triangular_part_examples() {
        let mut x = matrix![1.0, 2.0, 3.0;
                            4.0, 5.0, 6.0;
                            7.0, 8.0, 9.0];
        nullify_lower_triangular_part(&mut x);
        assert_matrix_eq!(x, matrix![
            1.0, 2.0, 3.0;
            0.0, 5.0, 6.0;
            0.0, 0.0, 9.0
        ]);
    }

    #[test]
    fn nullify_upper_triangular_part_examples() {
        let mut x = matrix![1.0, 2.0, 3.0;
                            4.0, 5.0, 6.0;
                            7.0, 8.0, 9.0];
        nullify_upper_triangular_part(&mut x);
        assert_matrix_eq!(x, matrix![
            1.0, 0.0, 0.0;
            4.0, 5.0, 0.0;
            7.0, 8.0, 9.0
        ]);
    }

    #[test]
    fn transpose_gemv_examples() {
        {
            let a = matrix![3.0, 4.0, 5.0;
                            2.0, 3.0, 1.0];
            let x = vec![2.0, 3.0];
            let mut y = vec![0.0; 3];
            transpose_gemv(&a, &x, &mut y);

            let y = Vector::new(y);
            assert_vector_eq!(y, vector![12.0, 17.0, 13.0]);
        }
    }

    #[test]
    fn ger_examples() {
        {
            let mut a = matrix![3.0, 4.0, 5.0;
                            2.0, 3.0, 1.0];
            let x = vec![3.0, 4.0];
            let y = vec![2.0, 1.0, 3.0];
            let alpha = 3.0;

            ger(&mut a, alpha, &x, &y);

            let expected = matrix![21.0, 13.0, 32.0;
                                   26.0, 15.0, 37.0];
            assert_matrix_eq!(a, expected, comp = float);
        }
    }
}