From f0a2248b9d9532e0814519408bff8b886bca2bbc Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 14 Nov 2008 16:16:04 +0000 Subject: subMatrix reimplemented and removed auxi.c --- lib/Data/Packed/Internal/Common.hs | 2 +- lib/Data/Packed/Internal/Matrix.hs | 43 ++++++++++++----------- lib/Data/Packed/Internal/auxi.c | 71 -------------------------------------- lib/Data/Packed/Internal/auxi.h | 19 ---------- 4 files changed, 24 insertions(+), 111 deletions(-) delete mode 100644 lib/Data/Packed/Internal/auxi.c delete mode 100644 lib/Data/Packed/Internal/auxi.h (limited to 'lib/Data/Packed/Internal') diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs index 1b63dd8..46e9612 100644 --- a/lib/Data/Packed/Internal/Common.hs +++ b/lib/Data/Packed/Internal/Common.hs @@ -82,7 +82,7 @@ errorCode n = "code "++show n -- | clear the fpu -foreign import ccall "auxi.h asm_finit" finit :: IO () +foreign import ccall "asm_finit" finit :: IO () -- | check the error code check :: String -> IO CInt -> IO () diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 44204a1..b1e7670 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -220,12 +220,12 @@ class (Storable a, Floating a) => Element a where constantD :: a -> Int -> Vector a instance Element Double where - subMatrixD = subMatrixR + subMatrixD = subMatrix' transdata = transdata' constantD = constant' instance Element (Complex Double) where - subMatrixD = subMatrixC + subMatrixD = subMatrix' transdata = transdata' constantD = constant' @@ -269,29 +269,32 @@ constant' v n = unsafePerformIO $ do ---------------------------------------------------------------------- --- | extraction of a submatrix from a real matrix -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' - app2 (c_submatrixR (fi r0) (fi $ r0+rt-1) (fi c0) (fi $ c0+ct-1)) mat x mat r "subMatrixR" - return r -foreign import ccall "auxi.h submatrixR" c_submatrixR :: CInt -> CInt -> CInt -> CInt -> TMM - --- | extraction of a submatrix from a complex matrix -subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double) -subMatrixC (r0,c0) (rt,ct) x = - reshape ct . asComplex . flatten . - subMatrixR (r0,2*c0) (rt,2*ct) . - reshape (2*cols x) . asReal . flatten $ x - -- | Extracts a submatrix from a matrix. subMatrix :: Element a => (Int,Int) -- ^ (r0,c0) starting position -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix -> Matrix a -- ^ input matrix -> Matrix a -- ^ result -subMatrix = subMatrixD +subMatrix (r0,c0) (rt,ct) m + | 0 <= r0 && 0 < rt && r0+rt <= (rows m) && + 0 <= c0 && 0 < ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m + | otherwise = error $ "wrong subMatrix "++ + show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m) + +subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do + w <- createVector (rt*ct) + withForeignPtr (fptr v) $ \p -> + withForeignPtr (fptr w) $ \q -> do + let go (-1) _ = return () + go !i (-1) = go (i-1) (ct-1) + go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0) + pokeElemOff q (i*ct+j) x + go i (j-1) + go (rt-1) (ct-1) + return w + +subMatrix' (r0,c0) (rt,ct) (MC _r c v) = MC rt ct $ subMatrix'' (r0,c0) (rt,ct) c v +subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m) -------------------------------------------------------------------------- @@ -316,4 +319,4 @@ fromFile filename (r,c) = do 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 +foreign import ccall "matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c deleted file mode 100644 index 5c06cb6..0000000 --- a/lib/Data/Packed/Internal/auxi.c +++ /dev/null @@ -1,71 +0,0 @@ -#include "auxi.h" -#include -#include -#include - -#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 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 BAD_SIZE 2000 -#define BAD_CODE 2001 -#define MEM 2002 -#define BAD_FILE 2003 - - -int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)) { - REQUIRES(0<=r1 && r1<=r2 && r2 - -#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 - -int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)); - -const char * gsl_strerror (const int gsl_errno); - -int matrix_fscanf(char*filename, RMAT(a)); - -void asm_finit(); -- cgit v1.2.3