9 files changed, 67 insertions, 108 deletions

diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index c3a733c..dc1c2b4 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -23,6 +23,8 @@ import Debug.Trace
 import Data.List(transpose,intersperse)
 import Data.Typeable
 import Data.Maybe(fromJust)
+import Foreign.C.String(peekCString)
+import Foreign.C.Types
 
 ----------------------------------------------------------------------
 instance (Storable a, RealFloat a) => Storable (Complex a) where    --
@@ -65,6 +67,13 @@ ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
 ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ \a1 -> ww2 w2 o2 w3 o3 (f a1)
 ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ \a1 -> ww3 w2 o2 w3 o3 w4 o4 (f a1)
 
+app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
+app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
+app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $
+     \a1 a2 a3 -> f // a1 // a2 // a3 // check s
+app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $ 
+     \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s
+
 -- GSL error codes are <= 1024
 -- | error codes for the auxiliary functions required by the wrappers
 errorCode :: Int -> String
@@ -78,6 +87,22 @@ errorCode 2006 = "the input matrix is not positive definite"
 errorCode 2007 = "not yet supported in this OS"
 errorCode n    = "code "++show n
 
+-- | check the error code
+check :: String -> IO Int -> IO ()
+check msg f = do
+    err <- f
+    when (err/=0) $ if err > 1024
+                      then (error (msg++": "++errorCode err)) -- our errors
+                      else do                                 -- GSL errors
+                        ps <- gsl_strerror err
+                        s <- peekCString ps
+                        error (msg++": "++s)
+    return ()
+
+-- | description of GSL error codes
+foreign import ccall "auxi.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar)
+
+
 {- | conversion of Haskell functions into function pointers that can be used in the C side
 -}
 foreign import ccall "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 90a96b5..0519603 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -79,8 +79,7 @@ cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transda
 fmat m@MF{} = m
 fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r}
 
---matc m f = f (rows m) (cols m) (ptr (cdat m))
---matf m f = f (rows m) (cols m) (ptr (fdat m))
+mat = withMatrix
 
 withMatrix MC {rows = r, cols = c, cdat = d } f =
     withForeignPtr (fptr d) $ \p -> do
@@ -308,8 +307,7 @@ subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double
 subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do
     r <- createMatrix RowMajor rt ct
     let x = cmat x'
-    ww2 withMatrix x withMatrix r $ \x r ->
-        c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // x // r // check "subMatrixR"
+    app2 (c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1)) mat x mat r "subMatrixR"
     return r
 foreign import ccall "auxi.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
 
@@ -333,8 +331,7 @@ subMatrix = subMatrixD
 
 diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
     m <- createMatrix RowMajor n n
-    ww2 withVector v withMatrix m $ \v m ->
-        fun // v // m // check msg
+    app2 fun vec v mat m msg
     return m
 
 -- | diagonal matrix from a real vector
@@ -356,19 +353,18 @@ diag = diagD
 constantAux fun x n = unsafePerformIO $ do
     v <- createVector n
     px <- newArray [x]
-    withVector v $ \v ->
-        fun px // v // check "constantAux"
+    app1 (fun px) vec v "constantAux"
     free px
     return v
 
 constantR :: Double -> Int -> Vector Double
 constantR = constantAux cconstantR
-foreign import ccall safe "auxi.h constantR"
+foreign import ccall "auxi.h constantR"
     cconstantR :: Ptr Double -> TV -- Double :> IO Int
 
 constantC :: Complex Double -> Int -> Vector (Complex Double)
 constantC = constantAux cconstantC
