restored C trans and constant for comparison

3 files changed, 91 insertions, 4 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index b1e7670..13ffc34 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -221,13 +221,13 @@ class (Storable a, Floating a) => Element a where
 
 instance Element Double where
     subMatrixD = subMatrix'
-    transdata  = transdata'
-    constantD   = constant'
+    transdata  = transdataAux ctransR     -- transdata'
+    constantD  = constantAux cconstantR   -- constant'
 
 instance Element (Complex Double) where
     subMatrixD = subMatrix'
-    transdata  = transdata'
-    constantD   = constant'
+    transdata  = transdataAux ctransC     -- transdata'
+    constantD  = constantAux cconstantC   -- constant'
 
 -------------------------------------------------------------------
 
@@ -257,6 +257,23 @@ transdata' c1 v c2 =
 -- The above pragmas only seem to work on top level defs
 -- Fortunately everything seems to work using the above class
 
+-- C versions, still a little faster:
+
+transdataAux fun c1 d c2 =
+    if noneed
+        then d
+        else unsafePerformIO $ do
+            v <- createVector (dim d)
+            withForeignPtr (fptr d) $ \pd ->
+                withForeignPtr (fptr v) $ \pv ->
+                    fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
+            return v
+  where r1 = dim d `div` c1
+        r2 = dim d `div` c2
+        noneed = r1 == 1 || c1 == 1
+
+foreign import ccall "transR" ctransR :: TMM
+foreign import ccall "transC" ctransC :: TCMCM
 ----------------------------------------------------------------------
 
 constant' v n = unsafePerformIO $ do
@@ -267,6 +284,22 @@ constant' v n = unsafePerformIO $ do
         go (n-1)
     return w
 
+-- C versions
+
+constantAux fun x n = unsafePerformIO $ do
+    v <- createVector n
+    px <- newArray [x]
+    app1 (fun px) vec v "constantAux"
+    free px
+    return v
+
+constantR :: Double -> Int -> Vector Double
+constantR = constantAux cconstantR
+foreign import ccall "constantR" cconstantR :: Ptr Double -> TV
+
+constantC :: Complex Double -> Int -> Vector (Complex Double)
+constantC = constantAux cconstantC
+foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV
 ----------------------------------------------------------------------
 
 -- | Extracts a submatrix from a matrix.
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index a8ccf5f..d7248d1 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -908,3 +908,51 @@ int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {
            (doublecomplex*)rp,&ldc);
     OK
 }
+
+//////////////////// transpose /////////////////////////
+
+int transR(KDMAT(x),DMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transR");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+int transC(KCMAT(x),CMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transC");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        ((doublecomplex*)tp)[i*tc+j] = ((doublecomplex*)xp)[j*xc+i];
+        }
+    }
+    OK
+}
+
+//////////////////// constant /////////////////////////
+
+int constantR(double * pval, DVEC(r)) {
+    DEBUGMSG("constantR")
+    int k;
+    double val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int constantC(doublecomplex* pval, CVEC(r)) {
+    DEBUGMSG("constantC")
+    int k;
+    doublecomplex val = *pval;
+    for(k=0;k<rn;k++) {
+        ((doublecomplex*)rp)[k]=val;
+    }
+    OK
+}
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index 3f58243..dc7a98f 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -55,6 +55,12 @@ typedef short ftnlen;
 int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r));
 int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r));
 
+int transR(KDMAT(x),DMAT(t));
+int transC(KCMAT(x),CMAT(t));
+
+int constantR(double * pval, DVEC(r));
+int constantC(doublecomplex* pval, CVEC(r));
+
 int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
 int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
 int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));