rust_blas/vector/
ops.rs

1// Copyright 2014 Michael Yang. All rights reserved.
2// Use of this source code is governed by a MIT-style
3// license that can be found in the LICENSE file.
4
5//! Wrappers for vector functions.
6
7use crate::default::Default;
8use crate::matrix::Matrix;
9use crate::pointer::CPtr;
10use crate::scalar::Scalar;
11use crate::vector::ll::*;
12use crate::vector::Vector;
13use num_complex::{Complex, Complex32, Complex64};
14use std::cmp;
15
16pub trait Copy: Sized {
17    /// Copies `src.len()` elements of `src` into `dst`.
18    fn copy<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(src: &V, dst: &mut W);
19    /// Copies the entire matrix `dst` into `src`.
20    fn copy_mat(src: &dyn Matrix<Self>, dst: &mut dyn Matrix<Self>);
21}
22
23macro_rules! copy_impl(($($t: ident), +) => (
24    $(
25        impl Copy for $t {
26            fn copy<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(src: &V, dst: &mut W) {
27                unsafe {
28                    prefix!($t, copy)(dst.len(),
29                        src.as_ptr().as_c_ptr(),  src.inc(),
30                        dst.as_mut_ptr().as_c_ptr(), dst.inc());
31                }
32            }
33
34            fn copy_mat(src: &dyn Matrix<Self>, dst: &mut dyn Matrix<Self>) {
35                let len = dst.rows() * dst.cols();
36
37                unsafe {
38                    prefix!($t, copy)(len,
39                        src.as_ptr().as_c_ptr(),  1,
40                        dst.as_mut_ptr().as_c_ptr(), 1);
41                }
42            }
43        }
44    )+
45));
46
47copy_impl!(f32, f64, Complex32, Complex64);
48
49/// Computes `a * x + y` and stores the result in `y`.
50pub trait Axpy: Sized {
51    fn axpy<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(alpha: &Self, x: &V, y: &mut W);
52    fn axpy_mat(alpha: &Self, x: &dyn Matrix<Self>, y: &mut dyn Matrix<Self>);
53}
54
55macro_rules! axpy_impl(($($t: ident), +) => (
56    $(
57        impl Axpy for $t {
58            fn axpy<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(alpha: &$t, x: &V, y: &mut W) {
59                unsafe {
60                    let n = cmp::min(x.len(), y.len());
61
62                    prefix!($t, axpy)(n,
63                        alpha.as_const(),
64                        x.as_ptr().as_c_ptr(), x.inc(),
65                        y.as_mut_ptr().as_c_ptr(), y.inc());
66                }
67            }
68
69            fn axpy_mat(alpha: &$t, x: &dyn Matrix<$t>, y: &mut dyn Matrix<$t>) {
70                unsafe {
71                    let x_len = x.rows() * x.cols();
72                    let y_len = y.rows() * y.cols();
73                    let n = cmp::min(x_len, y_len);
74
75                    prefix!($t, axpy)(n,
76                        alpha.as_const(),
77                        x.as_ptr().as_c_ptr(), 1,
78                        y.as_mut_ptr().as_c_ptr(), 1);
79                }
80            }
81        }
82    )+
83));
84
85axpy_impl!(f32, f64, Complex32, Complex64);
86
87#[cfg(test)]
88mod axpy_tests {
89    use crate::vector::ops::Axpy;
90    use num_complex::Complex;
91
92    #[test]
93    fn real() {
94        let x = vec![1f32, -2f32, 3f32, 4f32];
95        let y = vec![3f32, 7f32, -2f32, 2f32];
96        let mut z = y.clone();
97
98        Axpy::axpy(&1f32, &y, &mut z);
99        Axpy::axpy(&1f32, &x, &mut z);
100        assert_eq!(z, vec![7f32, 12f32, -1f32, 8f32]);
101    }
102
103    #[test]
104    fn slice() {
105        let x = vec![1f32, -2f32, 3f32, 4f32, 5f32];
106        let y = vec![3f32, 7f32, -2f32, 2f32];
107        let mut z = y.clone();
108
109        Axpy::axpy(&1f32, &y, &mut z);
110        Axpy::axpy(&1f32, &x[..4], &mut z);
111        assert_eq!(z, vec![7f32, 12f32, -1f32, 8f32]);
112    }
113
114    #[test]
115    fn complex() {
116        let x = vec![Complex::new(1f32, 1f32), Complex::new(1f32, 3f32)];
117        let y = vec![Complex::new(3f32, -2f32), Complex::new(2f32, 3f32)];
118        let mut z = x.clone();
119
120        Axpy::axpy(&Complex::new(-1f32, 1f32), &y, &mut z);
121        assert_eq!(z, vec![Complex::new(0f32, 6f32), Complex::new(-4f32, 2f32)]);
122    }
123}
124
125/// Computes `a * x` and stores the result in `x`.
126pub trait Scal: Sized {
127    fn scal<V: ?Sized + Vector<Self>>(alpha: &Self, x: &mut V);
128    fn scal_mat(alpha: &Self, x: &mut dyn Matrix<Self>);
129}
130
131macro_rules! scal_impl(($($t: ident), +) => (
132    $(
133        impl Scal for $t {
134            #[inline]
135            fn scal<V: ?Sized + Vector<Self>>(alpha: &$t, x: &mut V) {
136                unsafe {
137                    prefix!($t, scal)(x.len(),
138                        alpha.as_const(),
139                        x.as_mut_ptr().as_c_ptr(), x.inc());
140                }
141            }
142
143            fn scal_mat(alpha: &$t, x: &mut dyn Matrix<$t>) {
144                unsafe {
145                    prefix!($t, scal)(x.rows() * x.cols(),
146                        alpha.as_const(),
147                        x.as_mut_ptr().as_c_ptr(), 1);
148                }
149            }
150        }
151    )+
152));
153
154scal_impl!(f32, f64, Complex32, Complex64);
155
156#[cfg(test)]
157mod scal_tests {
158    use crate::vector::ops::Scal;
159    use num_complex::Complex;
160
161    #[test]
162    fn real() {
163        let mut x = vec![1f32, -2f32, 3f32, 4f32];
164
165        Scal::scal(&-2f32, &mut x);
166        assert_eq!(x, vec![-2f32, 4f32, -6f32, -8f32]);
167    }
168
169    #[test]
170    fn slice() {
171        let mut x = vec![1f32, -2f32, 3f32, 4f32];
172
173        Scal::scal(&-2f32, &mut x[..3]);
174        assert_eq!(x, vec![-2f32, 4f32, -6f32, 4f32]);
175    }
176
177    #[test]
178    fn complex() {
179        let mut x = vec![Complex::new(1f32, 1f32), Complex::new(1f32, 3f32)];
180
181        Scal::scal(&Complex::new(1f32, 1f32), &mut x);
182        assert_eq!(x, vec![Complex::new(0f32, 2f32), Complex::new(-2f32, 4f32)]);
183    }
184
185    #[test]
186    fn complex_real() {
187        let mut x = vec![Complex::new(1f32, 1f32), Complex::new(1f32, 3f32)];
188
189        Scal::scal(&Complex::new(2f32, 0f32), &mut x);
190        assert_eq!(x, vec![Complex::new(2f32, 2f32), Complex::new(2f32, 6f32)]);
191    }
192}
193
194/// Swaps the content of `x` and `y`.
195pub trait Swap: Sized {
196    /// If they are different lengths, the shorter length is used.
197    fn swap<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &mut V, y: &mut W);
198}
199
200macro_rules! swap_impl(($($t: ident), +) => (
201    $(
202        impl Swap for $t {
203            fn swap<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &mut V, y: &mut W) {
204                unsafe {
205                    let n = cmp::min(x.len(), y.len());
206
207                    prefix!($t, swap)(n,
208                        x.as_mut_ptr().as_c_ptr(), x.inc(),
209                        y.as_mut_ptr().as_c_ptr(), y.inc());
210                }
211            }
212        }
213    )+
214));
215
216swap_impl!(f32, f64, Complex32, Complex64);
217
218#[cfg(test)]
219mod swap_tests {
220    use crate::vector::ops::Swap;
221    use num_complex::Complex;
222
223    #[test]
224    fn real() {
225        let mut x = vec![1f32, -2f32, 3f32, 4f32];
226        let mut y = vec![2f32, -3f32, 4f32, 1f32];
227        let xr = y.clone();
228        let yr = x.clone();
229
230        Swap::swap(&mut x, &mut y);
231        assert_eq!(x, xr);
232        assert_eq!(y, yr);
233    }
234
235    #[test]
236    fn slice() {
237        let mut x = [1f32, -2f32, 3f32, 4f32];
238        let mut y = [2f32, -3f32, 4f32, 1f32];
239        let xr = [2f32, -3f32, 4f32, 1f32];
240        let yr = [1f32, -2f32, 3f32, 4f32];
241
242        Swap::swap(&mut x[..], &mut y[..]);
243        assert_eq!(x, xr);
244        assert_eq!(y, yr);
245    }
246
247    #[test]
248    fn complex() {
249        let mut x = vec![Complex::new(2f32, -3f32)];
250        let mut y = vec![Complex::new(-1f32, 4f32)];
251        let xr = y.clone();
252        let yr = x.clone();
253
254        Swap::swap(&mut x, &mut y);
255        assert_eq!(x, xr);
256        assert_eq!(y, yr);
257    }
258}
259
260/// Computes `x^T * y`.
261pub trait Dot: Sized {
262    fn dot<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &V, y: &W) -> Self;
263}
264
265macro_rules! real_dot_impl(($($t: ident), +) => (
266    $(
267        impl Dot for $t {
268            fn dot<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &V, y: &W) -> $t {
269                unsafe {
270                    let n = cmp::min(x.len(), y.len());
271
272                    prefix!($t, dot)(n,
273                        x.as_ptr().as_c_ptr(), x.inc(),
274                        y.as_ptr().as_c_ptr(), y.inc())
275                }
276            }
277        }
278    )+
279));
280
281macro_rules! complex_dot_impl(($($t: ident), +) => (
282    $(
283        impl Dot for $t {
284            fn dot<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &V, y: &W) -> $t {
285                let result: $t = Default::zero();
286
287                unsafe {
288                    let n = cmp::min(x.len(), y.len());
289
290                    prefix!($t, dotu_sub)(n,
291                        x.as_ptr().as_c_ptr(), x.inc(),
292                        y.as_ptr().as_c_ptr(), y.inc(),
293                        (&result).as_mut());
294                }
295
296                result
297            }
298        }
299    )+
300));
301
302real_dot_impl!(f32, f64);
303complex_dot_impl!(Complex32, Complex64);
304
305#[cfg(test)]
306mod dot_tests {
307    use crate::vector::ops::Dot;
308    use num_complex::Complex;
309
310    #[test]
311    fn real() {
312        let x = vec![1f32, -2f32, 3f32, 4f32];
313        let y = vec![1f32, 1f32, 1f32, 1f32];
314
315        let xr: f32 = Dot::dot(&x, &y);
316        assert_eq!(xr, 6f32);
317    }
318
319    #[test]
320    fn slice() {
321        let x = [1f32, -2f32, 3f32, 4f32];
322        let y = [1f32, 1f32, 1f32, 1f32];
323
324        let xr: f32 = Dot::dot(&x[..], &y[..]);
325        assert_eq!(xr, 6f32);
326    }
327
328    #[test]
329    fn complex() {
330        let x = vec![Complex::new(1f32, 1f32), Complex::new(1f32, 3f32)];
331        let y = vec![Complex::new(1f32, 1f32), Complex::new(1f32, 1f32)];
332
333        let xr: Complex<f32> = Dot::dot(&x, &y);
334        assert_eq!(xr, Complex::new(-2f32, 6f32));
335    }
336}
337
338/// Computes `x^H * y`.
339pub trait Dotc: Sized + Dot {
340    fn dotc<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &V, y: &W) -> Self {
341        Dot::dot(x, y)
342    }
343}
344
345macro_rules! dotc_impl(($($t: ident), +) => (
346    $(
347        impl Dotc for $t {
348            fn dotc<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &V, y: &W) -> $t {
349                let result: $t = Default::zero();
350
351                unsafe {
352                    let n = cmp::min(x.len(), y.len());
353
354                    prefix!($t, dotc_sub)(n,
355                        x.as_ptr().as_c_ptr(), x.inc(),
356                        y.as_ptr().as_c_ptr(), y.inc(),
357                        (&result).as_mut());
358                }
359
360                result
361            }
362        }
363    )+
364));
365
366impl Dotc for f32 {}
367impl Dotc for f64 {}
368dotc_impl!(Complex32, Complex64);
369
370#[cfg(test)]
371mod dotc_tests {
372    use crate::vector::ops::Dotc;
373    use num_complex::Complex;
374
375    #[test]
376    fn complex_conj() {
377        let x = vec![Complex::new(1f32, -1f32), Complex::new(1f32, -3f32)];
378        let y = vec![Complex::new(1f32, 2f32), Complex::new(1f32, 3f32)];
379
380        let xr: Complex<f32> = Dotc::dotc(&x, &y);
381        assert_eq!(xr, Complex::new(-9f32, 9f32));
382    }
383}
384
385/// Computes the sum of the absolute values of elements in a vector.
386///
387/// Complex vectors use `||Re(x)||_1 + ||Im(x)||_1`
388pub trait Asum: Sized {
389    fn asum<V: ?Sized + Vector<Self>>(x: &V) -> Self;
390}
391
392/// Computes the L2 norm (Euclidian length) of a vector.
393pub trait Nrm2: Sized {
394    fn nrm2<V: ?Sized + Vector<Self>>(x: &V) -> Self;
395}
396
397macro_rules! real_norm_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
398    $(
399        impl $trait_name for $t {
400            fn $fn_name<V: ?Sized + Vector<Self>>(x: &V) -> $t {
401                unsafe {
402                    prefix!($t, $fn_name)(x.len(),
403                        x.as_ptr().as_c_ptr(), x.inc())
404                }
405            }
406        }
407    )+
408));
409
410macro_rules! complex_norm_impl(
411    ($trait_name: ident, $fn_name: ident, $t: ty, $norm_fn: expr) => (
412        impl $trait_name for $t {
413            fn $fn_name<V: ?Sized + Vector<Self>>(x: &V) -> $t {
414                let re = unsafe {
415                    $norm_fn(x.len(),
416                        x.as_ptr().as_c_ptr(), x.inc())
417                };
418
419                Complex { im: 0.0, re: re }
420            }
421        }
422    );
423);
424
425real_norm_impl!(Asum, asum, f32, f64);
426real_norm_impl!(Nrm2, nrm2, f32, f64);
427complex_norm_impl!(Asum, asum, Complex32, cblas_s::casum);
428complex_norm_impl!(Asum, asum, Complex64, cblas_d::zasum);
429complex_norm_impl!(Nrm2, nrm2, Complex32, cblas_s::cnrm2);
430complex_norm_impl!(Nrm2, nrm2, Complex64, cblas_d::znrm2);
431
432#[cfg(test)]
433mod asum_tests {
434    use crate::vector::ops::Asum;
435    use num_complex::Complex;
436
437    #[test]
438    fn real() {
439        let x = vec![1f32, -2f32, 3f32, 4f32];
440
441        let r: f32 = Asum::asum(&x);
442        assert_eq!(r, 10f32);
443    }
444
445    #[test]
446    fn slice() {
447        let x = [1f32, -2f32, 3f32, 4f32];
448
449        let r: f32 = Asum::asum(&x[..]);
450        assert_eq!(r, 10f32);
451    }
452
453    #[test]
454    fn complex() {
455        let x = vec![Complex::new(3f32, 4f32)];
456
457        let r: Complex<f32> = Asum::asum(&x);
458        assert_eq!(r, Complex { im: 0.0, re: 7f32 });
459    }
460}
461
462#[cfg(test)]
463mod nrm2_tests {
464    use crate::vector::ops::Nrm2;
465    use num_complex::Complex;
466
467    #[test]
468    fn real() {
469        let x = vec![3f32, -4f32];
470
471        let xr: f32 = Nrm2::nrm2(&x);
472        assert_eq!(xr, 5f32);
473    }
474
475    #[test]
476    fn slice() {
477        let x = [3f32, -4f32];
478
479        let xr: f32 = Nrm2::nrm2(&x[..]);
480        assert_eq!(xr, 5f32);
481    }
482
483    #[test]
484    fn complex() {
485        let x = vec![Complex::new(3f32, 4f32)];
486
487        let xr: Complex<f32> = Nrm2::nrm2(&x);
488        assert_eq!(xr, Complex { im: 0.0, re: 5f32 });
489    }
490}
491
492/// Finds the index of the element with maximum absolute value in a vector.
493///
494/// Complex vectors maximize the value `|Re(x_k)| + |Im(x_k)|`.
495///
496/// The first index with a maximum is returned.
497pub trait Iamax: Sized {
498    fn iamax<V: ?Sized + Vector<Self>>(x: &V) -> usize;
499}
500
501macro_rules! iamax_impl(
502    ($t: ty, $iamax: expr) => (
503        impl Iamax for $t {
504            fn iamax<V: ?Sized + Vector<Self>>(x: &V) -> usize {
505                unsafe {
506                    $iamax(x.len(),
507                        x.as_ptr().as_c_ptr(), x.inc()) as usize
508                }
509            }
510        }
511    );
512);
513
514iamax_impl!(f32, cblas_i::samax);
515iamax_impl!(f64, cblas_i::damax);
516iamax_impl!(Complex32, cblas_i::camax);
517iamax_impl!(Complex64, cblas_i::zamax);
518
519#[cfg(test)]
520mod iamax_tests {
521    use crate::vector::ops::Iamax;
522    use num_complex::Complex;
523
524    #[test]
525    fn real() {
526        let x = vec![1f32, -2f32, 3f32, 4f32];
527
528        let xr = Iamax::iamax(&x);
529        assert_eq!(xr, 3usize);
530    }
531
532    #[test]
533    fn slice() {
534        let x = [1f32, -2f32, 3f32, 4f32];
535
536        let xr = Iamax::iamax(&x[..]);
537        assert_eq!(xr, 3usize);
538    }
539
540    #[test]
541    fn complex() {
542        let x = vec![Complex::new(3f32, 4f32), Complex::new(3f32, 5f32)];
543
544        let xr = Iamax::iamax(&x);
545        assert_eq!(xr, 1usize);
546    }
547}
548
549/// Applies a Givens rotation matrix to a pair of vectors, where `cos` is
550/// the value of the cosine of the angle in the Givens matrix, and `sin` is
551/// the sine.
552pub trait Rot: Sized {
553    fn rot<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
554        x: &mut V,
555        y: &mut W,
556        cos: &Self,
557        sin: &Self,
558    );
559}
560
561macro_rules! rot_impl(($($t: ident), +) => (
562    $(
563        impl Rot for $t {
564            fn rot<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(x: &mut V, y: &mut W, cos: &$t, sin: &$t) {
565                unsafe {
566                    prefix!($t, rot)(cmp::min(x.len(), y.len()),
567                        x.as_mut_ptr().as_c_ptr(), x.inc(),
568                        y.as_mut_ptr().as_c_ptr(), y.inc(),
569                        cos.as_const(), sin.as_const());
570                }
571            }
572        }
573    )+
574));
575
576rot_impl!(f32, f64);
577
578#[cfg(test)]
579mod rot_tests {
580    use crate::vector::ops::{Rot, Scal};
581
582    #[test]
583    fn real() {
584        let mut x = vec![1f32, -2f32, 3f32, 4f32];
585        let mut y = vec![3f32, 7f32, -2f32, 2f32];
586        let cos = 0f32;
587        let sin = 1f32;
588
589        let xr = y.clone();
590        let mut yr = x.clone();
591        Scal::scal(&-1f32, &mut yr);
592
593        Rot::rot(&mut x, &mut y, &cos, &sin);
594        assert_eq!(x, xr);
595        assert_eq!(y, yr);
596    }
597}
rust_blas/vector/ops.rs

rust_blas/vector/
ops.rs
