optimized conjugate

6 files changed, 71 insertions, 22 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 5892f1c..d39481d 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -255,27 +255,27 @@ class (Storable a, Floating a) => Element a where
     transdata = transdata'
     constantD  :: a -> Int -> Vector a
     constantD = constant'
-    ctrans' :: Matrix a -> Matrix a
+    conjugateD :: Vector a -> Vector a
 
 instance Element Float where
     transdata  = transdataAux ctransF
     constantD  = constantAux cconstantF
-    ctrans'    = trans
+    conjugateD = id
 
 instance Element Double where
     transdata  = transdataAux ctransR
     constantD  = constantAux cconstantR
-    ctrans'    = trans
+    conjugateD = id
 
 instance Element (Complex Float) where
     transdata  = transdataAux ctransQ
     constantD  = constantAux cconstantQ
-    ctrans'    = liftMatrix (mapVector conjugate) . trans
+    conjugateD = conjugateQ
 
 instance Element (Complex Double) where
     transdata  = transdataAux ctransC
     constantD  = constantAux cconstantC
-    ctrans'    = liftMatrix (mapVector conjugate) . trans
+    conjugateD = conjugateC
 
 -------------------------------------------------------------------
 
@@ -359,6 +359,21 @@ constantC :: Complex Double -> Int -> Vector (Complex Double)
 constantC = constantAux cconstantC
 foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV
 
+---------------------------------------
+
+conjugateAux fun x = unsafePerformIO $ do
+    v <- createVector (dim x)
+    app2 fun vec x vec v "conjugateAux"
+    return v
+
+conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
+conjugateQ = conjugateAux c_conjugateQ
+foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV
+
+conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
+conjugateC = conjugateAux c_conjugateC
+foreign import ccall "conjugateC" c_conjugateC :: TCVCV
+
 ----------------------------------------------------------------------
 
 -- | Extracts a submatrix from a matrix.
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index af937f4..0fc7876 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -270,7 +270,7 @@ asColumn v = reshape 1 v
 
 -- | conjugate transpose
 ctrans :: Element e => Matrix e -> Matrix e
-ctrans = ctrans'
+ctrans = liftMatrix conjugateD . trans
 
 
 {- | creates a Matrix of the specified size using the supplied function to
diff --git a/lib/Numeric/Conversion.hs b/lib/Numeric/Conversion.hs
index b05069c..809ac51 100644
--- a/lib/Numeric/Conversion.hs
+++ b/lib/Numeric/Conversion.hs
@@ -69,7 +69,6 @@ class ComplexContainer c where
     fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e)
     comp        :: (RealElement e) => c e -> c (Complex e)
     conj        :: (RealElement e) => c (Complex e) -> c (Complex e)
---    cmap        :: (Element a, Element b) => (a -> b) -> c a -> c b
     single'      :: Precision a b => c b -> c a
     double'      :: Precision a b => c a -> c b
 
@@ -78,16 +77,11 @@ instance ComplexContainer Vector where
     toComplex = toComplexV
     fromComplex = fromComplexV
     comp v = toComplex (v,constantD 0 (dim v))
-    conj = conjV
---    cmap = mapVector
+    conj = conjugateD
     single' = double2FloatG
     double' = float2DoubleG
 
 
--- | obtains the complex conjugate of a complex vector
-conjV :: (RealElement a) => Vector (Complex a) -> Vector (Complex a)
-conjV = mapVector conjugate
-
 -- | creates a complex vector from vectors with real and imaginary parts
 toComplexV :: (RealElement a) => (Vector a, Vector a) ->  Vector (Complex a)
 toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
@@ -104,7 +98,6 @@ instance ComplexContainer Matrix where
         where c = cols z
     comp = liftMatrix comp
     conj = liftMatrix conj
---    cmap f = liftMatrix (cmap f)
     single' = liftMatrix single'
     double' = liftMatrix double'
 
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 9e44431..2c4c647 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -1201,3 +1201,28 @@ int double2float(DVEC(x),FVEC(y)) {
     }
     OK
 }
+
+//////////////////// conjugate /////////////////////////
+
+int conjugateQ(KQVEC(x),QVEC(t)) {
+    REQUIRES(xn==tn,BAD_SIZE);
+    DEBUGMSG("conjugateQ");
+    int k;
+    for(k=0;k<xn;k++) {
+        ((complex*)tp)[k].r=((complex*)xp)[k].r;
+        ((complex*)tp)[k].i=-((complex*)xp)[k].i;
+    }
+    OK
+}
+
+int conjugateC(KCVEC(x),CVEC(t)) {
+    REQUIRES(xn==tn,BAD_SIZE);
+    DEBUGMSG("conjugateC");
+    int k;
+    for(k=0;k<xn;k++) {
+        ((doublecomplex*)tp)[k].r=((doublecomplex*)xp)[k].r;
+        ((doublecomplex*)tp)[k].i=-((doublecomplex*)xp)[k].i;
+    }
+    OK
+}
+
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 0b4e3bf..426700b 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -271,6 +271,10 @@ normsMTest = TestList [
 
 ---------------------------------------------------------------------
 
+conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m)
+
+---------------------------------------------------------------------
+
 
 -- | All tests must pass with a maximum dimension of about 20
 --  (some tests may fail with bigger sizes due to precision loss).
@@ -292,6 +296,9 @@ runTests n = do
     putStrLn "------ sub-trans"
     test (subProp . rM)
     test (subProp . cM)
+    putStrLn "------ ctrans"
+    test (conjuTest . cM)
+    test (conjuTest . zM)
     putStrLn "------ lu"
     test (luProp    . rM)
     test (luProp    . cM)
@@ -362,6 +369,9 @@ runTests n = do
     test (qrProp     . cM)
     test (rqProp     . rM)
     test (rqProp     . cM)
+    test (rqProp1     . cM)
+    test (rqProp2     . cM)
+    test (rqProp3     . cM)
     putStrLn "------ hess"
     test (hessProp   . rSq)
     test (hessProp   . cSq)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index b96f53e..e780c35 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -32,7 +32,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
     eigProp, eigSHProp, eigProp2, eigSHProp2,
-    qrProp, rqProp,
+    qrProp, rqProp, rqProp1, rqProp2, rqProp3,
     hessProp,
     schurProp1, schurProp2,
     cholProp,
@@ -210,15 +210,21 @@ eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
 qrProp m = q <> r |~| m && unitary q && upperTriang r
     where (q,r) = qr m
 
-rqProp m = r <> q |~| m && unitary q && utr
+rqProp m = r <> q |~| m && unitary q && upperTriang' r
     where (r,q) = rq m
-          upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
-              where t = f-c
-          utr = upptr (rows r) (cols r) * r |~| r
 
-upperTriang' m = rows m == 1 || down |~| z
-    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
-          z = constant 0 (dim down)
+rqProp1 m = r <> q |~| m
+    where (r,q) = rq m
+
+rqProp2 m = unitary q
+    where (r,q) = rq m
+
+rqProp3 m = upperTriang' r
+    where (r,q) = rq m
+
+upperTriang' r = upptr (rows r) (cols r) * r |~| r
+    where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
+              where t = f-c
 
 hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
     where (p,h) = hess m