summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--HSSL.cabal2
-rw-r--r--README47
-rw-r--r--examples/parallel.hs11
-rw-r--r--lib/Data/Packed/Internal/Common.hs1
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs20
-rw-r--r--lib/Data/Packed/Internal/Vector.hs2
-rw-r--r--lib/Data/Packed/Internal/auxi.c (renamed from lib/Data/Packed/Internal/aux.c)2
-rw-r--r--lib/Data/Packed/Internal/auxi.h (renamed from lib/Data/Packed/Internal/aux.h)0
-rw-r--r--lib/Numeric/GSL/gsl-aux.c6
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c10
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,
68 Numeric.LinearAlgebra.Interface, 68 Numeric.LinearAlgebra.Interface,
69 Numeric.LinearAlgebra.Algorithms, 69 Numeric.LinearAlgebra.Algorithms,
70 Graphics.Plot 70 Graphics.Plot
71C-sources: lib/Data/Packed/Internal/aux.c, 71C-sources: lib/Data/Packed/Internal/auxi.c,
72 lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c, 72 lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
73 lib/Numeric/GSL/gsl-aux.c 73 lib/Numeric/GSL/gsl-aux.c
74extra-libraries: gsl cblas lapack 74extra-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
46Prelude Numeric.GSL> :q 46Prelude Numeric.GSL> :q
47Leaving GHCi. 47Leaving GHCi.
48 48
49-------------------------------------------------------------------------------------- 49CHANGES
50
51This is a new version of GSLHaskell. The package is (provisionally)
52called \C{hssl} (a simple scientific library for Haskell) because only
53a small part of GSL is available and linear algebra is based on LAPACK.
54
55The code has been extensively refactored. There is a new internal representation
56which admits both C and Fortran matrices and avoids many transposes.
57
58There are only minor API changes:
59
60- the matrix product operator \C{(<>)} is now overloaded only for matrix-matrix,
61matrix-vector and vector-matrix, with the same base type. The dot product and the scaling
62of vectors or matrices is now denoted by `dot` and `scale`. Conversions from real to
63complex objects must be explicit.
64
65- Most linear algebra functions admit both real and complex objects. Utilities such as
66ident or constant are now polymorphic.
67
68- Runtime errors produced by GSL or LAPACK can be handled using \C{Control.Exeception.catch}.
69
70Old GSLHaskell code will work with small modifications.
71
50ACKNOWLEDGEMENTS 72ACKNOWLEDGEMENTS
51 73
52I thank Henning Thielemann and all the people in the Haskell mailing lists for their help. 74I 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
63- Pedro E. López de Teruel fixed the interface to lapack. 85- Pedro E. López de Teruel fixed the interface to lapack.
64 86
65- Antti Siira discovered a bug in the plotting functions. 87- Antti Siira discovered a bug in the plotting functions.
88
89INSTALLATION ON WINDOWS
90
911) Download the developer files gsl-1.8-lib.zip from
92 http://gnuwin32.sourceforge.net/packages/gsl.htm
93 and copy the gsl folder (under include) to the include folder of ghc:
94 C:\ghc\ghc.6.6.1\include
95
962) Install the package as usual from the command line in the hssl-0.1 folder:
97 runhaskell Setup.lhs configure
98 runhaskell Setup.lhs build
99 runhaskell Setup.lhs install
100
1013) Copy libgsl.dll, libcblas.dll (from the binaries package gsl-1.8.bin.zip)
102 and liblapack.dll (borrowed from the R system) to the folder where
103 hssl has been installed: C:\haskell\hss-0.1\ghc-6.6.1. They are needed to compile programs.
104 These three dlls are also available from http://perception.inf.um.es/darcs/HSSL/dll1.zip
105
1064) Copy the dlls at http://perception.inf.um.es/darcs/HSSL/dll2.zip to C:\windows\system
107 They are required to run the programs and ghci.
108
109Unfortunately the lapack dll supplied by the R system does not include zgels_ and zgelss_,
110so 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
5 5
6inParallel = parMap rwhnf id 6inParallel = parMap rwhnf id
7 7
8parMul x y = fromBlocks [[ay],[by]] 8parMul x y = fromBlocks [inParallel[x <> y1, x <> y2]]
9 where p = rows x `div` 2 9 where p = cols y `div` 2
10 a = takeRows p x 10 (y1,y2) = splitColumnsAt p y
11 b = dropRows p x
12 [ay,by] = inParallel [a<>y,b<>y]
13 11
14main = do 12main = do
15 n <- (read . head) `fmap` getArgs 13 n <- (read . head) `fmap` getArgs
@@ -20,6 +18,9 @@ main = do
20a = (2><3) [1..6::Double] 18a = (2><3) [1..6::Double]
21b = (3><4) [1..12::Double] 19b = (3><4) [1..12::Double]
22 20
21splitRowsAt p m = (takeRows p m, dropRows p m)
22splitColumnsAt p m = (takeColumns p m, dropColumns p m)
23
23time act = do 24time act = do
24 t0 <- getClockTime 25 t0 <- getClockTime
25 act 26 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"
70errorCode 2004 = "singular" 70errorCode 2004 = "singular"
71errorCode 2005 = "didn't converge" 71errorCode 2005 = "didn't converge"
72errorCode 2006 = "the input matrix is not positive definite" 72errorCode 2006 = "the input matrix is not positive definite"
73errorCode 2007 = "not yet supported in this OS"
73errorCode n = "code "++show n 74errorCode n = "code "++show n
74 75
75{- | conversion of Haskell functions into function pointers that can be used in the C side 76{- | 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 =
240 r2 = dim d `div` c2 240 r2 = dim d `div` c2
241 noneed = r1 == 1 || c1 == 1 241 noneed = r1 == 1 || c1 == 1
242 242
243foreign import ccall safe "aux.h transR" 243foreign import ccall safe "auxi.h transR"
244 ctransR :: TMM -- Double ::> Double ::> IO Int 244 ctransR :: TMM -- Double ::> Double ::> IO Int
245foreign import ccall safe "aux.h transC" 245foreign import ccall safe "auxi.h transC"
246 ctransC :: TCMCM -- Complex Double ::> Complex Double ::> IO Int 246 ctransC :: TCMCM -- Complex Double ::> Complex Double ::> IO Int
247 247
248------------------------------------------------------------------ 248------------------------------------------------------------------
@@ -258,14 +258,14 @@ multiplyAux fun a b = unsafePerformIO $ do
258 return r 258 return r
259 259
260multiplyR = multiplyAux cmultiplyR 260multiplyR = multiplyAux cmultiplyR
261foreign import ccall safe "aux.h multiplyR" 261foreign import ccall safe "auxi.h multiplyR"
262 cmultiplyR :: Int -> Int -> Int -> Ptr Double 262 cmultiplyR :: Int -> Int -> Int -> Ptr Double
263 -> Int -> Int -> Int -> Ptr Double 263 -> Int -> Int -> Int -> Ptr Double
264 -> Int -> Int -> Ptr Double 264 -> Int -> Int -> Ptr Double
265 -> IO Int 265 -> IO Int
266 266
267multiplyC = multiplyAux cmultiplyC 267multiplyC = multiplyAux cmultiplyC
268foreign import ccall safe "aux.h multiplyC" 268foreign import ccall safe "auxi.h multiplyC"
269 cmultiplyC :: Int -> Int -> Int -> Ptr (Complex Double) 269 cmultiplyC :: Int -> Int -> Int -> Ptr (Complex Double)
270 -> Int -> Int -> Int -> Ptr (Complex Double) 270 -> Int -> Int -> Int -> Ptr (Complex Double)
271 -> Int -> Int -> Ptr (Complex Double) 271 -> Int -> Int -> Ptr (Complex Double)
@@ -288,7 +288,7 @@ subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do
288 r <- createMatrix RowMajor rt ct 288 r <- createMatrix RowMajor rt ct
289 c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r] 289 c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r]
290 return r 290 return r
291foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM 291foreign import ccall "auxi.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
292 292
293-- | extraction of a submatrix from a complex matrix 293-- | extraction of a submatrix from a complex matrix
294subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double) 294subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double)
@@ -316,12 +316,12 @@ diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
316-- | diagonal matrix from a real vector 316-- | diagonal matrix from a real vector
317diagR :: Vector Double -> Matrix Double 317diagR :: Vector Double -> Matrix Double
318diagR = diagAux c_diagR "diagR" 318diagR = diagAux c_diagR "diagR"
319foreign import ccall "aux.h diagR" c_diagR :: TVM 319foreign import ccall "auxi.h diagR" c_diagR :: TVM
320 320
321-- | diagonal matrix from a real vector 321-- | diagonal matrix from a real vector
322diagC :: Vector (Complex Double) -> Matrix (Complex Double) 322diagC :: Vector (Complex Double) -> Matrix (Complex Double)
323diagC = diagAux c_diagC "diagC" 323diagC = diagAux c_diagC "diagC"
324foreign import ccall "aux.h diagC" c_diagC :: TCVCM 324foreign import ccall "auxi.h diagC" c_diagC :: TCVCM
325 325
326-- | creates a square matrix with the given diagonal 326-- | creates a square matrix with the given diagonal
327diag :: Field a => Vector a -> Matrix a 327diag :: Field a => Vector a -> Matrix a
@@ -338,12 +338,12 @@ constantAux fun x n = unsafePerformIO $ do
338 338
339constantR :: Double -> Int -> Vector Double 339constantR :: Double -> Int -> Vector Double
340constantR = constantAux cconstantR 340constantR = constantAux cconstantR
341foreign import ccall safe "aux.h constantR" 341foreign import ccall safe "auxi.h constantR"
342 cconstantR :: Ptr Double -> TV -- Double :> IO Int 342 cconstantR :: Ptr Double -> TV -- Double :> IO Int
343 343
344constantC :: Complex Double -> Int -> Vector (Complex Double) 344constantC :: Complex Double -> Int -> Vector (Complex Double)
345constantC = constantAux cconstantC 345constantC = constantAux cconstantC
346foreign import ccall safe "aux.h constantC" 346foreign import ccall safe "auxi.h constantC"
347 cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int 347 cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int
348 348
349{- | creates a vector with a given number of equal components: 349{- | creates a vector with a given number of equal components:
@@ -381,7 +381,7 @@ fromFile filename (r,c) = do
381 c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" [] 381 c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
382 --free charname -- TO DO: free the auxiliary CString 382 --free charname -- TO DO: free the auxiliary CString
383 return res 383 return res
384foreign import ccall "aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM 384foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
385 385
386------------------------------------------------------------------------- 386-------------------------------------------------------------------------
387 387
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
46 return () 46 return ()
47 47
48-- | description of GSL error codes 48-- | description of GSL error codes
49foreign import ccall "aux.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar) 49foreign import ccall "auxi.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar)
50 50
51-- | signature of foreign functions admitting C-style vectors 51-- | signature of foreign functions admitting C-style vectors
52type Vc t s = Int -> Ptr t -> s 52type 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 @@
1#include "aux.h" 1#include "auxi.h"
2#include <gsl/gsl_blas.h> 2#include <gsl/gsl_blas.h>
3#include <gsl/gsl_linalg.h> 3#include <gsl/gsl_linalg.h>
4#include <gsl/gsl_matrix.h> 4#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)) {
116 OP(7,sinh) 116 OP(7,sinh)
117 OP(8,cosh) 117 OP(8,cosh)
118 OP(9,tanh) 118 OP(9,tanh)
119 OP(10,asinh) 119 OP(10,gsl_asinh)
120 OP(11,acosh) 120 OP(11,gsl_acosh)
121 OP(12,atanh) 121 OP(12,gsl_atanh)
122 OP(13,exp) 122 OP(13,exp)
123 OP(14,log) 123 OP(14,log)
124 OP(15,sign) 124 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 @@
29#define SINGULAR 2004 29#define SINGULAR 2004
30#define NOCONVER 2005 30#define NOCONVER 2005
31#define NODEFPOS 2006 31#define NODEFPOS 2006
32#define NOSPRTD 2007
33
32 34
33//////////////////// real svd //////////////////////////////////// 35//////////////////// real svd ////////////////////////////////////
34 36
@@ -423,6 +425,9 @@ int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
423//////////////////// least squares complex linear system //////////// 425//////////////////// least squares complex linear system ////////////
424 426
425int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) { 427int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
428 #ifdef _WIN32
429 return NOSPRTD;
430 #else
426 integer m = ar; 431 integer m = ar;
427 integer n = ac; 432 integer n = ac;
428 integer nrhs = bc; 433 integer nrhs = bc;
@@ -463,6 +468,7 @@ int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
463 free(work); 468 free(work);
464 free(AC); 469 free(AC);
465 OK 470 OK
471 #endif
466} 472}
467 473
468//////////////////// least squares real linear system using SVD //////////// 474//////////////////// least squares real linear system using SVD ////////////
@@ -527,6 +533,9 @@ int zgelss_(integer *m, integer *n, integer *nhrs,
527 integer *info); 533 integer *info);
528 534
529int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) { 535int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
536 #ifdef _WIN32
537 return NOSPRTD;
538 #else
530 integer m = ar; 539 integer m = ar;
531 integer n = ac; 540 integer n = ac;
532 integer nrhs = bc; 541 integer nrhs = bc;
@@ -577,6 +586,7 @@ int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
577 free(S); 586 free(S);
578 free(AC); 587 free(AC);
579 OK 588 OK
589 #endif
580} 590}
581 591
582//////////////////// Cholesky factorization ///////////////////////// 592//////////////////// Cholesky factorization /////////////////////////