subMatrix reimplemented and removed auxi.c

5 files changed, 27 insertions, 111 deletions

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 <gsl/gsl_matrix.h>
-#include <string.h>
-#include <stdio.h>
-
-#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<xr && 0<=c1 && c1<=c2 && c2<xc &&
-            rr==r2-r1+1 && rc==c2-c1+1,BAD_SIZE);
-    DEBUGMSG("submatrixR");
-    KDMVIEW(x);
-    DMVIEW(r);
-    gsl_matrix_const_view S = gsl_matrix_const_submatrix(M(x),r1,c1,rr,rc);
-    int res = gsl_matrix_memcpy(M(r),M(S));
-    CHECK(res,res);
-    OK
-}
-
-//---------------------------------------
-void asm_finit() {
-    asm("finit");
-}
-//---------------------------------------
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
deleted file mode 100644
index c234ef5..0000000
--- a/lib/Data/Packed/Internal/auxi.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <gsl/gsl_complex.h>
-
-#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();
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index e85c1b7..6315041 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -34,6 +34,9 @@
 #define NODEFPOS 2006
 #define NOSPRTD  2007
 
+//---------------------------------------
+void asm_finit() { asm("finit");}
+//---------------------------------------
 
 //////////////////// real svd ////////////////////////////////////