5 files changed, 548 insertions, 0 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
new file mode 100644
index 0000000..2c57c07
--- /dev/null
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -0,0 +1,187 @@
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Fundamental types
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Internal.Matrix where
+
+import Data.Packed.Internal.Vector
+
+import Foreign hiding (xor)
+import Complex
+import Control.Monad(when)
+import Debug.Trace
+import Data.List(transpose,intersperse)
+import Data.Typeable
+import Data.Maybe(fromJust)
+
+debug x = trace (show x) x
+
+
+data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
+
+-- | 2D array
+data Matrix t = M { rows    :: Int
+                  , cols    :: Int
+                  , dat     :: Vector t
+                  , tdat    :: Vector t
+                  , isTrans :: Bool
+                  , order   :: MatrixOrder
+                  } deriving Typeable
+
+xor a b = a && not b || b && not a
+
+fortran m = order m == ColumnMajor
+
+cdat m = if fortran m `xor` isTrans m then tdat m else dat m
+fdat m = if fortran m `xor` isTrans m then dat m else tdat m
+
+trans m = m { rows = cols m
+            , cols = rows m
+            , isTrans = not (isTrans m)
+            }
+
+type Mt t s = Int -> Int -> Ptr t -> s
+infixr 6 ::>
+type t ::> s = Mt t s
+
+mat d m f = f (rows m) (cols m) (ptr (d m))
+
+instance (Show a, Storable a) => (Show (Matrix a)) where
+    show m = (sizes++) . dsp . map (map show) . toLists $ m
+        where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"
+
+partit :: Int -> [a] -> [[a]]
+partit _ [] = []
+partit n l  = take n l : partit n (drop n l)
+
+toLists m | fortran m = transpose $ partit (rows m) . toList . dat $ m
+          | otherwise = partit (cols m) . toList . dat $ m
+
+dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
+    where
+        mt = transpose as
+        longs = map (maximum . map length) mt
+        mtp = zipWith (\a b -> map (pad a) b) longs mt
+        pad n str = replicate (n - length str) ' ' ++ str
+        unwords' = concat . intersperse ", "
+
+matrixFromVector RowMajor c v =
+    M { rows = r
+      , cols = c
+      , dat  = v
+      , tdat = transdata c v r
+      , order = RowMajor
+      , isTrans = False
+      } where r = dim v `div` c -- TODO check mod=0
+
+matrixFromVector ColumnMajor c v =
+    M { rows = r
+      , cols = c
+      , dat  = v
+      , tdat = transdata r v c
+      , order = ColumnMajor
+      , isTrans = False
+      } where r = dim v `div` c -- TODO check mod=0
+
+createMatrix order r c = do
+    p <- createVector (r*c)
+    return (matrixFromVector order c p)
+
+transdataG :: Storable a => Int -> Vector a -> Int -> Vector a 
+transdataG c1 d c2 = fromList . concat . transpose . partit c1 . toList $ d
+
+transdataR :: Int -> Vector Double -> Int -> Vector Double
+transdataR = transdataAux ctransR
+
+transdataC :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double)
+transdataC = transdataAux ctransC
+
+transdataAux fun c1 d c2 = unsafePerformIO $ do
+    v <- createVector (dim d)
+    let r1 = dim d `div` c1
+        r2 = dim d `div` c2
+    fun r1 c1 (ptr d) r2 c2 (ptr v) // check "transdataAux" [d]
+    --putStrLn "---> transdataAux"
+    return v
+
+foreign import ccall safe "aux.h transR"
+    ctransR :: Double ::> Double ::> IO Int
+foreign import ccall safe "aux.h transC"
+    ctransC :: Complex Double ::> Complex Double ::> IO Int
+
+transdata :: Field a => Int -> Vector a -> Int -> Vector a
+transdata c1 d c2 | isReal baseOf d = scast $ transdataR c1 (scast d) c2
+                  | isComp baseOf d = scast $ transdataC c1 (scast d) c2
+                  | otherwise       = transdataG c1 d c2
+
+--transdata :: Storable a => Int -> Vector a -> Int -> Vector a 
+--transdata = transdataG
+--{-# RULES "transdataR" transdata=transdataR #-}
+--{-# RULES "transdataC" transdata=transdataC #-}
+
+-- | extracts the rows of a matrix as a list of vectors
+toRows :: Storable t => Matrix t -> [Vector t]
+toRows m = toRows' 0 where
+    v = cdat m
+    r = rows m
+    c = cols m
+    toRows' k | k == r*c  = []
+              | otherwise = subVector k c v : toRows' (k+c)
+
+------------------------------------------------------------------
+
+dotL a b = sum (zipWith (*) a b)
+
+multiplyL a b = [[dotL x y | y <- transpose b] | x <- a]
+
+transL m = matrixFromVector RowMajor (rows m) $ transdataG (cols m) (cdat m) (rows m)
+
+multiplyG a b = matrixFromVector RowMajor (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b)
+
+------------------------------------------------------------------
+
+gmatC m f | fortran m =
+                if (isTrans m)
+                    then f 0 (rows m) (cols m) (ptr (dat m))
+                    else f 1 (cols m) (rows m) (ptr (dat m))
+         | otherwise =
+                if isTrans m
+                    then f 1 (cols m) (rows m) (ptr (dat m))
+                    else f 0 (rows m) (cols m) (ptr (dat m))
+
+
+multiplyAux order fun a b = unsafePerformIO $ do
+    when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++
+                                      show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
+    r <- createMatrix order (rows a) (cols b)
+    fun // gmatC a // gmatC b // mat dat r // check "multiplyAux" [dat a, dat b]
+    return r
+
+foreign import ccall safe "aux.h multiplyR"
+    cmultiplyR :: Int -> Double ::> (Int -> Double ::> (Double ::> IO Int))
+
+foreign import ccall safe "aux.h multiplyC"
+    cmultiplyC :: Int -> Complex Double ::> (Int -> Complex Double ::> (Complex Double ::> IO Int))
+
+multiply :: (Num a, Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a
+multiply RowMajor a b    = multiplyD RowMajor a b
+multiply ColumnMajor a b = trans $ multiplyT ColumnMajor a b
+
+multiplyT order a b = multiplyD order (trans b) (trans a)
+
+multiplyD order a b 
+    | isReal (baseOf.dat) a = scast $ multiplyAux order cmultiplyR (scast a) (scast b)
+    | isComp (baseOf.dat) a = scast $ multiplyAux order cmultiplyC (scast a) (scast b)
+    | otherwise             = multiplyG a b
+
diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs
new file mode 100644
index 0000000..11101a9
--- /dev/null
+++ b/lib/Data/Packed/Internal/Tensor.hs
@@ -0,0 +1,32 @@
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Tensor
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Fundamental types
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Internal.Tensor where
+
+import Data.Packed.Internal.Vector
+import Data.Packed.Internal.Matrix
+
+
+data IdxTp = Covariant | Contravariant deriving Show
+
+data Tensor t = T { dims   :: [(Int,(IdxTp,String))]
+                  , ten    :: Vector t
+                  } deriving Show
+
+rank = length . dims
+
+outer u v = dat (multiply RowMajor r c)
+    where r = matrixFromVector RowMajor 1 u
+          c = matrixFromVector RowMajor (dim v) v
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
new file mode 100644
index 0000000..7dcefeb
--- /dev/null
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -0,0 +1,164 @@
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Fundamental types
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Internal.Vector where
+
+import Foreign
+import Complex
+import Control.Monad(when)
+import Debug.Trace
+import Data.List(transpose,intersperse)
+import Data.Typeable
+import Data.Maybe(fromJust)
+
+debug x = trace (show x) x
+
+----------------------------------------------------------------------
+instance (Storable a, RealFloat a) => Storable (Complex a) where    --
+    alignment x = alignment (realPart x)                            --
+    sizeOf x    = 2 * sizeOf (realPart x)                           --
+    peek p = do                                                     --
+        [re,im] <- peekArray 2 (castPtr p)                          --
+        return (re :+ im)                                           --
+    poke p (a :+ b) = pokeArray (castPtr p) [a,b]                   --
+----------------------------------------------------------------------
+
+(//) :: x -> (x -> y) -> y
+infixl 0 //
+(//) = flip ($)
+
+check msg ls f = do
+    err <- f
+    when (err/=0) (error msg)
+    mapM_ (touchForeignPtr . fptr) ls
+    return ()
+
+class (Storable a, Typeable a) => Field a where
+instance (Storable a, Typeable a) => Field a where
+
+isReal w x   = typeOf (undefined :: Double) == typeOf (w x)
+isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
+baseOf v = (v `at` 0)
+
+scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b
+scast = fromJust . cast
+
+
+
+----------------------------------------------------------------------
+
+data Vector t = V { dim  :: Int
+                  , fptr :: ForeignPtr t
+                  , ptr  :: Ptr t
+                  } deriving Typeable
+
+type Vc t s = Int -> Ptr t -> s
+infixr 5 :>
+type t :> s = Vc t s
+
+vec :: Vector t -> (Vc t s) -> s
+vec v f = f (dim v) (ptr v)
+
+createVector :: Storable a => Int -> IO (Vector a)
+createVector n = do
+    when (n <= 0) $ error ("trying to createVector of dim "++show n)
+    fp <- mallocForeignPtrArray n
+    let p = unsafeForeignPtrToPtr fp
+    --putStrLn ("\n---------> V"++show n)
+    return $ V n fp p
+
+fromList :: Storable a => [a] -> Vector a
+fromList l = unsafePerformIO $ do
+    v <- createVector (length l)
+    let f _ p = pokeArray p l >> return 0
+    f // vec v // check "fromList" []
+    return v
+
+toList :: Storable a => Vector a -> [a]
+toList v = unsafePerformIO $ peekArray (dim v) (ptr v)
+
+n # l = if length l == n then fromList l else error "# with wrong size"
+
+at' :: Storable a => Vector a -> Int -> a
+at' v n = unsafePerformIO $ peekElemOff (ptr v) n
+
+at :: Storable a => Vector a -> Int -> a
+at v n | n >= 0 && n < dim v = at' v n
+       | otherwise          = error "vector index out of range"
+
+instance (Show a, Storable a) => (Show (Vector a)) where
+    show v = (show (dim v))++" # " ++ show (toList v)
+
+-- | creates a Vector taking a number of consecutive toList from another Vector
+subVector :: Storable t => Int       -- ^ index of the starting element
+                        -> Int       -- ^ number of toList to extract
+                        -> Vector t  -- ^ source
+                        -> Vector t  -- ^ result
+subVector k l (v@V {dim=n, ptr=p, fptr=fp})
+    | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
+    | otherwise = unsafePerformIO $ do
+        r <- createVector l
+        let f = copyArray (ptr r) (advancePtr p k) l >> return 0
+        f // check "subVector" [v]
+        return r
+
+subVector' k l (v@V {dim=n, ptr=p, fptr=fp})
+    | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
+    | otherwise = v {dim=l, ptr=advancePtr p k}
+
+
+{-
+-- | creates a new Vector by joining a list of Vectors
+join :: Field t => [Vector t] -> Vector t
+join [] = error "joining an empty list"
+join as = unsafePerformIO $ do
+    let tot = sum (map size as)
+    p <- mallocForeignPtrArray tot
+    withForeignPtr p $ \p ->
+        joiner as tot p
+    return (V tot p)
+  where joiner [] _ _ = return ()
+        joiner (V n b : cs) _ p = do
+            withForeignPtr b  $ \b' -> copyArray p b' n
+            joiner cs 0 (advancePtr p n)
+-}
+
+
+constantG n x = fromList (replicate n x)
+
+constantR :: Int -> Double -> Vector Double
+constantR = constantAux cconstantR
+
+constantC :: Int -> Complex Double -> Vector (Complex Double)
+constantC = constantAux cconstantC
+
+constantAux fun n x = unsafePerformIO $ do
+    v <- createVector n
+    px <- newArray [x]
+    fun px // vec v // check "constantAux" []
+    free px
+    return v
+
+foreign import ccall safe "aux.h constantR"
+    cconstantR :: Ptr Double -> Double :> IO Int
+
+foreign import ccall safe "aux.h constantC"
+    cconstantC :: Ptr (Complex Double) -> Complex Double :> IO Int
+
+constant :: Field a => Int -> a -> Vector a
+constant n x | isReal id x = scast $ constantR n (scast x)
+             | isComp id x = scast $ constantC n (scast x)
+             | otherwise   = constantG n x
+
diff --git a/lib/Data/Packed/Internal/aux.c b/lib/Data/Packed/Internal/aux.c
new file mode 100644
index 0000000..da36035
--- /dev/null
+++ b/lib/Data/Packed/Internal/aux.c
@@ -0,0 +1,144 @@
+#include "aux.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_fft_complex.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_deriv.h>
+#include <gsl/gsl_poly.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <string.h>
+#include <stdio.h>
+
+#define MACRO(B) do {B} while (0)
+#define ERROR(CODE) MACRO(return CODE;)
+#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+#ifdef DBG
+#define DEBUGMSG(M) printf("GSL Wrapper "M": "); size_t t0 = time(NULL);
+#define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
+#else
+#define DEBUGMSG(M)
+#define OK return 0;
+#endif
+
+
+#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
+
+#ifdef DBG
+#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
+#else
+#define DEBUGMAT(MSG,X)
+#endif
+
+#ifdef DBG
+#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
+#else
+#define DEBUGVEC(MSG,X)
+#endif
+
+
+#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
+#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
+#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
+#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)
+#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
+#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
+#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n)
+#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)
+
+#define V(a) (&a.vector)
+#define M(a) (&a.matrix)
+
+#define GCVEC(A) int A##n, gsl_complex*A##p
+#define KGCVEC(A) int A##n, const gsl_complex*A##p
+
+#define BAD_SIZE 1000
+#define BAD_CODE 1001
+#define MEM      1002
+#define BAD_FILE 1003
+
+
+
+int transR(KRMAT(x),RMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transR");
+    KDMVIEW(x);
+    DMVIEW(t);
+    int res = gsl_matrix_transpose_memcpy(M(t),M(x));
+    CHECK(res,res);
+    OK
+}
+
+int transC(KCMAT(x),CMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transC");
+    KCMVIEW(x);
+    CMVIEW(t);
+    int res = gsl_matrix_complex_transpose_memcpy(M(t),M(x));
+    CHECK(res,res);
+    OK
+}
+
+
+int constantR(double * pval, RVEC(r)) {
+    DEBUGMSG("constantR")
+    int k;
+    double val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int constantC(gsl_complex* pval, CVEC(r)) {
+    DEBUGMSG("constantC")
+    int k;
+    gsl_complex val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int multiplyR(int ta, KRMAT(a), int tb, KRMAT(b),RMAT(r)) {
+    //printf("%d %d %d %d %d %d\n",ar,ac,br,bc,rr,rc);
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    DEBUGMSG("multiplyR (gsl_blas_dgemm)");
+    KDMVIEW(a);
+    KDMVIEW(b);
+    DMVIEW(r);
+    int res = gsl_blas_dgemm(
+         ta?CblasTrans:CblasNoTrans,
+         tb?CblasTrans:CblasNoTrans,
+         1.0, M(a), M(b),
+         0.0, M(r));
+    CHECK(res,res);
+    OK
+}
+
+int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    DEBUGMSG("multiplyC (gsl_blas_zgemm)");
+    KCMVIEW(a);
+    KCMVIEW(b);
+    CMVIEW(r);
+    gsl_complex alpha, beta;
+    GSL_SET_COMPLEX(&alpha,1.,0.);
+    GSL_SET_COMPLEX(&beta,0.,0.);
+    int res = gsl_blas_zgemm(
+         ta?CblasTrans:CblasNoTrans,
+         tb?CblasTrans:CblasNoTrans,
+         alpha, M(a), M(b),
+         beta, M(r));
+    CHECK(res,res);
+    OK
+}
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/aux.h
new file mode 100644
index 0000000..f45b55a
--- /dev/null
+++ b/lib/Data/Packed/Internal/aux.h
@@ -0,0 +1,21 @@
+#include <gsl/gsl_complex.h>
+
+#define RVEC(A) int A##n, double*A##p
+#define RMAT(A) int A##r, int A##c, double* A##p
+#define KRVEC(A) int A##n, const double*A##p
+#define KRMAT(A) int A##r, int A##c, const double* A##p
+
+#define CVEC(A) int A##n, gsl_complex*A##p
+#define CMAT(A) int A##r, int A##c, gsl_complex* A##p
+#define KCVEC(A) int A##n, const gsl_complex*A##p
+#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
+
+
+int transR(KRMAT(x),RMAT(t));
+int transC(KCMAT(x),CMAT(t));
+
+int constantR(double *val     , RVEC(r));
+int constantC(gsl_complex *val, CVEC(r));
+
+int multiplyR(int ta, KRMAT(a), int tb,  KRMAT(b),RMAT(r));
+int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r));