From f2cf177e93d4578b404909c68b24625a76466ee5 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Sat, 9 Jun 2007 16:29:46 +0000
Subject: eigC

---
 examples/tests.hs                  | 20 ++++++++++++++++++++
 lib/Data/Packed/Internal/Matrix.hs |  2 ++
 lib/LAPACK/Internal.hs             |  4 +++-
 3 files changed, 25 insertions(+), 1 deletion(-)

diff --git a/examples/tests.hs b/examples/tests.hs
index f167b92..5af33ba 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -66,6 +66,8 @@ aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length
 (|~|) = aprox abs
 (|~~|) = aprox magnitude
 
+v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2
+
 eps = 1E-8::Double
 
 asFortran m = (rows m >|< cols m) $ toList (fdat m)
@@ -116,6 +118,15 @@ instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where
         --return $ PairM ((a><b) l1) ((b><c) l2)
     coarbitrary = undefined
 
+data SqM a = SqM (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (SqM a) where
+    arbitrary = do
+        n <- choose (1,10)
+        l <- vector (n*n)
+        return $ SqM $ (n><n) l
+    coarbitrary = undefined
+
+
 type BaseType = Double
 
 
@@ -133,6 +144,14 @@ svdTestC fun prod m = u <> s' <> (trans v) |~~| m
           (<>) = prod
           s' = liftMatrix comp s
 
+eigTestC fun prod (SqM m) = (m <> v) |~~| (v <> diag s)
+                        && takeDiag ((liftMatrix conj (trans v)) `mulC` v) ~~ constant (rows m) 1
+    where (s,v) = fun m
+          (<>) = prod
+
+takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
+
+
 comp v = toComplex (v,constant (dim v) 0)
 
 main = do
@@ -147,3 +166,4 @@ main = do
     quickCheck (svdTestR svdR mulF)
     quickCheck (svdTestC svdC mulC)
     quickCheck (svdTestC svdC mulF)
+    quickCheck (eigTestC eigC mulC)
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 0664e9b..b8de245 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -108,6 +108,8 @@ createMatrix order r c = do
 
 reshape c v = matrixFromVector RowMajor c v
 
+singleton x = reshape 1 (fromList [x])
+
 transdataG :: Storable a => Int -> Vector a -> Int -> Vector a 
 transdataG c1 d c2 = fromList . concat . transpose . partit c1 . toList $ d
 
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs
index e39dd10..e3c9927 100644
--- a/lib/LAPACK/Internal.hs
+++ b/lib/LAPACK/Internal.hs
@@ -76,7 +76,9 @@ foreign import ccall "lapack-aux.h eig_l_C"
 -- The eigenvectors are the columns of v.
 -- The eigenvalues are not sorted.
 eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
-eigC (m@M {rows = r}) = unsafePerformIO $ do
+eigC (m@M {rows = r}) 
+    | r == 1 = (fromList [cdat m `at` 0], singleton 1)
+    | otherwise = unsafePerformIO $ do
     l <- createVector r
     v <- createMatrix ColumnMajor r r
     dummy <- createMatrix ColumnMajor 1 1
-- 
cgit v1.2.3