From 01a14ad32e0fd8586498ead61a426f20b724652b Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 22 Nov 2007 17:03:41 +0000 Subject: app1, app2, ... --- lib/Data/Packed/Internal/Common.hs | 25 +++++++++++++++++++++++++ lib/Data/Packed/Internal/Matrix.hs | 19 +++++++------------ lib/Data/Packed/Internal/Vector.hs | 30 ++++-------------------------- lib/Numeric/GSL/Fourier.hs | 3 +-- lib/Numeric/GSL/Matrix.hs | 36 ++++++++++++------------------------ lib/Numeric/GSL/Minimization.hs | 14 ++------------ lib/Numeric/GSL/Polynomials.hs | 3 +-- lib/Numeric/GSL/Vector.hs | 12 ++++-------- lib/Numeric/LinearAlgebra/LAPACK.hs | 33 +++++++++++---------------------- 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 -- cgit v1.2.3