10 files changed, 79 insertions, 22 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index 43569fe..6905a6a 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -68,7 +68,7 @@ Exposed-modules:    Data.Packed.Internal,
                     Numeric.LinearAlgebra.Interface,
                     Numeric.LinearAlgebra.Algorithms,
                     Graphics.Plot
-C-sources:          lib/Data/Packed/Internal/aux.c,
+C-sources:          lib/Data/Packed/Internal/auxi.c,
                     lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
                     lib/Numeric/GSL/gsl-aux.c
 extra-libraries:    gsl cblas lapack
diff --git a/README b/README
index 00188ce..cd54fcd 100644
--- a/README
+++ b/README
@@ -46,7 +46,29 @@ Prelude Numeric.LinearAlgebra> u <> d <> trans v
 Prelude Numeric.GSL> :q
 Leaving GHCi.
 
---------------------------------------------------------------------------------------
+CHANGES
+
+This is a new version of GSLHaskell. The package is (provisionally)
+called \C{hssl} (a simple scientific library for Haskell) because only
+a small part of GSL is available and linear algebra is based on LAPACK.
+
+The code has been extensively refactored. There is a new internal representation
+which admits both C and Fortran matrices and avoids many transposes.
+
+There are only minor API changes:
+
+- the matrix product operator \C{(<>)} is now overloaded only for matrix-matrix,
+matrix-vector and vector-matrix, with the same base type. The dot product and the scaling
+of vectors or matrices is now denoted by `dot` and `scale`. Conversions from real to
+complex objects must be explicit.
+
+- Most linear algebra functions admit both real and complex objects. Utilities such as
+ident or constant are now polymorphic.
+
+- Runtime errors produced by GSL or LAPACK can be handled using \C{Control.Exeception.catch}.
+
+Old GSLHaskell code will work with small modifications.
+
 ACKNOWLEDGEMENTS
 
 I thank Henning Thielemann and all the people in the Haskell mailing lists for their help.
@@ -63,3 +85,26 @@ I thank Henning Thielemann and all the people in the Haskell mailing lists for t
 - Pedro E. López de Teruel fixed the interface to lapack.
 
 - Antti Siira discovered a bug in the plotting functions.
+
+INSTALLATION ON WINDOWS
+
+1) Download the developer files gsl-1.8-lib.zip from
+   http://gnuwin32.sourceforge.net/packages/gsl.htm
+   and copy the gsl folder (under include) to the include folder of ghc:
+   C:\ghc\ghc.6.6.1\include
+
+2) Install the package as usual from the command line in the hssl-0.1 folder:
+       runhaskell Setup.lhs configure
+       runhaskell Setup.lhs build
+       runhaskell Setup.lhs install
+
+3) Copy libgsl.dll, libcblas.dll (from the binaries package gsl-1.8.bin.zip)
+   and liblapack.dll (borrowed from the R system) to the folder where 
+   hssl has been installed: C:\haskell\hss-0.1\ghc-6.6.1. They are needed to compile programs.
+   These three dlls are also available from http://perception.inf.um.es/darcs/HSSL/dll1.zip
+
+4) Copy the dlls at http://perception.inf.um.es/darcs/HSSL/dll2.zip to C:\windows\system
+   They are required to run the programs and ghci.
+
+Unfortunately the lapack dll supplied by the R system does not include zgels_ and zgelss_,
+so the functions depending on them will produced a "non supported" runtime error.
diff --git a/examples/parallel.hs b/examples/parallel.hs
index 8f67863..7256fb6 100644
--- a/examples/parallel.hs
+++ b/examples/parallel.hs
@@ -5,11 +5,9 @@ import System.Time
 
 inParallel = parMap rwhnf id
 
-parMul x y = fromBlocks [[ay],[by]]
-    where p = rows x `div` 2
-          a = takeRows p x
-          b = dropRows p x
-          [ay,by] = inParallel [a<>y,b<>y]
+parMul x y = fromBlocks [inParallel[x <> y1, x <> y2]]
+    where p = cols y `div` 2
+          (y1,y2) = splitColumnsAt p y
 
 main = do
     n <- (read . head) `fmap` getArgs
@@ -20,6 +18,9 @@ main = do
 a = (2><3) [1..6::Double]
 b = (3><4) [1..12::Double]
 
+splitRowsAt p m    = (takeRows p m, dropRows p m)
+splitColumnsAt p m = (takeColumns p m, dropColumns p m)
+
 time act = do
     t0 <- getClockTime
     act
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index 75a5b1e..2b3ec28 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -70,6 +70,7 @@ errorCode 2003 = "bad file"
 errorCode 2004 = "singular"
 errorCode 2005 = "didn't converge"
 errorCode 2006 = "the input matrix is not positive definite"
