C functions for Storable a => transdata, constantD

3 files changed, 57 insertions, 2 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 7a17ef0..090826d 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -252,10 +252,11 @@ class (Storable a, Floating a) => Element a where
                -> Matrix a -> Matrix a
     subMatrixD = subMatrix'
     transdata :: Int -> Vector a -> Int -> Vector a
-    transdata = transdata'
+    transdata = transdataP -- transdata'
     constantD  :: a -> Int -> Vector a
-    constantD = constant'
+    constantD = constantP -- constant'
     conjugateD :: Vector a -> Vector a
+    conjugateD = id
 
 instance Element Float where
     transdata  = transdataAux ctransF
@@ -320,10 +321,27 @@ transdataAux fun c1 d c2 =
         r2 = dim d `div` c2
         noneed = r1 == 1 || c1 == 1
 
+transdataP :: Storable a => Int -> Vector a -> Int -> Vector a
+transdataP c1 d c2 =
+    if noneed
+       then d
+       else unsafePerformIO $ do
+          v <- createVector (dim d)
+          unsafeWith d $ \pd ->
+              unsafeWith v $ \pv ->
+                  ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataStorable"
+          return v
+   where r1 = dim d `div` c1
+         r2 = dim d `div` c2
+         sz = sizeOf (d @> 0)
+         noneed = r1 == 1 || c1 == 1
+
 foreign import ccall "transF" ctransF :: TFMFM
 foreign import ccall "transR" ctransR :: TMM
 foreign import ccall "transQ" ctransQ :: TQMQM
 foreign import ccall "transC" ctransC :: TCMCM
+foreign import ccall "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
+
 ----------------------------------------------------------------------
 
 constant' v n = unsafePerformIO $ do
@@ -359,6 +377,17 @@ constantC :: Complex Double -> Int -> Vector (Complex Double)
 constantC = constantAux cconstantC
 foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV
 
+constantP :: Storable a => a -> Int -> Vector a
+constantP a n = unsafePerformIO $ do
+    let sz = sizeOf a
+    v <- createVector n
+    unsafeWith v $ \p -> do
+       alloca $ \k -> do
+                      poke k a
+                      cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"
+    return v
+foreign import ccall "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt
+
 ---------------------------------------
 
 conjugateAux fun x = unsafePerformIO $ do
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 09979cd..e8bbbdb 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -1139,6 +1139,19 @@ int transC(KCMAT(x),CMAT(t)) {
     OK
 }
 
+int transP(KPMAT(x), PMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    REQUIRES(xs==ts,NOCONVER);
+    DEBUGMSG("transP");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+          memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs);
+        }
+    }
+    OK
+}
+
 //////////////////// constant /////////////////////////
 
 int constantF(float * pval, FVEC(r)) {
@@ -1181,6 +1194,15 @@ int constantC(doublecomplex* pval, CVEC(r)) {
     OK
 }
 
+int constantP(void* pval, PVEC(r)) {
+    DEBUGMSG("constantP")
+    int k;
+    for(k=0;k<rn;k++) {
+      memcpy(rp+k*rs,pval,rs);
+    }
+    OK
+}
+
 //////////////////// float-double conversion /////////////////////////
 
 int float2double(FVEC(x),DVEC(y)) {
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index 3de5b19..12dcde4 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -44,19 +44,23 @@ typedef short ftnlen;
 #define DVEC(A) int A##n, double*A##p
 #define QVEC(A) int A##n, complex*A##p
 #define CVEC(A) int A##n, doublecomplex*A##p
+#define PVEC(A) int A##n, void* A##p, int A##s
 #define FMAT(A) int A##r, int A##c, float* A##p
 #define DMAT(A) int A##r, int A##c, double* A##p
 #define QMAT(A) int A##r, int A##c, complex* A##p
 #define CMAT(A) int A##r, int A##c, doublecomplex* A##p
+#define PMAT(A) int A##r, int A##c, void* A##p, int A##s
 
 #define KFVEC(A) int A##n, const float*A##p
 #define KDVEC(A) int A##n, const double*A##p
 #define KQVEC(A) int A##n, const complex*A##p
 #define KCVEC(A) int A##n, const doublecomplex*A##p
+#define KPVEC(A) int A##n, const void* A##p, int A##s
 #define KFMAT(A) int A##r, int A##c, const float* A##p
 #define KDMAT(A) int A##r, int A##c, const double* A##p
 #define KQMAT(A) int A##r, int A##c, const complex* A##p
 #define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p
+#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s
 
 /********************************************************/