1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

use crate::matrix::MatrixError;
use crate::{matrix::Matrix, number::c64};
use blas::dgemm;
use blas::zgemm;

impl Matrix {
  pub fn gemm(
    self,
    lhs: &Matrix,
    rhs: &Matrix,
    alpha: f64,
    beta: f64,
  ) -> Result<Matrix, MatrixError> {
    if self.rows != lhs.rows || self.cols != rhs.cols || lhs.cols != rhs.rows {
      return Err(MatrixError::DimensionMismatch);
    }

    let m = lhs.rows as i32;
    let k = lhs.cols as i32;
    let n = rhs.cols as i32;

    let mut slf = self;

    unsafe {
      dgemm(
        'N' as u8,
        'N' as u8,
        m,
        n,
        k,
        alpha,
        lhs.elems.as_slice(),
        m,
        rhs.elems.as_slice(),
        k,
        beta,
        &mut slf.elems,
        m,
      );
    }

    Ok(slf)
  }
}

impl Matrix<c64> {
  pub fn gemm(
    self,
    lhs: &Matrix<c64>,
    rhs: &Matrix<c64>,
    alpha: c64,
    beta: c64,
  ) -> Result<Matrix<c64>, MatrixError> {
    if self.rows != lhs.rows || self.cols != rhs.cols || lhs.cols != rhs.rows {
      return Err(MatrixError::DimensionMismatch);
    }

    let m = lhs.rows as i32;
    let k = lhs.cols as i32;
    let n = rhs.cols as i32;

    let mut slf = self;

    unsafe {
      zgemm(
        'N' as u8,
        'N' as u8,
        m,
        n,
        k,
        alpha,
        rhs.elems.as_slice(),
        m,
        lhs.elems.as_slice(),
        k,
        beta,
        &mut slf.elems,
        m,
      );
    }

    Ok(slf)
  }
}