use omat for multiplyI

2 files changed, 15 insertions, 20 deletions

diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c
index a977d5f..ac03120 100644
--- a/packages/base/src/C/lapack-aux.c
+++ b/packages/base/src/C/lapack-aux.c
@@ -1287,28 +1287,19 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
     OK
 }
 
-int multiplyI(int ta, int tb, KIMAT(a), KIMAT(b), IMAT(r)) {
-    int i,j,k;
-    int n;
-    int ai,ak,bk,bj;
-
-    n = ta ? ar : ac;
-    
-    if (ta==0) { ai =  1; ak = ar; } else { ai = ar; ak =  1; }
-    if (tb==0) { bk =  1; bj = br; } else { bk = br; bj =  1; }
-    
-    for (i=0;i<rr;i++) {
-        for (j=0;j<rc;j++) {
-            rp[i+rr*j] = 0;
-            for (k=0; k<n; k++) {
-                rp[i+rr*j] += ap[ai*i+ak*k] * bp[bk*k+bj*j];
-            }
+
+int multiplyI(KOIMAT(a), KOIMAT(b), OIMAT(r)) {
+    { TRAV(r,i,j) {
+        int k;
+        AT(r,i,j) = 0;
+        for (k=0;k<ac;k++) {
+            AT(r,i,j) += AT(a,i,k) * AT(b,k,j);
         }
+      }
     }
     OK
 }
 
-
 //////////////////// transpose /////////////////////////
 
 int transF(KFMAT(x),FMAT(t)) {
diff --git a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
index 219d996..6fb2b13 100644
--- a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
@@ -58,7 +58,7 @@ foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM
 foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM
 foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM
 foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM
-foreign import ccall unsafe "multiplyI" c_multiplyI :: CInt -> CInt -> (CM CInt (CM CInt (CM CInt (IO CInt))))
+foreign import ccall unsafe "multiplyI" c_multiplyI :: OM CInt (OM CInt (OM CInt (IO CInt)))
 
 isT Matrix{order = ColumnMajor} = 0
 isT Matrix{order = RowMajor} = 1
@@ -90,8 +90,12 @@ multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex
 multiplyQ a b = multiplyAux cgemmc "cgemmc" a b
 
 multiplyI :: Matrix CInt -> Matrix CInt -> Matrix CInt
-multiplyI = multiplyAux c_multiplyI "c_multiplyI"
-
+multiplyI a b = unsafePerformIO $ do
+    when (cols a /= rows b) $ error $
+        "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
+    s <- createMatrix ColumnMajor (rows a) (cols b)
+    app3 c_multiplyI omat a omat b omat s "c_multiplyI"
+    return s
 
 -----------------------------------------------------------------------------
 foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM