alternative multiply versions

3 files changed, 20 insertions, 113 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index caf3699..45a3955 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -212,7 +212,6 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
 class (Storable a, Floating a) => Element a where
     constantD :: a -> Int -> Vector a
     transdata :: Int -> Vector a -> Int -> Vector a
-    multiplyD :: Matrix a -> Matrix a -> Matrix a
     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
                -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
                -> Matrix a -> Matrix a
@@ -221,14 +220,12 @@ class (Storable a, Floating a) => Element a where
 instance Element Double where
     constantD  = constantR
     transdata = transdataR
-    multiplyD  = multiplyR
     subMatrixD = subMatrixR
     diagD      = diagR
 
 instance Element (Complex Double) where
     constantD  = constantC
     transdata  = transdataC
-    multiplyD  = multiplyC
     subMatrixD = subMatrixC
     diagD      = diagC
 
@@ -266,33 +263,6 @@ transdataAux fun c1 d c2 =
 foreign import ccall "auxi.h transR" ctransR :: TMM
 foreign import ccall "auxi.h transC" ctransC :: TCMCM
 
-------------------------------------------------------------------
-
-gmatC MF { rows = r, cols = c } p f = f 1 (fi c) (fi r) p
-gmatC MC { rows = r, cols = c } p f = f 0 (fi r) (fi c) p
-
-dtt MC { cdat = d } = d
-dtt MF { fdat = d } = d
-
-multiplyAux fun a b = unsafePerformIO $ do
-    when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++
-                                      show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
-    r <- createMatrix RowMajor (rows a) (cols b)
-    withForeignPtr (fptr (dtt a)) $ \pa -> withForeignPtr (fptr (dtt b)) $ \pb ->
-        withMatrix r $ \r' ->
-            fun // gmatC a pa // gmatC b pb // r' // check "multiplyAux"
-    return r
-
-multiplyR = multiplyAux cmultiplyR
-foreign import ccall "auxi.h multiplyR" cmultiplyR :: TauxMul Double
-
-multiplyC = multiplyAux cmultiplyC
-foreign import ccall "auxi.h multiplyC" cmultiplyC :: TauxMul (Complex Double)
-
--- | matrix product
-multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a
-multiply = multiplyD
-
 ----------------------------------------------------------------------
 
 -- | extraction of a submatrix from a real matrix
@@ -370,7 +340,12 @@ constant = constantD
 
 -- | obtains the complex conjugate of a complex vector
 conj :: Vector (Complex Double) -> Vector (Complex Double)
-conj v = asComplex $ flatten $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1])
+conj v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 cconjugate vec v vec r "cconjugate"
+    return r
+foreign import ccall "auxi.h conjugate" cconjugate :: TCVCV
+
 
 -- | creates a complex vector from vectors with real and imaginary parts
 toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index 7f83bcf..04dc7ad 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -4,14 +4,9 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>
-#include <gsl/gsl_fft_complex.h>
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_deriv.h>
-#include <gsl/gsl_poly.h>
-#include <gsl/gsl_multimin.h>
 #include <gsl/gsl_complex.h>
 #include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_cblas.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -118,78 +113,6 @@ int constantC(gsl_complex* pval, CVEC(r)) {
 }
 
 
-int multiplyR(int ta, KRMAT(a), int tb, KRMAT(b),RMAT(r)) {
-    //printf("%d %d %d %d %d %d\n",ar,ac,br,bc,rr,rc);
-    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("multiplyR (gsl_blas_dgemm)");
-    KDMVIEW(a);
-    KDMVIEW(b);
-    DMVIEW(r);
-    int k;
-    for(k=0;k<rr*rc;k++) rp[k]=0;
-    int debug = 0;
-    if(debug) {
-        printf("---------------------------\n");
-        printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
-    }
-    int res = gsl_blas_dgemm(
-         ta?CblasTrans:CblasNoTrans,
-         tb?CblasTrans:CblasNoTrans,
-         1.0, M(a), M(b),
-         0.0, M(r));
-    if(debug) {
-        printf("--------------\n");
-        printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
-    }
-    CHECK(res,res);
-    OK
-}
-
-int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r)) {
-    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("multiplyC (gsl_blas_zgemm)");
-    KCMVIEW(a);
-    KCMVIEW(b);
-    CMVIEW(r);
-    int k;
-    gsl_complex alpha, beta;
-    GSL_SET_COMPLEX(&alpha,1.,0.);
-    GSL_SET_COMPLEX(&beta,0.,0.);
-    //double *TEMP = (double*)malloc(rr*rc*2*sizeof(double));
-    //gsl_matrix_complex_view T = gsl_matrix_complex_view_array(TEMP,rr,rc);
-    for(k=0;k<rr*rc;k++) rp[k]=beta;
-    //for(k=0;k<2*rr*rc;k++) TEMP[k]=0;
-    int debug = 0;
-    if(debug) {
-        printf("---------------------------\n");
-        printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
-        //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
-    }
-    int res = gsl_blas_zgemm(
-         ta?CblasTrans:CblasNoTrans,
-         tb?CblasTrans:CblasNoTrans,
-         alpha, M(a), M(b),
-         beta, M(r)); 
-         //&T.matrix);
-    //memcpy(rp,TEMP,2*rr*rc*sizeof(double));
-    if(debug) {
-        printf("--------------\n");
-        printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
-        //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
-    }
-    CHECK(res,res);
-    OK
-}
-
-
 int diagR(KRVEC(d),RMAT(r)) {
     REQUIRES(dn==rr && rr==rc,BAD_SIZE);
     DEBUGMSG("diagR");
@@ -215,3 +138,14 @@ int diagC(KCVEC(d),CMAT(r)) {
     }
     OK
 }
+
+int conjugate(KCVEC(x),CVEC(t)) {
+    REQUIRES(xn==tn,BAD_SIZE);
+    DEBUGMSG("conjugate");
+    int k;
+    for (k=0; k<xn; k++) {
+        tp[k].dat[0] =   xp[k].dat[0];
+        tp[k].dat[1] = - xp[k].dat[1];
+    }
+    OK
+}
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
index 73334e3..377a4a1 100644
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -10,16 +10,12 @@
 #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 transR(KRMAT(x),RMAT(t));
 int transC(KCMAT(x),CMAT(t));
 
 int constantR(double *val     , RVEC(r));
 int constantC(gsl_complex *val, CVEC(r));
 
-int multiplyR(int ta, KRMAT(a), int tb,  KRMAT(b),RMAT(r));
-int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r));
-
 int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));
 
 int diagR(KRVEC(d),RMAT(r));
@@ -28,3 +24,5 @@ int diagC(KCVEC(d),CMAT(r));
 const char * gsl_strerror (const int gsl_errno);
 
 int matrix_fscanf(char*filename, RMAT(a));
+
+int conjugate(KCVEC(x),CVEC(t));