-foreign import ccall safe "auxi.h constantC"
+foreign import ccall "auxi.h constantC"
     cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int
 
 {- | creates a vector with a given number of equal components:
@@ -403,8 +399,7 @@ fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
 fromFile filename (r,c) = do
     charname <- newCString filename
     res <- createMatrix RowMajor r c
-    withMatrix res $ \res ->
-        c_gslReadMatrix charname // res // check "gslReadMatrix"
+    app1 (c_gslReadMatrix charname) mat res "gslReadMatrix"
     --free charname  -- TO DO: free the auxiliary CString
     return res
 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 386ebb5..7eee5fe 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -21,10 +21,6 @@ import Foreign
 import Complex
 import Control.Monad(when)
 import Data.List(transpose)
-import Debug.Trace(trace)
-import Foreign.C.String(peekCString)
-import Foreign.C.Types
-import Data.Monoid
 
 -- | A one-dimensional array of objects stored in a contiguous memory block.
 data Vector t = V { dim  :: Int              -- ^ number of elements
@@ -33,30 +29,14 @@ data Vector t = V { dim  :: Int              -- ^ number of elements
 
 --ptr (V _ fptr) = unsafeForeignPtrToPtr fptr
 
--- | check the error code
-check :: String -> IO Int -> IO ()
-check msg f = do
-    err <- f
-    when (err/=0) $ if err > 1024
-                      then (error (msg++": "++errorCode err)) -- our errors
-                      else do                                 -- GSL errors
-                        ps <- gsl_strerror err
-                        s <- peekCString ps
-                        error (msg++": "++s)
-    return ()
-
--- | description of GSL error codes
-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
 -- not yet admitted by my haddock version
 -- infixr 5 :>
 -- type t :> s = Vc t s
 
---- | adaptation of our vectors to be admitted by foreign functions: @f \/\/ vec v@
---vec :: Vector t -> (Vc t s) -> s
---vec v f = f (dim v) (ptr v)
+
+vec = withVector
 
 withVector (V n fp) f = withForeignPtr fp $ \p -> do
     let v f = do
@@ -80,8 +60,7 @@ fromList :: Storable a => [a] -> Vector a
 fromList l = unsafePerformIO $ do
     v <- createVector (length l)
     let f _ p = pokeArray p l >> return 0
-    withVector v $ \v ->
-        f // v // check "fromList"
+    app1 f vec v "fromList"
     return v
 
 safeRead v = unsafePerformIO . withForeignPtr (fptr v)
@@ -124,8 +103,7 @@ subVector k l (v@V {dim=n})
     | otherwise = unsafePerformIO $ do
         r <- createVector l
         let f _ s _ d = copyArray d (advancePtr s k) l >> return 0
-        ww2 withVector v withVector r $ \v r ->
-            f // v // r // check "subVector"
+        app2 f vec v vec r "subVector"
         return r
 
 {- | Reads a vector position:
diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs
index 4b08625..149c1e1 100644
--- a/lib/Numeric/GSL/Fourier.hs
+++ b/lib/Numeric/GSL/Fourier.hs
@@ -26,8 +26,7 @@ import Foreign
 
 genfft code v = unsafePerformIO $ do
     r <- createVector (dim v)
-    ww2 withVector v withVector r $ \ v r ->
-        c_fft code // v // r // check "fft"
+    app2 (c_fft code) vec v vec r "fft"
     return r
 
 foreign import ccall "gsl-aux.h fft" c_fft ::  Int -> TCVCV
diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
index 09a0be4..07d4660 100644
--- a/lib/Numeric/GSL/Matrix.hs
+++ b/lib/Numeric/GSL/Matrix.hs
@@ -51,8 +51,7 @@ eigSg' m
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix RowMajor r r
-        ww3 withMatrix m withVector l withMatrix v $ \m l v ->
-            c_eigS // m // l // v // check "eigSg"
+        app3 c_eigS mat m vec l mat v "eigSg"
         return (l,v)
   where r = rows m
 foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
@@ -85,8 +84,7 @@ eigHg' m
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix RowMajor r r
-        ww3 withMatrix m withVector l withMatrix v $ \m l v ->
-            c_eigH // m // l // v // check "eigHg"
+        app3 c_eigH mat m vec l mat v "eigHg"
         return (l,v)
   where r = rows m
 foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
@@ -122,8 +120,7 @@ svd' x = unsafePerformIO $ do
     u <- createMatrix RowMajor r c
     s <- createVector c
     v <- createMatrix RowMajor c c
-    ww4 withMatrix x withMatrix u withVector s withMatrix v $ \x u s v ->
-        c_svd // x // u // s // v // check "svdg"
+    app4 c_svd mat x mat u vec s mat v "svdg"
     return (u,s,v)
   where r = rows x
         c = cols x
@@ -152,8 +149,7 @@ qr = qr' . cmat
 qr' x = unsafePerformIO $ do
     q <- createMatrix RowMajor r r
     rot <- createMatrix RowMajor r c
-    ww3 withMatrix x withMatrix q withMatrix rot $ \x q rot ->
-        c_qr // x // q // rot // check "qr"
+    app3 c_qr mat x mat q mat rot  "qr"
     return (q,rot)
   where r = rows x
         c = cols x
@@ -165,8 +161,7 @@ qrPacked = qrPacked' . cmat
 qrPacked' x = unsafePerformIO $ do
     qr <- createMatrix RowMajor r c
     tau <- createVector (min r c)
-    ww3 withMatrix x withMatrix qr withVector tau $ \x qr tau ->
-        c_qrPacked // x // qr // tau // check "qrUnpacked"
+    app3 c_qrPacked mat x mat qr vec tau "qrUnpacked"
     return (qr,tau)
   where r = rows x
         c = cols x
@@ -178,8 +173,7 @@ unpackQR (qr,tau) = unpackQR' (cmat qr, tau)
 unpackQR' (qr,tau) = unsafePerformIO $ do
     q <- createMatrix RowMajor r r
     res <- createMatrix RowMajor r c
-    ww4 withMatrix qr withVector tau withMatrix q withMatrix res $ \qr tau q res ->
-        c_qrUnpack // qr // tau // q // res // check "qrUnpack"
+    app4 c_qrUnpack mat qr vec tau mat q mat res "qrUnpack"
     return (q,res)
   where r = rows qr
         c = cols qr
@@ -203,8 +197,7 @@ cholR = cholR' . cmat
 
 cholR' x = unsafePerformIO $ do
     r <- createMatrix RowMajor n n
-    ww2 withMatrix x withMatrix r $ \x r ->
-        c_cholR // x // r // check "cholR"
+    app2 c_cholR mat x mat r "cholR"
     return r
   where n = rows x
 foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
@@ -214,8 +207,7 @@ cholC = cholC' . cmat
 
 cholC' x = unsafePerformIO $ do
     r <- createMatrix RowMajor n n
-    ww2 withMatrix x withMatrix r $ \x r ->
-        c_cholC // x // r // check "cholC"
+    app2 c_cholC mat x mat r "cholC"
     return r
   where n = rows x
 foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
@@ -231,8 +223,7 @@ luSolveR a b = luSolveR' (cmat a) (cmat b)
 luSolveR' a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix RowMajor r c
-        ww3 withMatrix a withMatrix b withMatrix s $ \ a b s ->
-            c_luSolveR // a // b // s // check "luSolveR"
+        app3 c_luSolveR mat a mat b mat s "luSolveR"
         return s
     | otherwise = error "luSolveR of nonsquare matrix"
   where n1 = rows a
@@ -249,8 +240,7 @@ luSolveC a b = luSolveC' (cmat a) (cmat b)
 luSolveC' a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix RowMajor r c
-        ww3 withMatrix a withMatrix b withMatrix s $ \ a b s ->
-            c_luSolveC // a // b // s // check "luSolveC"
+        app3 c_luSolveC mat a mat b mat s "luSolveC"
         return s
     | otherwise = error "luSolveC of nonsquare matrix"
   where n1 = rows a
@@ -266,8 +256,7 @@ luRaux = luRaux' . cmat
 
 luRaux' x = unsafePerformIO $ do
     res <- createVector (r*r+r+1)
-    ww2 withMatrix x withVector res $ \x res ->
-        c_luRaux // x // res // check "luRaux"
+    app2 c_luRaux mat x vec res "luRaux"
     return res
   where r = rows x
         c = cols x
@@ -280,8 +269,7 @@ luCaux = luCaux' . cmat
 
 luCaux' x = unsafePerformIO $ do
     res <- createVector (r*r+r+1)
-    ww2 withMatrix x withVector res $ \x res ->
-        c_luCaux // x // res // check "luCaux"
+    app2 c_luCaux mat x vec res "luCaux"
     return res
   where r = rows x
         c = cols x
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index e44b2e5..235a88f 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -196,22 +196,12 @@ aux_vTov f n p r = g where
 
 --------------------------------------------------------------------
 
-
 createV n fun msg = unsafePerformIO $ do
     r <- createVector n
-    withVector r $ \ r ->
-        fun // r // check msg
+    app1 fun vec r msg
     return r
 
-{-
-createM r c fun msg ptrs = unsafePerformIO $ do
-    r <- createMatrix RowMajor r c
-    fun // matc r // check msg ptrs
-    return r
--}
-
 createMIO r c fun msg = do
     r <- createMatrix RowMajor r c
-    withMatrix r $ \ r ->
-        fun // r // check msg
+    app1 fun mat r msg
     return r
diff --git a/lib/Numeric/GSL/Polynomials.hs b/lib/Numeric/GSL/Polynomials.hs
index e663711..3e0cf81 100644
--- a/lib/Numeric/GSL/Polynomials.hs
+++ b/lib/Numeric/GSL/Polynomials.hs
@@ -47,8 +47,7 @@ polySolve = toList . polySolve' . fromList
 polySolve' :: Vector Double -> Vector (Complex Double)
 polySolve' v | dim v > 1 = unsafePerformIO $ do
     r <- createVector (dim v-1)
-    ww2 withVector v withVector r $ \ v r ->
-        c_polySolve // v // r // check "polySolve"
+    app2 c_polySolve vec v vec r "polySolve"
     return r
              | otherwise = error "polySolve on a polynomial of degree zero"
 
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 65f3a2e..86d8f10 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -73,28 +73,24 @@ data FunCodeS = Norm2
 
 toScalarAux fun code v = unsafePerformIO $ do
     r <- createVector 1
-    ww2 withVector v withVector r $ \v r ->
-        fun (fromEnum code) // v // r // check "toScalarAux"
+    app2 (fun (fromEnum code)) vec v vec r "toScalarAux"
     return (r `at` 0)
 
 vectorMapAux fun code v = unsafePerformIO $ do
     r <- createVector (dim v)
-    ww2 withVector v withVector r $ \v r ->
-        fun (fromEnum code) // v // r // check "vectorMapAux"
+    app2 (fun (fromEnum code)) vec v vec r "vectorMapAux"
     return r
 
 vectorMapValAux fun code val v = unsafePerformIO $ do
     r <- createVector (dim v)
     pval <- newArray [val]
-    ww2 withVector v withVector r $ \v r ->
-        fun (fromEnum code) pval // v // r // check "vectorMapValAux"
+    app2 (fun (fromEnum code) pval) vec v vec r "vectorMapValAux"
     free pval
     return r
 
 vectorZipAux fun code u v = unsafePerformIO $ do
     r <- createVector (dim u)
-    ww3 withVector u withVector v withVector r $ \u v r ->
-        fun (fromEnum code) // u // v // r // check "vectorZipAux"
+    app3 (fun (fromEnum code)) vec u vec v vec r "vectorZipAux"
     return r
 
 ---------------------------------------------------------------------
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 19516e3..936cb19 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -61,8 +61,7 @@ svdAux f st x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
-    ww4 withMatrix x withMatrix u withVector s withMatrix v $ \x u s v ->
-        f // x // u // s // v // check st
+    app4 f mat x mat u vec s mat v st
     return (u,s,trans v)
   where r = rows x
         c = cols x
@@ -74,8 +73,7 @@ eigAux f st m
         l <- createVector r
         v <- createMatrix ColumnMajor r r
         dummy <- createMatrix ColumnMajor 1 1
-        ww4 withMatrix m withMatrix dummy withVector l withMatrix v $ \m dummy l v ->
-            f // m // dummy // l // v // check st
+        app4 f mat m mat dummy vec l mat v st
         return (l,v)
   where r = rows m
 
@@ -117,8 +115,7 @@ eigRaux m
         l <- createVector r
         v <- createMatrix ColumnMajor r r
         dummy <- createMatrix ColumnMajor 1 1
-        ww4 withMatrix m withMatrix dummy withVector l withMatrix v $ \m dummy l v ->
-            dgeev // m // dummy // l // v // check "eigR"
+        app4 dgeev mat m mat dummy vec l mat v "eigR"
         return (l,v)
   where r = rows m
 
@@ -147,8 +144,7 @@ eigS' m
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
-        ww3 withMatrix m withVector l withMatrix v $ \m l v ->
-            dsyev // m // l // v // check "eigS"
+        app3 dsyev mat m vec l mat v "eigS"
         return (l,v)
   where r = rows m
 
@@ -170,8 +166,7 @@ eigH' m
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
-        ww3 withMatrix m withVector l withMatrix v $ \m l v ->
-            zheev // m // l // v // check "eigH"
+        app3 zheev mat m vec l mat v "eigH"
         return (l,v)
   where r = rows m
 
@@ -182,8 +177,7 @@ foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
 linearSolveSQAux f st a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix ColumnMajor r c
-        ww3 withMatrix a withMatrix b withMatrix s $ \a b s ->
-            f // a // b // s // check st
+        app3 f mat a mat b mat s st
         return s
     | otherwise = error $ st ++ " of nonsquare matrix"
   where n1 = rows a
@@ -207,8 +201,7 @@ foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double ->
 
 linearSolveAux f st a b = unsafePerformIO $ do
     r <- createMatrix ColumnMajor (max m n) nrhs
-    ww3 withMatrix a withMatrix b withMatrix r $ \a b r ->
-            f // a // b // r // check st
+    app3 f mat a mat b mat r st
     return r
   where m = rows a
         n = cols a
@@ -258,8 +251,7 @@ cholS = cholAux dpotrf "cholS" . fmat
 
 cholAux f st a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor n n
-    ww2 withMatrix a withMatrix r $ \a r ->
-        f // a // r // check st
+    app2 f mat a mat r st
     return r
   where n = rows a
 
@@ -278,8 +270,7 @@ qrC = qrAux zgeqr2 "qrC" . fmat
 qrAux f st a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor m n
     tau <- createVector mn
-    ww3 withMatrix a withMatrix r withVector tau $ \ a r tau ->
-        f // a // tau // r // check st
+    app3 f mat a vec tau mat r st
     return (r,tau)
   where m = rows a
         n = cols a
@@ -300,8 +291,7 @@ hessC = hessAux zgehrd "hessC" . fmat
 hessAux f st a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor m n
     tau <- createVector (mn-1)
-    ww3 withMatrix a withMatrix r withVector tau $ \ a r tau ->
-        f // a // tau // r // check st
+    app3 f mat a vec tau mat r st
     return (r,tau)
   where m = rows a
         n = cols a
@@ -322,8 +312,7 @@ schurC = schurAux zgees "schurC" . fmat
 schurAux f st a = unsafePerformIO $ do
     u <- createMatrix ColumnMajor n n
     s <- createMatrix ColumnMajor n n
-    ww3 withMatrix a withMatrix u withMatrix s $ \ a u s ->
-        f // a // u // s // check st
+    app3 f mat a mat u mat s st
     return (u,s)
   where n = rows a