update dependencies, move examples etc

author: Alberto Ruiz <aruiz@um.es> 2014-05-21 10:30:55 +0200
committer: Alberto Ruiz <aruiz@um.es> 2014-05-21 10:30:55 +0200
commit: 197e88c3b56d28840217010a2871c6ea3a4dd1a4 (patch)
tree: 825be9d6c9d87d23f7e5497c0133d11d52c63535 /packages/hmatrix
parent: e07c3dee7235496b71a89233106d93f6cc94ada1 (diff)
47 files changed, 0 insertions, 4803 deletions
diff --git a/packages/hmatrix/CHANGELOG b/packages/hmatrix/CHANGELOG
deleted file mode 100644
index 99a6845..0000000
--- a/packages/hmatrix/CHANGELOG
+++ /dev/null
@@ -1,219 +0,0 @@
-0.16.0.0
--------
-    * created hmatrix-base
-    * Added more organized reexport modules:
-      Numeric.HMatrix, Numeric.HMatrix.Data, Numeric.HMatrix.Devel
-      (The documentation is hidden for the other modules,
-       but they continue to be exposed and are not deprecated)
-    * added support for empty arrays
-    * join deprecated, use vjoin
-    * dot now conjugates the first input vector
-    * added udot (unconjugated dot product)
-    
-    * (<.>) overloaded to matrix and dot products
-    * added Monoid instance for Matrix using matrix product
-    
-    * improved build and konst
-    * improved linspace
-    
-    * improved error messages
-    * more usage examples
-    * simplified LSDiv
-    * Plot functions moved to Numeric.LinearAlgebra.Util
-    * removed (!) (use (¦)), added (——)
-0.15.2.0
--------
-    * general pinvTol and improved pinv
-0.15.1.0
--------
-    * One-dimensional minimization
-    * Doubly-adaptive quadrature for difficult integrands
-0.15.0.0
--------
-    * Data.Packed.Foreign (additional FFI helpers)
-    * NFData instance of Matrix
-    * Unidimensional root finding
-    * In Numeric.LinearAlgebra.Util:
-      pairwise2D, rowOuters, null1, null1sym, size, unitary, mt, (¦), (?), (¿)
-    * diagBlock
-    * meanCov moved to Container
-0.14.1.0
--------
-    * In Numeric.LinearAlgebra.Util:
-      convolution: corr, conv, corr2, conv2, separable, corrMin
-      kronecker: vec, vech, dup, vtrans
-0.14.0.0
--------
-    * integration over infinite intervals
-    * msadams and msbdf methods for ode
-    * Numeric.LinearAlgebra.Util
-    * (<\>) extended to multiple right-hand sides
-    * orth
-0.13.0.0
--------
-    * tests moved to new package hmatrix-tests
-0.11.2.0
--------
-    * geigSH' (symmetric generalized eigensystem)
-    * mapVectorWithIndex
-0.11.1.0
--------
-    * exported Mul
-    * mapMatrixWithIndex{,M,M_}
-0.11.0.0
--------
-    * flag -fvector default = True
-    * invlndet (inverse and log of determinant)
-    * step, cond
-    * find
-    * assoc, accum
-0.10.0.0
--------
-    * Module reorganization
-    * Support for Float and Complex Float elements (excluding LAPACK computations)
-    * Binary instances for Vector and Matrix
-    * optimiseMult
-    * mapVectorM, mapVectorWithIndexM, unzipVectorWith, and related functions.
-    * diagRect admits diagonal vectors of any length without producing an error,
-      and takes an additional argument for the off-diagonal elements.
-    * different signatures in some functions
-0.9.3.0
--------
-    * flag -fvector to optionally use Data.Vector.Storable.Vector
-  without any conversion.
-    * Simpler module structure.
-    * toBlocks, toBlocksEvery
-    * cholSolve, mbCholSH
-    * GSL Nonlinear Least-Squares fitting using Levenberg-Marquardt.
-    * GSL special functions moved to separate package hmatrix-special.
-    * Added offset of Vector, allowing fast, noncopy subVector (slice).
-  Vector is now identical to Roman Leshchinskiy's Data.Vector.Storable.Vector,
-  so we can convert from/to them in O(1).
-    * Removed Data.Packed.Convert, see examples/vector.hs
-0.8.3.0
--------
-    * odeSolve
-    * Matrix arithmetic automatically replicates matrix with single row/column
-    * latexFormat, dispcf
-0.8.2.0
--------
-    * fromRows/fromColumns now automatically expand vectors of dim 1
-      to match the common dimension.
-      fromBlocks also replicates single row/column matrices.
-      Previously all dimensions had to be exactly the same.
-    * display utilities: dispf, disps, vecdisp
-    * scalar
-    * minimizeV, minimizeVD, using Vector instead of lists.
-0.8.1.0
--------
-    * runBenchmarks
-0.8.0.0
--------
-    * singularValues, fullSVD, thinSVD, compactSVD, leftSV, rightSV
-      and complete interface to [d|z]gesdd.
-      Algorithms based on the SVD of large matrices can now be
-      significantly faster.
-    * eigenvalues, eigenvaluesSH
-    * linearSolveLS, rq
-0.7.2.0
--------
-    * ranksv
-0.7.1.0
--------
-    * buildVector/buildMatrix
-    * removed NFData instances
-0.6.0.0
--------
-    * added randomVector, gaussianSample, uniformSample, meanCov
-    * added rankSVD, nullspaceSVD
-    * rank, nullspacePrec, and economy svd defined in terms of ranksvd.
-    * economy svd now admits zero rank matrices and return a "degenerate
-      rank 1" decomposition with zero singular value.
-    * added NFData instances for Matrix and Vector.
-    * liftVector, liftVector2 replaced by mapVector, zipVector.
diff --git a/packages/hmatrix/INSTALL.md b/packages/hmatrix/INSTALL.md
deleted file mode 100644
index ef51167..0000000
--- a/packages/hmatrix/INSTALL.md
+++ /dev/null
@@ -1,117 +0,0 @@
-# [hmatrix][hmatrix2] installation
-This package requires the [Glasgow Haskell Compiler](http://www.haskell.org/ghc/index.html) ghc >= 6.10, and [cabal-install](http://www.haskell.org/haskellwiki/Cabal-Install), conveniently available in the [Haskell Platform](http://hackage.haskell.org/platform), and the development packages for [GSL](http://www.gnu.org/software/gsl) and BLAS/[LAPACK](http://www.netlib.org/lapack). (The graphical functions also require **gnuplot** and **imagemagick**.)
-[hmatrix]: http://code.haskell.org/hmatrix
-[hmatrix2]: http://perception.inf.um.es/hmatrix
-## Linux ##################################################
-Ubuntu/Debian:
-    $ sudo apt-get install libgsl0-dev liblapack-dev
-    $ cabal install hmatrix
-Arch Linux: If the automatic installation from Hackage fails, install atlas-lapack and gsl, unpack the source, change the build-type to Simple in hmatrix.cabal (line 28) and add extra-libraries: gsl lapack (line 194). 
-Other distributions may require additional libraries. They can be given in a **--configure-option**.
-## Mac OS/X ###############################################
-GSL must be installed via Homebrew or MacPorts.
-Via Homebrew:
-    $ brew install gsl
-    $ cabal install hmatrix
-Via MacPorts:
-    $ sudo port install gsl +universal
-    $ cabal install hmatrix
-(Contributed by Heinrich Apfelmus, Torsten Kemps-Benedix and Ted Fujimoto).
-## Windows ###############################################
-We use this [GSL binary](http://www.miscdebris.net/blog/2009/04/20/mingw-345-binaries-of-gnu-scientific-library-112-for-use-with-mingw-and-visual-c/), and blas/lapack dlls built with g77 (contributed by Gilberto Camara). All required files are in [gsl-lapack-windows.zip][winpack].
-(Due to [issue 21](https://github.com/albertoruiz/hmatrix/issues/21) we need hmatrix-0.13.1.0.)
-1) Install the Haskell Platform (tested on 2011.2.0.1)
-    > cabal update
-2) Download and unzip the following file into a stable folder %GSL%
-    http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip
-3.a) In a msys shell the installation should be fully automatic:
-    $ cabal install hmatrix-0.13.1.0 --extra-lib-dir=${GSL} --extra-include-dir=${GSL}
-3.b) Alternatively, in a normal windows cmd:
-     > cabal unpack hmatrix-0.13.1.0
-   Edit hmatrix.cabal, in line 28 change build-type to "Simple", and then
-     > cabal install --extra-lib-dir=%GSL% --extra-include-dir=%GSL%
-   It may be necessary to put the dlls in the search path.
-NOTE: The examples using graphics do not yet work in windows.
-[install]: http://code.haskell.org/hmatrix/INSTALL
-[install2]: http://patch-tag.com/r/aruiz/hmatrix/snapshot/current/content/pretty/INSTALL
-[winpack2]: http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip
-[winpack]: https://github.com/downloads/AlbertoRuiz/hmatrix/gsl-lapack-windows.zip
-## Tests ###############################################
-After installation we must verify that the library works as expected:
-    $ cabal install hmatrix-tests --enable-tests
-    $ ghci
-    > Numeric.LinearAlgebra.Tests.runTests 20
-    OK, passed 100 tests. 
-    OK, passed 100 tests.
-    ... etc...
-If you get any failure please run lapack's own tests to confirm that your version is not broken. For instance, in ubuntu 9.04, **libatlas-sse2** does not work (see this [bug report](https://bugs.launchpad.net/ubuntu/+source/atlas/+bug/368478)). If your lapack library is ok but hmatrix's tests fail please send a bug report!
-## Optimized BLAS/LAPACK ##########################################
-I have successfully tested ATLAS and MKL on Linux.
-### [ATLAS](http://math-atlas.sourceforge.net/)  ####################
-In Ubuntu >= 9.04 we need:
-    $ sudo apt-get install libatlas-base-dev
-In older Ubuntu/Debian versions we needed:
-    $ sudo apt-get install refblas3-dev lapack3-dev atlas3-base-dev
-We may use a version (sse2, 3dnow, etc.) optimized for the machine.
-### Intel's MKL  ###############################################
-There is a free noncommercial download available from Intel's website. To use it I have added the following lines in my .bashrc configuration file:
-    export LD_LIBRARY_PATH=/path/to/mkl/lib/arch
-    export LIBRARY_PATH=/path/to/mkl/lib/arch
-where arch = 32 or em64t.
-The library must be installed with the -fmkl flag:
-    $ cabal install hmatrix -fmkl
diff --git a/packages/hmatrix/LICENSE b/packages/hmatrix/LICENSE
deleted file mode 100644
index 3f67c2a..0000000
--- a/packages/hmatrix/LICENSE
+++ /dev/null
@@ -1,2 +0,0 @@
-Copyright Alberto Ruiz 2006-2007
-GPL license
diff --git a/packages/hmatrix/Setup.lhs b/packages/hmatrix/Setup.lhs
deleted file mode 100644
index 4b19c19..0000000
--- a/packages/hmatrix/Setup.lhs
+++ /dev/null
@@ -1,5 +0,0 @@
-#! /usr/bin/env runhaskell
-> import Distribution.Simple
-> main = defaultMain
diff --git a/packages/hmatrix/THANKS.md b/packages/hmatrix/THANKS.md
deleted file mode 100644
index b1417a6..0000000
--- a/packages/hmatrix/THANKS.md
+++ /dev/null
@@ -1,157 +0,0 @@
-I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus,
-and all the people in the Haskell mailing lists for their help.
-I am particularly grateful to Vivian McPhail for his excellent
-contributions: improved configure.hs, Binary instances for
-Vector and Matrix, support for Float and Complex Float elements,
-module reorganization, monadic mapVectorM, and many other improvements.
- Nico Mahlo discovered a bug in the eigendecomposition wrapper.
- Frederik Eaton discovered a bug in the design of the wrappers.
- Eric Kidd has created a wiki page explaining the installation on MacOS X:
-  http://www.haskell.org/haskellwiki/GSLHaskell_on_MacOS_X
- Fawzi Mohamed discovered a portability bug in the lapack wrappers.
- Pedro E. López de Teruel fixed the interface to lapack.
- Antti Siira discovered a bug in the plotting functions.
- Paulo Tanimoto helped to fix the configuration of the required libraries.
-  He also discovered the segfault of minimize.hs in ghci.
- Xiao-Yong Jin reported a bug on x86_64 caused by the assumptions in f2c.h,
-  which are wrong for this architecture.
- Jason Schroeder reported an error in the documentation.
- Bulat Ziganshin gave invaluable help for the ST monad interface to
-  in-place modifications.
- Don Stewart fixed the implementation of the internal data structures
-  to achieve excellent, C-like performance in Haskell functions which
-  explicitly work with the elements of vectors and matrices.
- Dylan Alex Simon improved the numeric instances to allow optimized
-  implementations of signum and abs on Vectors.
- Pedro E. López de Teruel discovered the need of asm("finit") to
-  avoid the wrong NaNs produced by foreign functions.
- Reiner Pope added support for luSolve, based on (d|z)getrs.
-  Made Matrix a product type and added changes to improve the code generated
-  by hmatrix-syntax.
- Simon Beaumont reported the need of QuickCheck<2 and the invalid
-  asm("finit") on ppc. He also contributed the configuration options
-  for the accelerate framework on OS X.
- Daniel Schüssler added compatibility with QuickCheck 2 as well
-  as QuickCheck 1 using the C preprocessor. He also added some
-  implementations for the new "shrink" method of class Arbitrary.
- Tracy Wadleigh improved the definitions of (|>) and (><), which now
-  apply an appropriate 'take' to the given lists so that they may be
-  safely used on lists that are too long (or infinite).
- Chris Waterson improved the configure.hs program for OS/X.
- Erik de Castro Lopo added buildVector and buildMatrix, which take a
-  size parameter(s) and a function that maps vector/matrix indices
-  to the values at that position.
- Jean-Francois Tremblay discovered an error in the tutorial.
- Gilberto Camara contributed improved blas and lapack dlls for Windows.
- Heinrich Apfelmus fixed hmatrix.cabal for OS/X. He also tested the package
-  on PPC discovering a problem in zgesdd.
- Felipe Lessa tested the performance of GSL special function bindings
-  and contributed the cabal flag "safe-cheap".
- Ozgur Akgun suggested better symbols for the Bound constructors in the
-  Linear Programming package.
- Tim Sears reported the zgesdd problem also in intel mac.
- Max Suica simplified the installation on Windows and improved the instructions.
- John Billings first reported an incompatibility with QuickCheck>=2.1.1
- Alexey Khudyakov cleaned up PRAGMAS and fixed some hlint suggestions.
- Torsten Kemps-Benedix reported an installation problem in OS/X.
- Stefan Kersten fixed hmatrix.cabal for 64-bit ghc-7 in OS/X
- Sacha Sokoloski reported an installation problem on Arch Linux and
-  helped with the configuration.
- Carter Schonwald helped with the configuration for Homebrew OS X and
-  found a tolerance problem in test "1E5 rots". He also discovered
-  a bug in the signature of cmap.
- Duncan Coutts reported a problem with configure.hs and contributed
-  a solution and a simplified Setup.lhs.
- Mark Wright fixed the import of vector >= 0.8.
- Bas van Dijk fixed the import of vector >= 0.8, got rid of some
-  deprecation warnings, used more explicit imports, and updated to ghc-7.4.
- Tom Nielsen discovered a problem in Config.hs, exposed by link problems
-  in Ubuntu 11.10 beta.
- Daniel Fischer reported some Haddock markup errors.
- Danny Chan added support for integration over infinite intervals, and fixed
-  Configure.hs using platform independent functions.
- Clark Gaebel removed superfluous thread safety.
- Jeffrey Burdges reported a glpk link problem on OS/X
- Jian Zhang reported the Windows installation problem due to new ODE interface.
- Mihaly Barasz and Ben Gamari fixed mapMatrix* and mapMatrixWithIndex
- Takano Akio fixed off-by-one errors in gsl-aux.c producing segfaults.
- Alex Lang implemented uniRoot and uniRootJ for one-dimensional root-finding.
- Mike Ledger contributed alternative FFI helpers for matrix interoperation with C
- Stephen J. Barr suggested flipping argument order in the double integral example
- Greg Horn fixed the bus error on ghci 64-bit.
- Kristof Bastiaensen added bindings for one-dimensional minimization.
- Matthew Peddie added bindings for gsl_integrate_cquad doubly-adaptive quadrature
-  for difficult integrands.
- Ben Gamari exposed matrixFromVector for Development.
- greg94301 reported tolerance issues in the tests.
- Clemens Lang updated the MacPort installation instructions.
- Henning Thielemann reported the pinv inefficient implementation.
- bdoering reported the problem of zero absolute tolerance in the integration functions.
- Alexei Uimanov replaced fromList by Vector.fromList.
- Adam Vogt updated the code for ghc-7.7
- Mike Meyer (mwm) added freeBSD library configuration information.
- tfgit updated the OSX installation instructions via Homebrew
- "yokto" and "ttylec" reported the problem with the dot product of complex vectors.
- Samium Gromoff reported a build failure caused by a size_t - int mismatch.
diff --git a/packages/hmatrix/examples/bool.hs b/packages/hmatrix/examples/bool.hs
deleted file mode 100644
index 679b8bf..0000000
--- a/packages/hmatrix/examples/bool.hs
+++ /dev/null
@@ -1,54 +0,0 @@
-- vectorized boolean operations defined in terms of step or cond
-import Numeric.LinearAlgebra
-infix  4  .==., ./=., .<., .<=., .>=., .>.
-infixr 3  .&&.
-infixr 2  .||.
-a .<.  b = step (b-a)
-a .<=. b = cond a b 1 1 0
-a .==. b = cond a b 0 1 0
-a ./=. b = cond a b 1 0 1
-a .>=. b = cond a b 0 1 1
-a .>.  b = step (a-b)
-a .&&. b  = step (a*b)
-a .||. b  = step (a+b)
-no a      = 1-a
-xor a b   = a ./=. b
-equiv a b = a .==. b
-imp a b   = no a .||. b
-taut x = minElement x == 1
-minEvery a b = cond a b a a b
-maxEvery a b = cond a b b b a
-- examples
-clip a b x = cond y b y y b where y = cond x a a x x
-disp = putStr . dispf 3
-eye n = ident n :: Matrix Double
-row = asRow . fromList    :: [Double] -> Matrix Double
-col = asColumn . fromList :: [Double] -> Matrix Double
-m = (3><4) [1..] :: Matrix Double
-p = row [0,0,1,1]
-q = row [0,1,0,1]
-main = do
-    print $ find (>6) m
-    disp $ assoc (6,8) 7 $ zip (find (/=0) (eye 5)) [10..]
-    disp $ accum (eye 5) (+) [((0,2),3), ((3,1),7), ((1,1),1)]
-    disp $ m .>=. 10  .||.  m .<. 4
-    (disp . fromColumns . map flatten) [p, q, p.&&.q, p .||.q, p `xor` q, p `equiv` q, p `imp` q]
-    print $ taut $ (p `imp` q ) `equiv` (no q `imp` no p)
-    print $ taut $ (xor p q) `equiv` (p .&&. no q .||. no p .&&. q)
-    disp $ clip 3 8 m
-    disp $ col [1..7] .<=. row [1..5]
-    disp $ cond (col [1..3]) (row [1..4]) m 50 (3*m)
diff --git a/packages/hmatrix/examples/data.txt b/packages/hmatrix/examples/data.txt
deleted file mode 100644
index 2b9bfea..0000000
--- a/packages/hmatrix/examples/data.txt
+++ /dev/null
@@ -1,10 +0,0 @@
- 0.9    1.1
- 2.1    3.9
- 3.1    9.2
- 4.0   51.8
- 4.9   25.3
- 6.1   35.7
- 7.0   49.4
- 7.9    3.6
- 9.1   81.5
-10.2   99.5
-\ No newline at end of file
diff --git a/packages/hmatrix/examples/deriv.hs b/packages/hmatrix/examples/deriv.hs
deleted file mode 100644
index c9456d1..0000000
--- a/packages/hmatrix/examples/deriv.hs
+++ /dev/null
@@ -1,8 +0,0 @@
-- Numerical differentiation
-import Numeric.GSL
-d :: (Double -> Double) -> (Double -> Double)
-d f x = fst $ derivCentral 0.01 f x
-main = print $ d (\x-> x * d (\y-> x+y) 1) 1
diff --git a/packages/hmatrix/examples/devel/ej1/functions.c b/packages/hmatrix/examples/devel/ej1/functions.c
deleted file mode 100644
index 02a4cdd..0000000
--- a/packages/hmatrix/examples/devel/ej1/functions.c
+++ /dev/null
@@ -1,35 +0,0 @@
-/* assuming row order */
-typedef struct { double r, i; } doublecomplex;
-#define DVEC(A) int A##n, double*A##p
-#define CVEC(A) int A##n, doublecomplex*A##p
-#define DMAT(A) int A##r, int A##c, double*A##p
-#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
-#define AT(M,row,col) (M##p[(row)*M##c + (col)])
-/*-----------------------------------------------------*/
-int c_scale_vector(double s, DVEC(x), DVEC(y)) {
-    int k;
-    for (k=0; k<=yn; k++) {
-        yp[k] = s*xp[k];
-    }
-    return 0;
-}
-/*-----------------------------------------------------*/
-int c_diag(DMAT(m),DVEC(y),DMAT(z)) {
-    int i,j;
-    for (j=0; j<yn; j++) {
-        yp[j] = AT(m,j,j);
-    }
-    for (i=0; i<mr; i++) {
-        for (j=0; j<mc; j++) {
-            AT(z,i,j) = i==j?yp[i]:0;
-        }
-    }
-    return 0;
-}
diff --git a/packages/hmatrix/examples/devel/ej1/wrappers.hs b/packages/hmatrix/examples/devel/ej1/wrappers.hs
deleted file mode 100644
index a88f74b..0000000
--- a/packages/hmatrix/examples/devel/ej1/wrappers.hs
+++ /dev/null
@@ -1,44 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-- $ ghc -O2 --make wrappers.hs functions.c
-import Numeric.LinearAlgebra
-import Data.Packed.Development
-import Foreign(Ptr,unsafePerformIO)
-import Foreign.C.Types(CInt)
-----------------------------------------------------
-main = do
-    print $ myScale 3.0 (fromList [1..10])
-    print $ myDiag $ (3><5) [1..]
-----------------------------------------------------
-foreign import ccall unsafe "c_scale_vector"
-    cScaleVector :: Double                -- scale
-                 -> CInt -> Ptr Double    -- argument
-                 -> CInt -> Ptr Double    -- result
-                 -> IO CInt               -- exit code
-myScale s x = unsafePerformIO $ do
-    y <- createVector (dim x)
-    app2 (cScaleVector s) vec x vec y "cScaleVector"
-    return y
-----------------------------------------------------
-- forcing row order
-foreign import ccall unsafe "c_diag"
-    cDiag :: CInt -> CInt -> Ptr Double  -- argument
-          -> CInt -> Ptr Double          -- result1
-          -> CInt -> CInt -> Ptr Double  -- result2
-          -> IO CInt                     -- exit code
-myDiag m = unsafePerformIO $ do
-    y <- createVector (min r c)
-    z <- createMatrix RowMajor r c
-    app3 cDiag mat (cmat m) vec y mat z "cDiag"
-    return (y,z)
-  where r = rows m
-        c = cols m
diff --git a/packages/hmatrix/examples/devel/ej2/functions.c b/packages/hmatrix/examples/devel/ej2/functions.c
deleted file mode 100644
index 4dcd377..0000000
--- a/packages/hmatrix/examples/devel/ej2/functions.c
+++ /dev/null
@@ -1,24 +0,0 @@
-/*  general element order   */
-typedef struct { double r, i; } doublecomplex;
-#define DVEC(A) int A##n, double*A##p
-#define CVEC(A) int A##n, doublecomplex*A##p
-#define DMAT(A) int A##r, int A##c, double*A##p
-#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
-#define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
-int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) {
-    int i,j,sr,sc;
-    if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
-    for (j=0; j<yn; j++) {
-        yp[j] = AT(m,j,j);
-    }
-    for (i=0; i<mr; i++) {
-        for (j=0; j<mc; j++) {
-            AT(z,i,j) = i==j?yp[i]:0;
-        }
-    }
-    return 0;
-}
diff --git a/packages/hmatrix/examples/devel/ej2/wrappers.hs b/packages/hmatrix/examples/devel/ej2/wrappers.hs
deleted file mode 100644
index 1c02a24..0000000
--- a/packages/hmatrix/examples/devel/ej2/wrappers.hs
+++ /dev/null
@@ -1,32 +0,0 @@
-{-# LANGUAGE ForeignFunctionInterface #-}
-- $ ghc -O2 --make wrappers.hs functions.c
-import Numeric.LinearAlgebra
-import Data.Packed.Development
-import Foreign(Ptr,unsafePerformIO)
-import Foreign.C.Types(CInt)
-----------------------------------------------------
-main = do
-    print $ myDiag $ (3><5) [1..]
-----------------------------------------------------
-- arbitrary data order
-foreign import ccall unsafe "c_diag"
-    cDiag :: CInt                        -- matrix order
-          -> CInt -> CInt -> Ptr Double  -- argument
-          -> CInt -> Ptr Double          -- result1
-          -> CInt -> CInt -> Ptr Double  -- result2
-          -> IO CInt                     -- exit code
-myDiag m = unsafePerformIO $ do
-    y <- createVector (min r c)
-    z <- createMatrix (orderOf m) r c
-    app3 (cDiag o) mat m vec y mat z "cDiag"
-    return (y,z)
-  where r = rows m
-        c = cols m
-        o = if orderOf m == RowMajor then 1 else 0
diff --git a/packages/hmatrix/examples/error.hs b/packages/hmatrix/examples/error.hs
deleted file mode 100644
index 5efae7c..0000000
--- a/packages/hmatrix/examples/error.hs
+++ /dev/null
@@ -1,21 +0,0 @@
-import Numeric.GSL
-import Numeric.GSL.Special
-import Numeric.LinearAlgebra
-import Prelude hiding (catch)
-import Control.Exception
-test x = catch
-       (print x)
-       (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]")
-main = do
-    setErrorHandlerOff
-    test $ log_e (-1)
-    test $ 5 + (fst.exp_e) 1000
-    test $ bessel_zero_Jnu_e (-0.3) 2
-    test $ (linearSolve 0 4 :: Matrix Double)
-    test $ (linearSolve 5 (sqrt (-1)) :: Matrix Double)
-    putStrLn "Bye"
-\ No newline at end of file
diff --git a/packages/hmatrix/examples/fitting.hs b/packages/hmatrix/examples/fitting.hs
deleted file mode 100644
index a8f6b1c..0000000
--- a/packages/hmatrix/examples/fitting.hs
+++ /dev/null
@@ -1,24 +0,0 @@
-- nonlinear least-squares fitting
-import Numeric.GSL.Fitting
-import Numeric.LinearAlgebra
-xs = map return [0 .. 39]
-sigma = 0.1
-ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs)
-            + scalar sigma * (randomVector 0 Gaussian 40)
-dat :: [([Double],([Double],Double))]
-dat = zip xs (zip ys (repeat sigma))
-expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
-expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
-(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
-main = do
-    print dat
-    print path
-    print sol
diff --git a/packages/hmatrix/examples/inplace.hs b/packages/hmatrix/examples/inplace.hs
deleted file mode 100644
index dcfff56..0000000
--- a/packages/hmatrix/examples/inplace.hs
+++ /dev/null
@@ -1,152 +0,0 @@
-- some tests of the interface for pure
-- computations with inplace updates
-import Numeric.LinearAlgebra
-import Data.Packed.ST
-import Data.Packed.Convert
-import Data.Array.Unboxed
-import Data.Array.ST
-import Control.Monad.ST
-import Control.Monad
-main = sequence_[
-    print test1,
-    print test2,
-    print test3,
-    print test4,
-    test5,
-    test6,
-    print test7,
-    test8,
-    test0]
-- helper functions
-vector l = fromList l :: Vector Double
-norm v = pnorm PNorm2 v
-- hmatrix vector and matrix
-v = vector [1..10]
-m = (5><10) [1..50::Double]
----------------------------------------------------------------------
-- vector creation by in-place updates on a copy of the argument
-test1 = fun v
-fun :: Element t => Vector t -> Vector t
-fun x = runSTVector $ do
-    a <- thawVector x
-    mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]
-    return a
-- another example: creation of an antidiagonal matrix from a list
-test2 = antiDiag 5 8 [1..] :: Matrix Double
-antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
-antiDiag r c l = runSTMatrix $ do
-    m <- newMatrix 0 r c
-    let d = min r c - 1
-    sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l
-    return m
-- using vector or matrix functions on mutable objects requires freezing:
-test3 = g1 v
-g1 x = runST $ do
-    a <- thawVector x
-    writeVector a (dim x -1) 0
-    b <- freezeVector a
-    return (norm b)
--  another possibility:
-test4 = g2 v
-g2 x = runST $ do
-    a <- thawVector x
-    writeVector a (dim x -1) 0
-    t <- liftSTVector norm a
-    return t
--------------------------------------------------------------
-- haskell arrays
-hv = listArray (0,9) [1..10::Double]
-hm = listArray ((0,0),(4,9)) [1..50::Double]
-- conversion from standard Haskell arrays
-test5 = do
-    print $ norm (vectorFromArray hv)
-    print $ norm v
-    print $ rcond (matrixFromArray hm)
-    print $ rcond m
-- conversion to mutable ST arrays
-test6 = do
-    let y = clearColumn m 1
-    print y
-    print (matrixFromArray y)
-clearColumn x c = runSTUArray $ do
-    a <- mArrayFromMatrix x
-    forM_ [0..rows x-1] $ \i->
-        writeArray a (i,c) (0::Double)
-    return a
-- hmatrix functions applied to mutable ST arrays
-test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)
-unitary v = runSTUArray $ do
-    a <- thaw v
-    n <- norm `fmap` vectorFromMArray a
-    b <- mapArray (/n) a
-    return b
-------------------------------------------------
-- (just to check that they are not affected)
-test0 = do
-    print v
-    print m
-    --print hv
-    --print hm
-------------------------------------------------
-histogram n ds = runSTVector $ do
-    h <- newVector (0::Double) n -- number of bins
-    let inc k = modifyVector h k (+1)
-    mapM_ inc ds
-    return h
-- check that newVector is really called with a fresh new array
-histoCheck ds = runSTVector $ do
-    h <- newVector (0::Double) 15 -- > constant for this test
-    let inc k = modifyVector h k (+1)
-    mapM_ inc ds
-    return h
-hc = fromList [1 .. 15::Double]
-- check that thawVector creates a new array
-histoCheck2 ds = runSTVector $ do
-    h <- thawVector hc
-    let inc k = modifyVector h k (+1)
-    mapM_ inc ds
-    return h
-test8 = do
-    let ds = [0..14]
-    print $ histogram 15 ds
-    print $ histogram 15 ds
-    print $ histogram 15 ds
-    print $ histoCheck ds
-    print $ histoCheck ds
-    print $ histoCheck ds
-    print $ histoCheck2 ds
-    print $ histoCheck2 ds
-    print $ histoCheck2 ds
-    putStrLn "----------------------"
diff --git a/packages/hmatrix/examples/integrate.hs b/packages/hmatrix/examples/integrate.hs
deleted file mode 100644
index 3a03da6..0000000
--- a/packages/hmatrix/examples/integrate.hs
+++ /dev/null
@@ -1,24 +0,0 @@
-- Numerical integration
-import Numeric.GSL
-quad f a b = fst $ integrateQAGS 1E-9 100 f a b  
-- A multiple integral can be easily defined using partial application
-quad2 f y1 y2 g1 g2 = quad h y1 y2
-  where
-    h y = quad (flip f y) (g1 y) (g2 y)
-volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
-                        0 r (const 0) (\x->sqrt (r*r-x*x))
-- wikipedia example
-exw = quad2 f 7 10 (const 11) (const 14)
-  where
-    f x y = x**2 + 4*y
-main = do
-    print $ quad (\x -> 4/(x^2+1)) 0 1
-    print pi
-    print $ volSphere 2.5
-    print $ 4/3*pi*2.5**3
-    print $ exw
diff --git a/packages/hmatrix/examples/kalman.hs b/packages/hmatrix/examples/kalman.hs
deleted file mode 100644
index 7fbe3d2..0000000
--- a/packages/hmatrix/examples/kalman.hs
+++ /dev/null
@@ -1,51 +0,0 @@
-import Numeric.LinearAlgebra
-import Graphics.Plot
-vector l = fromList l :: Vector Double
-matrix ls = fromLists ls :: Matrix Double
-diagl = diag . vector
-f = matrix [[1,0,0,0],
-            [1,1,0,0],
-            [0,0,1,0],
-            [0,0,0,1]]
-h = matrix [[0,-1,1,0],
-            [0,-1,0,1]]
-q = diagl [1,1,0,0]
-r = diagl [2,2]
-s0 = State (vector [0, 0, 10, -10]) (diagl [10,0, 100, 100])
-data System = System {kF, kH, kQ, kR :: Matrix Double}
-data State = State {sX :: Vector Double , sP :: Matrix Double}
-type Measurement = Vector Double
-kalman :: System -> State -> Measurement -> State
-kalman (System f h q r) (State x p) z = State x' p' where
-    px = f <> x                            -- prediction
-    pq = f <> p <> trans f + q             -- its covariance
-    y  = z - h <> px                       -- residue
-    cy = h <> pq <> trans h + r            -- its covariance
-    k  = pq <> trans h <> inv cy           -- kalman gain
-    x' = px + k <> y                       -- new state
-    p' = (ident (dim x) - k <> h) <> pq    -- its covariance
-sys = System f h q r
-zs = [vector [15-k,-20-k] | k <- [0..]]
-xs = s0 : zipWith (kalman sys) xs zs
-des = map (sqrt.takeDiag.sP) xs
-evolution n (xs,des) =
-    vector [1.. fromIntegral n]:(toColumns $ fromRows $ take n (zipWith (-) (map sX xs) des)) ++
-    (toColumns $ fromRows $ take n (zipWith (+) (map sX xs) des))
-main = do
-    print $ fromRows $ take 10 (map sX xs)
-    mapM_ (print . sP) $ take 10 xs
-    mplot (evolution 20 (xs,des))
diff --git a/packages/hmatrix/examples/lie.hs b/packages/hmatrix/examples/lie.hs
deleted file mode 100644
index db21ea8..0000000
--- a/packages/hmatrix/examples/lie.hs
+++ /dev/null
@@ -1,65 +0,0 @@
-- The magic of Lie Algebra
-import Numeric.LinearAlgebra
-disp = putStrLn . dispf 5
-rot1 :: Double -> Matrix Double
-rot1 a = (3><3)
- [ 1, 0, 0
- , 0, c, s
- , 0,-s, c ]
-    where c = cos a
-          s = sin a
-g1,g2,g3 :: Matrix Double
-g1 = (3><3) [0, 0,0
-            ,0, 0,1
-            ,0,-1,0]
-rot2 :: Double -> Matrix Double
-rot2 a = (3><3)
- [ c, 0, s
- , 0, 1, 0
- ,-s, 0, c ]
-    where c = cos a
-          s = sin a
-g2 = (3><3) [ 0,0,1
-            , 0,0,0
-            ,-1,0,0]
-rot3 :: Double -> Matrix Double
-rot3 a = (3><3)
- [ c, s, 0
- ,-s, c, 0
- , 0, 0, 1 ]
-    where c = cos a
-          s = sin a
-g3 = (3><3) [ 0,1,0
-            ,-1,0,0
-            , 0,0,0]
-deg=pi/180
-- commutator
-infix 8 &
-a & b = a <> b - b <> a
-infixl 6 |+|
-a |+| b = a + b + a&b /2  + (a-b)&(a & b) /12
-main = do
-   let a =  45*deg
-       b =  50*deg
-       c = -30*deg
-       exact = rot3 a <> rot1 b <> rot2 c
-       lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2
-   putStrLn "position in the tangent space:"
-   disp lie
-   putStrLn "exponential map back to the group (2 terms):"
-   disp (expm lie)
-   putStrLn "exact position:"
-   disp exact
diff --git a/packages/hmatrix/examples/minimize.hs b/packages/hmatrix/examples/minimize.hs
deleted file mode 100644
index 19b2cb3..0000000
--- a/packages/hmatrix/examples/minimize.hs
+++ /dev/null
@@ -1,50 +0,0 @@
-- the multidimensional minimization example in the GSL manual
-import Numeric.GSL
-import Numeric.LinearAlgebra
-import Graphics.Plot
-import Text.Printf(printf)
-- the function to be minimized
-f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
-- exact gradient
-df [x,y] = [20*(x-1), 40*(y-2)]
-- a minimization algorithm which does not require the gradient
-minimizeS f xi = minimize NMSimplex2 1E-2 100 (replicate (length xi) 1) f xi
-- Numerical estimation of the gradient
-gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]]
-partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where
-    g x = f (concat [a,x:b])
-    (a,_:b) = splitAt n v
-disp = putStrLn . format "  " (printf "%.3f")
-allMethods :: (Enum a, Bounded a) => [a]
-allMethods = [minBound .. maxBound]
-test method = do
-    print method
-    let (s,p) = minimize method 1E-2 30 [1,1] f [5,7]
-    print s
-    disp p
-testD method = do
-    print method
-    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7]
-    print s
-    disp p
-testD' method = do
-    putStrLn $ show method ++ " with estimated gradient"
-    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7]
-    print s
-    disp p
-main = do
-    mapM_ test [NMSimplex, NMSimplex2]
-    mapM_ testD allMethods
-    testD' ConjugateFR
-    mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7]
diff --git a/packages/hmatrix/examples/monadic.hs b/packages/hmatrix/examples/monadic.hs
deleted file mode 100644
index 7c6e0f4..0000000
--- a/packages/hmatrix/examples/monadic.hs
+++ /dev/null
@@ -1,118 +0,0 @@
-- monadic computations
-- (contributed by Vivian McPhail)
-import Numeric.LinearAlgebra
-import Control.Monad.State.Strict
-import Control.Monad.Maybe
-import Foreign.Storable(Storable)
-import System.Random(randomIO)
-------------------------------------------
-- an instance of MonadIO, a monad transformer
-type VectorMonadT = StateT Int IO
-test1 :: Vector Int -> IO (Vector Int)
-test1 = mapVectorM $ \x -> do
-    putStr $ (show x) ++ " "
-    return (x + 1)
-- we can have an arbitrary monad AND do IO
-addInitialM :: Vector Int -> VectorMonadT ()
-addInitialM = mapVectorM_ $ \x -> do
-    i <- get
-    liftIO $ putStr $ (show $ x + i) ++ " "
-    put $ x + i
-- sum the values of the even indiced elements
-sumEvens :: Vector Int -> Int
-sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0
-- sum and print running total of evens
-sumEvensAndPrint :: Vector Int -> VectorMonadT ()
-sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do
-    when (i `mod` 2 == 0) $ do
-        v <- get
-        put $ v + x
-        v' <- get
-        liftIO $ putStr $ (show v') ++ " "
-indexPlusSum :: Vector Int -> VectorMonadT ()
-indexPlusSum v' = do
-    let f i x = do
-            s <- get
-            let inc = x+s
-            liftIO $ putStr $ show (i,inc) ++ " "
-            put inc
-            return inc
-    v <- mapVectorWithIndexM f v'
-    liftIO $ do
-        putStrLn ""
-        putStrLn $ show v
-------------------------------------------
-- short circuit
-monoStep :: Double -> MaybeT (State Double) ()
-monoStep d = do
-    dp <- get
-    when (d < dp) (fail "negative difference")
-    put d
-{-# INLINE monoStep #-}
-isMonotoneIncreasing :: Vector Double -> Bool
-isMonotoneIncreasing v =
-    let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0)
-     in case res of
-        Nothing -> False
-        Just _  -> True
-------------------------------------------
-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
-successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
-successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)
-   where step e = do
-                  ep <- lift $ get
-                  if t e ep
-                     then lift $ put e
-                     else (fail "successive_ test failed")
-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
-successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
-successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)
-   where step e = do
-                  ep <- get
-                  put e
-                  return $ f ep e
-------------------------------------------
-v :: Vector Int
-v = 10 |> [0..]
-w = fromList ([1..10]++[10,9..1]) :: Vector Double
-main = do
-    v' <- test1 v
-    putStrLn ""
-    putStrLn $ show v'
-    evalStateT (addInitialM v) 0
-    putStrLn ""
-    putStrLn $ show (sumEvens v)
-    evalStateT (sumEvensAndPrint v) 0
-    putStrLn ""
-    evalStateT (indexPlusSum v) 0
-    putStrLn "-----------------------"
-    mapVectorM_ print v
-    print =<< (mapVectorM (const randomIO) v :: IO (Vector Double))
-    print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double))
-    putStrLn "-----------------------"
-    print $ isMonotoneIncreasing w
-    print $ isMonotoneIncreasing (subVector 0 7 w)
-    print $ successive_ (>) v
-    print $ successive_ (>) w
-    print $ successive (+) v
diff --git a/packages/hmatrix/examples/multiply.hs b/packages/hmatrix/examples/multiply.hs
deleted file mode 100644
index 572961c..0000000
--- a/packages/hmatrix/examples/multiply.hs
+++ /dev/null
@@ -1,104 +0,0 @@
-{-# LANGUAGE UnicodeSyntax
-           , MultiParamTypeClasses
-           , FunctionalDependencies
-           , FlexibleInstances
-           , FlexibleContexts
--           , OverlappingInstances
-           , UndecidableInstances #-}
-import Numeric.LinearAlgebra
-class Scaling a b c | a b -> c where
- -- ^ 0x22C5    8901    DOT OPERATOR, scaling
- infixl 7 ⋅
- (⋅) :: a -> b -> c
-instance (Num t) => Scaling t t t where
-    (⋅) = (*)
-instance Container Vector t => Scaling t (Vector t) (Vector t) where
-    (⋅) = scale
-instance Container Vector t => Scaling (Vector t) t (Vector t) where
-    (⋅) = flip scale
-instance Container Vector t => Scaling t (Matrix t) (Matrix t) where
-    (⋅) = scale
-instance Container Vector t => Scaling (Matrix t) t (Matrix t) where
-    (⋅) = flip scale
-class Mul a b c | a b -> c, a c -> b, b c -> a where
- -- ^ 0x00D7    215     MULTIPLICATION SIGN     ×, contraction
- infixl 7 ×
- (×) :: a -> b -> c
-------
-instance Product t => Mul (Vector t) (Vector t) t where
-    (×) = udot
-instance Product t => Mul (Matrix t) (Vector t) (Vector t) where
-    (×) = mXv
-instance Product t => Mul (Vector t) (Matrix t) (Vector t) where
-    (×) = vXm
-instance Product t => Mul (Matrix t) (Matrix t) (Matrix t) where
-    (×) = mXm
--instance Scaling a b c => Contraction a b c where
--    (×) = (⋅)
--------------------------------------------------------------------------------
-class Outer a
-  where
-    infixl 7 ⊗
-    -- | unicode 0x2297 8855 CIRCLED TIMES      ⊗
-    --
-    -- vector outer product and matrix Kronecker product
-    (⊗) :: Product t => a t -> a t -> Matrix t
-instance Outer Vector where
-    (⊗) = outer
-instance Outer Matrix where
-    (⊗) = kronecker
--------------------------------------------------------------------------------
-v = 3 |> [1..] :: Vector Double
-m = (3 >< 3) [1..] :: Matrix Double
-s = 3 :: Double
-a = s ⋅ v × m × m × v ⋅ s
--b = (v ⊗ m) ⊗ (v ⊗ m)
--c = v ⊗ m ⊗ v ⊗ m
-d = s ⋅ (3 |> [10,20..] :: Vector Double)
-u = fromList [3,0,5]
-w = konst 1 (2,3) :: Matrix Double
-main = do
-    print $ (scale s v <> m) `udot` v
-    print $ scale s v `udot` (m <> v)
-    print $ s * ((v <> m) `udot` v)
-    print $ s ⋅ v × m × v
-    print a
--    print (b == c)
-    print d
-    print $ asColumn u ⊗ w
-    print $ w ⊗ asColumn u
diff --git a/packages/hmatrix/examples/ode.hs b/packages/hmatrix/examples/ode.hs
deleted file mode 100644
index dc6e0ec..0000000
--- a/packages/hmatrix/examples/ode.hs
+++ /dev/null
@@ -1,50 +0,0 @@
-{-# LANGUAGE ViewPatterns #-}
-import Numeric.GSL.ODE
-import Numeric.LinearAlgebra
-import Graphics.Plot
-import Debug.Trace(trace)
-debug x = trace (show x) x
-vanderpol mu = do
-    let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
-        ts = linspace 1000 (0,50)
-        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
-    mplot (ts : sol)
-    mplot sol
-harmonic w d = do
-    let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w
-        ts = linspace 100 (0,20)
-        sol = odeSolve (xdot w d) [1,0] ts
-    mplot (ts : toColumns sol)
-kepler v a = mplot (take 2 $ toColumns sol) where
-    xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]
-        where g=1
-              k=(-g)*(x*x+y*y)**(-1.5)
-    ts = linspace 100 (0,30)
-    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
-    degree = pi/180
-main = do
-    vanderpol 2
-    harmonic 1 0
-    harmonic 1 0.1
-    kepler 0.3 60
-    kepler 0.4 70
-    vanderpol' 2
-- example of odeSolveV with jacobian
-vanderpol' mu = do
-    let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)]
-        jac t (toList->[x,v]) = (2><2) [ 0          ,          1
-                                       , -1-2*x*v*mu, mu*(1-x**2) ]
-        ts = linspace 1000 (0,50)
-        hi = (ts@>1 - ts@>0)/100
-        sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
-    mplot sol
diff --git a/packages/hmatrix/examples/parallel.hs b/packages/hmatrix/examples/parallel.hs
deleted file mode 100644
index e875407..0000000
--- a/packages/hmatrix/examples/parallel.hs
+++ /dev/null
@@ -1,28 +0,0 @@
-- $ ghc --make -O -rtsopts -threaded parallel.hs
-- $ ./parallel 3000 +RTS -N4 -s -A200M
-import System.Environment(getArgs)
-import Numeric.LinearAlgebra
-import Control.Parallel.Strategies
-import System.Time
-inParallel = parMap rwhnf id
-- matrix product decomposed into p parallel subtasks
-parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]
-    where [ys] = toBlocksEvery (rows y) (cols y `div` p) y
-main = do
-    n <- (read . head) `fmap` getArgs
-    let m = ident n :: Matrix Double
-    time $ print $ maxElement $ takeDiag $ m <> m
-    time $ print $ maxElement $ takeDiag $ parMul 2 m m
-    time $ print $ maxElement $ takeDiag $ parMul 4 m m
-    time $ print $ maxElement $ takeDiag $ parMul 8 m m
-time act = do
-    t0 <- getClockTime
-    act
-    t1 <- getClockTime
-    print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0
diff --git a/packages/hmatrix/examples/pca1.hs b/packages/hmatrix/examples/pca1.hs
deleted file mode 100644
index a11eba9..0000000
--- a/packages/hmatrix/examples/pca1.hs
+++ /dev/null
@@ -1,46 +0,0 @@
-- Principal component analysis
-import Numeric.LinearAlgebra
-import System.Directory(doesFileExist)
-import System.Process(system)
-import Control.Monad(when)
-type Vec = Vector Double
-type Mat = Matrix Double
-- Vector with the mean value of the columns of a matrix
-mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
-- covariance matrix of a list of observations stored as rows
-cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
-    where xc = x - asRow (mean x)
-- creates the compression and decompression functions from the desired number of components
-pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec)
-pca n dataSet = (encode,decode)
-  where
-    encode x = vp <> (x - m)
-    decode x = x <> vp + m
-    m = mean dataSet
-    c = cov dataSet
-    (_,v) = eigSH' c
-    vp = takeRows n (trans v)
-norm = pnorm PNorm2
-main = do
-    ok <- doesFileExist ("mnist.txt")
-    when (not ok)  $ do
-        putStrLn "\nTrying to download test datafile..."
-        system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
-        system("gunzip mnist.txt.gz")
-        return ()
-    m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)
-    let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)
-    let x = toRows xs !! 4  -- an arbitrary test Vec
-    let (pe,pd) = pca 10 xs
-    let y = pe x
-    print y  -- compressed version
-    print $ norm (x - pd y) / norm x --reconstruction quality
diff --git a/packages/hmatrix/examples/pca2.hs b/packages/hmatrix/examples/pca2.hs
deleted file mode 100644
index e7ea95f..0000000
--- a/packages/hmatrix/examples/pca2.hs
+++ /dev/null
@@ -1,65 +0,0 @@
-- Improved PCA, including illustrative graphics
-import Numeric.LinearAlgebra
-import Graphics.Plot
-import System.Directory(doesFileExist)
-import System.Process(system)
-import Control.Monad(when)
-type Vec = Vector Double
-type Mat = Matrix Double
-- Vector with the mean value of the columns of a matrix
-mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
-- covariance matrix of a list of observations stored as rows
-cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
-    where xc = x - asRow (mean x)
-type Stat = (Vec, [Double], Mat)
-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)
-stat :: Mat -> Stat
-stat x = (m, toList s, trans v) where
-    m = mean x
-    (s,v) = eigSH' (cov x)
-- creates the compression and decompression functions from the desired reconstruction
-- quality and the statistics of a data set
-pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec)
-pca prec (m,s,v) = (encode,decode)    
-  where    
-    encode x = vp <> (x - m)
-    decode x = x <> vp + m
-    vp = takeRows n v
-    n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)
-    cumSum = tail . scanl (+) 0.0
-    prec' = if prec <=0.0 || prec >= 1.0
-                then error "the precision in pca must be 0<prec<1"
-                else prec
-shdigit :: Vec -> IO ()
-shdigit v = imshow (reshape 28 (-v))
-- shows the effect of a given reconstruction quality on a test vector
-test :: Stat -> Double -> Vec -> IO ()
-test st prec x = do
-    let (pe,pd) = pca prec st
-    let y = pe x
-    print $ dim y
-    shdigit (pd y)
-main = do
-    ok <- doesFileExist ("mnist.txt")
-    when (not ok)  $ do
-        putStrLn "\nTrying to download test datafile..."
-        system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
-        system("gunzip mnist.txt.gz")
-        return ()
-    m <- loadMatrix "mnist.txt"
-    let xs = takeColumns (cols m -1) m
-    let x = toRows xs !! 4  -- an arbitrary test vector
-    shdigit x
-    let st = stat xs
-    test st 0.90 x
-    test st 0.50 x
diff --git a/packages/hmatrix/examples/pinv.hs b/packages/hmatrix/examples/pinv.hs
deleted file mode 100644
index 7de50b8..0000000
--- a/packages/hmatrix/examples/pinv.hs
+++ /dev/null
@@ -1,20 +0,0 @@
-import Numeric.LinearAlgebra
-import Graphics.Plot
-import Text.Printf(printf)
-expand :: Int -> Vector Double -> Matrix Double
-expand n x = fromColumns $ map (x^) [0 .. n]
-polynomialModel :: Vector Double -> Vector Double -> Int
-                -> (Vector Double -> Vector Double)
-polynomialModel x y n = f where
-    f z = expand n z <> ws
-    ws  = expand n x <\> y
-main = do
-    [x,y] <- (toColumns . readMatrix) `fmap` readFile "data.txt"
-    let pol = polynomialModel x y
-    let view = [x, y, pol 1 x, pol 2 x, pol 3 x]
-    putStrLn $ "  x      y      p 1    p 2    p 3"
-    putStrLn $ format "  " (printf "%.2f") $ fromColumns view
-    mplot view
diff --git a/packages/hmatrix/examples/plot.hs b/packages/hmatrix/examples/plot.hs
deleted file mode 100644
index f950aa5..0000000
--- a/packages/hmatrix/examples/plot.hs
+++ /dev/null
@@ -1,20 +0,0 @@
-import Numeric.LinearAlgebra
-import Graphics.Plot
-import Numeric.GSL.Special(erf_Z, erf)
-sombrero n = f x y where 
-    (x,y) = meshdom range range
-    range = linspace n (-2,2)
-    f x y = exp (-r2) * cos (2*r2) where 
-        r2 = x*x+y*y
-f x =  sin x + 0.5 * sin (5*x)
-gaussianPDF = erf_Z
-cumdist x = 0.5 * (1+ erf (x/sqrt 2))
-main = do
-    let x = linspace 1000 (-4,4)
-    mplot [f x]
-    mplot [x, mapVector cumdist x, mapVector gaussianPDF x]
-    mesh (sombrero 40)
-\ No newline at end of file
diff --git a/packages/hmatrix/examples/root.hs b/packages/hmatrix/examples/root.hs
deleted file mode 100644
index 8546ff5..0000000
--- a/packages/hmatrix/examples/root.hs
+++ /dev/null
@@ -1,31 +0,0 @@
-- root finding examples
-import Numeric.GSL
-import Numeric.LinearAlgebra
-import Text.Printf(printf)
-rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
-test method = do
-    print method
-    let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5]
-    print s -- solution
-    disp p -- evolution of the algorithm
-jacobian a b [x,y] = [ [-a    , 0]
-                     , [-2*b*x, b] ]
-testJ method = do
-    print method
-    let (s,p) = rootJ method 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5]
-    print s
-    disp p
-disp = putStrLn . format "  " (printf "%.3f")
-main = do
-    test Hybrids
-    test Hybrid
-    test DNewton
-    test Broyden
-    mapM_ testJ [HybridsJ .. GNewton]
diff --git a/packages/hmatrix/examples/vector.hs b/packages/hmatrix/examples/vector.hs
deleted file mode 100644
index f531cbd..0000000
--- a/packages/hmatrix/examples/vector.hs
+++ /dev/null
@@ -1,31 +0,0 @@
-- conversion to/from Data.Vector.Storable
-- from Roman Leshchinskiy "vector" package
--
-- In the future Data.Packed.Vector will be replaced by Data.Vector.Storable
-------------------------------------------
-import Numeric.LinearAlgebra as H
-import Data.Packed.Development(unsafeFromForeignPtr, unsafeToForeignPtr)
-import Foreign.Storable
-import qualified Data.Vector.Storable as V
-fromVector :: Storable t => V.Vector t -> H.Vector t
-fromVector v = unsafeFromForeignPtr p i n where
-    (p,i,n) = V.unsafeToForeignPtr v
-toVector :: Storable t => H.Vector t -> V.Vector t
-toVector v = V.unsafeFromForeignPtr p i n where
-    (p,i,n) = unsafeToForeignPtr v
-------------------------------------------
-v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7)
-w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double
-main = do
-    print v
-    print $ fromVector v
-    print w
-    print $ toVector w
diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal
deleted file mode 100644
index 83a2188..0000000
--- a/packages/hmatrix/hmatrix.cabal
+++ /dev/null
@@ -1,175 +0,0 @@
-Name:               hmatrix-gsl
-Version:            0.16.0.0
-License:            GPL
-License-file:       LICENSE
-Author:             Alberto Ruiz
-Maintainer:         Alberto Ruiz <aruiz@um.es>
-Stability:          provisional
-Homepage:           https://github.com/albertoruiz/hmatrix
-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
-cabal-version:      >=1.8
-build-type:         Simple
-extra-source-files: THANKS.md INSTALL.md CHANGELOG
-extra-source-files: examples/deriv.hs
-                    examples/integrate.hs
-                    examples/minimize.hs
-                    examples/root.hs
-                    examples/ode.hs
-                    examples/pca1.hs
-                    examples/pca2.hs
-                    examples/pinv.hs
-                    examples/data.txt
-                    examples/lie.hs
-                    examples/kalman.hs
-                    examples/parallel.hs
-                    examples/plot.hs
-                    examples/inplace.hs
-                    examples/error.hs
-                    examples/fitting.hs
-                    examples/devel/ej1/wrappers.hs
-                    examples/devel/ej1/functions.c
-                    examples/devel/ej2/wrappers.hs
-                    examples/devel/ej2/functions.c
-                    examples/vector.hs
-                    examples/monadic.hs
-                    examples/bool.hs
-                    examples/multiply.hs
-extra-source-files: src/Numeric/GSL/gsl-ode.c
-flag dd
-    description:    svd = zgesdd
-    default:        True
-flag mkl
-    description:    Link with Intel's MKL optimized libraries.
-    default:        False
-flag unsafe
-    description:    Compile the library with bound checking disabled.
-    default:        False
-flag finit
-    description:    Force FPU initialization in foreing calls
-    default:        False
-flag debugfpu
-    description:    Check FPU stack
-    default:        False
-flag debugnan
-    description:    Check NaN
-    default:        False
-library
-    Build-Depends:      base, hmatrix, array, vector,
-                        process, random
-    Extensions:         ForeignFunctionInterface,
-                        CPP
-    hs-source-dirs:     src
-    Exposed-modules:    Numeric.GSL.Differentiation,
-                        Numeric.GSL.Integration,
-                        Numeric.GSL.Fourier,
-                        Numeric.GSL.Polynomials,
-                        Numeric.GSL.Minimization,
-                        Numeric.GSL.Root,
-                        Numeric.GSL.Fitting,
-                        Numeric.GSL.ODE,
-                        Numeric.GSL,
-                        Numeric.GSL.LinearAlgebra,
-                        Graphics.Plot
-    other-modules:      Numeric.GSL.Internal,
-                        Numeric.GSL.Vector,
-                        Numeric.GSL.IO,
-                        Numeric.GSL.Random
-    C-sources:          src/Numeric/GSL/gsl-aux.c
-    cc-options:         -O4 -msse2 -Wall
-    cpp-options:        -DBINARY
- -- ghc-prof-options:   -auto
-    ghc-options:  -Wall -fno-warn-missing-signatures
-                        -fno-warn-orphans
-                        -fno-warn-unused-binds
-    if flag(unsafe)
-        cpp-options: -DUNSAFE
-    if !flag(dd)
-        cpp-options: -DNOZGESDD
-    if impl(ghc < 6.10.2)
-        cpp-options: -DFINIT
-    if impl(ghc == 7.0.1)
-        cpp-options: -DFINIT
-    if impl(ghc == 7.0.2)
-        cpp-options: -DFINIT
-    if flag(finit)
-        cpp-options: -DFINIT
-    if flag(debugfpu)
-        cc-options: -DFPUDEBUG
-    if flag(debugnan)
-        cc-options: -DNANDEBUG
-    if impl(ghc == 7.0.1)
-        cpp-options: -DNONORMVTEST
-    if flag(mkl)
-      if arch(x86_64)
-        extra-libraries:   gsl mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core
-      else
-        extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core
-    if os(OSX)
-        extra-lib-dirs: /opt/local/lib/
-        include-dirs: /opt/local/include/
-        extra-lib-dirs: /usr/local/lib/
-        include-dirs: /usr/local/include/
-        extra-libraries: gsl
-        if arch(i386)
-            cc-options: -arch i386
-        frameworks: Accelerate
-    if os(freebsd)
-       extra-lib-dirs: /usr/local/lib
-       include-dirs: /usr/local/include
-       extra-libraries: gsl blas lapack
-    if os(windows)
-        extra-libraries: gsl-0 blas lapack
-    if os(linux)
-        if arch(x86_64)
-            cc-options: -fPIC
-    extra-libraries: gsl
-source-repository head
-    type:     git
-    location: https://github.com/albertoruiz/hmatrix
-- The tests are in package hmatrix-tests
diff --git a/packages/hmatrix/src/Graphics/Plot.hs b/packages/hmatrix/src/Graphics/Plot.hs
deleted file mode 100644
index 0ea41ac..0000000
--- a/packages/hmatrix/src/Graphics/Plot.hs
+++ /dev/null
@@ -1,184 +0,0 @@
-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Plot
-- Copyright   :  (c) Alberto Ruiz 2005-8
-- License     :  GPL-style
--
-- Maintainer  :  Alberto Ruiz (aruiz at um dot es)
-- Stability   :  provisional
-- Portability :  uses gnuplot and ImageMagick
--
-- This module is deprecated. It can be replaced by improved drawing tools
-- available in the plot\\plot-gtk packages by Vivian McPhail or Gnuplot by Henning Thielemann.
-----------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
-module Graphics.Plot(
-    mplot,
-    plot, parametricPlot, 
-    splot, mesh, meshdom,
-    matrixToPGM, imshow,
-    gnuplotX, gnuplotpdf, gnuplotWin
-) where
-import Numeric.Container
-import Data.List(intersperse)
-import System.Process (system)
-- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions.
-meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double)
-meshdom r1 r2 = (outer r1 (constant 1 (dim r2)), outer (constant 1 (dim r1)) r2)
-{- | Draws a 3D surface representation of a real matrix.
-> > mesh $ build (10,10) (\\i j -> i + (j-5)^2)
-In certain versions you can interactively rotate the graphic using the mouse.
-}
-mesh :: Matrix Double -> IO ()
-mesh m = gnuplotX (command++dat) where
-    command = "splot "++datafollows++" matrix with lines\n"
-    dat = prep $ toLists m
-{- | Draws the surface represented by the function f in the desired ranges and number of points, internally using 'mesh'.
-> > let f x y = cos (x + y) 
-> > splot f (0,pi) (0,2*pi) 50    
-}
-splot :: (Matrix Double->Matrix Double->Matrix Double) -> (Double,Double) -> (Double,Double) -> Int -> IO () 
-splot f rx ry n = mesh z where
-    (x,y) = meshdom (linspace n rx) (linspace n ry)
-    z = f x y
-{- | plots several vectors against the first one 
-> > let t = linspace 100 (-3,3) in mplot [t, sin t, exp (-t^2)]
-}
-mplot :: [Vector Double] -> IO ()
-mplot m = gnuplotX (commands++dats) where
-    commands = if length m == 1 then command1 else commandmore
-    command1 = "plot "++datafollows++" with lines\n" ++ dat
-    commandmore = "plot " ++ plots ++ "\n"
-    plots = concat $ intersperse ", " (map cmd [2 .. length m])
-    cmd k = datafollows++" using 1:"++show k++" with lines"
-    dat = prep $ toLists $ fromColumns m
-    dats = concat (replicate (length m-1) dat)
-{- | Draws a list of functions over a desired range and with a desired number of points 
-> > plot [sin, cos, sin.(3*)] (0,2*pi) 1000
-}
-plot :: [Vector Double->Vector Double] -> (Double,Double) -> Int -> IO ()
-plot fs rx n = mplot (x: mapf fs x)
-    where x = linspace n rx
-          mapf gs y = map ($ y) gs
-{- | Draws a parametric curve. For instance, to draw a spiral we can do something like:
-> > parametricPlot (\t->(t * sin t, t * cos t)) (0,10*pi) 1000
-}
-parametricPlot :: (Vector Double->(Vector Double,Vector Double)) -> (Double, Double) -> Int -> IO ()
-parametricPlot f rt n = mplot [fx, fy]
-    where t = linspace n rt
-          (fx,fy) = f t
-- | writes a matrix to pgm image file
-matrixToPGM :: Matrix Double -> String
-matrixToPGM m = header ++ unlines (map unwords ll) where
-    c = cols m
-    r = rows m
-    header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"
-    maxgray = 255.0
-    maxval = maxElement m
-    minval = minElement m
-    scale' = if maxval == minval
-        then 0.0
-        else maxgray / (maxval - minval)
-    f x = show ( round ( scale' *(x - minval) ) :: Int )
-    ll = map (map f) (toLists m)
-- | imshow shows a representation of a matrix as a gray level image using ImageMagick's display.
-imshow :: Matrix Double -> IO ()
-imshow m = do
-    _ <- system $ "echo \""++ matrixToPGM m ++"\"| display -antialias -resize 300 - &"
-    return ()
----------------------------------------------------
-gnuplotX :: String -> IO ()
-gnuplotX command = do { _ <- system cmdstr; return()} where
-    cmdstr = "echo \""++command++"\" | gnuplot -persist"
-datafollows = "\\\"-\\\""
-prep = (++"e\n\n") . unlines . map (unwords . map show)
-gnuplotpdf :: String -> String -> [([[Double]], String)] -> IO ()
-gnuplotpdf title command ds = gnuplot (prelude ++ command ++" "++ draw) >> postproc where
-    prelude = "set terminal epslatex color; set output '"++title++".tex';"
-    (dats,defs) = unzip ds
-    draw = concat (intersperse ", " (map ("\"-\" "++) defs)) ++ "\n" ++
-           concatMap pr dats
-    postproc = do
-        _ <- system $ "epstopdf "++title++".eps"
-        mklatex
-        _ <- system $ "pdflatex "++title++"aux.tex > /dev/null"
-        _ <- system $ "pdfcrop "++title++"aux.pdf > /dev/null"
-        _ <- system $ "mv "++title++"aux-crop.pdf "++title++".pdf"
-        _ <- system $ "rm "++title++"aux.* "++title++".eps "++title++".tex"
-        return ()
-    mklatex = writeFile (title++"aux.tex") $
-       "\\documentclass{article}\n"++
-       "\\usepackage{graphics}\n"++
-       "\\usepackage{nopageno}\n"++
-       "\\usepackage{txfonts}\n"++
-       "\\renewcommand{\\familydefault}{phv}\n"++
-       "\\usepackage[usenames]{color}\n"++
-       "\\begin{document}\n"++
-       "\\begin{center}\n"++
-       "  \\input{./"++title++".tex}\n"++
-       "\\end{center}\n"++
-       "\\end{document}"
-    pr = (++"e\n") . unlines . map (unwords . map show)
-    gnuplot cmd = do
-        writeFile "gnuplotcommand" cmd
-        _ <- system "gnuplot gnuplotcommand"
-        _ <- system "rm gnuplotcommand"
-        return ()
-gnuplotWin :: String -> String -> [([[Double]], String)] -> IO ()
-gnuplotWin title command ds = gnuplot (prelude ++ command ++" "++ draw) where
-    (dats,defs) = unzip ds
-    draw = concat (intersperse ", " (map ("\"-\" "++) defs)) ++ "\n" ++
-           concatMap pr dats
-    pr = (++"e\n") . unlines . map (unwords . map show)
-    prelude = "set title \""++title++"\";"
-    gnuplot cmd = do
-        writeFile "gnuplotcommand" cmd
-        _ <- system "gnuplot -persist gnuplotcommand"
-        _ <- system "rm gnuplotcommand"
-        return ()
diff --git a/packages/hmatrix/src/Numeric/GSL.hs b/packages/hmatrix/src/Numeric/GSL.hs
deleted file mode 100644
index 61b8d7b..0000000
--- a/packages/hmatrix/src/Numeric/GSL.hs
+++ /dev/null
@@ -1,43 +0,0 @@
-{- |
-Module      :  Numeric.GSL
-Copyright   :  (c) Alberto Ruiz 2006-14
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-This module reexports all available GSL functions.
-The GSL special functions are in the separate package hmatrix-special.
-}
-module Numeric.GSL (
-  module Numeric.GSL.Integration
-, module Numeric.GSL.Differentiation
-, module Numeric.GSL.Fourier
-, module Numeric.GSL.Polynomials
-, module Numeric.GSL.Minimization
-, module Numeric.GSL.Root
-, module Numeric.GSL.ODE
-, module Numeric.GSL.Fitting
-, module Data.Complex
-, setErrorHandlerOff
-) where
-import Numeric.GSL.Integration
-import Numeric.GSL.Differentiation
-import Numeric.GSL.Fourier
-import Numeric.GSL.Polynomials
-import Numeric.GSL.Minimization
-import Numeric.GSL.Root
-import Numeric.GSL.ODE
-import Numeric.GSL.Fitting
-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
deleted file mode 100644
index 0fb58ef..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
+++ /dev/null
@@ -1,85 +0,0 @@
-{- |
-Module      :  Numeric.GSL.Differentiation
-Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-Numerical differentiation.
-<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#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,
-    derivBackward
-) where
-import Foreign.C.Types
-import Foreign.Marshal.Alloc(malloc, free)
-import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
-import Foreign.Storable(peek)
-import Numeric.GSL.Internal
-import System.IO.Unsafe(unsafePerformIO)
-derivGen ::
-    CInt                   -- ^ type: 0 central, 1 forward, 2 backward
-    -> Double             -- ^ initial step size
-    -> (Double -> Double) -- ^ function
-    -> Double             -- ^ point where the derivative is taken
-    -> (Double, Double)   -- ^ result and error
-derivGen c h f x = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    fp <- mkfun (\y _ -> f y) 
-    c_deriv c fp x h r e // check "deriv"
-    vr <- peek r
-    ve <- peek e
-    let result = (vr,ve)
-    free r
-    free e
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "gsl-aux.h deriv" 
- c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double 
-                    -> Ptr Double -> Ptr Double -> IO CInt
-{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
->>> let deriv = derivCentral 0.01
->>> deriv sin (pi/4)
-(0.7071067812000676,1.0600063101654055e-10)
->>> cos (pi/4)
-0.7071067811865476
-}
-derivCentral :: Double                  -- ^ initial step size
-                -> (Double -> Double)   -- ^ function 
-                -> Double               -- ^ point where the derivative is taken
-                -> (Double, Double)     -- ^ result and absolute error
-derivCentral = derivGen 0
-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step.
-derivForward :: Double                  -- ^ initial step size
-                -> (Double -> Double)   -- ^ function 
-                -> Double               -- ^ point where the derivative is taken
-                -> (Double, Double)     -- ^ result and absolute error
-derivForward = derivGen 1
-- | Adaptive backward difference algorithm, /gsl_deriv_backward/. 
-derivBackward ::Double                  -- ^ initial step size
-                -> (Double -> Double)   -- ^ function 
-                -> Double               -- ^ point where the derivative is taken
-                -> (Double, Double)     -- ^ result and absolute error
-derivBackward = derivGen 2
-{- | conversion of Haskell functions into function pointers that can be used in the C side
-}
-foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 
diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
deleted file mode 100644
index 93fb281..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs
+++ /dev/null
@@ -1,180 +0,0 @@
-{- |
-Module      :  Numeric.GSL.Fitting
-Copyright   :  (c) Alberto Ruiz 2010
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-Nonlinear Least-Squares Fitting
-<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>
-The example program in the GSL manual (see examples/fitting.hs):
-@
-dat = [
- ([0.0],([6.0133918608118675],0.1)),
- ([1.0],([5.5153769909966535],0.1)),
- ([2.0],([5.261094606015287],0.1)),
- ...
- ([39.0],([1.0619821710802808],0.1))]
-expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
-expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
-(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
-@
->>> path
-(6><5)
- [ 1.0,  76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
- , 2.0, 37.683816318260355,  2.858760367632973,  8.092094813253975e-2, 1.4479636296208662
- , 3.0,    9.5807893736187,  4.948995119561291,   0.11942927999921617, 1.0945766509238248
- , 4.0,  5.630494933603935,  5.021755718065913,   0.10287787128056883, 1.0338835440862608
- , 5.0,  5.443976278682909,  5.045204331329302,   0.10405523433131504,  1.019416067207375
- , 6.0, 5.4439736648994685,  5.045357818922331,   0.10404905846029407, 1.0192487112786812 ]
->>> sol
-[(5.045357818922331,6.027976702418132e-2),
-(0.10404905846029407,3.157045047172834e-3),
-(1.0192487112786812,3.782067731353722e-2)]
-}
-module Numeric.GSL.Fitting (
-    -- * Levenberg-Marquardt
-    nlFitting, FittingMethod(..),
-    -- * Utilities
-    fitModelScaled, fitModel
-) where
-import Numeric.LinearAlgebra
-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 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)
-- | Nonlinear multidimensional least-squares fitting.
-nlFitting :: FittingMethod
-      -> Double                     -- ^ absolute tolerance
-      -> Double                     -- ^ relative tolerance
-      -> Int                        -- ^ maximum number of iterations allowed
-      -> (Vector Double -> Vector Double)     -- ^ function to be minimized
-      -> (Vector Double -> Matrix Double)   -- ^ Jacobian
-      -> Vector Double                   -- ^ starting point
-      -> (Vector Double, Matrix Double)  -- ^ solution vector and optimization path
-nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
-nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
-    let p   = dim xiv
-        n   = dim (f xiv)
-    fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
-    jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
-    rawpath <- createMatrix RowMajor maxit (2+p)
-    app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toRows $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    freeHaskellFunPtr jp
-    return (subVector 2 p sol, path)
-foreign import ccall safe "nlfit"
-    c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
-------------------------------------------------------
-checkdim1 n _p v
-    | dim v == n = v
-    | otherwise = error $ "Error: "++ show n
-                        ++ " components expected in the result of the function supplied to nlFitting"
-checkdim2 n p m
-    | rows m == n && cols m == p = m
-    | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
-                        ++ " Jacobian expected in nlFitting"
------------------------------------------------------------
-err (model,deriv) dat vsol = zip sol errs where
-    sol = toList vsol
-    c = max 1 (chi/sqrt (fromIntegral dof))
-    dof = length dat - (rows cov)
-    chi = norm2 (fromList $ cost (resMs model) dat sol)
-    js = fromLists $ jacobian (resDs deriv) dat sol
-    cov = inv $ trans js <.> js
-    errs = toList $ scalar c * sqrt (takeDiag cov)
-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
-- instances (x, (y,sigma)) to be fitted.
-fitModelScaled
-         :: Double -- ^ absolute tolerance
-         -> Double -- ^ relative tolerance
-         -> Int    -- ^ maximum number of iterations allowed
-         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
-         -> [(x, ([Double], Double))] -- ^ instances
-         -> [Double] -- ^ starting point
-         -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
-fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
-    (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
-                (fromList . cost (resMs model) dt . toList)
-                (fromLists . jacobian (resDs deriv) dt . toList)
-                (fromList xin)
-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
-- instances (x,y) to be fitted.
-fitModel :: Double -- ^ absolute tolerance
-         -> Double -- ^ relative tolerance
-         -> Int    -- ^ maximum number of iterations allowed
-         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
-         -> [(x, [Double])]  -- ^ instances
-         -> [Double] -- ^ starting point
-         -> ([Double], Matrix Double) -- ^ solution and optimization path
-fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
-    (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
-                (fromList . cost (resM model) dt . toList)
-                (fromLists . jacobian (resD deriv) dt . toList)
-                (fromList xin)
-cost model ds vs = concatMap (model vs) ds
-jacobian modelDer ds vs = concatMap (modelDer vs) ds
-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
-resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
-resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
-- | Associated derivative for 'resMs'.
-resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
-resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
-- to 'resMs' with all sigmas = 1.
-resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
-resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
-- | Associated derivative for 'resM'.
-resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
-resD m v = \(x,_) -> m v x
diff --git a/packages/hmatrix/src/Numeric/GSL/Fourier.hs b/packages/hmatrix/src/Numeric/GSL/Fourier.hs
deleted file mode 100644
index 734325b..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Fourier.hs
+++ /dev/null
@@ -1,44 +0,0 @@
-{- |
-Module      : Numeric.GSL.Fourier
-Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-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
-import Numeric.GSL.Internal
-import Data.Complex
-import Foreign.C.Types
-import System.IO.Unsafe (unsafePerformIO)
-genfft code v = unsafePerformIO $ do
-    r <- createVector (dim v)
-    app2 (c_fft code) vec v vec r "fft"
-    return r
-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.
->>> fft (fromList [1,2,3,4])
-fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
-}
-fft :: Vector (Complex Double) -> Vector (Complex Double)
-fft = genfft 0
-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
-ifft :: Vector (Complex Double) -> Vector (Complex Double)
-ifft = genfft 1
diff --git a/packages/hmatrix/src/Numeric/GSL/IO.hs b/packages/hmatrix/src/Numeric/GSL/IO.hs
deleted file mode 100644
index 0d6031a..0000000
--- a/packages/hmatrix/src/Numeric/GSL/IO.hs
+++ /dev/null
@@ -1,42 +0,0 @@
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.GSL.IO
-- Copyright   :  (c) Alberto Ruiz 2007-14
-- License     :  GPL
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-----------------------------------------------------------------------------
-module Numeric.GSL.IO (
-    saveMatrix,
-    fwriteVector, freadVector, fprintfVector, fscanfVector,
-    fileDimensions, loadMatrix, fromFile
-) where
-import Data.Packed
-import Numeric.GSL.Vector
-import System.Process(readProcess)
-{- |  obtains the number of rows and columns in an ASCII data file
-      (provisionally using unix's wc).
-}
-fileDimensions :: FilePath -> IO (Int,Int)
-fileDimensions fname = do
-    wcres <- readProcess "wc" ["-w",fname] ""
-    contents <- readFile fname
-    let tot = read . head . words $ wcres
-        c   = length . head . dropWhile null . map words . lines $ contents
-    if tot > 0
-        then return (tot `div` c, c)
-        else return (0,0)
-- | Loads a matrix from an ASCII file formatted as a 2D table.
-loadMatrix :: FilePath -> IO (Matrix Double)
-loadMatrix file = fromFile file =<< fileDimensions file
-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
-fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
-fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
diff --git a/packages/hmatrix/src/Numeric/GSL/Integration.hs b/packages/hmatrix/src/Numeric/GSL/Integration.hs
deleted file mode 100644
index 9c1d43a..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Integration.hs
+++ /dev/null
@@ -1,246 +0,0 @@
-{- |
-Module      :  Numeric.GSL.Integration
-Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-Numerical integration routines.
-<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
-}
-module Numeric.GSL.Integration (
-    integrateQNG,
-    integrateQAGS,
-    integrateQAGI,
-    integrateQAGIU,
-    integrateQAGIL,
-    integrateCQUAD
-) where
-import Foreign.C.Types
-import Foreign.Marshal.Alloc(malloc, free)
-import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
-import Foreign.Storable(peek)
-import Numeric.GSL.Internal
-import System.IO.Unsafe(unsafePerformIO)
-eps = 1e-12
-{- | conversion of Haskell functions into function pointers that can be used in the C side
-}
-foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
--------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
->>> let quad = integrateQAGS 1E-9 1000
->>> let f a x = x**(-0.5) * log (a*x)
->>> quad (f 1) 0 1
-(-3.999999999999974,4.871658632055187e-13)
-}
-integrateQAGS :: Double               -- ^ precision (e.g. 1E-9)
-                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
-                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
-                 -> Double           -- ^ a
-                 -> Double           -- ^ b
-                 -> (Double, Double) -- ^ result of the integration and error
-integrateQAGS prec n f a b = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    fp <- mkfun (\x _ -> f x)
-    c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
-    vr <- peek r
-    ve <- peek e
-    let result = (vr,ve)
-    free r
-    free e
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "integrate_qags" c_integrate_qags
-    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
-    -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
-----------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
->>> let quad = integrateQNG 1E-6
->>> quad (\x -> 4/(1+x*x)) 0 1
-(3.141592653589793,3.487868498008632e-14)
-}
-integrateQNG :: Double               -- ^ precision (e.g. 1E-9)
-                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
-                 -> Double           -- ^ a
-                 -> Double           -- ^ b
-                 -> (Double, Double) -- ^ result of the integration and error
-integrateQNG prec f a b = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    fp <- mkfun (\x _ -> f x)
-    c_integrate_qng fp a b eps prec r e  // check "integrate_qng"
-    vr <- peek r
-    ve <- peek e
-    let result = (vr,ve)
-    free r
-    free e
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "integrate_qng" c_integrate_qng
-    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
-    -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
--------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
-For example:
->>> let quad = integrateQAGI 1E-9 1000
->>> let f a x = exp(-a * x^2)
->>> quad (f 0.5)
-(2.5066282746310002,6.229215880648858e-11)
-}
-integrateQAGI :: Double               -- ^ precision (e.g. 1E-9)
-                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
-                 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
-                 -> (Double, Double) -- ^ result of the integration and error
-integrateQAGI prec n f = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    fp <- mkfun (\x _ -> f x)
-    c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
-    vr <- peek r
-    ve <- peek e
-    let result = (vr,ve)
-    free r
-    free e
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "integrate_qagi" c_integrate_qagi
-    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
-    -> CInt -> Ptr Double -> Ptr Double -> IO CInt
--------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
-For example:
->>> let quad = integrateQAGIU 1E-9 1000
->>> let f a x = exp(-a * x^2)
->>> quad (f 0.5) 0
-(1.2533141373155001,3.114607940324429e-11)
-}
-integrateQAGIU :: Double               -- ^ precision (e.g. 1E-9)
-                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
-                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
-                 -> Double             -- ^ a
-                 -> (Double, Double) -- ^ result of the integration and error
-integrateQAGIU prec n f a = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    fp <- mkfun (\x _ -> f x)
-    c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
-    vr <- peek r
-    ve <- peek e
-    let result = (vr,ve)
-    free r
-    free e
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
-    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
-    -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
--------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
-For example:
->>> let quad = integrateQAGIL 1E-9 1000
->>> let f a x = exp(-a * x^2)
->>> quad (f 0.5) 0
-(1.2533141373155001,3.114607940324429e-11)
-}
-integrateQAGIL :: Double               -- ^ precision (e.g. 1E-9)
-                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
-                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
-                 -> Double             -- ^ b
-                 -> (Double, Double) -- ^ result of the integration and error
-integrateQAGIL prec n f b = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    fp <- mkfun (\x _ -> f x)
-    c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
-    vr <- peek r
-    ve <- peek e
-    let result = (vr,ve)
-    free r
-    free e
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil
-    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
-    -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
--------------------------------------------------------------------
-{- | Numerical integration using /gsl_integration_cquad/ (quadrature
-for general integrands).  From the GSL manual:
-@CQUAD is a new doubly-adaptive general-purpose quadrature routine
-which can handle most types of singularities, non-numerical function
-values such as Inf or NaN, as well as some divergent integrals. It
-generally requires more function evaluations than the integration
-routines in QUADPACK, yet fails less often for difficult integrands.@
-For example:
->>> let quad = integrateCQUAD 1E-12 1000
->>> let f a x = exp(-a * x^2)
->>> quad (f 0.5) 2 5
-(5.7025405463957006e-2,9.678874441303705e-16,95)
-Unlike other quadrature methods, integrateCQUAD also returns the
-number of function evaluations required.
-}
-integrateCQUAD :: Double               -- ^ precision (e.g. 1E-9)
-                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
-                 -> (Double -> Double) -- ^ function to be integrated on the interval (a, b)
-                 -> Double             -- ^ a
-                 -> Double             -- ^ b
-                 -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed
-integrateCQUAD prec n f a b = unsafePerformIO $ do
-    r <- malloc
-    e <- malloc
-    neval <- malloc
-    fp <- mkfun (\x _ -> f x)
-    c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad"
-    vr <- peek r
-    ve <- peek e
-    vneval <- peek neval
-    let result = (vr,ve,vneval)
-    free r
-    free e
-    free neval
-    freeHaskellFunPtr fp
-    return result
-foreign import ccall safe "integrate_cquad" c_integrate_cquad
-    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
-    -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt
diff --git a/packages/hmatrix/src/Numeric/GSL/Internal.hs b/packages/hmatrix/src/Numeric/GSL/Internal.hs
deleted file mode 100644
index a1c4e0c..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Internal.hs
+++ /dev/null
@@ -1,126 +0,0 @@
-- |
-- Module      :  Numeric.GSL.Internal
-- Copyright   :  (c) Alberto Ruiz 2009
-- License     :  GPL
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
--
-- Auxiliary functions.
--
-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
-    copy n' q = do
-        copyArray q p (fromIntegral n')
-        return 0
-- | conversion of Haskell functions into function pointers that can be used in the C side
-foreign import ccall safe "wrapper"
-    mkVecfun :: (CInt -> Ptr Double -> Double)
-             -> IO( FunPtr (CInt -> Ptr Double -> Double))
-foreign import ccall safe "wrapper"
-    mkVecVecfun :: TVV -> IO (FunPtr TVV)
-foreign import ccall safe "wrapper"
-    mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
-foreign import ccall safe "wrapper"
-    mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double))
-aux_vTov :: (Vector Double -> Vector Double) -> TVV
-aux_vTov f n p nr r = g where
-    v = f x
-    x = createV (fromIntegral n) copy "aux_vTov"
-    copy n' q = do
-        copyArray q p (fromIntegral n')
-        return 0
-    g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
-           return 0
-foreign import ccall safe "wrapper"
-    mkVecMatfun :: TVM -> IO (FunPtr TVM)
-foreign import ccall safe "wrapper"
-    mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
-aux_vTom :: (Vector Double -> Matrix Double) -> TVM
-aux_vTom f n p rr cr r = g where
-    v = flatten $ f x
-    x = createV (fromIntegral n) copy "aux_vTov"
-    copy n' q = do
-        copyArray q p (fromIntegral n')
-        return 0
-    g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
-           return 0
-createV n fun msg = unsafePerformIO $ do
-    r <- createVector n
-    app1 fun vec r msg
-    return r
-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/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
deleted file mode 100644
index 17e2258..0000000
--- a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
+++ /dev/null
@@ -1,135 +0,0 @@
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.GSL.LinearAlgebra
-- Copyright   :  (c) Alberto Ruiz 2007-14
-- License     :  GPL
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-----------------------------------------------------------------------------
-module Numeric.GSL.LinearAlgebra (
-    RandDist(..), randomVector,
-    saveMatrix,
-    fwriteVector, freadVector, fprintfVector, fscanfVector,
-    fileDimensions, loadMatrix, fromFile
-) where
-import Data.Packed
-import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
-import Foreign.Marshal.Alloc(free)
-import Foreign.Ptr(Ptr)
-import Foreign.C.Types
-import Foreign.C.String(newCString)
-import System.IO.Unsafe(unsafePerformIO)
-import System.Process(readProcess)
-fromei x = fromIntegral (fromEnum x) :: CInt
-----------------------------------------------------------------------
-data RandDist = Uniform  -- ^ uniform distribution in [0,1)
-              | Gaussian -- ^ normal distribution with mean zero and standard deviation one
-              deriving Enum
-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
-randomVector :: Int      -- ^ seed
-             -> RandDist -- ^ distribution
-             -> Int      -- ^ vector size
-             -> Vector Double
-randomVector seed dist n = unsafePerformIO $ do
-    r <- createVector n
-    app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
-    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 PD = Ptr Double                            --
-type TV = CInt -> PD -> IO CInt                 --
-type TM = CInt -> CInt -> PD -> IO CInt         --
--------------------------------------------------------------------------------
-{- |  obtains the number of rows and columns in an ASCII data file
-      (provisionally using unix's wc).
-}
-fileDimensions :: FilePath -> IO (Int,Int)
-fileDimensions fname = do
-    wcres <- readProcess "wc" ["-w",fname] ""
-    contents <- readFile fname
-    let tot = read . head . words $ wcres
-        c   = length . head . dropWhile null . map words . lines $ contents
-    if tot > 0
-        then return (tot `div` c, c)
-        else return (0,0)
-- | Loads a matrix from an ASCII file formatted as a 2D table.
-loadMatrix :: FilePath -> IO (Matrix Double)
-loadMatrix file = fromFile file =<< fileDimensions file
-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
-fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
-fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
diff --git a/packages/hmatrix/src/Numeric/GSL/Minimization.hs b/packages/hmatrix/src/Numeric/GSL/Minimization.hs
deleted file mode 100644
index 056d463..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Minimization.hs
+++ /dev/null
@@ -1,222 +0,0 @@
-{- |
-Module      :  Numeric.GSL.Minimization
-Copyright   :  (c) Alberto Ruiz 2006-9
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-Minimization of a multidimensional function using some of the algorithms described in:
-<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
-The example in the GSL manual:
-@
-f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
-main = do
-    let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7]
-    print s
-    print p
-@
->>> main
-[0.9920430849306288,1.9969168063253182]
- 0.000  512.500  1.130  6.500  5.000
- 1.000  290.625  1.409  5.250  4.000
- 2.000  290.625  1.409  5.250  4.000
- 3.000  252.500  1.409  5.500  1.000
- ...
-22.000   30.001  0.013  0.992  1.997
-23.000   30.001  0.008  0.992  1.997
-The path to the solution can be graphically shown by means of:
-@'Graphics.Plot.mplot' $ drop 3 ('toColumns' p)@
-Taken from the GSL manual:
-The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.
-The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's Practical Methods of Optimization, Algorithms 2.6.2 and 2.6.4. It supercedes the original bfgs routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance tol corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches).
-The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimplex minimiser. It calculates the size of simplex as the rms distance of each vertex from the center rather than the mean distance, which has the advantage of allowing a linear update.
-}
-module Numeric.GSL.Minimization (
-    minimize, minimizeV, MinimizeMethod(..),
-    minimizeD, minimizeVD, MinimizeMethodD(..),
-    uniMinimize, UniMinimizeMethod(..),
-    minimizeNMSimplex,
-    minimizeConjugateGradient,
-    minimizeVectorBFGS2
-) where
-import Data.Packed
-import Numeric.GSL.Internal
-import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
-import Foreign.C.Types
-import System.IO.Unsafe(unsafePerformIO)
------------------------------------------------------------------------
-{-# DEPRECATED minimizeNMSimplex "use minimize NMSimplex2 eps maxit sizes f xi" #-}
-minimizeNMSimplex f xi szs eps maxit = minimize NMSimplex eps maxit szs f xi
-{-# DEPRECATED minimizeConjugateGradient "use minimizeD ConjugateFR eps maxit step tol f g xi" #-}
-minimizeConjugateGradient step tol eps maxit f g xi = minimizeD ConjugateFR eps maxit step tol f g xi
-{-# DEPRECATED minimizeVectorBFGS2 "use minimizeD VectorBFGS2 eps maxit step tol f g xi" #-}
-minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit step tol f g xi
-------------------------------------------------------------------------
-data UniMinimizeMethod = GoldenSection
-                       | BrentMini
-                       | QuadGolden
-                       deriving (Enum, Eq, Show, Bounded)
-- | Onedimensional minimization.
-uniMinimize :: UniMinimizeMethod -- ^ The method used.
-        -> Double               -- ^ desired precision of the solution
-        -> Int                  -- ^ maximum number of iterations allowed
-        -> (Double -> Double)   -- ^ function to minimize
-        -> Double               -- ^ guess for the location of the minimum
-        -> Double               -- ^ lower bound of search interval
-        -> Double               -- ^ upper bound of search interval
-        -> (Double, Matrix Double) -- ^ solution and optimization path
-uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit
-uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
-    fp <- mkDoublefun f
-    rawpath <- createMIO maxit 4
-                         (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
-                         "uniMinimize"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toLists $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    return (sol !! 1, path)
-foreign import ccall safe "uniMinimize"
-    c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM Res
-data MinimizeMethod = NMSimplex
-                    | NMSimplex2
-                    deriving (Enum,Eq,Show,Bounded)
-- | Minimization without derivatives
-minimize :: MinimizeMethod
-         -> Double              -- ^ desired precision of the solution (size test)
-         -> Int                 -- ^ maximum number of iterations allowed
-         -> [Double]            -- ^ sizes of the initial search box
-         -> ([Double] -> Double) -- ^ function to minimize
-         -> [Double]            -- ^ starting point
-         -> ([Double], Matrix Double) -- ^ solution vector and optimization path
-- | Minimization without derivatives (vector version)
-minimizeV :: MinimizeMethod
-         -> Double              -- ^ desired precision of the solution (size test)
-         -> Int                 -- ^ maximum number of iterations allowed
-         -> Vector Double       -- ^ sizes of the initial search box
-         -> (Vector Double -> Double) -- ^ function to minimize
-         -> Vector Double            -- ^ starting point
-         -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
-minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
-    where v2l (v,m) = (toList v, m)
-ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
-minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
-    let n   = dim xiv
-    fp <- mkVecfun (iv f)
-    rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
-                   createMIO maxit (n+3)
-                         (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
-                         "minimize"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        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 -> TV(TV(TM Res))
----------------------------------------------------------------------------------
-data MinimizeMethodD = ConjugateFR
-                     | ConjugatePR
-                     | VectorBFGS
-                     | VectorBFGS2
-                     | SteepestDescent
-                     deriving (Enum,Eq,Show,Bounded)
-- | Minimization with derivatives.
-minimizeD :: MinimizeMethodD
-    -> Double                 -- ^ desired precision of the solution (gradient test)
-    -> Int                    -- ^ maximum number of iterations allowed
-    -> Double                 -- ^ size of the first trial step
-    -> Double                 -- ^ tol (precise meaning depends on method)
-    -> ([Double] -> Double)   -- ^ function to minimize
-    -> ([Double] -> [Double]) -- ^ gradient
-    -> [Double]               -- ^ starting point
-    -> ([Double], Matrix Double) -- ^ solution vector and optimization path
-- | Minimization with derivatives (vector version)
-minimizeVD :: MinimizeMethodD
-    -> Double                 -- ^ desired precision of the solution (gradient test)
-    -> Int                    -- ^ maximum number of iterations allowed
-    -> Double                 -- ^ size of the first trial step
-    -> Double                 -- ^ tol (precise meaning depends on method)
-    -> (Vector Double -> Double)   -- ^ function to minimize
-    -> (Vector Double -> Vector Double) -- ^ gradient
-    -> Vector Double               -- ^ starting point
-    -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
-minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
-          method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)
-    where v2l (v,m) = (toList v, m)
-minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
-    let n = dim xiv
-        f' = f
-        df' = (checkdim1 n . df)
-    fp <- mkVecfun (iv f')
-    dfp <- mkVecVecfun (aux_vTov df')
-    rawpath <- vec xiv $ \xiv' ->
-                    createMIO maxit (n+2)
-                         (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
-                         "minimizeD"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        sol = flatten $ dropColumns 2 $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    freeHaskellFunPtr dfp
-    return (sol,path)
-foreign import ccall safe "gsl-aux.h minimizeD"
-    c_minimizeD :: CInt
-                -> FunPtr (CInt -> Ptr Double -> Double)
-                -> FunPtr (TV (TV Res))
-                -> Double -> Double -> Double -> CInt
-                -> TV (TM Res)
---------------------------------------------------------------------
-checkdim1 n v
-    | dim v == n = v
-    | otherwise = error $ "Error: "++ show n
-                        ++ " components expected in the result of the gradient supplied to minimizeD"
diff --git a/packages/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs
deleted file mode 100644
index 7549a65..0000000
--- a/packages/hmatrix/src/Numeric/GSL/ODE.hs
+++ /dev/null
@@ -1,140 +0,0 @@
-{- |
-Module      :  Numeric.GSL.ODE
-Copyright   :  (c) Alberto Ruiz 2010
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-Solution of ordinary differential equation (ODE) initial value problems.
-<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
-A simple example:
-@
-import Numeric.GSL.ODE
-import Numeric.LinearAlgebra
-import Numeric.LinearAlgebra.Util(mplot)
-xdot t [x,v] = [v, -0.95*x - 0.1*v]
-ts = linspace 100 (0,20 :: Double)
-sol = odeSolve xdot [10,0] ts
-main = mplot (ts : toColumns sol)
-@
-}
-----------------------------------------------------------------------------
-module Numeric.GSL.ODE (
-    odeSolve, odeSolveV, ODEMethod(..), Jacobian
-) where
-import Data.Packed
-import Numeric.GSL.Internal
-import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
-import Foreign.C.Types
-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
-data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
-               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods.
-               | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
-               | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
-               | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
-               | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
-               | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
-               | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
-               | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.
-               | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. 
-               | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.
-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
-odeSolve
-    :: (Double -> [Double] -> [Double])        -- ^ xdot(t,x)
-    -> [Double]        -- ^ initial conditions
-    -> Vector Double   -- ^ desired solution times
-    -> Matrix Double   -- ^ solution
-odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
-    where hi = (ts@>1 - ts@>0)/100
-          epsAbs = 1.49012e-08
-          epsRel = 1.49012e-08
-          l2v f = \t -> fromList  . f t . toList
-- | Evolution of the system with adaptive step-size control.
-odeSolveV
-    :: ODEMethod
-    -> Double -- ^ initial step size
-    -> Double -- ^ absolute tolerance for the state vector
-    -> Double -- ^ relative tolerance for the state vector
-    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
-    -> Vector Double     -- ^ initial conditions
-    -> Vector Double     -- ^ desired solution times
-    -> Matrix Double     -- ^ solution
-odeSolveV RK2 = odeSolveV' 0 Nothing
-odeSolveV RK4 = odeSolveV' 1 Nothing
-odeSolveV RKf45 = odeSolveV' 2 Nothing
-odeSolveV RKck = odeSolveV' 3 Nothing
-odeSolveV RK8pd = odeSolveV' 4 Nothing
-odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
-odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
-odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
-odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
-odeSolveV MSAdams = odeSolveV' 9 Nothing
-odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
-odeSolveV'
-    :: CInt
-    -> Maybe (Double -> Vector Double -> Matrix Double)   -- ^ optional jacobian
-    -> Double -- ^ initial step size
-    -> Double -- ^ absolute tolerance for the state vector
-    -> Double -- ^ relative tolerance for the state vector
-    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
-    -> Vector Double     -- ^ initial conditions
-    -> Vector Double     -- ^ desired solution times
-    -> Matrix Double     -- ^ solution
-odeSolveV' method mbjac h epsAbs epsRel f  xiv ts = unsafePerformIO $ do
-    let n   = dim xiv
-    fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
-    jp <- case mbjac of
-        Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
-        Nothing  -> return nullFunPtr
-    sol <- vec xiv $ \xiv' ->
-            vec (checkTimes ts) $ \ts' ->
-             createMIO (dim ts) n
-              (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
-              "ode"
-    freeHaskellFunPtr fp
-    return sol
-foreign import ccall safe "ode"
-    ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
-------------------------------------------------------
-checkdim1 n v
-    | dim v == n = v
-    | otherwise = error $ "Error: "++ show n
-                        ++ " components expected in the result of the function supplied to odeSolve"
-checkdim2 n m
-    | rows m == n && cols m == n = m
-    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
-                        ++ " Jacobian expected in odeSolve"
-checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
-              | otherwise = error "odeSolve requires increasing times"
-    where ts' = toList ts
diff --git a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
deleted file mode 100644
index b1be85d..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
+++ /dev/null
@@ -1,57 +0,0 @@
-{- |
-Module      :  Numeric.GSL.Polynomials
-Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-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
-import Numeric.GSL.Internal
-import Data.Complex
-import System.IO.Unsafe (unsafePerformIO)
-#if __GLASGOW_HASKELL__ >= 704
-import Foreign.C.Types (CInt(..))
-#endif
-{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. 
-For example, the three solutions of x^3 + 8 = 0
->>> polySolve [8,0,0,1]
-[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]
-The example in the GSL manual: To find the roots of x^5 -1 = 0:
->>> polySolve [-1, 0, 0, 0, 0, 1]
-[(-0.8090169943749472) :+ 0.5877852522924731,
-(-0.8090169943749472) :+ (-0.5877852522924731),
-0.30901699437494756 :+ 0.9510565162951535,
-0.30901699437494756 :+ (-0.9510565162951535),
-1.0000000000000002 :+ 0.0]
-}  
-polySolve :: [Double] -> [Complex Double]
-polySolve = toList . polySolve' . fromList
-polySolve' :: Vector Double -> Vector (Complex Double)
-polySolve' v | dim v > 1 = unsafePerformIO $ do
-    r <- createVector (dim v-1)
-    app2 c_polySolve vec v vec r "polySolve"
-    return r
-             | otherwise = error "polySolve on a polynomial of degree zero"
-foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TV (TCV Res)
diff --git a/packages/hmatrix/src/Numeric/GSL/Random.hs b/packages/hmatrix/src/Numeric/GSL/Random.hs
deleted file mode 100644
index 2872b17..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Random.hs
+++ /dev/null
@@ -1,84 +0,0 @@
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.GSL.Random
-- Copyright   :  (c) Alberto Ruiz 2009-14
-- License     :  GPL
--
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-- Random vectors and matrices.
--
-----------------------------------------------------------------------------
-module Numeric.GSL.Random (
-    Seed,
-    RandDist(..),
-    randomVector,
-    gaussianSample,
-    uniformSample,
-    rand, randn
-) where
-import Numeric.GSL.Vector
-import Numeric.Container
-import Numeric.LinearAlgebra(Seed,RandDist(..),cholSH)
-import System.Random(randomIO)
-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
-- Gaussian distribution.
-gaussianSample :: Seed
-               -> Int -- ^ number of rows
-               -> Vector Double -- ^ mean vector
-               -> Matrix Double -- ^ covariance matrix
-               -> Matrix Double -- ^ result
-gaussianSample seed n med cov = m where
-    c = dim med
-    meds = konst' 1 n `outer` med
-    rs = reshape c $ randomVector seed Gaussian (c * n)
-    m = rs `mXm` cholSH cov `add` meds
-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
-- uniform distribution.
-uniformSample :: Seed
-               -> Int -- ^ number of rows
-               -> [(Double,Double)] -- ^ ranges for each column
-               -> Matrix Double -- ^ result
-uniformSample seed n rgs = m where
-    (as,bs) = unzip rgs
-    a = fromList as
-    cs = zipWith subtract as bs
-    d = dim a
-    dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
-    am = konst' 1 n `outer` a
-    m = fromColumns (zipWith scale cs dat) `add` am
-- | pseudorandom matrix with uniform elements between 0 and 1
-randm :: RandDist
-     -> Int -- ^ rows
-     -> Int -- ^ columns
-     -> IO (Matrix Double)
-randm d r c = do
-    seed <- randomIO
-    return (reshape c $ randomVector seed d (r*c))
-- | pseudorandom matrix with uniform elements between 0 and 1
-rand :: Int -> Int -> IO (Matrix Double)
-rand = randm Uniform
-{- | pseudorandom matrix with normal elements
->>> x <- randn 3 5
->>> disp 3 x
-3x5
-0.386  -1.141   0.491  -0.510   1.512
-0.069  -0.919   1.022  -0.181   0.745
-0.313  -0.670  -0.097  -1.575  -0.583
-}
-randn :: Int -> Int -> IO (Matrix Double)
-randn = randm Gaussian
diff --git a/packages/hmatrix/src/Numeric/GSL/Root.hs b/packages/hmatrix/src/Numeric/GSL/Root.hs
deleted file mode 100644
index b9f3b94..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Root.hs
+++ /dev/null
@@ -1,199 +0,0 @@
-{- |
-Module      :  Numeric.GSL.Root
-Copyright   :  (c) Alberto Ruiz 2009
-License     :  GPL
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-Multidimensional root finding.
-<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
-The example in the GSL manual:
->>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
->>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
->>> sol
-[1.0,1.0]
->>> disp 3 path
-11x5
- 1.000  -10.000  -5.000  11.000  -1050.000
- 2.000   -3.976  24.827   4.976     90.203
- 3.000   -3.976  24.827   4.976     90.203
- 4.000   -3.976  24.827   4.976     90.203
- 5.000   -1.274  -5.680   2.274    -73.018
- 6.000   -1.274  -5.680   2.274    -73.018
- 7.000    0.249   0.298   0.751      2.359
- 8.000    0.249   0.298   0.751      2.359
- 9.000    1.000   0.878  -0.000     -1.218
-10.000    1.000   0.989  -0.000     -0.108
-11.000    1.000   1.000   0.000      0.000
-}
-----------------------------------------------------------------------------
-module Numeric.GSL.Root (
-    uniRoot, UniRootMethod(..),
-    uniRootJ, UniRootMethodJ(..),
-    root, RootMethod(..),
-    rootJ, RootMethodJ(..),
-) where
-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
-                   | Brent
-                   deriving (Enum, Eq, Show, Bounded)
-uniRoot :: UniRootMethod
-        -> Double
-        -> Int
-        -> (Double -> Double)
-        -> Double
-        -> Double
-        -> (Double, Matrix Double)
-uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
-uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
-    fp <- mkDoublefun f
-    rawpath <- createMIO maxit 4
-                         (c_root m fp epsrel (fi maxit) xl xu)
-                         "root"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toLists $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    return (sol !! 1, path)
-foreign import ccall safe "root"
-    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM Res
-------------------------------------------------------------------------
-data UniRootMethodJ = UNewton
-                    | Secant
-                    | Steffenson
-                    deriving (Enum, Eq, Show, Bounded)
-uniRootJ :: UniRootMethodJ
-        -> Double
-        -> Int
-        -> (Double -> Double)
-        -> (Double -> Double)
-        -> Double
-        -> (Double, Matrix Double)
-uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun
-    dfun x epsrel maxit
-uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
-    fp <- mkDoublefun f
-    dfp <- mkDoublefun df
-    rawpath <- createMIO maxit 2
-                         (c_rootj m fp dfp epsrel (fi maxit) x)
-                         "rootj"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toLists $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    return (sol !! 1, path)
-foreign import ccall safe "rootj"
-    c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
-            -> Double -> CInt -> Double -> TM Res
-------------------------------------------------------------------------
-data RootMethod = Hybrids
-                | Hybrid
-                | DNewton
-                | Broyden
-                deriving (Enum,Eq,Show,Bounded)
-- | Nonlinear multidimensional root finding using algorithms that do not require 
-- any derivative information to be supplied by the user.
-- Any derivatives needed are approximated by finite differences.
-root :: RootMethod
-     -> Double                     -- ^ maximum residual
-     -> Int                        -- ^ maximum number of iterations allowed
-     -> ([Double] -> [Double])     -- ^ function to minimize
-     -> [Double]                   -- ^ starting point
-     -> ([Double], Matrix Double)  -- ^ solution vector and optimization path
-root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
-rootGen m f xi epsabs maxit = unsafePerformIO $ do
-    let xiv = fromList xi
-        n   = dim xiv
-    fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
-    rawpath <- vec xiv $ \xiv' ->
-                   createMIO maxit (2*n+1)
-                         (c_multiroot m fp epsabs (fi maxit) // xiv')
-                         "multiroot"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toLists $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    return (take n $ drop 1 sol, path)
-foreign import ccall safe "multiroot"
-    c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
-------------------------------------------------------------------------
-data RootMethodJ = HybridsJ
-                 | HybridJ
-                 | Newton
-                 | GNewton
-                deriving (Enum,Eq,Show,Bounded)
-- | Nonlinear multidimensional root finding using both the function and its derivatives.
-rootJ :: RootMethodJ
-      -> Double                     -- ^ maximum residual
-      -> Int                        -- ^ maximum number of iterations allowed
-      -> ([Double] -> [Double])     -- ^ function to minimize
-      -> ([Double] -> [[Double]])   -- ^ Jacobian
-      -> [Double]                   -- ^ starting point
-      -> ([Double], Matrix Double)  -- ^ solution vector and optimization path
-rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit
-rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
-    let xiv = fromList xi
-        n   = dim xiv
-    fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
-    jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
-    rawpath <- vec xiv $ \xiv' ->
-                   createMIO maxit (2*n+1)
-                         (c_multirootj m fp jp epsabs (fi maxit) // xiv')
-                         "multiroot"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toLists $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    freeHaskellFunPtr jp
-    return (take n $ drop 1 sol, path)
-foreign import ccall safe "multirootj"
-    c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
-------------------------------------------------------
-checkdim1 n v
-    | dim v == n = v
-    | otherwise = error $ "Error: "++ show n
-                        ++ " components expected in the result of the function supplied to root"
-checkdim2 n m
-    | rows m == n && cols m == n = m
-    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
-                        ++ " Jacobian expected in rootJ"
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
deleted file mode 100644
index af79f32..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ /dev/null
@@ -1,107 +0,0 @@
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.GSL.Vector
-- Copyright   :  (c) Alberto Ruiz 2007-14
-- License     :  GPL
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-----------------------------------------------------------------------------
-module Numeric.GSL.Vector (
-    randomVector,
-    saveMatrix,
-    fwriteVector, freadVector, fprintfVector, fscanfVector
-) where
-import Data.Packed
-import Numeric.LinearAlgebra(RandDist(..))
-import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
-import Foreign.Marshal.Alloc(free)
-import Foreign.Ptr(Ptr)
-import Foreign.C.Types
-import Foreign.C.String(newCString)
-import System.IO.Unsafe(unsafePerformIO)
-fromei x = fromIntegral (fromEnum x) :: CInt
-----------------------------------------------------------------------
-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
-randomVector :: Int      -- ^ seed
-             -> RandDist -- ^ distribution
-             -> Int      -- ^ vector size
-             -> Vector Double
-randomVector seed dist n = unsafePerformIO $ do
-    r <- createVector n
-    app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL"
-    return r
-foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: 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 PD = Ptr Double                            --
-type TV = CInt -> PD -> IO CInt                 --
-type TM = CInt -> CInt -> PD -> IO CInt         --
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
deleted file mode 100644
index ffc5c20..0000000
--- a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
+++ /dev/null
@@ -1,945 +0,0 @@
-#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
-#define FVEC(A) int A##n, float*A##p
-#define FMAT(A) int A##r, int A##c, float* A##p
-#define KFVEC(A) int A##n, const float*A##p
-#define KFMAT(A) int A##r, int A##c, const float* A##p
-#define QVEC(A) int A##n, gsl_complex_float*A##p
-#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p
-#define KQVEC(A) int A##n, const gsl_complex_float*A##p
-#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_fft_complex.h>
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_deriv.h>
-#include <gsl/gsl_poly.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_multiroots.h>
-#include <gsl/gsl_min.h>
-#include <gsl/gsl_complex_math.h>
-#include <gsl/gsl_rng.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_roots.h>
-#include <gsl/gsl_multifit_nlin.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 OK return 0;
-#define MIN(A,B) ((A)<(B)?(A):(B))
-#define MAX(A,B) ((A)>(B)?(A):(B))
-#ifdef DBG
-#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
-#else
-#define DEBUGMSG(M)
-#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 FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)
-#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)
-#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
-#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)
-#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)
-#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)
-#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)
-#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)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 GQVEC(A) int A##n, gsl_complex_float*A##p
-#define KGQVEC(A) int A##n, const gsl_complex_float*A##p
-#define BAD_SIZE 2000
-#define BAD_CODE 2001
-#define MEM      2002
-#define BAD_FILE 2003
-void no_abort_on_error() {
-    gsl_set_error_handler_off();
-}
-int fft(int code, KCVEC(X), CVEC(R)) {
-    REQUIRES(Xn == Rn,BAD_SIZE);
-    DEBUGMSG("fft");
-    int s = Xn;
-    gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
-    gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
-    gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn);
-    gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn);
-    gsl_blas_dcopy(&X.vector,&R.vector);
-    if(code==0) {
-        gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace);
-    } else {
-        gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace);
-    }
-    gsl_fft_complex_wavetable_free (wavetable);
-    gsl_fft_complex_workspace_free (workspace);
-    OK
-}
-int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
-{
-    gsl_function F;
-    F.function = f;
-    F.params = 0;
-    if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
-    if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
-    if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
-    return 0;
-}
-int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec,
-                   double *result, double*error) {
-    DEBUGMSG("integrate_qng");
-    gsl_function F;
-    F.function = f;
-    F.params = NULL;
-    size_t neval;
-    int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval);
-    CHECK(res,res);
-    OK
-}
-int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w,
-               double *result, double* error) {
-    DEBUGMSG("integrate_qags");
-    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
-    gsl_function F;
-    F.function = f;
-    F.params = NULL;
-    int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error);
-    CHECK(res,res);
-    gsl_integration_workspace_free (wk); 
-    OK
-}
-int integrate_qagi(double f(double,void*), double aprec, double prec, int w,
-               double *result, double* error) {
-    DEBUGMSG("integrate_qagi");
-    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
-    gsl_function F;
-    F.function = f;
-    F.params = NULL;
-    int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error);
-    CHECK(res,res);
-    gsl_integration_workspace_free (wk); 
-    OK
-}
-int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w,
-               double *result, double* error) {
-    DEBUGMSG("integrate_qagiu");
-    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
-    gsl_function F;
-    F.function = f;
-    F.params = NULL;
-    int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error);
-    CHECK(res,res);
-    gsl_integration_workspace_free (wk); 
-    OK
-}
-int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w,
-               double *result, double* error) {
-    DEBUGMSG("integrate_qagil");
-    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
-    gsl_function F;
-    F.function = f;
-    F.params = NULL;
-    int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error);
-    CHECK(res,res);
-    gsl_integration_workspace_free (wk); 
-    OK
-}
-int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec,
-                    int w, double *result, double* error, int *neval) {
-    DEBUGMSG("integrate_cquad");
-    gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w);
-    gsl_function F;
-    F.function = f;
-    F.params = NULL;
-    size_t * sneval = NULL;
-    int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval);
-    *neval = *sneval;
-    CHECK(res,res);
-    gsl_integration_cquad_workspace_free (wk); 
-    OK
-}
-int polySolve(KRVEC(a), CVEC(z)) {
-    DEBUGMSG("polySolve");
-    REQUIRES(an>1,BAD_SIZE);
-    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an);
-    int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp);
-    CHECK(res,res);
-    gsl_poly_complex_workspace_free (w);
-    OK;
-}
-int vector_fscanf(char*filename, RVEC(a)) {
-    DEBUGMSG("gsl_vector_fscanf");
-    DVVIEW(a);
-    FILE * f = fopen(filename,"r");
-    CHECK(!f,BAD_FILE);
-    int res = gsl_vector_fscanf(f,V(a));
-    CHECK(res,res);
-    fclose (f);
-    OK
-}
-int vector_fprintf(char*filename, char*fmt, RVEC(a)) {
-    DEBUGMSG("gsl_vector_fprintf");
-    DVVIEW(a);
-    FILE * f = fopen(filename,"w");
-    CHECK(!f,BAD_FILE);
-    int res = gsl_vector_fprintf(f,V(a),fmt);
-    CHECK(res,res);
-    fclose (f);
-    OK
-}
-int vector_fread(char*filename, RVEC(a)) {
-    DEBUGMSG("gsl_vector_fread");
-    DVVIEW(a);
-    FILE * f = fopen(filename,"r");
-    CHECK(!f,BAD_FILE);
-    int res = gsl_vector_fread(f,V(a));
-    CHECK(res,res);
-    fclose (f);
-    OK
-}
-int vector_fwrite(char*filename, RVEC(a)) {
-    DEBUGMSG("gsl_vector_fwrite");
-    DVVIEW(a);
-    FILE * f = fopen(filename,"w");
-    CHECK(!f,BAD_FILE);
-    int res = gsl_vector_fwrite(f,V(a));
-    CHECK(res,res);
-    fclose (f);
-    OK
-}
-int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) {
-    DEBUGMSG("matrix_fprintf");
-    FILE * f = fopen(filename,"w");
-    CHECK(!f,BAD_FILE);
-    int i,j,sr,sc;
-    if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
-    #define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
-    for (i=0; i<mr; i++) {
-        for (j=0; j<mc-1; j++) {
-            fprintf(f,fmt,AT(m,i,j));
-            fprintf(f," ");
-        }
-        fprintf(f,fmt,AT(m,i,j));
-        fprintf(f,"\n");
-    }
-    fclose (f);
-    OK
-}
-//---------------------------------------------------------------
-typedef double Trawfun(int, double*);
-double only_f_aux_min(const gsl_vector*x, void *pars) {
-    Trawfun * f = (Trawfun*) pars;  
-    double* p = (double*)calloc(x->size,sizeof(double));
-    int k;
-    for(k=0;k<x->size;k++) {
-        p[k] = gsl_vector_get(x,k);
-    }  
-    double res = f(x->size,p);
-    free(p);
-    return res;
-}
-double only_f_aux_root(double x, void *pars);
-int uniMinimize(int method, double f(double),
-                double epsrel, int maxit, double min,
-                double xl, double xu, RMAT(sol)) {
-   REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
-   DEBUGMSG("minimize_only_f");
-   gsl_function my_func;
-   my_func.function = only_f_aux_root;
-   my_func.params = f;
-   size_t iter = 0;
-   int status;
-   const gsl_min_fminimizer_type *T;
-   gsl_min_fminimizer *s;
-   // Starting point
-   switch(method) {
-     case 0 : {T = gsl_min_fminimizer_goldensection; break; }
-     case 1 : {T = gsl_min_fminimizer_brent; break; }
-     case 2 : {T = gsl_min_fminimizer_quad_golden; break; }
-     default: ERROR(BAD_CODE);
-   }
-   s = gsl_min_fminimizer_alloc (T);
-   gsl_min_fminimizer_set (s, &my_func, min, xl, xu);
-   do {
-       double current_min, current_lo, current_hi;
-       status = gsl_min_fminimizer_iterate (s);
-       current_min = gsl_min_fminimizer_x_minimum (s);
-       current_lo = gsl_min_fminimizer_x_lower (s);
-       current_hi = gsl_min_fminimizer_x_upper (s);
-       solp[iter*solc] = iter + 1;
-       solp[iter*solc+1] = current_min;
-       solp[iter*solc+2] = current_lo;
-       solp[iter*solc+3] = current_hi;
-       iter++;
-       if (status)   /* check if solver is stuck */
-          break;
-       
-       status =
-         gsl_min_test_interval (current_lo, current_hi, 0, epsrel);
-   }
-   while (status == GSL_CONTINUE && iter < maxit);
-   int i;
-   for (i=iter; i<solr; i++) {
-       solp[i*solc+0] = iter;
-       solp[i*solc+1]=0.;
-       solp[i*solc+2]=0.;
-       solp[i*solc+3]=0.;
-   }
-   gsl_min_fminimizer_free(s);
-   OK
-}
-   
-// this version returns info about intermediate steps
-int minimize(int method, double f(int, double*), double tolsize, int maxit, 
-                 KRVEC(xi), KRVEC(sz), RMAT(sol)) {
-    REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE);
-    DEBUGMSG("minimizeList (nmsimplex)");
-    gsl_multimin_function my_func;
-    // extract function from pars
-    my_func.f = only_f_aux_min;
-    my_func.n = xin; 
-    my_func.params = f;
-    size_t iter = 0;
-    int status;
-    double size;
-    const gsl_multimin_fminimizer_type *T;
-    gsl_multimin_fminimizer *s = NULL;
-    // Initial vertex size vector 
-    KDVVIEW(sz);
-    // Starting point
-    KDVVIEW(xi);
-    // Minimizer nmsimplex, without derivatives
-    switch(method) {
-        case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
-#ifdef GSL110
-        case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
-#else
-        case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; }
-#endif
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_multimin_fminimizer_alloc (T, my_func.n);
-    gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz));
-    do {
-        status = gsl_multimin_fminimizer_iterate (s);
-        size = gsl_multimin_fminimizer_size (s);
-        solp[iter*solc+0] = iter+1;
-        solp[iter*solc+1] = s->fval;
-        solp[iter*solc+2] = size;
-        int k;
-        for(k=0;k<xin;k++) {
-            solp[iter*solc+k+3] = gsl_vector_get(s->x,k);
-        }
-        iter++;
-        if (status) break;
-        status = gsl_multimin_test_size (size, tolsize);
-    } while (status == GSL_CONTINUE && iter < maxit);
-    int i,j;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        for(j=1;j<solc;j++) {
-            solp[i*solc+j]=0.;
-        }
-    }
-    gsl_multimin_fminimizer_free(s);
-    OK
-}
-// working with the gradient
-typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;
-double f_aux_min(const gsl_vector*x, void *pars) {
-    Tfdf * fdf = ((Tfdf*) pars);
-    double* p = (double*)calloc(x->size,sizeof(double));
-    int k;
-    for(k=0;k<x->size;k++) {
-        p[k] = gsl_vector_get(x,k);
-    }
-    double res = fdf->f(x->size,p);
-    free(p);
-    return res;
-}
-void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
-    Tfdf * fdf = ((Tfdf*) pars);
-    double* p = (double*)calloc(x->size,sizeof(double));
-    double* q = (double*)calloc(g->size,sizeof(double));
-    int k;
-    for(k=0;k<x->size;k++) {
-        p[k] = gsl_vector_get(x,k);
-    }
-    fdf->df(x->size,p,g->size,q);
-    for(k=0;k<x->size;k++) {
-        gsl_vector_set(g,k,q[k]);
-    }
-    free(p);
-    free(q);
-}
-void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) {
-    *f = f_aux_min(x,pars);
-    df_aux_min(x,pars,g);
-}
-int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*),
-              double initstep, double minimpar, double tolgrad, int maxit, 
-              KRVEC(xi), RMAT(sol)) {
-    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
-    DEBUGMSG("minimizeWithDeriv (conjugate_fr)");
-    gsl_multimin_function_fdf my_func;
-    // extract function from pars
-    my_func.f = f_aux_min;
-    my_func.df = df_aux_min;
-    my_func.fdf = fdf_aux_min;
-    my_func.n = xin; 
-    Tfdf stfdf;
-    stfdf.f = f;
-    stfdf.df = df;
-    my_func.params = &stfdf;
-    size_t iter = 0;
-    int status;
-    const gsl_multimin_fdfminimizer_type *T;
-    gsl_multimin_fdfminimizer *s = NULL;
-    // Starting point
-    KDVVIEW(xi);
-    // conjugate gradient fr
-    switch(method) {
-        case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }
-        case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; }
-        case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; }
-        case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }
-        case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; }
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
-    gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
-    do {
-        status = gsl_multimin_fdfminimizer_iterate (s);
-        solp[iter*solc+0] = iter+1;
-        solp[iter*solc+1] = s->f;
-        int k;
-        for(k=0;k<xin;k++) {
-            solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
-        }
-        iter++;
-        if (status) break;
-        status = gsl_multimin_test_gradient (s->gradient, tolgrad);
-    } while (status == GSL_CONTINUE && iter < maxit);
-    int i,j;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        for(j=1;j<solc;j++) {
-            solp[i*solc+j]=0.;
-        }
-    }
-    gsl_multimin_fdfminimizer_free(s);
-    OK
-}
-//---------------------------------------------------------------
-double only_f_aux_root(double x, void *pars) {
-    double (*f)(double) = (double (*)(double)) pars;
-    return f(x);
-}
-int root(int method, double f(double),
-         double epsrel, int maxit,
-         double xl, double xu, RMAT(sol)) {
-    REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
-    DEBUGMSG("root_only_f");
-    gsl_function my_func;
-    // extract function from pars
-    my_func.function = only_f_aux_root;
-    my_func.params = f;
-    size_t iter = 0;
-    int status;
-    const gsl_root_fsolver_type *T;
-    gsl_root_fsolver *s;
-    // Starting point
-    switch(method) {
-        case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
-        case 1 : {T = gsl_root_fsolver_falsepos; break; }
-        case 2 : {T = gsl_root_fsolver_brent; break; }
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_root_fsolver_alloc (T);
-    gsl_root_fsolver_set (s, &my_func, xl, xu);
-    do {
-           double best, current_lo, current_hi;
-           status = gsl_root_fsolver_iterate (s);
-           best = gsl_root_fsolver_root (s);
-           current_lo = gsl_root_fsolver_x_lower (s);
-           current_hi = gsl_root_fsolver_x_upper (s);
-           solp[iter*solc] = iter + 1;
-           solp[iter*solc+1] = best;
-           solp[iter*solc+2] = current_lo;
-           solp[iter*solc+3] = current_hi;
-           iter++;
-           if (status)   /* check if solver is stuck */
-             break;
-           status =
-               gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
-        }
-        while (status == GSL_CONTINUE && iter < maxit);
-    int i;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        solp[i*solc+1]=0.;
-        solp[i*solc+2]=0.;
-        solp[i*solc+3]=0.;
-    }
-    gsl_root_fsolver_free(s);
-    OK
-}
-typedef struct {
-    double (*f)(double);
-    double (*jf)(double);
-} uniTfjf;
-double f_aux_uni(double x, void *pars) {
-    uniTfjf * fjf = ((uniTfjf*) pars);
-    return (fjf->f)(x);
-}
-double jf_aux_uni(double x, void * pars) {
-    uniTfjf * fjf = ((uniTfjf*) pars);
-    return (fjf->jf)(x);
-}
-void fjf_aux_uni(double x, void * pars, double * f, double * g) {
-    *f = f_aux_uni(x,pars);
-    *g = jf_aux_uni(x,pars);
-}
-int rootj(int method, double f(double),
-          double df(double),
-         double epsrel, int maxit,
-         double x, RMAT(sol)) {
-    REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
-    DEBUGMSG("root_fjf");
-    gsl_function_fdf my_func;
-    // extract function from pars
-    my_func.f = f_aux_uni;
-    my_func.df = jf_aux_uni;
-    my_func.fdf = fjf_aux_uni;
-    uniTfjf stfjf;
-    stfjf.f = f;
-    stfjf.jf = df;
-    my_func.params = &stfjf;
-    size_t iter = 0;
-    int status;
-    const gsl_root_fdfsolver_type *T;
-    gsl_root_fdfsolver *s;
-    // Starting point
-    switch(method) {
-        case 0 : {T = gsl_root_fdfsolver_newton;; break; }
-        case 1 : {T = gsl_root_fdfsolver_secant; break; }
-        case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_root_fdfsolver_alloc (T);
-    gsl_root_fdfsolver_set (s, &my_func, x);
-    do {
-           double x0;
-           status = gsl_root_fdfsolver_iterate (s);
-           x0 = x;
-           x = gsl_root_fdfsolver_root(s);
-           solp[iter*solc+0] = iter+1;
-           solp[iter*solc+1] = x;
-           iter++;
-           if (status)   /* check if solver is stuck */
-             break;
-           status =
-               gsl_root_test_delta (x, x0, 0, epsrel);
-        }
-        while (status == GSL_CONTINUE && iter < maxit);
-    int i;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        solp[i*solc+1]=0.;
-    }
-    gsl_root_fdfsolver_free(s);
-    OK
-}
-//---------------------------------------------------------------
-typedef void TrawfunV(int, double*, int, double*);
-int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) {
-    TrawfunV * f = (TrawfunV*) pars;
-    double* p = (double*)calloc(x->size,sizeof(double));
-    double* q = (double*)calloc(y->size,sizeof(double));
-    int k;
-    for(k=0;k<x->size;k++) {
-        p[k] = gsl_vector_get(x,k);
-    }
-    f(x->size,p,y->size,q);
-    for(k=0;k<y->size;k++) {
-        gsl_vector_set(y,k,q[k]);
-    }
-    free(p);
-    free(q);
-    return 0; //hmmm
-}
-int multiroot(int method, void f(int, double*, int, double*),
-         double epsabs, int maxit,
-         KRVEC(xi), RMAT(sol)) {
-    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
-    DEBUGMSG("root_only_f");
-    gsl_multiroot_function my_func;
-    // extract function from pars
-    my_func.f = only_f_aux_multiroot;
-    my_func.n = xin;
-    my_func.params = f;
-    size_t iter = 0;
-    int status;
-    const gsl_multiroot_fsolver_type *T;
-    gsl_multiroot_fsolver *s;
-    // Starting point
-    KDVVIEW(xi);
-    switch(method) {
-        case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; }
-        case 1 : {T = gsl_multiroot_fsolver_hybrid; break; }
-        case 2 : {T = gsl_multiroot_fsolver_dnewton; break; }
-        case 3 : {T = gsl_multiroot_fsolver_broyden; break; }
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_multiroot_fsolver_alloc (T, my_func.n);
-    gsl_multiroot_fsolver_set (s, &my_func, V(xi));
-    do {
-           status = gsl_multiroot_fsolver_iterate (s);
-           solp[iter*solc+0] = iter+1;
-           int k;
-           for(k=0;k<xin;k++) {
-               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
-           }
-           for(k=xin;k<2*xin;k++) {
-               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
-           }
-           iter++;
-           if (status)   /* check if solver is stuck */
-             break;
-           status =
-             gsl_multiroot_test_residual (s->f, epsabs);
-        }
-        while (status == GSL_CONTINUE && iter < maxit);
-    int i,j;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        for(j=1;j<solc;j++) {
-            solp[i*solc+j]=0.;
-        }
-    }
-    gsl_multiroot_fsolver_free(s);
-    OK
-}
-// working with the jacobian
-typedef struct {int (*f)(int, double*, int, double *);
-                int (*jf)(int, double*, int, int, double*);} Tfjf;
-int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) {
-    Tfjf * fjf = ((Tfjf*) pars);
-    double* p = (double*)calloc(x->size,sizeof(double));
-    double* q = (double*)calloc(y->size,sizeof(double));
-    int k;
-    for(k=0;k<x->size;k++) {
-        p[k] = gsl_vector_get(x,k);
-    }
-    (fjf->f)(x->size,p,y->size,q);
-    for(k=0;k<y->size;k++) {
-        gsl_vector_set(y,k,q[k]);
-    }
-    free(p);
-    free(q);
-    return 0;
-}
-int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) {
-    Tfjf * fjf = ((Tfjf*) pars);
-    double* p = (double*)calloc(x->size,sizeof(double));
-    double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double));
-    int i,j,k;
-    for(k=0;k<x->size;k++) {
-        p[k] = gsl_vector_get(x,k);
-    }
-    (fjf->jf)(x->size,p,jac->size1,jac->size2,q);
-    k=0;
-    for(i=0;i<jac->size1;i++) {
-        for(j=0;j<jac->size2;j++){
-            gsl_matrix_set(jac,i,j,q[k++]);
-        }
-    }
-    free(p);
-    free(q);
-    return 0;
-}
-int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
-    f_aux(x,pars,f);
-    jf_aux(x,pars,g);
-    return 0;
-}
-int multirootj(int method, int f(int, double*, int, double*),
-                      int jac(int, double*, int, int, double*),
-         double epsabs, int maxit,
-         KRVEC(xi), RMAT(sol)) {
-    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
-    DEBUGMSG("root_fjf");
-    gsl_multiroot_function_fdf my_func;
-    // extract function from pars
-    my_func.f = f_aux;
-    my_func.df = jf_aux;
-    my_func.fdf = fjf_aux;
-    my_func.n = xin;
-    Tfjf stfjf;
-    stfjf.f = f;
-    stfjf.jf = jac;
-    my_func.params = &stfjf;
-    size_t iter = 0;
-    int status;
-    const gsl_multiroot_fdfsolver_type *T;
-    gsl_multiroot_fdfsolver *s;
-    // Starting point
-    KDVVIEW(xi);
-    switch(method) {
-        case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; }
-        case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; }
-        case 2 : {T = gsl_multiroot_fdfsolver_newton; break; }
-        case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; }
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_multiroot_fdfsolver_alloc (T, my_func.n);
-    gsl_multiroot_fdfsolver_set (s, &my_func, V(xi));
-    do {
-           status = gsl_multiroot_fdfsolver_iterate (s);
-           solp[iter*solc+0] = iter+1;
-           int k;
-           for(k=0;k<xin;k++) {
-               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
-           }
-           for(k=xin;k<2*xin;k++) {
-               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
-           }
-           iter++;
-           if (status)   /* check if solver is stuck */
-             break;
-           status =
-             gsl_multiroot_test_residual (s->f, epsabs);
-        }
-        while (status == GSL_CONTINUE && iter < maxit);
-    int i,j;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        for(j=1;j<solc;j++) {
-            solp[i*solc+j]=0.;
-        }
-    }
-    gsl_multiroot_fdfsolver_free(s);
-    OK
-}
-//-------------- non linear least squares fitting -------------------
-int nlfit(int method, int f(int, double*, int, double*),
-                      int jac(int, double*, int, int, double*),
-         double epsabs, double epsrel, int maxit, int p,
-         KRVEC(xi), RMAT(sol)) {
-    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
-    DEBUGMSG("nlfit");
-    const gsl_multifit_fdfsolver_type *T;
-    gsl_multifit_fdfsolver *s;
-    gsl_multifit_function_fdf my_f;
-    // extract function from pars
-    my_f.f = f_aux;
-    my_f.df = jf_aux;
-    my_f.fdf = fjf_aux;
-    my_f.n = p;
-    my_f.p = xin;  // !!!!
-    Tfjf stfjf;
-    stfjf.f = f;
-    stfjf.jf = jac;
-    my_f.params = &stfjf;
-    size_t iter = 0;
-    int status;
-    KDVVIEW(xi);
-    //DMVIEW(cov);
-    switch(method) {
-        case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; }
-        case 1 : { T = gsl_multifit_fdfsolver_lmder; break; }
-        default: ERROR(BAD_CODE);
-    }
-    s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p);
-    gsl_multifit_fdfsolver_set (s, &my_f, V(xi));
-    do {   status = gsl_multifit_fdfsolver_iterate (s);
-           solp[iter*solc+0] = iter+1;
-           solp[iter*solc+1] = gsl_blas_dnrm2 (s->f);
-           int k;
-           for(k=0;k<xin;k++) {
-               solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
-           }
-           iter++;
-           if (status)   /* check if solver is stuck */
-             break;
-           status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
-        }
-        while (status == GSL_CONTINUE && iter < maxit);
-    int i,j;
-    for (i=iter; i<solr; i++) {
-        solp[i*solc+0] = iter;
-        for(j=1;j<solc;j++) {
-            solp[i*solc+j]=0.;
-        }
-    }
-    //gsl_multifit_covar (s->J, 0.0, M(cov));
-    gsl_multifit_fdfsolver_free (s);
-    OK
-}
-//////////////////////////////////////////////////////
-#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
-int random_vector_GSL(int seed, int code, RVEC(r)) {
-    DEBUGMSG("random_vector_GSL")
-    static gsl_rng * gen = NULL;
-    if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
-    gsl_rng_set (gen, seed);
-    int k;
-    switch (code) {
-        RAN(0,gsl_rng_uniform)
-        RAN(1,gsl_ran_ugaussian)
-        default: ERROR(BAD_CODE);
-    }
-}
-#undef RAN
-//////////////////////////////////////////////////////
-#include "gsl-ode.c"
-//////////////////////////////////////////////////////
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
deleted file mode 100644
index 3f2771b..0000000
--- a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
+++ /dev/null
@@ -1,182 +0,0 @@
-#ifdef GSLODE1
-////////////////////////////// ODE V1 //////////////////////////////////////////
-#include <gsl/gsl_odeiv.h>
-typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
-int odefunc (double t, const double y[], double f[], void *params) { 
-    Tode * P = (Tode*) params;
-    (P->f)(t,P->n,y,P->n,f);
-    return GSL_SUCCESS;
-}
-int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
-     Tode * P = ((Tode*) params);
-     (P->j)(t,P->n,y,P->n,P->n,dfdy);
-     int j;
-     for (j=0; j< P->n; j++)
-        dfdt[j] = 0.0;
-     return GSL_SUCCESS;
-}
-int ode(int method, double h, double eps_abs, double eps_rel,
-        int f(double, int, const double*, int, double*),
-        int jac(double, int, const double*, int, int, double*),
-        KRVEC(xi), KRVEC(ts), RMAT(sol)) {
-    const gsl_odeiv_step_type * T;
-    switch(method) {
-        case 0 : {T = gsl_odeiv_step_rk2; break; }
-        case 1 : {T = gsl_odeiv_step_rk4; break; }
-        case 2 : {T = gsl_odeiv_step_rkf45; break; }
-        case 3 : {T = gsl_odeiv_step_rkck; break; }
-        case 4 : {T = gsl_odeiv_step_rk8pd; break; }
-        case 5 : {T = gsl_odeiv_step_rk2imp; break; }
-        case 6 : {T = gsl_odeiv_step_rk4imp; break; }
-        case 7 : {T = gsl_odeiv_step_bsimp; break; }
-        case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); }
-        case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); }
-        case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); }
-        default: ERROR(BAD_CODE);
-    }
-    gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
-    gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
-    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
-    Tode P;
-    P.f = f;
-    P.j = jac;
-    P.n = xin;
-    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
-    double t = tsp[0];
-    double* y = (double*)calloc(xin,sizeof(double));
-    int i,j;
-    for(i=0; i< xin; i++) {
-        y[i] = xip[i];
-        solp[i] = xip[i];
-    }
-       for (i = 1; i < tsn ; i++)
-         {
-           double ti = tsp[i];
-           while (t < ti)
-             {
-               gsl_odeiv_evolve_apply (e, c, s,
-                                       &sys,
-                                       &t, ti, &h,
-                                       y);
-               // if (h < hmin) h = hmin;
-             }
-           for(j=0; j<xin; j++) {
-               solp[i*xin + j] = y[j];
-           }
-         }
-    free(y);
-    gsl_odeiv_evolve_free (e);
-    gsl_odeiv_control_free (c);
-    gsl_odeiv_step_free (s);
-    return 0;
-}
-#else
-///////////////////// ODE V2 ///////////////////////////////////////////////////
-                   
-#include <gsl/gsl_odeiv2.h>
-typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
-int odefunc (double t, const double y[], double f[], void *params) { 
-    Tode * P = (Tode*) params;
-    (P->f)(t,P->n,y,P->n,f);
-    return GSL_SUCCESS;
-}
-int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
-     Tode * P = ((Tode*) params);
-     (P->j)(t,P->n,y,P->n,P->n,dfdy);
-     int j;
-     for (j=0; j< P->n; j++)
-        dfdt[j] = 0.0;
-     return GSL_SUCCESS;
-}
-int ode(int method, double h, double eps_abs, double eps_rel,
-        int f(double, int, const double*, int, double*),
-        int jac(double, int, const double*, int, int, double*),
-        KRVEC(xi), KRVEC(ts), RMAT(sol)) {
-    const gsl_odeiv2_step_type * T;
-    switch(method) {
-        case 0 : {T = gsl_odeiv2_step_rk2; break; }
-        case 1 : {T = gsl_odeiv2_step_rk4; break; }
-        case 2 : {T = gsl_odeiv2_step_rkf45; break; }
-        case 3 : {T = gsl_odeiv2_step_rkck; break; }
-        case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
-        case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
-        case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
-        case 7 : {T = gsl_odeiv2_step_bsimp; break; }
-        case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
-        case 9 : {T = gsl_odeiv2_step_msadams; break; }
-        case 10: {T = gsl_odeiv2_step_msbdf; break; }
-        default: ERROR(BAD_CODE);
-    }
-    Tode P;
-    P.f = f;
-    P.j = jac;
-    P.n = xin;
-    gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
-    gsl_odeiv2_driver * d =
-         gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
-    double t = tsp[0];
-    double* y = (double*)calloc(xin,sizeof(double));
-    int i,j;
-    int status=0;
-    for(i=0; i< xin; i++) {
-        y[i] = xip[i];
-        solp[i] = xip[i];
-    }
-       for (i = 1; i < tsn ; i++)
-         {
-           double ti = tsp[i];
-           
-           status = gsl_odeiv2_driver_apply (d, &t, ti, y);
-     
-           if (status != GSL_SUCCESS) {
-                  printf ("error in ode, return value=%d\n", status);
-                  break;
-                }
-//           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
-           
-           for(j=0; j<xin; j++) {
-               solp[i*xin + j] = y[j];
-           }
-         }
-    free(y);
-    gsl_odeiv2_driver_free (d);
-    
-    return status;
-}
-#endif
author	Alberto Ruiz <aruiz@um.es>	2014-05-21 10:30:55 +0200
committer	Alberto Ruiz <aruiz@um.es>	2014-05-21 10:30:55 +0200
commit	197e88c3b56d28840217010a2871c6ea3a4dd1a4 (patch)
tree	825be9d6c9d87d23f7e5497c0133d11d52c63535 /packages/hmatrix
parent	e07c3dee7235496b71a89233106d93f6cc94ada1 (diff)