1use format::{Compressed, Conventional, Diagonal};
2use operation::{Multiply, MultiplyInto, MultiplySelf, Transpose};
3use {Element, Number};
4
5impl<T> Multiply<Diagonal<T>, Compressed<T>> for Compressed<T>
6where
7 T: Element + Number,
8{
9 #[inline]
10 fn multiply(&self, right: &Diagonal<T>) -> Self {
11 let mut result = self.clone();
12 result.multiply_self(right);
13 result
14 }
15}
16
17impl<'l, T> MultiplyInto<[T], [T]> for Compressed<T>
18where
19 T: Element + Number,
20{
21 #[inline]
22 fn multiply_into(&self, right: &[T], result: &mut [T]) {
23 let (m, p) = (self.rows, self.columns);
24 let n = right.len() / p;
25 multiply_matrix_left(self, right, result, m, p, n)
26 }
27}
28
29impl<'l, T> MultiplyInto<Compressed<T>, [T]> for Conventional<T>
30where
31 T: Element + Number,
32{
33 #[inline]
34 fn multiply_into(&self, right: &Compressed<T>, result: &mut [T]) {
35 let (m, p, n) = (self.rows, self.columns, right.columns);
36 multiply_matrix_right(&self.values, right, result, m, p, n)
37 }
38}
39
40impl<T> MultiplySelf<Diagonal<T>> for Compressed<T>
41where
42 T: Element + Number,
43{
44 fn multiply_self(&mut self, right: &Diagonal<T>) {
45 let (m, n) = (self.rows, right.columns);
46 debug_assert_eq!(self.columns, right.rows);
47 self.resize((m, n));
48 for (_, j, value) in self.iter_mut() {
49 *value = *value * right[j];
50 }
51 }
52}
53
54impl<T: Element> Transpose for Compressed<T> {
55 fn transpose(&self) -> Self {
56 let &Compressed {
57 rows,
58 columns,
59 nonzeros,
60 variant,
61 ..
62 } = self;
63 let mut matrix = Compressed::with_capacity((columns, rows), variant, nonzeros);
64 for (i, j, &value) in self.iter() {
65 matrix.set((j, i), value);
66 }
67 matrix
68 }
69}
70
71fn multiply_matrix_left<T>(a: &Compressed<T>, b: &[T], c: &mut [T], m: usize, p: usize, n: usize)
72where
73 T: Element + Number,
74{
75 debug_assert_eq!(a.rows * a.columns, m * p);
76 debug_assert_eq!(b.len(), p * n);
77 debug_assert_eq!(c.len(), m * n);
78 let &Compressed {
79 ref values,
80 ref indices,
81 ref offsets,
82 ..
83 } = a;
84 for j in 0..n {
85 let bo = j * p;
86 let co = j * m;
87 for l in 0..p {
88 let bi = bo + l;
89 for k in offsets[l]..offsets[l + 1] {
90 let i = co + indices[k];
91 c[i] = c[i] + values[k] * b[bi];
92 }
93 }
94 }
95}
96
97fn multiply_matrix_right<T>(a: &[T], b: &Compressed<T>, c: &mut [T], m: usize, p: usize, n: usize)
98where
99 T: Element + Number,
100{
101 debug_assert_eq!(a.len(), m * p);
102 debug_assert_eq!(b.rows * b.columns, p * n);
103 debug_assert_eq!(c.len(), m * n);
104 let &Compressed {
105 ref values,
106 ref indices,
107 ref offsets,
108 ..
109 } = b;
110 for j in 0..n {
111 let co = j * m;
112 for k in offsets[j]..offsets[j + 1] {
113 let ao = indices[k] * m;
114 for i in 0..m {
115 c[co + i] = c[co + i] + values[k] * a[ao + i];
116 }
117 }
118 }
119}
120
121#[cfg(test)]
122mod tests {
123 use format::compressed::Variant;
124 use prelude::*;
125
126 #[test]
127 fn multiply_self() {
128 let mut matrix = new!(
129 3,
130 2,
131 3,
132 Variant::Column,
133 vec![1.0, 2.0, 3.0],
134 vec![1, 0, 2],
135 vec![0, 1, 3]
136 );
137 let right = Diagonal::from_vec((2, 4), vec![4.0, 5.0]);
138 matrix.multiply_self(&right);
139 assert_eq!(
140 matrix,
141 new!(
142 3,
143 4,
144 3,
145 Variant::Column,
146 vec![4.0, 10.0, 15.0],
147 vec![1, 0, 2],
148 vec![0, 1, 3, 3, 3]
149 )
150 );
151 }
152
153 #[test]
154 fn multiply_into_left() {
155 let matrix = Compressed::from(Conventional::from_vec(
156 (4, 3),
157 matrix![
158 1.0, 5.0, 4.0;
159 2.0, 6.0, 3.0;
160 3.0, 6.0, 2.0;
161 4.0, 5.0, 1.0;
162 ],
163 ));
164 let right = Conventional::from_vec(
165 (3, 2),
166 matrix![
167 1.0, 4.0;
168 2.0, 5.0;
169 3.0, 6.0;
170 ],
171 );
172 let mut result = Conventional::from_vec(
173 (4, 2),
174 matrix![
175 1.0, 1.0;
176 1.0, 1.0;
177 1.0, 1.0;
178 1.0, 1.0;
179 ],
180 );
181 matrix.multiply_into(&right, &mut result);
182 assert_eq!(
183 &result.values,
184 &matrix![
185 24.0, 54.0;
186 24.0, 57.0;
187 22.0, 55.0;
188 18.0, 48.0;
189 ]
190 );
191 }
192
193 #[test]
194 fn multiply_into_right() {
195 let matrix = Conventional::from_vec(
196 (4, 3),
197 matrix![
198 1.0, 5.0, 4.0;
199 2.0, 6.0, 3.0;
200 3.0, 6.0, 2.0;
201 4.0, 5.0, 1.0;
202 ],
203 );
204 let right = Compressed::from(Conventional::from_vec(
205 (3, 2),
206 matrix![
207 1.0, 4.0;
208 2.0, 5.0;
209 3.0, 6.0;
210 ],
211 ));
212 let mut result = Conventional::from_vec(
213 (4, 2),
214 matrix![
215 1.0, 1.0;
216 1.0, 1.0;
217 1.0, 1.0;
218 1.0, 1.0;
219 ],
220 );
221 matrix.multiply_into(&right, &mut result);
222 assert_eq!(
223 &result.values,
224 &matrix![
225 24.0, 54.0;
226 24.0, 57.0;
227 22.0, 55.0;
228 18.0, 48.0;
229 ]
230 );
231 }
232
233 #[test]
234 fn transpose() {
235 let matrix = new!(
236 5,
237 7,
238 5,
239 Variant::Column,
240 vec![1.0, 2.0, 3.0, 4.0, 5.0],
241 vec![1, 0, 3, 1, 4],
242 vec![0, 0, 0, 1, 2, 2, 3, 5]
243 );
244 let matrix = matrix.transpose();
245 assert_eq!(
246 matrix,
247 new!(
248 7,
249 5,
250 5,
251 Variant::Column,
252 vec![2.0, 1.0, 4.0, 3.0, 5.0],
253 vec![3, 2, 6, 5, 6],
254 vec![0, 1, 3, 3, 4, 5]
255 )
256 );
257 }
258}