8 files changed, 279 insertions, 57 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index bbc5986..c7118c1 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -20,6 +20,7 @@ imported from "Numeric.LinearAlgebra.LAPACK".
 
 module Numeric.LinearAlgebra.Algorithms (
 -- * Linear Systems
+    multiply, dot,
     linearSolve,
     inv, pinv,
     pinvTol, det, rank, rcond,
@@ -51,6 +52,8 @@ module Numeric.LinearAlgebra.Algorithms (
 -- * Misc
     ctrans,
     eps, i,
+    outer, kronecker,
+    mulH,
 -- * Util
     haussholder,
     unpackQR, unpackHess,
@@ -60,13 +63,14 @@ module Numeric.LinearAlgebra.Algorithms (
 
 import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj, (//))
 import Data.Packed
-import qualified Numeric.GSL.Matrix as GSL
 import Numeric.GSL.Vector
 import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Complex
 import Numeric.LinearAlgebra.Linear
 import Data.List(foldl1')
 import Data.Array
+import Foreign
+import Foreign.C.Types
 
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
@@ -105,6 +109,7 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     schur       :: Matrix t -> (Matrix t, Matrix t)
     -- | Conjugate transpose.
     ctrans :: Matrix t -> Matrix t
+    multiply :: Matrix t -> Matrix t -> Matrix t
 
 
 instance Field Double where
@@ -116,9 +121,10 @@ instance Field Double where
     eig = eigR
     eigSH' = eigS
     cholSH = cholS
-    qr = GSL.unpackQR . qrR
+    qr = unpackQR . qrR
     hess = unpackHess hessR
     schur = schurR
+    multiply = multiplyR3
 
 instance Field (Complex Double) where
     svd = svdC
@@ -132,6 +138,8 @@ instance Field (Complex Double) where
     qr = unpackQR . qrC
     hess = unpackHess hessC
     schur = schurC
+    multiply = mulCW -- workaround
+               -- multiplyC3
 
 -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
 --
@@ -501,3 +509,162 @@ luFact (lu,perm) | r <= c    = (l ,u ,p, s)
     u' = takeRows c (lu |*| tu)
     (|+|) = add
     (|*|) = mul
+
+--------------------------------------------------
+
+-- | euclidean inner product
+dot :: (Field t) => Vector t -> Vector t -> t
+dot u v = multiply r c  @@> (0,0)
+    where r = asRow u
+          c = asColumn v
+
+
+{- | Outer product of two vectors.
+
+@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
+(3><3)
+ [  5.0, 2.0, 3.0
+ , 10.0, 4.0, 6.0
+ , 15.0, 6.0, 9.0 ]@
+-}
+outer :: (Field t) => Vector t -> Vector t -> Matrix t
+outer u v = asColumn u `multiply` asRow v
+
+{- | Kronecker product of two matrices.
+
+@m1=(2><3)
+ [ 1.0,  2.0, 0.0
+ , 0.0, -1.0, 3.0 ]
+m2=(4><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0 ]@
+
+@\> kronecker m1 m2
+(8><9)
+ [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
+ ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
+ ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
+ , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
+ ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
+ ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
+ ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
+ ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
+-}
+kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t
+kronecker a b = fromBlocks
+              . partit (cols a)
+              . map (reshape (cols b))
+              . toRows
+              $ flatten a `outer` flatten b
+
+---------------------------------------------------------------------
+-- reference multiply
+---------------------------------------------------------------------
+
+mulH a b = fromLists [[ dot ai bj | bj <- toColumns b] | ai <- toRows a ]
+    where dot u v = sum $ zipWith (*) (toList u) (toList v)
+
+-----------------------------------------------------------------------------------
+-- workaround
+-----------------------------------------------------------------------------------
+
+mulCW a b = toComplex (rr,ri)
+    where rr = multiply ar br `sub` multiply ai bi
+          ri = multiply ar bi `add` multiply ai br
+          (ar,ai) = fromComplex a
+          (br,bi) = fromComplex b
+
+-----------------------------------------------------------------------------------
+-- Direct CBLAS
+-----------------------------------------------------------------------------------
+
+newtype CBLASOrder = CBLASOrder CInt deriving (Eq, Show)
+newtype CBLASTrans = CBLASTrans CInt deriving (Eq, Show)
+
+rowMajor, colMajor :: CBLASOrder
+rowMajor = CBLASOrder 101
+colMajor = CBLASOrder 102
+
+noTrans, trans', conjTrans :: CBLASTrans
+noTrans   = CBLASTrans 111
+trans'     = CBLASTrans 112
+conjTrans = CBLASTrans 113
+
+foreign import ccall "cblas.h cblas_dgemm"
+    dgemm  :: CBLASOrder -> CBLASTrans -> CBLASTrans -> CInt -> CInt -> CInt -> Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Double -> Ptr Double -> CInt -> IO ()
+
+
+multiplyR3 :: Matrix Double -> Matrix Double -> Matrix Double
+multiplyR3 a b = multiply3 dgemm "cblas_dgemm" (fmat a) (fmat b)
+    where
+        multiply3 f st a b
+            | cols a == rows b = unsafePerformIO $ do
+                s <- createMatrix ColumnMajor (rows a) (cols b)
+                let g ar ac ap br bc bp rr rc rp = f colMajor noTrans noTrans ar bc ac 1 ap ar bp br 0 rp rr >> return 0
+                app3 g mat a mat b mat s st
+                return s
+            | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
+
+
+foreign import ccall "cblas.h cblas_zgemm"
+    zgemm  :: CBLASOrder -> CBLASTrans -> CBLASTrans -> CInt -> CInt -> CInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> CInt -> Ptr (Complex Double) -> CInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> CInt -> IO ()
+
+multiplyC3 :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+multiplyC3 a b = unsafePerformIO $ multiply3 zgemm "cblas_zgemm" (fmat a) (fmat b)
+    where
+        multiply3 f st a b
+            | cols a == rows b = do
+                s <- createMatrix ColumnMajor (rows a) (cols b)
+                palpha <- new 1
+                pbeta  <- new 0
+                let g ar ac ap br bc bp rr rc rp = f colMajor noTrans noTrans ar bc ac palpha ap ar bp br pbeta rp rr >> return 0
+                app3 g mat a mat b mat s st
+                free palpha
+                free pbeta
+                return s
+                -- if toLists s== toLists s then return s else error $ "HORROR " ++ (show (toLists s))
+            | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
+
+-----------------------------------------------------------------------------------
+-- BLAS via auxiliary C
+-----------------------------------------------------------------------------------
+
+foreign import ccall "multiply.h multiplyR2" dgemmc :: TMMM
+foreign import ccall "multiply.h multiplyC2" zgemmc :: TCMCMCM
+
+multiply2 f st a b
+    | cols a == rows b = unsafePerformIO $ do
+        s <- createMatrix ColumnMajor (rows a) (cols b)
+        app3 f mat a mat b mat s st 
+        if toLists s== toLists s then return s else error $ "AYYY " ++ (show (toLists s))
+    | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
+
+multiplyR2 :: Matrix Double -> Matrix Double -> Matrix Double
+multiplyR2 a b = multiply2 dgemmc "dgemmc" (fmat a) (fmat b)
+
+multiplyC2 :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+multiplyC2 a b = multiply2 zgemmc "zgemmc" (fmat a) (fmat b)
+
+-----------------------------------------------------------------------------------
+-- direct C multiplication
+-----------------------------------------------------------------------------------
+
+foreign import ccall "multiply.h multiplyR" cmultiplyR :: TMMM
+foreign import ccall "multiply.h multiplyC" cmultiplyC :: TCMCMCM
+
+cmultiply f st a b
+--    | cols a == rows b = 
+      = unsafePerformIO $ do
+        s <- createMatrix RowMajor (rows a) (cols b)
+        app3 f mat a mat b mat s st
+        if toLists s== toLists s then return s else error $ "BRUTAL " ++ (show (toLists s))
+        -- return s
+--    | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
+
+multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
+multiplyR a b = cmultiply cmultiplyR "cmultiplyR" (cmat a) (cmat b)
+
+multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+multiplyC a b = cmultiply cmultiplyC "cmultiplyR" (cmat a) (cmat b)
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
index 4a9b309..0ae9698 100644
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -29,7 +29,7 @@ import Numeric.LinearAlgebra.Algorithms
 class Mul a b c | a b -> c where
  infixl 7 <>
  -- | matrix product
- (<>) :: Element t => a t -> b t -> c t
+ (<>) :: Field t => a t -> b t -> c t
 
 instance Mul Matrix Matrix Matrix where
     (<>) = multiply
@@ -43,7 +43,7 @@ instance Mul Vector Matrix Vector where
 ---------------------------------------------------
 
 -- | @u \<.\> v = dot u v@
-(<.>) :: (Element t) => Vector t -> Vector t -> t
+(<.>) :: (Field t) => Vector t -> Vector t -> t
 infixl 7 <.>
 (<.>) = dot
 
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 310f6ee..0dccea2 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -814,3 +814,77 @@ int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
     free(auxipiv);
     OK
 }
+
+////////////////////////////////////////////////////////////
+
+int multiplyR(KDMAT(a),KDMAT(b),DMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    int i,j,k;
+    for (i=0;i<ar;i++) {
+        for(j=0;j<bc;j++) {
+            double temp = 0;
+            for(k=0;k<ac;k++) {
+                temp += ap[i*ac+k]*bp[k*bc+j];
+            }
+            rp[i*rc+j] = temp;
+        }
+    }
+    OK
+}
+
+int multiplyC(KCMAT(a),KCMAT(b),CMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    int i,j,k;
+    for (i=0;i<ar;i++) {
+        for(j=0;j<bc;j++) {
+            doublecomplex temp = {0,0};
+            for(k=0;k<ac;k++) {
+                doublecomplex aik = ((doublecomplex*)ap)[i*ac+k];
+                doublecomplex bkj = ((doublecomplex*)bp)[k*bc+j];
+                //double w = aik.r+aik.i+bkj.r+bkj.i;
+                //if (w>w) exit(1);
+                //printf("%d",w>w);
+                temp.r += aik.r * bkj.r - aik.i * bkj.i;
+                temp.i += aik.r * bkj.i + aik.i * bkj.r;
+                //printf("%f %f %f %f \n",aik.r,aik.i,bkj.r,bkj.i);
+                //printf("%f %f %f \n",w,temp.r,temp.i);
+
+            }
+            ((doublecomplex*)rp)[i*rc+j] = temp;
+            //printf("%f %f\n",temp.r,temp.i);
+        }
+    }
+    OK
+}
+
+void dgemm_(char *, char *, integer *, integer *, integer *,
+           double *, const double *, integer *, const double *,
+           integer *, double *, double *, integer *);
+
+void zgemm_(char *, char *, integer *, integer *, integer *,
+           doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
+           integer *, doublecomplex *, doublecomplex *, integer *);
+
+
+int multiplyR2(KDMAT(a),KDMAT(b),DMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    double alpha = 1;
+    double beta = 0;
+    integer m = ar;
+    integer n = bc;
+    integer k = ac;
+    int i,j;
+    dgemm_("N","N",&m,&n,&k,&alpha,ap,&m,bp,&k,&beta,rp,&m);
+    OK
+}
+
+int multiplyC2(KCMAT(a),KCMAT(b),CMAT(r)) {
+    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    integer m = ar;
+    integer n = bc;
+    integer k = ac;
+    doublecomplex alpha = {1,0};
+    doublecomplex beta = {0,0};
+    zgemm_("N","N",&m,&n,&k,&alpha,(doublecomplex*)ap,&m,(doublecomplex*)bp,&k,&beta,(doublecomplex*)rp,&m);
+    OK
+}
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index 79e52be..c0361a6 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -84,3 +84,9 @@ int schur_l_C(KCMAT(a), CMAT(u), CMAT(s));
 
 int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r));
 int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r));
+
+int multiplyR(KDMAT(a),KDMAT(b),DMAT(r));
+int multiplyC(KCMAT(a),KCMAT(b),CMAT(r));
+
+int multiplyR2(KDMAT(a),KDMAT(b),DMAT(r));
+int multiplyC2(KCMAT(a),KCMAT(b),CMAT(r));
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 0ddbb55..1bf8b04 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -15,12 +15,11 @@ Basic optimized operations on vectors and matrices.
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Linear (
-    Linear(..),
-    multiply, dot, outer, kronecker
+    Linear(..)
 ) where
 
 
-import Data.Packed.Internal(multiply,partit)
+import Data.Packed.Internal(partit)
 import Data.Packed
 import Numeric.GSL.Vector
 import Complex
@@ -69,52 +68,3 @@ instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
     equal a b = cols a == cols b && flatten a `equal` flatten b
-
---------------------------------------------------
-
--- | euclidean inner product
-dot :: (Element t) => Vector t -> Vector t -> t
-dot u v = multiply r c  @@> (0,0)
-    where r = asRow u
-          c = asColumn v
-
-
-{- | Outer product of two vectors.
-
-@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
-(3><3)
- [  5.0, 2.0, 3.0
- , 10.0, 4.0, 6.0
- , 15.0, 6.0, 9.0 ]@
--}
-outer :: (Element t) => Vector t -> Vector t -> Matrix t
-outer u v = asColumn u `multiply` asRow v
-
-{- | Kronecker product of two matrices.
-
-@m1=(2><3)
- [ 1.0,  2.0, 0.0
- , 0.0, -1.0, 3.0 ]
-m2=(4><3)
- [  1.0,  2.0,  3.0
- ,  4.0,  5.0,  6.0
- ,  7.0,  8.0,  9.0
- , 10.0, 11.0, 12.0 ]@
-
-@\> kronecker m1 m2
-(8><9)
- [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
- ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
- ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
- , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
- ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
- ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
- ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
- ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
--}
-kronecker :: (Element t) => Matrix t -> Matrix t -> Matrix t
-kronecker a b = fromBlocks
-              . partit (cols a)
-              . map (reshape (cols b))
-              . toRows
-              $ flatten a `outer` flatten b
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 7b28075..07b9f63 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -123,6 +123,11 @@ runTests :: Int  -- ^ maximum dimension
 runTests n = do
     setErrorHandlerOff
     let test p = qCheck n p
+    putStrLn "------ mult"
+    test (multProp1  . rConsist)
+    test (multProp1  . cConsist)
+    test (multProp2  . rConsist)
+    test (multProp2  . cConsist)
     putStrLn "------ lu"
     test (luProp    . rM)
     test (luProp    . cM)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
index af486c8..e7fecf2 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -20,6 +20,7 @@ module Numeric.LinearAlgebra.Tests.Instances(
     WC(..),     rWC,cWC,
     SqWC(..),   rSqWC, cSqWC,
     PosDef(..), rPosDef, cPosDef,
+    Consistent(..), rConsist, cConsist,
     RM,CM, rM,cM
 ) where
 
@@ -116,6 +117,19 @@ instance (Field a, Arbitrary a) => Arbitrary (PosDef a) where
         return $ PosDef (0.5 .* p + 0.5 .* ctrans p)
     coarbitrary = undefined
 
+-- a pair of matrices that can be multiplied
+newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
+    arbitrary = do
+        n <- chooseDim
+        k <- chooseDim
+        m <- chooseDim
+        la <- vector (n*k)
+        lb <- vector (k*m)
+        return $ Consistent ((n><k) la, (k><m) lb)
+    coarbitrary = undefined
+
+
 type RM = Matrix Double
 type CM = Matrix (Complex Double)
 
@@ -140,3 +154,5 @@ cSqWC (SqWC m) = m :: CM
 rPosDef (PosDef m) = m :: RM
 cPosDef (PosDef m) = m :: CM
 
+rConsist (Consistent (a,b)) = (a,b::RM)
+cConsist (Consistent (a,b)) = (a,b::CM)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 55e9a1b..5663b86 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -34,7 +34,8 @@ module Numeric.LinearAlgebra.Tests.Properties (
     hessProp,
     schurProp1, schurProp2,
     cholProp,
-    expmDiagProp
+    expmDiagProp,
+    multProp1, multProp2
 ) where
 
 import Numeric.LinearAlgebra
@@ -151,3 +152,6 @@ cholProp m = m |~| ctrans c <> c && upperTriang c
 expmDiagProp m = expm (logm m) :~ 7 ~: complex m
     where logm m = matFunc log m
 
+multProp1 (a,b) = a <> b |~| mulH a b
+
+multProp2 (a,b) = trans (a <> b) |~| trans b <> trans a