3 files changed, 101 insertions, 14 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index b19c0ec..79cc64d 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -69,7 +69,7 @@ import Data.List(foldl1')
 import Data.Array
 
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
-class (Normed (Matrix t), Linear Matrix t) => Field t where
+class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     -- | Singular value decomposition using lapack's dgesvd or zgesvd.
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     luPacked    :: Matrix t -> (Matrix t, [Int])
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index a39df50..0ddbb55 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Linear
@@ -60,18 +60,7 @@ instance Linear Vector (Complex Double) where
     divide = vectorZipC Div
     equal u v = dim u == dim v && vectorMax (liftVector magnitude (sub u v)) == 0.0
 
-instance Linear Matrix Double where
-    scale x = liftMatrix (scale x)
-    scaleRecip x = liftMatrix (scaleRecip x)
-    addConstant x = liftMatrix (addConstant x)
-    add = liftMatrix2 add
-    sub = liftMatrix2 sub
-    mul = liftMatrix2 mul
-    divide = liftMatrix2 divide
-    equal a b = cols a == cols b && flatten a `equal` flatten b
-
-
-instance Linear Matrix (Complex Double) where
+instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where
     scale x = liftMatrix (scale x)
     scaleRecip x = liftMatrix (scaleRecip x)
     addConstant x = liftMatrix (addConstant x)
diff --git a/lib/Numeric/LinearAlgebra/Testing.hs b/lib/Numeric/LinearAlgebra/Testing.hs
new file mode 100644
index 0000000..dcf1d8e
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Testing.hs
@@ -0,0 +1,98 @@
+{-# OPTIONS_GHC -XPatternSignatures #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Testing
+Copyright   :  (c) Alberto Ruiz 2007
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Some consistency tests.
+
+-}
+
+module Numeric.LinearAlgebra.Testing(
+    runTests, runBigTests
+) where
+
+import Numeric.LinearAlgebra
+import Test.QuickCheck
+import Debug.Trace
+
+qCheck n = check defaultConfig {configSize = const n}
+
+debug x = trace (show x) x
+
+type RM = Matrix Double
+type CM = Matrix (Complex Double)
+
+instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
+    arbitrary = do
+        r <- arbitrary
+        i <- arbitrary
+        return (r:+i)
+    coarbitrary = undefined
+
+chooseDim = sized $ \m -> choose (1,max 1 m)
+
+instance (Field a, Arbitrary a) => Arbitrary (Vector a) where 
+   arbitrary = do m <- chooseDim
+                  l <- vector m
+                  return $ fromList l
+   coarbitrary = undefined
+
+instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 
+    arbitrary = do
+        m <- chooseDim
+        n <- chooseDim
+        l <- vector (m*n)
+        return $ (m><n) l
+    coarbitrary = undefined
+
+
+newtype (Sq a) = Sq (Matrix a) deriving Show
+sq (Sq m) = m
+
+instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
+    arbitrary = do
+        n <- chooseDim
+        l <- vector (n*n)
+        return $ Sq $ (n><n) l
+    coarbitrary = undefined
+
+newtype (Rot a) = Rot (Matrix a) deriving Show
+rot (Rot a) = a
+
+instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let (q,_) = qr m
+        return (Rot q)
+    coarbitrary = undefined
+
+newtype (Her a) = Her (Matrix a) deriving Show
+her (Her a) = a
+
+instance (Field a, Arbitrary a) => Arbitrary (Her a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let m' = m/2
+        return $ Her (m' + ctrans m')
+    coarbitrary = undefined
+
+-------------------------------------------------------------------
+
+herR x = her x :: RM
+
+-- | It runs all the tests.
+runTests :: Int  -- ^ maximum dimension
+         -> IO ()
+runTests n = do
+    qCheck n (\(Her (m::CM))-> m == ctrans m)
+    qCheck n $ (\m->m==ctrans m) . herR
+
+-- | Some additional tests on big matrices. They take a few minutes.
+runBigTests :: IO ()
+runBigTests = undefined