+errorCode 2007 = "not yet supported in this OS"
 errorCode n    = "code "++show n
 
 {- | conversion of Haskell functions into function pointers that can be used in the C side
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index e76500b..bf7f0ec 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -240,9 +240,9 @@ transdataAux fun c1 d c2 =
         r2 = dim d `div` c2
         noneed = r1 == 1 || c1 == 1
 
-foreign import ccall safe "aux.h transR"
+foreign import ccall safe "auxi.h transR"
     ctransR :: TMM -- Double ::> Double ::> IO Int
-foreign import ccall safe "aux.h transC"
+foreign import ccall safe "auxi.h transC"
     ctransC :: TCMCM -- Complex Double ::> Complex Double ::> IO Int
 
 ------------------------------------------------------------------
@@ -258,14 +258,14 @@ multiplyAux fun a b = unsafePerformIO $ do
     return r
 
 multiplyR = multiplyAux cmultiplyR
-foreign import ccall safe "aux.h multiplyR"
+foreign import ccall safe "auxi.h multiplyR"
     cmultiplyR :: Int -> Int -> Int -> Ptr Double
                -> Int -> Int -> Int -> Ptr Double
                -> Int -> Int -> Ptr Double
                -> IO Int
 
 multiplyC = multiplyAux cmultiplyC
-foreign import ccall safe "aux.h multiplyC"
+foreign import ccall safe "auxi.h multiplyC"
     cmultiplyC :: Int -> Int -> Int -> Ptr (Complex Double)
                -> Int -> Int -> Int -> Ptr (Complex Double)
                -> Int -> Int -> Ptr (Complex Double)
@@ -288,7 +288,7 @@ subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do
     r <- createMatrix RowMajor rt ct
     c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r]
     return r
-foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
+foreign import ccall "auxi.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
 
 -- | extraction of a submatrix from a complex matrix
 subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double)
@@ -316,12 +316,12 @@ diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
 -- | diagonal matrix from a real vector
 diagR :: Vector Double -> Matrix Double
 diagR = diagAux c_diagR "diagR"
-foreign import ccall "aux.h diagR" c_diagR :: TVM
+foreign import ccall "auxi.h diagR" c_diagR :: TVM
 
 -- | diagonal matrix from a real vector
 diagC :: Vector (Complex Double) -> Matrix (Complex Double)
 diagC = diagAux c_diagC "diagC"
-foreign import ccall "aux.h diagC" c_diagC :: TCVCM
+foreign import ccall "auxi.h diagC" c_diagC :: TCVCM
 
 -- | creates a square matrix with the given diagonal
 diag :: Field a => Vector a -> Matrix a
@@ -338,12 +338,12 @@ constantAux fun x n = unsafePerformIO $ do
 
 constantR :: Double -> Int -> Vector Double
 constantR = constantAux cconstantR
-foreign import ccall safe "aux.h constantR"
+foreign import ccall safe "auxi.h constantR"
     cconstantR :: Ptr Double -> TV -- Double :> IO Int
 
 constantC :: Complex Double -> Int -> Vector (Complex Double)
 constantC = constantAux cconstantC
-foreign import ccall safe "aux.h constantC"
+foreign import ccall safe "auxi.h constantC"
     cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int
 
 {- | creates a vector with a given number of equal components:
@@ -381,7 +381,7 @@ fromFile filename (r,c) = do
     c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
     --free charname  -- TO DO: free the auxiliary CString
     return res
-foreign import ccall "aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
+foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
 
 -------------------------------------------------------------------------
 
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index ebe6371..9557206 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -46,7 +46,7 @@ check msg ls f = do
     return ()
 
 -- | description of GSL error codes
-foreign import ccall "aux.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar)
+foreign import ccall "auxi.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar)
 
 -- | signature of foreign functions admitting C-style vectors
 type Vc t s = Int -> Ptr t -> s
diff --git a/lib/Data/Packed/Internal/aux.c b/lib/Data/Packed/Internal/auxi.c
index 3db3535..b53d9b7 100644
--- a/lib/Data/Packed/Internal/aux.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -1,4 +1,4 @@
-#include "aux.h"
+#include "auxi.h"
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/auxi.h
index 73334e3..73334e3 100644
--- a/lib/Data/Packed/Internal/aux.h
+++ b/lib/Data/Packed/Internal/auxi.h
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 83c3bf8..bd0a6bd 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -116,9 +116,9 @@ int mapR(int code, KRVEC(x), RVEC(r)) {
         OP(7,sinh)
         OP(8,cosh)
         OP(9,tanh)
-        OP(10,asinh)
-        OP(11,acosh)
-        OP(12,atanh)
+        OP(10,gsl_asinh)
+        OP(11,gsl_acosh)
+        OP(12,gsl_atanh)
         OP(13,exp)
         OP(14,log)
         OP(15,sign)
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index bf0cfba..9b6c1db 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -29,6 +29,8 @@
 #define SINGULAR 2004
 #define NOCONVER 2005
 #define NODEFPOS 2006
+#define NOSPRTD  2007
+
 
 //////////////////// real svd ////////////////////////////////////
 
@@ -423,6 +425,9 @@ int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
 //////////////////// least squares complex linear system ////////////
 
 int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+    #ifdef _WIN32
+    return NOSPRTD;
+    #else
     integer m = ar;
     integer n = ac;
     integer nrhs = bc;
@@ -463,6 +468,7 @@ int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
     free(work);
     free(AC);
     OK
+    #endif
 }
 
 //////////////////// least squares real linear system using SVD ////////////
@@ -527,6 +533,9 @@ int zgelss_(integer *m, integer *n, integer *nhrs,
     integer *info);
 
 int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
+    #ifdef _WIN32
+    return NOSPRTD;
+    #else
     integer m = ar;
     integer n = ac;
     integer nrhs = bc;
@@ -577,6 +586,7 @@ int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
     free(S);
     free(AC);
     OK
+    #endif
 }
 
 //////////////////// Cholesky factorization /////////////////////////