algorithms refactoring

8 files changed, 137 insertions, 634 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index b841588..22ac658 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -14,7 +14,7 @@ tested-with:        GHC ==6.6.1
 Build-Depends:      base, haskell98
 Extensions:         ForeignFunctionInterface
 --ghc-options:        -Wall
-ghc-options:        -O0
+ghc-options:        -O
 hs-source-dirs:     lib
 Exposed-modules:    Data.Packed.Internal,
                     Data.Packed.Internal.Common,
@@ -51,6 +51,7 @@ Exposed-modules:    Data.Packed.Internal,
                     LinearAlgebra.Instances,
                     LinearAlgebra.Interface,
                     LinearAlgebra.Algorithms
+--                    , LinearAlgebra.Tests
                     , Graphics.Plot
 --                    , GSLHaskell
 Other-modules:
diff --git a/examples/oldtests.hs b/examples/oldtests.hs
deleted file mode 100644
index 7d4701c..0000000
--- a/examples/oldtests.hs
+++ /dev/null
@@ -1,120 +0,0 @@
-import Test.HUnit
-import LinearAlgebra
-import GSL hiding (exp)
-import System.Random(randomRs,mkStdGen)
-
-realMatrix = fromLists :: [[Double]] -> Matrix Double
-realVector = fromList ::  [Double] -> Vector Double
-
-
-
-infixl 2 =~=
-a =~= b = pnorm PNorm1 (flatten (a - b)) < 1E-6
-
-randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed 
-
-randomMatrixC seed (n,m) = toComplex (randomMatrix seed (n,m), randomMatrix (seed+1) (n,m))
-
-besselTest = do
-    let (r,e) = bessel_J0_e 5.0
-    let expected = -0.17759677131433830434739701
-    assertBool "bessel_J0_e" ( abs (r-expected) < e ) 
-
-exponentialTest = do
-    let (v,e,err) = exp_e10_e 30.0
-    let expected = exp 30.0
-    assertBool "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) 
-
-disp m = putStrLn (format " " show m)
-
-ms = realMatrix [[1,2,3]
-                ,[-4,1,7]]
-
-ms' = randomMatrix 27 (50,100)
-
-ms'' = toComplex (randomMatrix 100 (50,100),randomMatrix 101 (50,100))
-
-fullsvdTest method mat msg = do
-    let (u,s,vt) = method mat
-    assertBool msg (u <> s <> trans vt =~= mat)
-
-svdg' m = (u, diag s, v) where (u,s,v) = svdg m
-
-full_svd_Rd = svdRdd
-
---------------------------------------------------------------------
-
-mcu = toComplex (randomMatrix 33 (20,20),randomMatrix 34 (20,20))
-
-mcur = randomMatrix 35 (40,40)
-
--- eigenvectors are columns
-eigTest method m msg = do
-    let (s,v) = method m
-    assertBool msg $ m <> v =~= v <> diag s
-
-bigmat = m + trans m where m = randomMatrix 18 (1000,1000)
-bigmatc = mc + conjTrans mc where mc = toComplex(m,m)
-                                  m = randomMatrix 19 (1000,1000)
-
---------------------------------------------------------------------
-
-invTest msg m = do
-    assertBool msg $ m <> inv m =~= ident (rows m)
-
-invComplexTest msg m = do
-    assertBool msg $ m <> invC m =~= identC (rows m)
-
-invC m = linearSolveC m (identC (rows m))
-
-identC = comp . ident
-
---------------------------------------------------------------------
-
-pinvTest f msg m = do
-    assertBool msg $ m <> f m <> m =~= m
-
-pinvC m = linearSolveLSC m (identC (rows m))
-
-pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
-
-pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m))
-
---------------------------------------------------------------------
-
-
-tests = TestList [
-      TestCase $ besselTest
-    , TestCase $ exponentialTest
-    , TestCase $ invTest "inv 100x100" (randomMatrix 18 (100,100))
-    , TestCase $ invComplexTest "complex inv 100x100" (randomMatrixC 18 (100,100))
-    , TestCase $ pinvTest (pinvTolg 1) "pinvg 100x50" (randomMatrix 18 (100,50))
-    , TestCase $ pinvTest pinv "pinv 100x50" (randomMatrix 18 (100,50))
-    , TestCase $ pinvTest pinv "pinv 50x100" (randomMatrix 18 (50,100))
-    , TestCase $ pinvTest pinvSVDR "pinvSVDR 100x50" (randomMatrix 18 (100,50))
-    , TestCase $ pinvTest pinvSVDR "pinvSVDR 50x100" (randomMatrix 18 (50,100))
-    , TestCase $ pinvTest pinvC "pinvC 100x50" (randomMatrixC 18 (100,50))
-    , TestCase $ pinvTest pinvC "pinvC 50x100" (randomMatrixC 18 (50,100))
-    , TestCase $ pinvTest pinvSVDC "pinvSVDC 100x50" (randomMatrixC 18 (100,50))
-    , TestCase $ pinvTest pinvSVDC "pinvSVDC 50x100" (randomMatrixC 18 (50,100))
-    , TestCase $ eigTest eigC mcu "eigC" 
-    , TestCase $ eigTest eigR mcur "eigR"
-    , TestCase $ eigTest eigS (mcur+trans mcur) "eigS"
-    , TestCase $ eigTest eigSg (mcur+trans mcur) "eigSg"
-    , TestCase $ eigTest eigH (mcu+ (conjTrans) mcu) "eigH"
-    , TestCase $ eigTest eigHg (mcu+ (conjTrans) mcu) "eigHg"
-    , TestCase $ fullsvdTest svdg' ms "GSL svd small"
-    , TestCase $ fullsvdTest svdR ms "fullsvdR small"
-    , TestCase $ fullsvdTest svdR (trans ms) "fullsvdR small"
-    , TestCase $ fullsvdTest svdR ms' "fullsvdR"
-    , TestCase $ fullsvdTest svdR (trans ms') "fullsvdR"
-    , TestCase $ fullsvdTest full_svd_Rd ms' "fullsvdRd"
-    , TestCase $ fullsvdTest full_svd_Rd (trans ms') "fullsvdRd"
-    , TestCase $ fullsvdTest svdC ms'' "fullsvdC"
-    , TestCase $ fullsvdTest svdC (trans ms'') "fullsvdC"
-    , TestCase $ eigTest eigS bigmat "big eigS"
-    , TestCase $ eigTest eigH bigmatc "big eigH"
-    , TestCase $ eigTest eigR bigmat "big eigR"
-    ]
-
-main = runTestTT tests
diff --git a/examples/tests.hs b/examples/tests.hs
deleted file mode 100644
index dcc3cbf..0000000
--- a/examples/tests.hs
+++ /dev/null
@@ -1,357 +0,0 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
-
---
--- QuickCheck tests
---
-
------------------------------------------------------------------------------
-
-import Data.Packed.Internal((>|<), fdat, cdat, multiply', multiplyG, MatrixOrder(..))
-import GSL hiding (sin,cos,exp,choose)
-import LinearAlgebra hiding ((<>))
-import Test.QuickCheck
-import Test.HUnit hiding ((~:))
-
-
-dist :: (Normed t, Num t) => t -> t -> Double
-dist a b = pnorm Infinity (a-b)
-
-infixl 4 |~|
-a |~| b = a :~8~: b
-
-data Aprox a = (:~) a Int
-
-(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
-a :~n~: b = dist a b < 10^^(-n)
-
-
-{-
--- Bravo por quickCheck!
-
-pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m  ~~ ident (cols m)
-    where infix 2 ~~
-          (~~) = approxEqual tol 
-
-pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m  ~~ m)
-    where r = rank m
-          c = cols m
-          infix 2 ~~
-          (~~) = approxEqual tol 
-
-nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros
-    where nt    = trans (nullspace m)
-          cr    = corank m
-          r     = rows m
-          zeros = create [r,cr] $ replicate (r*cr) 0  
-
--}
-
-ac = (2><3) [1 .. 6::Double]
-bc = (3><4) [7 .. 18::Double]
-
-mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)]
-
-af = (2>|<3) [1,4,2,5,3,6::Double]
-bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
-
-
-{-
-aprox fun a b = rows a == rows b &&
-          cols a == cols b &&
-          epsTol > aproxL fun (toList (t a)) (toList (t b))
-    where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat
-
-aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1)
-
-normVR a b = toScalarR AbsSum (vectorZipR Sub a b)
-
-a |~| b = rows a == rows b && cols a == cols b && epsTol > normVR (t a) (t b)
-    where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat
-
-(|~~|) = aprox magnitude
-
-v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2
-
-u ~|~ v = normVR u v < epsTol
--}
-
-epsTol = 1E-8::Double
-
-asFortran m = (rows m >|< cols m) $ toList (fdat m)
-asC m = (rows m >< cols m) $ toList (cdat m)
-
-mulC a b = multiply' RowMajor a b
-mulF a b = multiply' ColumnMajor a b
-
-identC = comp . ident
-
-infixl 7 <>
-a <> b = mulF a b
-
-cc = mulC ac bf
-cf = mulF af bc
-
-r = mulC cc (trans cf)
-
-rd = (2><2)
- [ 27736.0,  65356.0
- , 65356.0, 154006.0 ::Double]
-
-instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
-    arbitrary = do
-        r <- arbitrary
-        i <- arbitrary
-        return (r:+i)
-    coarbitrary = undefined
-
-instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where 
-   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
-                  m <- choose (1,10)
-                  n <- choose (1,10)
-                  l <- vector (m*n)
-                  ctype <- arbitrary
-                  let h = if ctype then (m><n) else (m>|<n)
-                  trMode <- arbitrary
-                  let tr = if trMode then trans else id
-                  return $ tr (h l)
-   coarbitrary = undefined
-
-data PairM a = PairM (Matrix a) (Matrix a) deriving Show
-instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where
-    arbitrary = do
-        a <- choose (1,10)
-        b <- choose (1,10)
-        c <- choose (1,10)
-        l1 <- vector (a*b)
-        l2 <- vector (b*c)
-        return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int])))
-        --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
-
-data Sym a = Sym (Matrix a) deriving Show
-instance (Linear Vector a, Arbitrary a) => Arbitrary (Sym a) where
-    arbitrary = do
-        SqM m <- arbitrary
-        return $ Sym (m + trans m)
-    coarbitrary = undefined
-
-data Her = Her (Matrix (Complex Double)) deriving Show
-instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where
-    arbitrary = do
-        SqM m <- arbitrary
-        return $ Her (m + conjTrans m)
-    coarbitrary = undefined
-
-data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show
-instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where
-    arbitrary = do
-        a <- choose (1,10)
-        c <- choose (1,10)
-        l1 <- vector (a*a)
-        l2 <- vector (a*c)
-        return $ PairSM ((a><a) (map fromIntegral (l1::[Int]))) ((a><c) (map fromIntegral (l2::[Int])))
-        --return $ PairSM ((a><a) l1) ((a><c) l2)
-    coarbitrary = undefined
-
-instance (Field a, Arbitrary a) => Arbitrary (Vector a) where 
-   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
-                  m <- choose (1,100)
-                  l <- vector m
-                  return $ fromList l
-   coarbitrary = undefined
-
-data PairV a = PairV (Vector a) (Vector a)
-instance (Field a, Arbitrary a) => Arbitrary (PairV a) where 
-   arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)
-                  m <- choose (1,100)
-                  l1 <- vector m
-                  l2 <- vector m
-                  return $ PairV (fromList l1) (fromList l2)
-   coarbitrary = undefined
-
-
-
-type BaseType = Complex Double
-
-svdTestR fun m = u <> s <> trans v |~| m
-               && u <> trans u |~| ident (rows m)
-               && v <> trans v |~| ident (cols m)
-    where (u,s,v) = fun m
-
-
-svdTestC m = u <> s' <> (trans v) |~| m
-                  && u <> conjTrans u |~| identC (rows m)
-                  && v <> conjTrans v |~| identC (cols m)
-    where (u,s,v) = svdC m
-          s' = liftMatrix comp s
-
---svdg' m = (u,s',v) where 
-
-eigTestC (SqM m) = (m <> v) |~| (v <> diag s)
-                        && takeDiag (conjTrans v <> v) |~| comp (constant 1 (rows m)) --normalized
-    where (s,v) = eigC m
-
-eigTestR (SqM m) = (liftMatrix comp m <> v) |~| (v <> diag s)
-                      -- && takeDiag ((liftMatrix conj (trans v)) <> v) |~| constant 1 (rows m) --normalized ???
-    where (s,v) = eigR m
-
-eigTestS (Sym m) = (m <> v) |~| (v <> diag s)
-                       && v <> trans v |~| ident (cols m)
-    where (s,v) = eigS m
-
-eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s))
-                        && v <> conjTrans v |~| identC (cols m)
-    where (s,v) = eigH m
-
-linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b
-
-prec = 1E-15
-
-singular fun m = s1 < prec || s2/s1 < prec
-    where (_,ss,v) = fun m
-          s = toList ss
-          s1 = maximum s
-          s2 = minimum s
-
-{-
-invTest msg m = do
-    assertBool msg $ m <> inv m =~= ident (rows m)
-
-invComplexTest msg m = do
-    assertBool msg $ m <> invC m =~= identC (rows m)
-
-invC m = linearSolveC m (identC (rows m))
-
-identC n = toComplex(ident n, (0::Double) <>ident n)
--}
-
---------------------------------------------------------------------
-
-pinvTest f m = (m <> f m <> m) |~| m
-
-pinvR m = linearSolveLSR m (ident (rows m))
-pinvC m = linearSolveLSC m (identC (rows m))
-
-pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
-
-pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m))
-
---------------------------------------------------------------------
-
-polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
-
-polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
-
-
-polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4])
-
----------------------------------------------------------------------
-
-quad f a b = fst $ integrateQAGS 1E-9 100 f a b  
-
--- A multiple integral can be easily defined using partial application
-quad2 f a b g1 g2 = quad h a b
-    where h x = quad (f x) (g1 x) (g2 x)
-
-volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
-                        0 r (const 0) (\x->sqrt (r*r-x*x))
-
-integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)
-
-
----------------------------------------------------------------------
-
-arit1 u = sin u ^ 2 + cos u ^ 2 |~| 1
-    where _ = u :: Vector Double
-
-arit2 u = sin u ** 2 + cos u ** 2 |~| 1
-    where _ = u :: Vector Double
-
-arit3 u = cos u * tan u |~| sin u
-    where _ = u :: Vector Double
-
-arit4 u = (cos u * tan u) :~6~: sin u
-    where _ = u :: Vector (Complex Double)
-
----------------------------------------------------------------------
-
-besselTest = do
-    let (r,e) = bessel_J0_e 5.0
-    let expected = -0.17759677131433830434739701
-    assertBool "bessel_J0_e" ( abs (r-expected) < e ) 
-
-exponentialTest = do
-    let (v,e,err) = exp_e10_e 30.0
-    let expected = exp 30.0
-    assertBool "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) 
-
-gammaTest = do
-    assertBool "gamma" (gamma 5 == 24.0)
-
-tests = TestList
-    [ TestCase $ besselTest
-    , TestCase $ exponentialTest
-    , TestCase $ gammaTest
-    , TestCase $ polySolveTest
-    , TestCase $ integrateTest
-    ]
-
-----------------------------------------------------------------------
-
-main = do
-    putStrLn "--------- general -----"
-    quickCheck (\(Sym m) -> m == (trans m:: Matrix BaseType))    
-    quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType])
-
-    quickCheck $ \m -> m == asC (m :: Matrix BaseType)
-    quickCheck $ \m -> m == asFortran (m :: Matrix BaseType)
-    quickCheck $ \m -> m == (asC . asFortran) (m :: Matrix BaseType)
-    putStrLn "--------- MULTIPLY ----"
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: Matrix BaseType)
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: Matrix BaseType))
-    quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: Matrix BaseType)
-    putStrLn "--------- SVD ---------"
-    quickCheck (svdTestR svdR)
-    quickCheck (svdTestR svdRdd)
---    quickCheck (svdTestR svdg)
-    quickCheck svdTestC
-    putStrLn "--------- EIG ---------"
-    quickCheck eigTestC
-    quickCheck eigTestR
-    quickCheck eigTestS
-    quickCheck eigTestH
-    putStrLn "--------- SOLVE ---------"
-    quickCheck (linearSolveSQTest linearSolveR (singular svdR'))
-    quickCheck (linearSolveSQTest linearSolveC (singular svdC'))
-    quickCheck (pinvTest pinvR)
-    quickCheck (pinvTest pinvC)
-    quickCheck (pinvTest pinvSVDR)
-    quickCheck (pinvTest pinvSVDC)
-    putStrLn "--------- VEC OPER ------"
-    quickCheck arit1
-    quickCheck arit2
-    quickCheck arit3
-    quickCheck arit4
-    putStrLn "--------- GSL ------"
-    runTestTT tests
-    quickCheck $ \v -> ifft (fft v) |~| v
-
-kk = (2><2)
- [  1.0, 0.0
- , -1.5, 1.0 ::Double]
-
-v = 11 |> [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double]
-
-pol =  [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]
-
-mm = (2><2)
- [ 0.5, 0.0
- , 0.0, 0.0 ] :: Matrix Double
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index f0ef8b6..ebe6371 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -24,6 +24,7 @@ import Data.List(transpose)
 import Debug.Trace(trace)
 import Foreign.C.String(peekCString)
 import Foreign.C.Types
+import Data.Monoid
 
 -- | A one-dimensional array of objects stored in a contiguous memory block.
 data Vector t = V { dim  :: Int            -- ^ number of elements
@@ -177,3 +178,4 @@ liftVector  f = fromList . map f . toList
 liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
 liftVector2 f u v = fromList $ zipWith f (toList u) (toList v)
 
+-----------------------------------------------------------------
diff --git a/lib/GSL.hs b/lib/GSL.hs
index 8b6365b..d65f8ff 100644
--- a/lib/GSL.hs
+++ b/lib/GSL.hs
@@ -8,7 +8,7 @@ Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses -fffi and -fglasgow-exts
 
-This module reexports all the GSL functions (except those in "LinearAlgebra").
+This module reexports all the available GSL functions (except those in "LinearAlgebra").
 
 -}
 
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index 2b92a2a..54eea8a 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -14,7 +14,7 @@
 -----------------------------------------------------------------------------
 
 module LAPACK (
-    svdR, svdR', svdRdd, svdRdd', svdC, svdC',
+    svdR, svdRdd, svdC,
     eigC, eigR, eigS, eigH,
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
@@ -26,7 +26,6 @@ import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
 import Data.Packed.Vector
 import Data.Packed.Matrix
---import LinearAlgebra.Linear(scale)
 import GSL.Vector(vectorMapValR, FunCodeSV(Scale))
 import Complex
 import Foreign
@@ -36,14 +35,9 @@ foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM
 
 -- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
 --
--- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
-svdR :: Matrix Double -> (Matrix Double, Matrix Double, Matrix Double)
-svdR x = (u, diagRect s r c, v) where (u,s,v) = svdR' x
-                                      r = rows x
-                                      c = cols x
-
-svdR' :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdR' x = unsafePerformIO $ do
+-- @(u,s,v)=full svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
+svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+svdR x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
@@ -56,14 +50,9 @@ foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM
 
 -- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
 --
--- @(u,s,v)=svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@.
-svdRdd :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
-svdRdd x = (u, diagRect s r c, v) where (u,s,v) = svdRdd' x
-                                        r = rows x
-                                        c = cols x
-
-svdRdd' :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdRdd' x = unsafePerformIO $ do
+-- @(u,s,v)=full svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@.
+svdRdd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+svdRdd x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
@@ -77,14 +66,9 @@ foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
 
 -- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
 --
--- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@.
-svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double))
-svdC x = (u, diagRect s r c, v) where (u,s,v) = svdC' x
-                                      r = rows x
-                                      c = cols x
-
-svdC' :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-svdC' x = unsafePerformIO $ do
+-- @(u,s,v)=full svdC m@ so that @m=u \<\> comp s \<\> 'trans' v@.
+svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
+svdC x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 3112ce6..7953386 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -14,15 +14,13 @@ Portability :  uses ffi
 -----------------------------------------------------------------------------
 
 module LinearAlgebra.Algorithms (
-    -- mXv, vXm,
-    inv,
-    pinv,
+    GMatrix(..),
+    Normed(..), NormType(..),
+    det,inv,pinv,full,economy,
     pinvTol,
-    pinvTolg,
+--    pinvTolg,
     nullspacePrec,
     nullVector,
-    Normed(..), NormType(..),
-    det,
     eps, i
 ) where
 
@@ -33,6 +31,69 @@ import GSL.Matrix
 import GSL.Vector
 import LAPACK
 import Complex
+import LinearAlgebra.Linear
+
+class (Linear Matrix t) => GMatrix t where
+    svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
+    lu          :: Matrix t -> (Matrix t, Matrix t, [Int], t)
+    linearSolve :: Matrix t -> Matrix t -> Matrix t
+    linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
+    ctrans :: Matrix t -> Matrix t
+    eig         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
+    eigSH       :: Matrix t -> (Vector Double, Matrix t)
+
+instance GMatrix Double where
+    svd = svdR
+    lu  = luR
+    linearSolve = linearSolveR
+    linearSolveSVD = linearSolveSVDR Nothing
+    ctrans = trans
+    eig = eigR
+    eigSH = eigS
+
+instance GMatrix (Complex Double) where
+    svd = svdC
+    lu  = luC
+    linearSolve = linearSolveC
+    linearSolveSVD = linearSolveSVDC Nothing
+    ctrans = conjTrans
+    eig = eigC
+    eigSH = eigH
+
+square m = rows m == cols m
+
+det :: GMatrix t => Matrix t -> t
+det m | square m = s * (product $ toList $ takeDiag $ u)
+      | otherwise = error "det of nonsquare matrix"
+    where (_,u,_,s) = lu m
+
+inv :: GMatrix t => Matrix t -> Matrix t
+inv m | square m = m `linearSolve` ident (rows m)
+      | otherwise = error "inv of nonsquare matrix"
+
+pinv :: GMatrix t => Matrix t -> Matrix t
+pinv m = linearSolveSVD m (ident (rows m))
+
+
+full svd m = (u, d ,v) where
+    (u,s,v) = svd m
+    d = diagRect s r c
+    r = rows m
+    c = cols m
+
+economy svd m = (u', subVector 0 d s, v') where
+    (u,s,v) = svd m
+    sl@(g:_) = toList (complex s)
+    s' = fromList . filter rec $ sl
+    rec x = magnitude x > magnitude g*tol
+    t = 1
+    tol = (fromIntegral (max (rows m) (cols m)) * magnitude g * t * eps)
+    r = rows m
+    c = cols m
+    d = dim s'
+    u' = takeColumns d u
+    v' = takeColumns d v
+
 
 {- | Machine precision of a Double.
 
@@ -70,119 +131,6 @@ vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t
 vXm v m = flatten $ (asRow v) `mXm` m
 
 
-
--- | Pseudoinverse of a real matrix
---
--- @dispR 3 $ pinv (fromLists [[1,2],
---                           [3,4],
---                           [5,6]])
---matrix (2x3)
--- -1.333 | -0.333 |  0.667
---  1.083 |  0.333 | -0.417@
---
-
-pinv :: Matrix Double -> Matrix Double
-pinv m = pinvTol 1 m
---pinv m = linearSolveSVDR Nothing m (ident (rows m))
-
-{- -| Pseudoinverse of a real matrix with the default tolerance used by GNU-Octave: the singular values less than max (rows, colums) * greatest singular value * 'eps' are ignored. See 'pinvTol'.
-
-@\> let m = 'fromLists' [[ 1, 2]
-                    ,[ 5, 8]
-                    ,[10,-5]]
-\> pinv m
-9.353e-3 4.539e-2  7.637e-2
-2.231e-2 8.993e-2 -4.719e-2
-\ 
-\> m \<\> pinv m \<\> m
- 1.  2.
- 5.  8.
-10. -5.@
-
--}
---pinvg :: Matrix Double -> Matrix Double
-pinvg m = pinvTolg 1 m
-
-{- | Pseudoinverse of a real matrix with the desired tolerance, expressed as a
-multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
-
-@\> let m = 'fromLists' [[1,0,    0]
-                    ,[0,1,    0]
-                    ,[0,0,1e-10]]
-\ 
-\> 'pinv' m 
-1. 0.           0.
-0. 1.           0.
-0. 0. 10000000000.
-\ 
-\> pinvTol 1E8 m
-1. 0. 0.
-0. 1. 0.
-0. 0. 1.@
-
--}
-pinvTol :: Double -> Matrix Double -> Matrix Double
-pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
-    (u,s,v) = svdR' m
-    sl@(g:_) = toList s
-    s' = fromList . map rec $ sl
-    rec x = if x < g*tol then 1 else 1/x
-    tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps)
-    r = rows m
-    c = cols m
-    d = dim s
-    u' = takeColumns d u
-    v' = takeColumns d v
-
-
-pinvTolg :: Double -> Matrix Double -> Matrix Double
-pinvTolg t m = v `mXm` diag s' `mXm` trans u where
-    (u,s,v) = svdg m
-    sl@(g:_) = toList s
-    s' = fromList . map rec $ sl
-    rec x = if x < g*tol then 1 else 1/x
-    tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps)
-
-
-
-{- | Inverse of a square matrix.
-
-inv m = 'linearSolveR' m ('ident' ('rows' m))
-
-@\>inv ('fromLists' [[1,4]
-                ,[0,2]])
-1.   -2.
-0. 0.500@
--}
-inv :: Matrix Double -> Matrix Double
-inv m = if rows m == cols m
-    then m `linearSolveR` ident (rows m)
-    else error "inv of nonsquare matrix"
-
-
-{- - | Shortcut for the 2-norm ('pnorm' 2)
-
-@ > norm $ 'hilb' 5
-1.5670506910982311
-@
-
-@\> norm $ 'fromList' [1,-1,'i',-'i']
-2.0@
-
--}
-
-
-
-{- | Determinant of a square matrix, computed from the LU decomposition.
-
-@\> det ('fromLists' [[7,2],[3,8]])
-50.0@
-
--}
-det :: Matrix Double -> Double
-det m = s * (product $ toList $ takeDiag $ u)
-    where (_,u,_,s) = luR m
-
 ---------------------------------------------------------------------------
 
 norm2 :: Vector Double -> Double
@@ -212,12 +160,12 @@ pnormCV PNorm1 = norm1 . liftVector magnitude
 pnormCV Infinity = vectorMax . liftVector magnitude
 --pnormCV _ = error "pnormCV not yet defined"
 
-pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR' m
+pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR m
 pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m
 pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)
 --pnormRM _ _ = error "p norm not yet defined"
 
-pnormCM PNorm2 m = head (toList s) where (_,s,_) = svdC' m
+pnormCM PNorm2 m = head (toList s) where (_,s,_) = svdC m
 pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (liftVector magnitude) m
 pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` constant 1 (cols m)
 --pnormCM _ _ = error "p norm not yet defined"
@@ -245,17 +193,52 @@ instance Normed (Matrix (Complex Double)) where
 
 -----------------------------------------------------------------------
 
--- | The nullspace of a real matrix from its SVD decomposition.
-nullspacePrec :: Double          -- ^ relative tolerance in 'eps' units
-              -> Matrix Double   -- ^ input matrix
-              -> [Vector Double] -- ^ list of unitary vectors spanning the nullspace
+-- | The nullspace of a matrix from its SVD decomposition.
+nullspacePrec :: GMatrix t
+              => Double          -- ^ relative tolerance in 'eps' units
+              -> Matrix t   -- ^ input matrix
+              -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
 nullspacePrec t m = ns where
-    (_,s,v) = svdR' m
+    (_,s,v) = svd m
     sl@(g:_) = toList s
     tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps)
     rank = length (filter (> g*tol) sl)
-    ns = drop rank (toColumns v)
+--    ns = drop rank (toColumns v)
+    ns = drop rank $ toRows $ ctrans v
 
--- | The nullspace of a real matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
-nullVector :: Matrix Double -> Vector Double
+-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
+nullVector :: GMatrix t => Matrix t -> Vector t
 nullVector = last . nullspacePrec 1
+
+------------------------------------------------------------------------
+
+{- | Pseudoinverse of a real matrix with the desired tolerance, expressed as a
+multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
+
+@\> let m = 'fromLists' [[1,0,    0]
+                    ,[0,1,    0]
+                    ,[0,0,1e-10]]
+\ 
+\> 'pinv' m 
+1. 0.           0.
+0. 1.           0.
+0. 0. 10000000000.
+\ 
+\> pinvTol 1E8 m
+1. 0. 0.
+0. 1. 0.
+0. 0. 1.@
+
+-}
+pinvTol :: Double -> Matrix Double -> Matrix Double
+pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
+    (u,s,v) = svdR m
+    sl@(g:_) = toList s
+    s' = fromList . map rec $ sl
+    rec x = if x < g*tol then 1 else 1/x
+    tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps)
+    r = rows m
+    c = cols m
+    d = dim s
+    u' = takeColumns d u
+    v' = takeColumns d v
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index c12e30b..2f1bc6f 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -37,6 +37,8 @@ class (Field e) => Linear c e where
     fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
     comp        :: RealFloat e => c e -> c (Complex e)
     conj        :: RealFloat e => c (Complex e) -> c (Complex e)
+    real        :: c Double -> c e
+    complex     :: c e -> c (Complex Double)
 
 instance Linear Vector Double where
     scale = vectorMapValR Scale
@@ -50,6 +52,8 @@ instance Linear Vector Double where
     fromComplex = Data.Packed.Internal.fromComplex
     comp = Data.Packed.Internal.comp
     conj = Data.Packed.Internal.conj
+    real = id
+    complex = LinearAlgebra.Linear.comp
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
@@ -63,6 +67,8 @@ instance Linear Vector (Complex Double) where
     fromComplex = undefined
     comp = undefined
     conj = undefined
+    real = LinearAlgebra.Linear.comp
+    complex = id
 
 instance Linear Matrix Double where
     scale x = liftMatrix (scale x)
@@ -78,6 +84,8 @@ instance Linear Matrix Double where
               c = cols z
     comp = liftMatrix Data.Packed.Internal.comp
     conj = liftMatrix Data.Packed.Internal.conj
+    real = id
+    complex = LinearAlgebra.Linear.comp
 
 instance Linear Matrix (Complex Double) where
     scale x = liftMatrix (scale x)
@@ -91,6 +99,8 @@ instance Linear Matrix (Complex Double) where
     fromComplex = undefined
     comp = undefined
     conj = undefined
+    real = LinearAlgebra.Linear.comp
+    complex = id
 
 --------------------------------------------------