From 7430630fa0504296b796223e01cbd417b88650ef Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 4 Jun 2007 19:10:28 +0000 Subject: separation of Internal --- lib/Data/Packed/Internal.hs | 293 ------------------------------------- lib/Data/Packed/Internal/Matrix.hs | 187 +++++++++++++++++++++++ lib/Data/Packed/Internal/Tensor.hs | 32 ++++ lib/Data/Packed/Internal/Vector.hs | 164 +++++++++++++++++++++ lib/Data/Packed/Internal/aux.c | 144 ++++++++++++++++++ lib/Data/Packed/Internal/aux.h | 21 +++ lib/Data/Packed/aux.c | 144 ------------------ lib/Data/Packed/aux.h | 21 --- 8 files changed, 548 insertions(+), 458 deletions(-) delete mode 100644 lib/Data/Packed/Internal.hs create mode 100644 lib/Data/Packed/Internal/Matrix.hs create mode 100644 lib/Data/Packed/Internal/Tensor.hs create mode 100644 lib/Data/Packed/Internal/Vector.hs create mode 100644 lib/Data/Packed/Internal/aux.c create mode 100644 lib/Data/Packed/Internal/aux.h delete mode 100644 lib/Data/Packed/aux.c delete mode 100644 lib/Data/Packed/aux.h (limited to 'lib/Data/Packed') diff --git a/lib/Data/Packed/Internal.hs b/lib/Data/Packed/Internal.hs deleted file mode 100644 index b06f044..0000000 --- a/lib/Data/Packed/Internal.hs +++ /dev/null @@ -1,293 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable (uses FFI) --- --- Fundamental types --- ------------------------------------------------------------------------------ - -module Data.Packed.Internal where - -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 - ----------------------------------------------------------------------- -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 () - ----------------------------------------------------------------------- - -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) - ------------------------------------------------------------------------- - -data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) - --- | 2D array -data Matrix t = M { rows :: Int - , cols :: Int - , cmat :: Vector t - , fmat :: Vector t - , isTrans :: Bool - , order :: MatrixOrder - } deriving Typeable - -xor a b = a && not b || b && not a - -fortran m = order m == ColumnMajor - -dat m = if fortran m `xor` isTrans m then fmat m else cmat m - -pref m = if fortran m then fmat m else cmat 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 :: Matrix t -> (Mt t s) -> s -mat m f = f (rows m) (cols m) (ptr (dat m)) - -gmat m f | fortran m = - if (isTrans m) - then f 0 (rows m) (cols m) (ptr (fmat m)) - else f 1 (cols m) (rows m) (ptr (fmat m)) - | otherwise = - if isTrans m - then f 1 (cols m) (rows m) (ptr (cmat m)) - else f 0 (rows m) (cols m) (ptr (cmat 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 = partit (cols m) . toList . cmat $ 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 - , cmat = v - , fmat = 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 - , fmat = v - , cmat = transdata c v r - , 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 - - -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 - -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 #-} - ------------------------------------------------------------------- - -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 - ------------------------------------------------------------------- - -dotL a b = sum (zipWith (*) a b) - -multiplyL a b = [[dotL x y | y <- transpose b] | x <- a] - -transL m = m {rows = cols m, cols = rows m, cmat = v, fmat = cmat m} - where v = transdataG (cols m) (cmat m) (rows m) - ------------------------------------------------------------------- - -multiplyG a b = matrixFromVector RowMajor (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b) - -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 // gmat a // gmat b // mat r // check "multiplyAux" [pref a, pref 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 - --------------------------------------------------------------------- - -data IdxTp = Covariant | Contravariant - --- | multidimensional array -data Tensor t = T { rank :: Int - , dims :: [Int] - , idxNm :: [String] - , idxTp :: [IdxTp] - , ten :: Vector t - } - 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 +-- 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 +-- 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 +-- 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 + +#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)); diff --git a/lib/Data/Packed/aux.c b/lib/Data/Packed/aux.c deleted file mode 100644 index da36035..0000000 --- a/lib/Data/Packed/aux.c +++ /dev/null @@ -1,144 +0,0 @@ -#include "aux.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#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 - -#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)); -- cgit v1.2.3