clarabel/algebra/dense/blas/
syr2k.rs

1#![allow(non_snake_case)]
2
3use crate::algebra::*;
4
5impl<S, T> MultiplySYR2K<T> for DenseStorageMatrix<S, T>
6where
7    T: FloatT,
8    S: AsMut<[T]> + AsRef<[T]>,
9{
10    // implements self = C = α(A*B' + B*A') + βC
11    fn syr2k<S1, S2>(
12        &mut self,
13        A: &DenseStorageMatrix<S1, T>,
14        B: &DenseStorageMatrix<S2, T>,
15        α: T,
16        β: T,
17    ) where
18        S1: AsRef<[T]>,
19        S2: AsRef<[T]>,
20    {
21        assert!(self.nrows() == A.nrows());
22        assert!(self.nrows() == B.nrows());
23        assert!(self.ncols() == B.nrows());
24        assert!(A.ncols() == B.ncols());
25
26        if self.nrows() == 0 {
27            return;
28        }
29
30        let n = self.nrows().try_into().unwrap();
31        let k = A.ncols().try_into().unwrap();
32        let a = A.data();
33        let lda = n;
34        let b = B.data();
35        let ldb = n;
36        let c = self.data_mut();
37        let ldc = n;
38
39        #[rustfmt::skip]
40        T::xsyr2k(
41            MatrixTriangle::Triu.as_blas_char(),
42            MatrixShape::N.as_blas_char(),
43            n,k,α,a,lda,b,ldb,β,c,ldc,
44        );
45    }
46}
47
48macro_rules! generate_test_syr2k {
49    ($fxx:ty, $test_name:ident) => {
50        #[test]
51        fn $test_name() {
52
53            #[rustfmt::skip]
54            let A = Matrix::<$fxx>::from(&[
55                [ 1., -5.], 
56                [-4.,  3.], 
57                [ 2.,  6.],
58            ]);
59
60            let B = Matrix::<$fxx>::from(&[
61                [ 4.,  5.], 
62                [ 2., -2.], 
63                [-3., -2.],
64            ]);
65
66            let mut C = Matrix::<$fxx>::identity(3);
67
68            //NB: modifies upper triangle only
69            C.syr2k(&A, &B, 2., 1.);
70
71            let Ctest = Matrix::<$fxx>::from(&[
72                [-83.,  22.,  90.], 
73                [  0., -55.,  -4.], 
74                [  0.,   0., -71.],
75            ]);
76
77            assert_eq!(C, Ctest);
78        }
79    };
80}
81
82generate_test_syr2k!(f32, test_syr2k_f32);
83generate_test_syr2k!(f64, test_syr2k_f64);