rust_blas/matrix_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 matrix-vector functions.
6
7use crate::attribute::{Diagonal, Symmetry, Transpose};
8use crate::matrix::{BandMatrix, Matrix};
9use crate::matrix_vector::ll::*;
10use crate::pointer::CPtr;
11use crate::scalar::Scalar;
12use crate::vector::Vector;
13use num_complex::{Complex, Complex32, Complex64};
14
15/// General multiply with vector
16///
17/// A ← αA<sup>OP</sup>x + βy
18pub trait Gemv: Sized {
19    fn gemv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
20        trans: Transpose,
21        alpha: &Self,
22        a: &dyn Matrix<Self>,
23        x: &V,
24        beta: &Self,
25        y: &mut W,
26    );
27}
28
29macro_rules! gemv_impl(($($t: ident), +) => (
30    $(
31        impl Gemv for $t {
32            fn gemv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(trans: Transpose, alpha: &$t, a: &dyn Matrix<$t>, x: &V, beta: &$t, y: &mut W){
33                unsafe {
34                    prefix!($t, gemv)(a.order(), trans,
35                        a.rows(), a.cols(),
36                        alpha.as_const(),
37                        a.as_ptr().as_c_ptr(), a.lead_dim(),
38                        x.as_ptr().as_c_ptr(), x.inc(),
39                        beta.as_const(),
40                        y.as_mut_ptr().as_c_ptr(), y.inc());
41                }
42            }
43        }
44    )+
45));
46
47gemv_impl!(f32, f64, Complex32, Complex64);
48
49#[cfg(test)]
50mod gemv_tests {
51    use crate::attribute::Transpose;
52    use crate::matrix::tests::M;
53    use crate::matrix_vector::ops::Gemv;
54
55    #[test]
56    fn real() {
57        let a = M(2, 2, vec![1.0, -2.0, 2.0, -4.0]);
58        let x = vec![2.0, 1.0];
59        let mut y = vec![1.0, 2.0];
60        let t = Transpose::NoTrans;
61
62        Gemv::gemv(t, &1f32, &a, &x, &0f32, &mut y);
63
64        assert_eq!(y, vec![0.0, 0.0]);
65    }
66
67    #[test]
68    fn non_square() {
69        let a = M(2, 3, vec![1.0, -3.0, 1.0, 2.0, -6.0, 2.0]);
70        let x = vec![2.0, 1.0, 1.0];
71        let mut y = vec![1.0, 2.0];
72        let t = Transpose::NoTrans;
73
74        Gemv::gemv(t, &1f32, &a, &x, &0f32, &mut y);
75        assert_eq!(y, vec![0.0, 0.0]);
76    }
77
78    #[test]
79    fn transpose() {
80        let a = M(3, 2, vec![1.0, 2.0, -3.0, -6.0, 1.0, 2.0]);
81
82        let x = vec![2.0, 1.0, 1.0];
83        let mut y = vec![1.0, 2.0];
84        let t = Transpose::Trans;
85
86        Gemv::gemv(t, &1f32, &a, &x, &0f32, &mut y);
87
88        assert_eq!(y, vec![0.0, 0.0]);
89    }
90}
91
92/// Symmetric multiply with vector
93///
94/// A ← αAx + βy
95pub trait Symv: Sized {
96    fn symv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
97        symmetry: Symmetry,
98        alpha: &Self,
99        a: &dyn Matrix<Self>,
100        x: &V,
101        beta: &Self,
102        y: &mut W,
103    );
104}
105
106/// Hermitian multiply with vector
107///
108/// A ← αAx + βy
109pub trait Hemv: Sized {
110    fn hemv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
111        symmetry: Symmetry,
112        alpha: &Self,
113        a: &dyn Matrix<Self>,
114        x: &V,
115        beta: &Self,
116        y: &mut W,
117    );
118}
119
120macro_rules! symv_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
121    $(
122        impl $trait_name for $t {
123            fn $fn_name<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, a: &dyn Matrix<$t>, x: &V, beta: &$t, y: &mut W){
124                unsafe {
125                    prefix!($t, $fn_name)(a.order(), symmetry,
126                        a.rows(),
127                        alpha.as_const(),
128                        a.as_ptr().as_c_ptr(), a.lead_dim(),
129                        x.as_ptr().as_c_ptr(), x.inc(),
130                        beta.as_const(),
131                        y.as_mut_ptr().as_c_ptr(), y.inc());
132                }
133            }
134        }
135    )+
136));
137
138symv_impl!(Symv, symv, f32, f64, Complex32, Complex64);
139symv_impl!(Hemv, hemv, Complex32, Complex64);
140
141#[cfg(test)]
142mod symv_tests {
143    use crate::attribute::{Symmetry, Transpose};
144    use crate::matrix::tests::M;
145    use crate::matrix_vector::ops::{Gemv, Symv};
146
147    #[test]
148    fn real() {
149        let x = vec![2.0, 1.0];
150        let gemv = {
151            let a = M(2, 2, vec![1.0, -2.0, -2.0, 1.0]);
152            let mut y = vec![1.0, 2.0];
153            let t = Transpose::NoTrans;
154
155            Gemv::gemv(t, &1.0, &a, &x, &0.0, &mut y);
156            y
157        };
158
159        let symv_upper = {
160            // symv shouldn't look at some elements
161            let a = M(2, 2, vec![1.0, -2.0, 0.0, 1.0]);
162            let mut y = vec![1.0, 2.0];
163            let s = Symmetry::Upper;
164
165            Symv::symv(s, &1.0, &a, &x, &0.0, &mut y);
166            y
167        };
168
169        let symv_lower = {
170            // symv shouldn't look at some elements
171            let a = M(2, 2, vec![1.0, 0.0, -2.0, 1.0]);
172            let mut y = vec![1.0, 2.0];
173            let s = Symmetry::Lower;
174
175            Symv::symv(s, &1.0, &a, &x, &0.0, &mut y);
176            y
177        };
178
179        assert_eq!(gemv, symv_upper);
180        assert_eq!(gemv, symv_lower);
181    }
182}
183
184/// General rank-1 update
185///
186/// A ← A + αxy<sup>T</sup>
187pub trait Ger: Sized {
188    fn ger<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
189        alpha: &Self,
190        x: &V,
191        y: &W,
192        a: &mut dyn Matrix<Self>,
193    );
194}
195
196/// General rank-1 update (using hermitian conjugate)
197///
198/// A ← A + αxy<sup>H</sup>
199pub trait Gerc: Ger {
200    fn gerc<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
201        alpha: &Self,
202        x: &V,
203        y: &W,
204        a: &mut dyn Matrix<Self>,
205    ) {
206        Ger::ger(alpha, x, y, a);
207    }
208}
209
210macro_rules! ger_impl(
211    ($trait_name: ident, $fn_name: ident, $t: ty, $ger_fn: expr) => (
212        impl $trait_name for $t {
213            fn $fn_name<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(alpha: &$t, x: &V, y: &W, a: &mut dyn Matrix<$t>) {
214                unsafe {
215                    $ger_fn(a.order(),
216                        a.rows(), a.cols(),
217                        alpha.as_const(),
218                        x.as_ptr().as_c_ptr(), x.inc(),
219                        y.as_ptr().as_c_ptr(), y.inc(),
220                        a.as_mut_ptr().as_c_ptr(), a.lead_dim());
221                }
222            }
223        }
224    );
225);
226
227ger_impl!(Ger, ger, f32, cblas_s::ger);
228ger_impl!(Ger, ger, f64, cblas_d::ger);
229ger_impl!(Ger, ger, Complex32, cblas_c::geru);
230ger_impl!(Ger, ger, Complex64, cblas_z::geru);
231
232impl Gerc for f32 {}
233impl Gerc for f64 {}
234ger_impl!(Gerc, gerc, Complex32, cblas_c::gerc);
235ger_impl!(Gerc, gerc, Complex64, cblas_z::gerc);
236
237#[cfg(test)]
238mod ger_tests {
239    use crate::matrix::tests::M;
240    use crate::matrix_vector::ops::Ger;
241    use std::iter::repeat;
242
243    #[test]
244    fn real() {
245        let mut a = M(3, 3, repeat(0.0).take(9).collect());
246        let x = vec![2.0, 1.0, 4.0];
247        let y = vec![3.0, 6.0, -1.0];
248
249        Ger::ger(&1f32, &x, &y, &mut a);
250
251        let result = vec![6.0, 12.0, -2.0, 3.0, 6.0, -1.0, 12.0, 24.0, -4.0];
252        assert_eq!(a.2, result);
253    }
254}
255
256/// Symmetric rank-1 update
257///
258/// A ← A + αxx<sup>T</sup>
259pub trait Syr: Sized {
260    fn syr<V: ?Sized + Vector<Self>>(
261        symmetry: Symmetry,
262        alpha: &Self,
263        x: &V,
264        a: &mut dyn Matrix<Self>,
265    );
266}
267
268/// Hermitian rank-1 update
269///
270/// A ← A + αxx<sup>H</sup>
271pub trait Her: Sized {
272    fn her<V: ?Sized + Vector<Complex<Self>>>(
273        symmetry: Symmetry,
274        alpha: &Self,
275        x: &V,
276        a: &mut dyn Matrix<Complex<Self>>,
277    );
278}
279
280macro_rules! her_impl(($($t: ident), +) => (
281    $(
282        impl Her for $t {
283            fn her<V: ?Sized + Vector<Complex<Self>>>(symmetry: Symmetry, alpha: &$t, x: &V, a: &mut dyn Matrix<Complex<$t>>) {
284                unsafe {
285                    prefix!(Complex<$t>, her)(a.order(), symmetry,
286                        a.rows(),
287                        *alpha,
288                        x.as_ptr().as_c_ptr(), x.inc(),
289                        a.as_mut_ptr().as_c_ptr(), a.lead_dim());
290                }
291            }
292        }
293    )+
294));
295
296her_impl!(f32, f64);
297
298macro_rules! syr_impl(($($t: ident), +) => (
299    $(
300        impl Syr for $t {
301            fn syr<V: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, x: &V, a: &mut dyn Matrix<$t>) {
302                unsafe {
303                    prefix!($t, syr)(a.order(), symmetry,
304                        a.rows(),
305                        *alpha,
306                        x.as_ptr().as_c_ptr(), x.inc(),
307                        a.as_mut_ptr().as_c_ptr(), a.lead_dim());
308                }
309            }
310        }
311    )+
312));
313
314syr_impl!(f32, f64);
315
316/// Symmetric rank-2 update
317///
318/// A ← A + αxy<sup>T</sup> + αyx<sup>T</sup>
319pub trait Syr2: Sized {
320    fn syr2<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
321        symmetry: Symmetry,
322        alpha: &Self,
323        x: &V,
324        y: &W,
325        a: &mut dyn Matrix<Self>,
326    );
327}
328
329/// Hermitian rank-2 update
330///
331/// A ← A + αxy<sup>H</sup> + y(αx)<sup>H</sup>
332pub trait Her2: Sized {
333    fn her2<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
334        symmetry: Symmetry,
335        alpha: &Self,
336        x: &V,
337        y: &W,
338        a: &mut dyn Matrix<Self>,
339    );
340}
341
342macro_rules! syr2_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
343    $(
344        impl $trait_name for $t {
345            fn $fn_name<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, x: &V, y: &W, a: &mut dyn Matrix<$t>) {
346                unsafe {
347                    prefix!($t, $fn_name)(a.order(), symmetry,
348                        a.rows(),
349                        alpha.as_const(),
350                        x.as_ptr().as_c_ptr(), x.inc(),
351                        y.as_ptr().as_c_ptr(), y.inc(),
352                        a.as_mut_ptr().as_c_ptr(), a.lead_dim());
353                }
354            }
355        }
356    )+
357));
358
359syr2_impl!(Syr2, syr2, f32, f64);
360syr2_impl!(Her2, her2, Complex32, Complex64);
361
362/// General band matrix multiply with vector.
363///
364/// A ← αA<sup>OP</sup>x + βy
365pub trait Gbmv: Sized {
366    fn gbmv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
367        trans: Transpose,
368        alpha: &Self,
369        a: &dyn BandMatrix<Self>,
370        x: &V,
371        beta: &Self,
372        y: &mut W,
373    );
374}
375
376macro_rules! gbmv_impl(($($t: ident), +) => (
377    $(
378        impl Gbmv for $t {
379            fn gbmv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(trans: Transpose, alpha: &$t, a: &dyn BandMatrix<$t>, x: &V, beta: &$t, y: &mut W){
380                unsafe {
381                    prefix!($t, gbmv)(a.order(), trans,
382                        a.rows(), a.cols(),
383                        a.sub_diagonals(), a.sup_diagonals(),
384                        alpha.as_const(),
385                        a.as_ptr().as_c_ptr(), a.lead_dim(),
386                        x.as_ptr().as_c_ptr(), x.inc(),
387                        beta.as_const(),
388                        y.as_mut_ptr().as_c_ptr(), y.inc());
389                }
390            }
391        }
392    )+
393));
394
395gbmv_impl!(f32, f64, Complex32, Complex64);
396
397/// Symmetric band matrix multiply with vector
398///
399/// A ← αAx + βy
400pub trait Sbmv: Sized {
401    fn sbmv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
402        symmetry: Symmetry,
403        alpha: &Self,
404        a: &dyn BandMatrix<Self>,
405        x: &V,
406        beta: &Self,
407        y: &mut W,
408    );
409}
410
411/// Hermitian band matrix multiply with vector
412///
413/// A ← αAx + βy
414pub trait Hbmv: Sized {
415    fn hbmv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
416        symmetry: Symmetry,
417        alpha: &Self,
418        a: &dyn BandMatrix<Self>,
419        x: &V,
420        beta: &Self,
421        y: &mut W,
422    );
423}
424
425macro_rules! sbmv_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
426    $(
427        impl $trait_name for $t {
428            fn $fn_name<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, a: &dyn BandMatrix<$t>, x: &V, beta: &$t, y: &mut W) {
429                unsafe {
430                    prefix!($t, $fn_name)(a.order(), symmetry,
431                        a.rows(), a.sub_diagonals(),
432                        alpha.as_const(),
433                        a.as_ptr().as_c_ptr(), a.lead_dim(),
434                        x.as_ptr().as_c_ptr(), x.inc(),
435                        beta.as_const(),
436                        y.as_mut_ptr().as_c_ptr(), y.inc());
437                }
438            }
439        }
440    )+
441));
442
443sbmv_impl!(Sbmv, sbmv, f32, f64);
444sbmv_impl!(Hbmv, hbmv, Complex32, Complex64);
445
446/// Triangular band matrix multiply with vector
447///
448/// A ← A<sup>OP</sup>x
449pub trait Tbmv: Sized {
450    fn tbmv<V: ?Sized + Vector<Self>>(
451        symmetry: Symmetry,
452        trans: Transpose,
453        diagonal: Diagonal,
454        a: &dyn BandMatrix<Self>,
455        x: &mut V,
456    );
457}
458
459/// Solve triangular band matrix system
460///
461/// A ← A<sup>-1 OP</sup>x
462pub trait Tbsv: Sized {
463    fn tbsv<V: ?Sized + Vector<Self>>(
464        symmetry: Symmetry,
465        trans: Transpose,
466        diagonal: Diagonal,
467        a: &dyn BandMatrix<Self>,
468        x: &mut V,
469    );
470}
471
472macro_rules! tbmv_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
473    $(
474        impl $trait_name for $t {
475            fn $fn_name<V: ?Sized + Vector<Self>>(symmetry: Symmetry, trans: Transpose, diagonal: Diagonal, a: &dyn BandMatrix<$t>, x: &mut V) {
476                unsafe {
477                    prefix!($t, $fn_name)(a.order(), symmetry,
478                        trans, diagonal,
479                        a.rows(), a.sub_diagonals(),
480                        a.as_ptr().as_c_ptr(),
481                        x.as_mut_ptr().as_c_ptr(), x.inc());
482                }
483            }
484        }
485    )+
486));
487
488tbmv_impl!(Tbmv, tbmv, f32, f64, Complex32, Complex64);
489tbmv_impl!(Tbsv, tbsv, f32, f64, Complex32, Complex64);
490
491/// Symmetric packed matrix multiply with vector
492///
493/// A ← αAx + βy
494pub trait Spmv: Sized {
495    fn spmv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
496        symmetry: Symmetry,
497        alpha: &Self,
498        a: &dyn Matrix<Self>,
499        x: &V,
500        beta: &Self,
501        y: &mut W,
502    );
503}
504
505/// Hermitian packed matrix multiply with vector
506///
507/// A ← αAx + βy
508pub trait Hpmv: Sized {
509    fn hpmv<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
510        symmetry: Symmetry,
511        alpha: &Self,
512        a: &dyn Matrix<Self>,
513        x: &V,
514        beta: &Self,
515        y: &mut W,
516    );
517}
518
519macro_rules! spmv_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
520    $(
521        impl $trait_name for $t {
522            fn $fn_name<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, a: &dyn Matrix<$t>, x: &V, beta: &$t, y: &mut W) {
523                unsafe {
524                    prefix!($t, $fn_name)(a.order(), symmetry,
525                        a.rows(),
526                        alpha.as_const(),
527                        a.as_ptr().as_c_ptr(),
528                        x.as_ptr().as_c_ptr(), x.inc(),
529                        beta.as_const(),
530                        y.as_mut_ptr().as_c_ptr(), y.inc());
531                }
532            }
533        }
534    )+
535));
536
537spmv_impl!(Spmv, spmv, f32, f64);
538spmv_impl!(Hpmv, hpmv, Complex32, Complex64);
539
540/// Triangular packed matrix multiply with vector
541///
542/// A ← A<sup>OP</sup>x
543pub trait Tpmv: Sized {
544    fn tpmv<V: ?Sized + Vector<Self>>(
545        symmetry: Symmetry,
546        trans: Transpose,
547        diagonal: Diagonal,
548        a: &dyn Matrix<Self>,
549        x: &mut V,
550    );
551}
552
553/// Solve triangular packed matrix system
554///
555/// A ← A<sup>-1 OP</sup>x
556pub trait Tpsv: Sized {
557    fn tpsv<V: ?Sized + Vector<Self>>(
558        symmetry: Symmetry,
559        trans: Transpose,
560        diagonal: Diagonal,
561        a: &dyn Matrix<Self>,
562        x: &mut V,
563    );
564}
565
566macro_rules! tpmv_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
567    $(
568        impl $trait_name for $t {
569            fn $fn_name<V: ?Sized + Vector<Self>>(symmetry: Symmetry, trans: Transpose, diagonal: Diagonal, a: &dyn Matrix<$t>, x: &mut V) {
570                unsafe {
571                    prefix!($t, $fn_name)(a.order(), symmetry,
572                        trans, diagonal,
573                        a.rows(),
574                        a.as_ptr().as_c_ptr(),
575                        x.as_mut_ptr().as_c_ptr(), x.inc());
576                }
577            }
578        }
579    )+
580));
581
582tpmv_impl!(Tpmv, tpmv, f32, f64, Complex32, Complex64);
583tpmv_impl!(Tpsv, tpsv, f32, f64, Complex32, Complex64);
584
585/// Hermitian packed matrix rank-1 update
586///
587/// A ← A + αxx<sup>H</sup>
588pub trait Hpr: Sized {
589    fn hpr<V: ?Sized + Vector<Complex<Self>>>(
590        symmetry: Symmetry,
591        alpha: &Self,
592        x: &V,
593        a: &mut dyn Matrix<Complex<Self>>,
594    );
595}
596
597macro_rules! hpr_impl(($($t: ident), +) => (
598    $(
599        impl Hpr for $t {
600            fn hpr<V: ?Sized + Vector<Complex<Self>>>(symmetry: Symmetry, alpha: &$t, x: &V, a: &mut dyn Matrix<Complex<$t>>) {
601                unsafe {
602                    prefix!(Complex<$t>, hpr)(a.order(), symmetry,
603                        a.rows(),
604                        *alpha,
605                        x.as_ptr().as_c_ptr(), x.inc(),
606                        a.as_mut_ptr().as_c_ptr());
607                }
608            }
609        }
610    )+
611));
612
613hpr_impl!(f32, f64);
614
615/// Symmetric packed matrix rank-1 update
616///
617/// A ← A + αxx<sup>T</sup>
618pub trait Spr: Sized {
619    fn spr<V: ?Sized + Vector<Self>>(
620        symmetry: Symmetry,
621        alpha: &Self,
622        x: &V,
623        a: &mut dyn Matrix<Self>,
624    );
625}
626
627macro_rules! spr_impl(($($t: ident), +) => (
628    $(
629        impl Spr for $t {
630            fn spr<V: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, x: &V, a: &mut dyn Matrix<$t>) {
631                unsafe {
632                    prefix!($t, spr)(a.order(), symmetry,
633                        a.rows(),
634                        *alpha,
635                        x.as_ptr().as_c_ptr(), x.inc(),
636                        a.as_mut_ptr().as_c_ptr());
637                }
638            }
639        }
640    )+
641));
642
643spr_impl!(f32, f64);
644
645/// Symmetric packed matrix rank-2 update
646///
647/// A ← A + αxy<sup>T</sup> + αyx<sup>T</sup>
648pub trait Spr2: Sized {
649    fn spr2<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
650        symmetry: Symmetry,
651        alpha: &Self,
652        x: &V,
653        y: &W,
654        a: &mut dyn Matrix<Self>,
655    );
656}
657
658/// Hermitian packed matrix rank-2 update
659///
660/// A ← A + αxy<sup>H</sup> + y(αx)<sup>H</sup>
661pub trait Hpr2: Sized {
662    fn hpr2<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(
663        symmetry: Symmetry,
664        alpha: &Self,
665        x: &V,
666        y: &W,
667        a: &mut dyn Matrix<Self>,
668    );
669}
670
671macro_rules! spr2_impl(($trait_name: ident, $fn_name: ident, $($t: ident), +) => (
672    $(
673        impl $trait_name for $t {
674            fn $fn_name<V: ?Sized + Vector<Self>, W: ?Sized + Vector<Self>>(symmetry: Symmetry, alpha: &$t, x: &V, y: &W, a: &mut dyn Matrix<$t>) {
675                unsafe {
676                    prefix!($t, $fn_name)(a.order(), symmetry,
677                        a.rows(),
678                        alpha.as_const(),
679                        x.as_ptr().as_c_ptr(), x.inc(),
680                        y.as_ptr().as_c_ptr(), y.inc(),
681                        a.as_mut_ptr().as_c_ptr());
682                }
683            }
684        }
685    )+
686));
687
688spr2_impl!(Spr2, spr2, f32, f64);
689spr2_impl!(Hpr2, hpr2, Complex32, Complex64);