2 files changed, 60 insertions, 0 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 2843ab5..4d48594 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1398,6 +1398,65 @@ ROWOP(int64_t)
 ROWOP_MOD(int32_t,mod)
 ROWOP_MOD(int64_t,mod_l)
 
+/////////////////////////////// inplace GEMM ////////////////////////////////
+
+#define GEMM(T) int gemm_##T(VECG(T,c),VECG(int,p),MATG(T,a),MATG(T,b),MATG(T,r)) {  \
+    T a = cp[0], b = cp[1];                                    \
+    int r1a = pp[0],              c1a = pp[2],  c2a = pp[3] ;  \
+    int r1b = pp[4],              c1b = pp[6]               ;  \
+    int r1r = pp[8], r2r = pp[9], c1r = pp[10], c2r = pp[11];  \
+    int dra = r1a - r1r;                                       \
+    int dcb = c1b-c1r;                                         \
+    int nk  = c2a-c1a+1;                                       \
+    int i,j,k;                                                 \
+    T t;                                                       \
+    for (i=r1r; i<=r2r; i++) {                                 \
+        for (j=c1r; j<=c2r; j++) {                             \
+            t = 0;                                             \
+            for(k=0; k<nk; k++) {                              \
+                t += AT(a,i+dra,k+c1a) * AT(b,k+r1b,j+dcb);    \
+            }                                                  \
+            AT(r,i,j) = b*AT(r,i,j) + a*t;                     \
+        }                                                      \
+    }                                                          \
+    OK                                                         \
+}
+
+GEMM(double)
+GEMM(float)
+GEMM(TCD)
+GEMM(TCF)
+GEMM(int32_t)
+GEMM(int64_t)
+
+#define GEMM_MOD(T,M) int gemm_mod_##T(T m, VECG(T,c),VECG(int,p),MATG(T,a),MATG(T,b),MATG(T,r)) {  \
+    T a = cp[0], b = cp[1];                                    \
+    int r1a = pp[0],              c1a = pp[2],  c2a = pp[3] ;  \
+    int r1b = pp[4],              c1b = pp[6]               ;  \
+    int r1r = pp[8], r2r = pp[9], c1r = pp[10], c2r = pp[11];  \
+    int dra = r1a - r1r;                                       \
+    int dcb = c1b-c1r;                                         \
+    int nk  = c2a-c1a+1;                                       \
+    int i,j,k;                                                 \
+    T t;                                                       \
+    for (i=r1r; i<=r2r; i++) {                                 \
+        for (j=c1r; j<=c2r; j++) {                             \
+            t = 0;                                             \
+            for(k=0; k<nk; k++) {                                     \
+                t = M(t+M(AT(a,i+dra,k+c1a) * AT(b,k+r1b,j+dcb)));    \
+            }                                                         \
+            AT(r,i,j) = M(M(b*AT(r,i,j)) + M(a*t));                   \
+        }                                                      \
+    }                                                          \
+    OK                                                         \
+}
+
+#define MOD32(X) mod(X,m)
+#define MOD64(X) mod_l(X,m)
+
+GEMM_MOD(int32_t,MOD32)
+GEMM_MOD(int64_t,MOD64)
+
 ////////////////// sparse matrix-product ///////////////////////////////////////
 
 
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
index e4d95bc..bf8c5e9 100644
--- a/packages/base/src/Internal/C/lapack-aux.h
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -59,6 +59,7 @@ typedef short ftnlen;
 #define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p
 #define OCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, doublecomplex* A##p
 
+#define VECG(T,A) int A##n, T* A##p
 #define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
 
 #define KIVEC(A) int A##n, const int*A##p