1use 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
15pub 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
92pub 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
106pub 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 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 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
184pub 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
196pub 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
256pub 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
268pub 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
316pub 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
329pub 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
362pub 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
397pub 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
411pub 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
446pub 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
459pub 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
491pub 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
505pub 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
540pub 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
553pub 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
585pub 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
615pub 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
645pub 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
658pub 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);