4 files changed, 91 insertions, 1 deletions

diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
new file mode 100644
index 0000000..668d2f7
--- /dev/null
+++ b/lib/Data/Packed.hs
@@ -0,0 +1,70 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Data.Packed
+Copyright   :  (c) Alberto Ruiz 2006-7
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+The Vector and Matrix types and some utilities.
+
+-}
+-----------------------------------------------------------------------------
+
+module Data.Packed (
+    module Data.Packed.Vector,
+    module Data.Packed.Matrix,
+    module Data.Complex,
+    Container(..)
+) where
+
+import Data.Packed.Vector
+import Data.Packed.Matrix
+import Data.Complex
+import Data.Packed.Internal
+
+-- | conversion utilities
+class (Field e) => Container c e where
+    toComplex   :: RealFloat e => (c e, c e) -> c (Complex e)
+    fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
+    comp        :: RealFloat e => c e -> c (Complex e)
+    conj        :: RealFloat e => c (Complex e) -> c (Complex e)
+    real        :: c Double -> c e
+    complex     :: c e -> c (Complex Double)
+
+instance Container Vector Double where
+    toComplex = Data.Packed.Internal.toComplex
+    fromComplex = Data.Packed.Internal.fromComplex
+    comp = Data.Packed.Internal.comp
+    conj = Data.Packed.Internal.conj
+    real = id
+    complex = Data.Packed.comp
+
+instance Container Vector (Complex Double) where
+    toComplex = undefined -- can't match
+    fromComplex = undefined
+    comp = undefined
+    conj = undefined
+    real = Data.Packed.comp
+    complex = id
+
+instance Container Matrix Double where
+    toComplex = uncurry $ liftMatrix2 $ curry Data.Packed.toComplex
+    fromComplex z = (reshape c r, reshape c i)
+        where (r,i) = Data.Packed.fromComplex (cdat z)
+              c = cols z
+    comp = liftMatrix Data.Packed.Internal.comp
+    conj = liftMatrix Data.Packed.Internal.conj
+    real = id
+    complex = Data.Packed.comp
+
+instance Container Matrix (Complex Double) where
+    toComplex = undefined
+    fromComplex = undefined
+    comp = undefined
+    conj = undefined
+    real = Data.Packed.comp
+    complex = id
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 63ebddf..e76500b 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -24,6 +24,8 @@ import Complex
 import Control.Monad(when)
 import Data.List(transpose,intersperse)
 import Data.Maybe(fromJust)
+import Foreign.C.String
+import Foreign.C.Types
 
 -----------------------------------------------------------------
 
@@ -371,6 +373,16 @@ fromComplex z = (r,i) where
 comp :: Vector Double -> Vector (Complex Double)
 comp v = toComplex (v,constant 0 (dim v))
 
+-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = do
+    charname <- newCString filename
+    res <- createMatrix RowMajor r c
+    c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
+    --free charname  -- TO DO: free the auxiliary CString
+    return res
+foreign import ccall "aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
+
 -------------------------------------------------------------------------
 
 -- Generic definitions
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/aux.h
index 83111e5..73334e3 100644
--- a/lib/Data/Packed/Internal/aux.h
+++ b/lib/Data/Packed/Internal/aux.h
@@ -26,3 +26,5 @@ int diagR(KRVEC(d),RMAT(r));
 int diagC(KCVEC(d),CMAT(r));
 
 const char * gsl_strerror (const int gsl_errno);
+
+int matrix_fscanf(char*filename, RMAT(a));
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 404fde7..2b93348 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -25,9 +25,10 @@ module Data.Packed.Matrix (
     fromBlocks,
     flipud, fliprl,
     subMatrix, takeRows, dropRows, takeColumns, dropColumns,
+    extractRows,
     ident, diag, diagRect, takeDiag,
     liftMatrix, liftMatrix2,
-    format, dispR, readMatrix, fromArray2D
+    format, dispR, readMatrix, fromFile, fromArray2D
 ) where
 
 import Data.Packed.Internal
@@ -209,3 +210,8 @@ dispC d m = disp m (shfc d)
 -- | creates a matrix from a table of numbers.
 readMatrix :: String -> Matrix Double
 readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
+
+-- | rearranges the rows of a matrix according to the order given in a list of integers. 
+extractRows :: Field t => [Int] -> Matrix t -> Matrix t
+extractRows l m = fromRows $ extract (toRows $ m) l
+    where extract l is = [l!!i |i<-is]