separation ok

16 files changed, 256 insertions, 131 deletions

diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal
index afd678f..8496663 100644
--- a/packages/hmatrix/hmatrix.cabal
+++ b/packages/hmatrix/hmatrix.cabal
@@ -6,11 +6,11 @@ Author:             Alberto Ruiz
 Maintainer:         Alberto Ruiz <aruiz@um.es>
 Stability:          provisional
 Homepage:           https://github.com/albertoruiz/hmatrix
-Synopsis:           Linear algebra and numerical computation
+Synopsis:           Numerical computation
 Description:        Purely functional interface to basic linear algebra
                     and other numerical computations, internally implemented using
                     GSL, BLAS and LAPACK.
-                    
+
 Category:           Math
 tested-with:        GHC ==7.8
 
@@ -45,8 +45,7 @@ extra-source-files: examples/deriv.hs
                     examples/bool.hs
                     examples/multiply.hs
 
-extra-source-files: src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,
-                    src/Numeric/GSL/gsl-ode.c
+extra-source-files: src/Numeric/GSL/gsl-ode.c
 
 flag dd
     description:    svd = zgesdd
@@ -74,23 +73,15 @@ flag debugnan
 
 library
 
-    Build-Depends:      base, hmatrix-base,
-                        array,
-                        storable-complex,
-                        process, random,
-                        vector >= 0.8,
-                        binary,
-                        deepseq
+    Build-Depends:      base, hmatrix-base, array, vector,
+                        process, random
+
 
     Extensions:         ForeignFunctionInterface,
                         CPP
 
     hs-source-dirs:     src
-    Exposed-modules:    Data.Packed,
-                        Data.Packed.Vector,
-                        Data.Packed.Matrix,
-                        Data.Packed.Foreign,
-                        Numeric.GSL.Differentiation,
+    Exposed-modules:    Numeric.GSL.Differentiation,
                         Numeric.GSL.Integration,
                         Numeric.GSL.Fourier,
                         Numeric.GSL.Polynomials,
@@ -101,24 +92,15 @@ library
                         Numeric.GSL,
                         Numeric.Container,
                         Numeric.LinearAlgebra,
-                        Numeric.LinearAlgebra.LAPACK,
                         Numeric.LinearAlgebra.Algorithms,
                         Numeric.LinearAlgebra.Util,
-                        Data.Packed.ST,
-                        Data.Packed.Development
                         Graphics.Plot,
                         Numeric.HMatrix,
                         Numeric.HMatrix.Data,
                         Numeric.HMatrix.Devel
-    other-modules:      Data.Packed.Internal,
-                        Data.Packed.Internal.Common,
-                        Data.Packed.Internal.Signatures,
-                        Data.Packed.Internal.Vector,
-                        Data.Packed.Internal.Matrix,
-                        Numeric.Random,
+    other-modules:      Numeric.Random,
                         Numeric.GSL.Internal,
                         Numeric.GSL.Vector,
-                        Numeric.Conversion,
                         Numeric.ContainerBoot,
                         Numeric.IO,
                         Numeric.Chain,
@@ -127,8 +109,7 @@ library
                         Numeric.LinearAlgebra.Util.Convolution
                         
 
-    C-sources:          src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
-                        src/Numeric/GSL/gsl-aux.c
+    C-sources:          src/Numeric/GSL/gsl-aux.c
 
     cc-options:         -O4 -msse2 -Wall
 
diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs
index ff649b7..3a8dd94 100644
--- a/packages/hmatrix/src/Numeric/Container.hs
+++ b/packages/hmatrix/src/Numeric/Container.hs
@@ -30,7 +30,7 @@ module Numeric.Container (
     -- * Basic functions
     module Data.Packed,
     konst, build,
-    constant, linspace,
+    linspace,
     diag, ident,
     ctrans,
     -- * Generic operations
@@ -64,8 +64,7 @@ module Numeric.Container (
     fscanfVector, fprintfVector, freadVector, fwriteVector,
 ) where
 
-import Data.Packed
-import Data.Packed.Internal(constantD)
+import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
 import Numeric.ContainerBoot
 import Numeric.Chain
 import Numeric.IO
@@ -75,15 +74,6 @@ import Numeric.Random
 
 ------------------------------------------------------------------
 
-{- | creates a vector with a given number of equal components:
-
-@> constant 2 7
-7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
--}
-constant :: Element a => a -> Int -> Vector a
--- constant x n = runSTVector (newVector x n)
-constant = constantD-- about 2x faster
-
 {- | Creates a real vector containing a range of values:
 
 >>> linspace 5 (-3,7::Double)
diff --git a/packages/hmatrix/src/Numeric/ContainerBoot.hs b/packages/hmatrix/src/Numeric/ContainerBoot.hs
index ea4262c..ef21763 100644
--- a/packages/hmatrix/src/Numeric/ContainerBoot.hs
+++ b/packages/hmatrix/src/Numeric/ContainerBoot.hs
@@ -42,7 +42,7 @@ module Numeric.ContainerBoot (
 import Data.Packed
 import Data.Packed.ST as ST
 import Numeric.Conversion
-import Data.Packed.Internal
+import Data.Packed.Development
 import Numeric.GSL.Vector
 import Data.Complex
 import Control.Applicative((<*>))
@@ -201,7 +201,7 @@ instance Container Vector Float where
     equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
     arctan2 = vectorZipF ATan2
     scalar x = fromList [x]
-    konst' = constantD
+    konst' = constant
     build' = buildV
     conj = id
     cmap = mapVector
@@ -229,7 +229,7 @@ instance Container Vector Double where
     equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
     arctan2 = vectorZipR ATan2
     scalar x = fromList [x]
-    konst' = constantD
+    konst' = constant
     build' = buildV
     conj = id
     cmap = mapVector
@@ -257,7 +257,7 @@ instance Container Vector (Complex Double) where
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
     arctan2 = vectorZipC ATan2
     scalar x = fromList [x]
-    konst' = constantD
+    konst' = constant
     build' = buildV
     conj = conjugateC
     cmap = mapVector
@@ -285,7 +285,7 @@ instance Container Vector (Complex Float) where
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
     arctan2 = vectorZipQ ATan2
     scalar x = fromList [x]
-    konst' = constantD
+    konst' = constant
     build' = buildV
     conj = conjugateQ
     cmap = mapVector
@@ -569,7 +569,7 @@ diag v = diagRect 0 v n n where n = dim v
 
 -- | creates the identity matrix of given dimension
 ident :: (Num a, Element a) => Int -> Matrix a
-ident n = diag (constantD 1 n)
+ident n = diag (constant 1 n)
 
 --------------------------------------------------------
 
diff --git a/packages/hmatrix/src/Numeric/GSL.hs b/packages/hmatrix/src/Numeric/GSL.hs
index 5f39a3e..61b8d7b 100644
--- a/packages/hmatrix/src/Numeric/GSL.hs
+++ b/packages/hmatrix/src/Numeric/GSL.hs
@@ -1,12 +1,11 @@
 {- |
 
 Module      :  Numeric.GSL
-Copyright   :  (c) Alberto Ruiz 2006-7
-License     :  GPL-style
+Copyright   :  (c) Alberto Ruiz 2006-14
+License     :  GPL
 
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses -fffi and -fglasgow-exts
 
 This module reexports all available GSL functions.
 
@@ -41,3 +40,4 @@ import Data.Complex
 -- | This action removes the GSL default error handler (which aborts the program), so that
 -- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
 foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
+
diff --git a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
index 93c5007..0fb58ef 100644
--- a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
@@ -1,13 +1,10 @@
-{-# OPTIONS  #-}
------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.GSL.Differentiation
 Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL-style
+License     :  GPL
 
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Numerical differentiation.
 
@@ -15,7 +12,8 @@ Numerical differentiation.
 
 From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\"
 -}
------------------------------------------------------------------------------
+
+
 module Numeric.GSL.Differentiation (
     derivCentral,
     derivForward,
@@ -26,7 +24,7 @@ import Foreign.C.Types
 import Foreign.Marshal.Alloc(malloc, free)
 import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
 import Foreign.Storable(peek)
-import Data.Packed.Internal(check,(//))
+import Numeric.GSL.Internal
 import System.IO.Unsafe(unsafePerformIO)
 
 derivGen ::
diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
index c4f3a91..0a92373 100644
--- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
@@ -3,9 +3,8 @@ Module      :  Numeric.GSL.Fitting
 Copyright   :  (c) Alberto Ruiz 2010
 License     :  GPL
 
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Nonlinear Least-Squares Fitting
 
@@ -42,7 +41,7 @@ expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]
 (1.0192487112786812,3.782067731353722e-2)]
 
 -}
------------------------------------------------------------------------------
+
 
 module Numeric.GSL.Fitting (
     -- * Levenberg-Marquardt
@@ -51,7 +50,6 @@ module Numeric.GSL.Fitting (
     fitModelScaled, fitModel
 ) where
 
-import Data.Packed.Internal
 import Numeric.LinearAlgebra
 import Numeric.GSL.Internal
 
@@ -61,6 +59,9 @@ import System.IO.Unsafe(unsafePerformIO)
 
 -------------------------------------------------------------------------
 
+type TVV = TV (TV Res)
+type TVM = TV (TM Res)
+
 data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom.
                    | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately.
         deriving (Enum,Eq,Show,Bounded)
diff --git a/packages/hmatrix/src/Numeric/GSL/Fourier.hs b/packages/hmatrix/src/Numeric/GSL/Fourier.hs
index 86aedd6..734325b 100644
--- a/packages/hmatrix/src/Numeric/GSL/Fourier.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Fourier.hs
@@ -1,26 +1,23 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
------------------------------------------------------------------------------
 {- |
 Module      : Numeric.GSL.Fourier
 Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+License     :  GPL
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Fourier Transform.
 
 <http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html#Fast-Fourier-Transforms>
 
 -}
------------------------------------------------------------------------------
+
 module Numeric.GSL.Fourier (
     fft,
     ifft
 ) where
 
-import Data.Packed.Internal
+import Data.Packed
+import Numeric.GSL.Internal
 import Data.Complex
 import Foreign.C.Types
 import System.IO.Unsafe (unsafePerformIO)
@@ -30,7 +27,7 @@ genfft code v = unsafePerformIO $ do
     app2 (c_fft code) vec v vec r "fft"
     return r
 
-foreign import ccall unsafe "gsl-aux.h fft" c_fft ::  CInt -> TCVCV
+foreign import ccall unsafe "gsl-aux.h fft" c_fft ::  CInt -> TCV (TCV Res)
 
 
 {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
diff --git a/packages/hmatrix/src/Numeric/GSL/Integration.hs b/packages/hmatrix/src/Numeric/GSL/Integration.hs
index 5f0a415..9c1d43a 100644
--- a/packages/hmatrix/src/Numeric/GSL/Integration.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Integration.hs
@@ -1,19 +1,15 @@
-{-# OPTIONS #-}
------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.GSL.Integration
 Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+License     :  GPL
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Numerical integration routines.
 
 <http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
 -}
------------------------------------------------------------------------------
+
 
 module Numeric.GSL.Integration (
     integrateQNG,
@@ -28,7 +24,7 @@ import Foreign.C.Types
 import Foreign.Marshal.Alloc(malloc, free)
 import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
 import Foreign.Storable(peek)
-import Data.Packed.Internal(check,(//))
+import Numeric.GSL.Internal
 import System.IO.Unsafe(unsafePerformIO)
 
 eps = 1e-12
diff --git a/packages/hmatrix/src/Numeric/GSL/Internal.hs b/packages/hmatrix/src/Numeric/GSL/Internal.hs
index 69a9750..a1c4e0c 100644
--- a/packages/hmatrix/src/Numeric/GSL/Internal.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Internal.hs
@@ -1,23 +1,43 @@
+-- |
 -- Module      :  Numeric.GSL.Internal
 -- Copyright   :  (c) Alberto Ruiz 2009
 -- License     :  GPL
---
--- Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  uses ffi
+--
 --
 -- Auxiliary functions.
 --
--- #hide
 
-module Numeric.GSL.Internal where
 
-import Data.Packed.Internal
+module Numeric.GSL.Internal(
+    iv,
+    mkVecfun,
+    mkVecVecfun,
+    mkDoubleVecVecfun,
+    mkDoublefun,
+    aux_vTov,
+    mkVecMatfun,
+    mkDoubleVecMatfun,
+    aux_vTom,
+    createV,
+    createMIO,
+    module Data.Packed.Development,
+    check,
+    Res,TV,TM,TCV,TCM
+) where
+
+import Data.Packed
+import Data.Packed.Development hiding (check)
+import Data.Complex
 
 import Foreign.Marshal.Array(copyArray)
 import Foreign.Ptr(Ptr, FunPtr)
 import Foreign.C.Types
+import Foreign.C.String(peekCString)
 import System.IO.Unsafe(unsafePerformIO)
+import Data.Vector.Storable(unsafeWith)
+import Control.Monad(when)
 
 iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
 iv f n p = f (createV (fromIntegral n) copy "iv") where
@@ -74,3 +94,33 @@ createMIO r c fun msg = do
     res <- createMatrix RowMajor r c
     app1 fun mat res msg
     return res
+
+--------------------------------------------------------------------------------
+
+-- | check the error code
+check :: String -> IO CInt -> IO ()
+check msg f = do
+    err <- f
+    when (err/=0) $ do
+        ps <- gsl_strerror err
+        s <- peekCString ps
+        error (msg++": "++s)
+    return ()
+
+-- | description of GSL error codes
+foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
+
+type PF = Ptr Float
+type PD = Ptr Double
+type PQ = Ptr (Complex Float)
+type PC = Ptr (Complex Double)
+
+type Res = IO CInt
+type TV x  = CInt -> PD -> x
+type TM x  = CInt -> CInt -> PD -> x
+type TCV x = CInt -> PC -> x
+type TCM x = CInt -> CInt -> PC -> x
+
+type TVV = TV (TV Res)
+type TVM = TV (TM Res)
+
diff --git a/packages/hmatrix/src/Numeric/GSL/Minimization.hs b/packages/hmatrix/src/Numeric/GSL/Minimization.hs
index 1879dab..056d463 100644
--- a/packages/hmatrix/src/Numeric/GSL/Minimization.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Minimization.hs
@@ -1,13 +1,9 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.GSL.Minimization
 Copyright   :  (c) Alberto Ruiz 2006-9
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+License     :  GPL
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Minimization of a multidimensional function using some of the algorithms described in:
 
@@ -48,7 +44,7 @@ The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimp
 
 -}
 
------------------------------------------------------------------------------
+
 module Numeric.GSL.Minimization (
     minimize, minimizeV, MinimizeMethod(..),
     minimizeD, minimizeVD, MinimizeMethodD(..),
@@ -60,8 +56,7 @@ module Numeric.GSL.Minimization (
 ) where
 
 
-import Data.Packed.Internal
-import Data.Packed.Matrix
+import Data.Packed
 import Numeric.GSL.Internal
 
 import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
@@ -112,7 +107,7 @@ uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
 
 
 foreign import ccall safe "uniMinimize"
-    c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM
+    c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM Res
 
 data MinimizeMethod = NMSimplex
                     | NMSimplex2
@@ -150,13 +145,13 @@ minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
                          "minimize"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
-        sol = cdat $ dropColumns 3 $ dropRows (it-1) path
+        sol = flatten $ dropColumns 3 $ dropRows (it-1) path
     freeHaskellFunPtr fp
     return (sol, path)
 
 
 foreign import ccall safe "gsl-aux.h minimize"
-    c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM
+    c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TV(TV(TM Res))
 
 ----------------------------------------------------------------------------------
 
@@ -207,7 +202,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
                          "minimizeD"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
-        sol = cdat $ dropColumns 2 $ dropRows (it-1) path
+        sol = flatten $ dropColumns 2 $ dropRows (it-1) path
     freeHaskellFunPtr fp
     freeHaskellFunPtr dfp
     return (sol,path)
@@ -215,9 +210,9 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
 foreign import ccall safe "gsl-aux.h minimizeD"
     c_minimizeD :: CInt
                 -> FunPtr (CInt -> Ptr Double -> Double)
-                -> FunPtr TVV
+                -> FunPtr (TV (TV Res))
                 -> Double -> Double -> Double -> CInt
-                -> TVM
+                -> TV (TM Res)
 
 ---------------------------------------------------------------------
 
diff --git a/packages/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs
index 9a29085..7549a65 100644
--- a/packages/hmatrix/src/Numeric/GSL/ODE.hs
+++ b/packages/hmatrix/src/Numeric/GSL/ODE.hs
@@ -2,10 +2,8 @@
 Module      :  Numeric.GSL.ODE
 Copyright   :  (c) Alberto Ruiz 2010
 License     :  GPL
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Solution of ordinary differential equation (ODE) initial value problems.
 
@@ -34,7 +32,7 @@ module Numeric.GSL.ODE (
     odeSolve, odeSolveV, ODEMethod(..), Jacobian
 ) where
 
-import Data.Packed.Internal
+import Data.Packed
 import Numeric.GSL.Internal
 
 import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
@@ -43,6 +41,10 @@ import System.IO.Unsafe(unsafePerformIO)
 
 -------------------------------------------------------------------------
 
+type TVV  = TV (TV Res)
+type TVM  = TV (TM Res)
+type TVVM = TV (TV (TM Res))
+
 type Jacobian = Double -> Vector Double -> Matrix Double
 
 -- | Stepping functions
diff --git a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
index 290c615..b1be85d 100644
--- a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
@@ -1,25 +1,23 @@
-{-# LANGUAGE CPP, ForeignFunctionInterface #-}
------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.GSL.Polynomials
 Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+License     :  GPL
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Polynomials.
 
 <http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>
 
 -}
------------------------------------------------------------------------------
+
+
 module Numeric.GSL.Polynomials (
     polySolve
 ) where
 
-import Data.Packed.Internal
+import Data.Packed
+import Numeric.GSL.Internal
 import Data.Complex
 import System.IO.Unsafe (unsafePerformIO)
 
@@ -55,4 +53,5 @@ polySolve' v | dim v > 1 = unsafePerformIO $ do
     return r
              | otherwise = error "polySolve on a polynomial of degree zero"
 
-foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TVCV
+foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TV (TCV Res)
+
diff --git a/packages/hmatrix/src/Numeric/GSL/Root.hs b/packages/hmatrix/src/Numeric/GSL/Root.hs
index 9d561c4..b9f3b94 100644
--- a/packages/hmatrix/src/Numeric/GSL/Root.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Root.hs
@@ -2,10 +2,8 @@
 Module      :  Numeric.GSL.Root
 Copyright   :  (c) Alberto Ruiz 2009
 License     :  GPL
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Multidimensional root finding.
 
@@ -41,14 +39,16 @@ module Numeric.GSL.Root (
     rootJ, RootMethodJ(..),
 ) where
 
-import Data.Packed.Internal
-import Data.Packed.Matrix
+import Data.Packed
 import Numeric.GSL.Internal
 import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
 import Foreign.C.Types
 import System.IO.Unsafe(unsafePerformIO)
 
 -------------------------------------------------------------------------
+type TVV = TV (TV Res)
+type TVM = TV (TM Res)
+
 
 data UniRootMethod = Bisection
                    | FalsePos
@@ -76,7 +76,7 @@ uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
     return (sol !! 1, path)
 
 foreign import ccall safe "root"
-    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM
+    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM Res
 
 -------------------------------------------------------------------------
 data UniRootMethodJ = UNewton
@@ -108,7 +108,7 @@ uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
 
 foreign import ccall safe "rootj"
     c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
-            -> Double -> CInt -> Double -> TM
+            -> Double -> CInt -> Double -> TM Res
 
 -------------------------------------------------------------------------
 
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
index 6204b8e..29c8bb7 100644
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs
@@ -2,11 +2,9 @@
 -- |
 -- Module      :  Numeric.GSL.Vector
 -- Copyright   :  (c) Alberto Ruiz 2007
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  portable (uses FFI)
 --
 -- Low level interface to vector operations.
 --
@@ -20,18 +18,20 @@ module Numeric.GSL.Vector (
     FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
     FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
     FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,
-    RandDist(..), randomVector
+    RandDist(..), randomVector,
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector
 ) where
 
-import Data.Packed.Internal.Common
-import Data.Packed.Internal.Signatures
-import Data.Packed.Internal.Vector
+import Data.Packed
+import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
 
 import Data.Complex
 import Foreign.Marshal.Alloc(free)
 import Foreign.Marshal.Array(newArray)
 import Foreign.Ptr(Ptr)
 import Foreign.C.Types
+import Foreign.C.String(newCString)
 import System.IO.Unsafe(unsafePerformIO)
 import Control.Monad(when)
 
@@ -186,7 +186,7 @@ foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV
 toScalarAux fun code v = unsafePerformIO $ do
     r <- createVector 1
     app2 (fun (fromei code)) vec v vec r "toScalarAux"
-    return (r `at` 0)
+    return (r @> 0)
 
 vectorMapAux fun code v = unsafePerformIO $ do
     r <- createVector (dim v)
@@ -326,3 +326,118 @@ randomVector seed dist n = unsafePerformIO $ do
     return r
 
 foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
+
+--------------------------------------------------------------------------------
+
+-- | Saves a matrix as 2D ASCII table.
+saveMatrix :: FilePath
+           -> String     -- ^ format (%f, %g, %e)
+           -> Matrix Double
+           -> IO ()
+saveMatrix filename fmt m = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    let o = if orderOf m == RowMajor then 1 else 0
+    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
+
+--------------------------------------------------------------------------------
+
+-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
+fscanfVector :: FilePath -> Int -> IO (Vector Double)
+fscanfVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
+fprintfVector :: FilePath -> String -> Vector Double -> IO ()
+fprintfVector filename fmt v = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
+
+-- | Loads a vector from a binary file (the number of elements must be known in advance).
+freadVector :: FilePath -> Int -> IO (Vector Double)
+freadVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector to a binary file.
+fwriteVector :: FilePath -> Vector Double -> IO ()
+fwriteVector filename v = do
+    charname <- newCString filename
+    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
+    free charname
+
+foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
+
+type PF = Ptr Float                             --
+type PD = Ptr Double                            --
+type PQ = Ptr (Complex Float)                   --
+type PC = Ptr (Complex Double)                  --
+type TF = CInt -> PF -> IO CInt                 --
+type TFF = CInt -> PF -> TF                     --
+type TFV = CInt -> PF -> TV                     --
+type TVF = CInt -> PD -> TF                     --
+type TFFF = CInt -> PF -> TFF                   --
+type TV = CInt -> PD -> IO CInt                 --
+type TVV = CInt -> PD -> TV                     --
+type TVVV = CInt -> PD -> TVV                   --
+type TFM = CInt -> CInt -> PF -> IO CInt        --
+type TFMFM =  CInt -> CInt -> PF -> TFM         --
+type TFMFMFM =  CInt -> CInt -> PF -> TFMFM     --
+type TM = CInt -> CInt -> PD -> IO CInt         --
+type TMM =  CInt -> CInt -> PD -> TM            --
+type TVMM = CInt -> PD -> TMM                   --
+type TMVMM = CInt -> CInt -> PD -> TVMM         --
+type TMMM =  CInt -> CInt -> PD -> TMM          --
+type TVM = CInt -> PD -> TM                     --
+type TVVM = CInt -> PD -> TVM                   --
+type TMV = CInt -> CInt -> PD -> TV             --
+type TMMV = CInt -> CInt -> PD -> TMV           --
+type TMVM = CInt -> CInt -> PD -> TVM           --
+type TMMVM = CInt -> CInt -> PD -> TMVM         --
+type TCM = CInt -> CInt -> PC -> IO CInt        --
+type TCVCM = CInt -> PC -> TCM                  --
+type TCMCVCM = CInt -> CInt -> PC -> TCVCM      --
+type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM   --
+type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM  --
+type TCMCM = CInt -> CInt -> PC -> TCM          --
+type TVCM = CInt -> PD -> TCM                   --
+type TCMVCM = CInt -> CInt -> PC -> TVCM        --
+type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM    --
+type TCMCMCM = CInt -> CInt -> PC -> TCMCM      --
+type TCV = CInt -> PC -> IO CInt                --
+type TCVCV = CInt -> PC -> TCV                  --
+type TCVCVCV = CInt -> PC -> TCVCV              --
+type TCVV = CInt -> PC -> TV                    --
+type TQV = CInt -> PQ -> IO CInt                --
+type TQVQV = CInt -> PQ -> TQV                  --
+type TQVQVQV = CInt -> PQ -> TQVQV              --
+type TQVF = CInt -> PQ -> TF                    --
+type TQM = CInt -> CInt -> PQ -> IO CInt        --
+type TQMQM = CInt -> CInt -> PQ -> TQM          --
+type TQMQMQM = CInt -> CInt -> PQ -> TQMQM      --
+type TCMCV = CInt -> CInt -> PC -> TCV          --
+type TVCV = CInt -> PD -> TCV                   --
+type TCVM = CInt -> PC -> TM                    --
+type TMCVM = CInt -> CInt -> PD -> TCVM         --
+type TMMCVM = CInt -> CInt -> PD -> TMCVM       --
+
diff --git a/packages/hmatrix/src/Numeric/IO.hs b/packages/hmatrix/src/Numeric/IO.hs
index 836f352..58fa2b4 100644
--- a/packages/hmatrix/src/Numeric/IO.hs
+++ b/packages/hmatrix/src/Numeric/IO.hs
@@ -20,11 +20,12 @@ module Numeric.IO (
 ) where
 
 import Data.Packed
-import Data.Packed.Internal
+import Data.Packed.Development
 import System.Process(readProcess)
 import Text.Printf(printf)
 import Data.List(intersperse)
 import Data.Complex
+import Numeric.GSL.Vector
 
 {- | Creates a string from a matrix given a separator and a function to show each entry. Using
 this function the user can easily define any desired display function:
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs
index 8c4b610..0a6eaa1 100644
--- a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -79,8 +79,8 @@ module Numeric.LinearAlgebra.Algorithms (
 ) where
 
 
-import Data.Packed.Internal hiding ((//))
-import Data.Packed.Matrix
+import Data.Packed.Development hiding ((//))
+import Data.Packed
 import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Data.List(foldl1')
 import Data.Array