l.a. algorithms, etc.

12 files changed, 1119 insertions, 89 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index fc9db60..84ca5f5 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -24,7 +24,7 @@ Exposed-modules:    Data.Packed.Internal,
                     Data.Packed.Internal.Matrix, Data.Packed.Matrix,
                     Data.Packed.Internal.Tensor, Data.Packed.Tensor,
                     Data.Packed.Plot
-                    -- Data.Packed.Instances
+                    Data.Packed.Instances
                     LAPACK,
                     GSL.Vector, 
                     GSL.Matrix
diff --git a/examples/tests.hs b/examples/tests.hs
index 37800cd..9de9339 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -16,12 +16,26 @@ import GSL.Fourier
 import GSL.Polynomials
 import LAPACK
 import Test.QuickCheck
-import Test.HUnit
+import Test.HUnit hiding ((~:))
 import Complex
+import LinearAlgebra.Algorithms
+import GSL.Matrix
+import Data.Packed.Instances hiding ((<>))
+
+dist :: (Normed t, Num t) => t -> t -> Double
+dist a b = norm (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!    
+-- Bravo por quickCheck!
 
 pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m  ~~ ident (cols m)
     where infix 2 ~~
@@ -32,7 +46,7 @@ pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m  ~~ m)
           c = cols m
           infix 2 ~~
           (~~) = approxEqual tol 
-        
+
 nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros
     where nt    = trans (nullspace m)
           cr    = corank m
@@ -49,31 +63,28 @@ 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]
 
-a |=| b = rows a == rows b &&
-          cols a == cols b &&
-          toList (cdat a) == toList (cdat b) &&
-          toList (fdat a) == toList (fdat b)
 
+{-
 aprox fun a b = rows a == rows b &&
           cols a == cols b &&
-          eps > aproxL fun (toList (t a)) (toList (t 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 && eps > normVR (t a) (t 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 < eps
-
+u ~|~ v = normVR u v < epsTol
+-}
 
-eps = 1E-8::Double
+epsTol = 1E-8::Double
 
 asFortran m = (rows m >|< cols m) $ toList (fdat m)
 asC m = (rows m >< cols m) $ toList (cdat m)
@@ -81,6 +92,9 @@ asC m = (rows m >< cols m) $ toList (cdat m)
 mulC a b = multiply RowMajor a b
 mulF a b = multiply ColumnMajor a b
 
+infixl 7 <>
+a <> b = mulF a b
+
 cc = mulC ac bf
 cf = mulF af bc
 
@@ -133,14 +147,14 @@ data Sym a = Sym (Matrix a) deriving Show
 instance (Field a, Arbitrary a, Num a) => Arbitrary (Sym a) where
     arbitrary = do
         SqM m <- arbitrary
-        return $ Sym (m `addM` trans m)
+        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 `addM` (liftMatrix conj) (trans m))
+        return $ Her (m + conjTrans m)
     coarbitrary = undefined
 
 data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show
@@ -171,49 +185,40 @@ instance (Field a, Arbitrary a) => Arbitrary (PairV a) where
    coarbitrary = undefined
 
 
-addM m1 m2 = liftMatrix2 add m1 m2
-
 
 type BaseType = Double
 
-svdTestR fun prod m = u <> s <> trans v |~| m
-                  && u <> trans u |~| ident (rows m)
-                  && v <> trans v |~| ident (cols m)
+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
-          (<>) = prod
 
 
-svdTestC prod m = u <> s' <> (trans v) |~~| m
-                  && u <> (liftMatrix conj) (trans u) |~~| ident (rows m)
-                  && v <> (liftMatrix conj) (trans v) |~~| ident (cols m)
+svdTestC m = u <> s' <> (trans v) |~| m
+                  && u <> conjTrans u |~| ident (rows m)
+                  && v <> conjTrans v |~| ident (cols m)
     where (u,s,v) = svdC m
-          (<>) = prod
           s' = liftMatrix comp s
 
-eigTestC prod (SqM m) = (m <> v) |~~| (v <> diag s)
-                        && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized
+--svdg' m = (u,s',v) where 
+
+eigTestC (SqM m) = (m <> v) |~| (v <> diag s)
+                        && takeDiag (conjTrans v <> v) |~| constant 1 (rows m) --normalized
     where (s,v) = eigC m
-          (<>) = prod
 
-eigTestR prod (SqM m) = (liftMatrix comp m <> v) |~~| (v <> diag s)
-                      -- && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized ???
+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
-          (<>) = prod
 
-eigTestS prod (Sym m) = (m <> v) |~| (v <> diag s)
+eigTestS (Sym m) = (m <> v) |~| (v <> diag s)
                        && v <> trans v |~| ident (cols m)
     where (s,v) = eigS m
-          (<>) = prod
 
-eigTestH prod (Her m) = (m <> v) |~~| (v <> diag (comp s))
-                        && v <> (liftMatrix conj) (trans v) |~~| ident (cols m)
+eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s))
+                        && v <> conjTrans v |~| ident (cols m)
     where (s,v) = eigH m
-          (<>) = prod
-
-linearSolveSQTest fun eqfun singu prod (PairSM a b) = singu a || (a <> fun a b) ==== b
-    where (<>) = prod
-          (====) = eqfun
 
+linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b
 
 prec = 1E-15
 
@@ -237,8 +242,7 @@ identC n = toComplex(ident n, (0::Double) <>ident n)
 
 --------------------------------------------------------------------
 
-pinvTest f feq m = (m <> f m <> m) `feq` m
-    where (<>) = mulF
+pinvTest f m = (m <> f m <> m) |~| m
 
 pinvR m = linearSolveLSR m (ident (rows m))
 pinvC m = linearSolveLSC m (ident (rows m))
@@ -252,7 +256,7 @@ pinvSVDC m = linearSolveSVDC Nothing m (ident (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))
-    where l1 |~~| l2 = eps > aproxL magnitude l1 l2
+
 
 polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4])
 
@@ -267,17 +271,17 @@ quad2 f a b g1 g2 = quad h a b
 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) < eps)
+integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)
 
 
 ---------------------------------------------------------------------
 
 arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u)
           `add` vectorMapValR PowVS 2 (vectorMapR Cos u)
-          ~|~ constant 1 (dim u)
+          |~| constant 1 (dim u)
 
 arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u)
-          ~|~ vectorMapR Sin u
+          |~| vectorMapR Sin u
 
 
 -- arit3 (PairV u v) =
@@ -305,50 +309,48 @@ tests = TestList
 
 main = do
     putStrLn "--------- general -----"
-    quickCheck (\(Sym m) -> m |=| (trans m:: Matrix BaseType))    
+    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)
+    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)
+    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 mulC)
-    quickCheck (svdTestR svdR mulF)
-    quickCheck (svdTestR svdRdd mulC)
-    quickCheck (svdTestR svdRdd mulF)
-    quickCheck (svdTestC mulC)
-    quickCheck (svdTestC mulF)
+    quickCheck (svdTestR svdR)
+    quickCheck (svdTestR svdRdd)
+--    quickCheck (svdTestR svdg)
+    quickCheck svdTestC
     putStrLn "--------- EIG ---------"
-    quickCheck (eigTestC mulC)
-    quickCheck (eigTestC mulF)
-    quickCheck (eigTestR mulC)
-    quickCheck (eigTestR mulF)
-    quickCheck (eigTestS mulC)
-    quickCheck (eigTestS mulF)
-    quickCheck (eigTestH mulC)
-    quickCheck (eigTestH mulF)
+    quickCheck eigTestC
+    quickCheck eigTestR
+    quickCheck eigTestS
+    quickCheck eigTestH
     putStrLn "--------- SOLVE ---------"
-    quickCheck (linearSolveSQTest linearSolveR (|~|) (singular svdR') mulC)
-    quickCheck (linearSolveSQTest linearSolveC (|~~|) (singular svdC') mulF)
-    quickCheck (pinvTest pinvR (|~|))
-    quickCheck (pinvTest pinvC (|~~|))
-    quickCheck (pinvTest pinvSVDR (|~|))
-    quickCheck (pinvTest pinvSVDC (|~~|))
+    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
     putStrLn "--------- GSL ------"
     runTestTT tests
-    quickCheck $ \v -> ifft (fft v) ~~ v
+    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]
+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]
 
-pol =  [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]
\ No newline at end of file
+mm = (2><2)
+ [ 0.5, 0.0
+ , 0.0, 0.0 ] :: Matrix Double
diff --git a/lib/Data/Packed/Instances.hs b/lib/Data/Packed/Instances.hs
new file mode 100644
index 0000000..4478469
--- /dev/null
+++ b/lib/Data/Packed/Instances.hs
@@ -0,0 +1,391 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Data.Packed.Instances
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses -fffi and -fglasgow-exts
+
+Creates reasonable numeric instances for Vectors and Matrices. In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
+
+-}
+-----------------------------------------------------------------------------
+
+module Data.Packed.Instances(
+    Contractible(..)
+) where
+
+import Data.Packed.Internal
+import Data.Packed.Vector
+import Data.Packed.Matrix
+import GSL.Vector
+import GSL.Matrix
+import LinearAlgebra.Algorithms
+import Complex
+
+instance (Eq a, Field a) => Eq (Vector a) where
+    a == b = dim a == dim b && toList a == toList b
+
+instance (Num a, Field a) => Num (Vector a) where
+    (+) = add
+    (-) = sub
+    (*) = mul
+    signum = liftVector signum
+    abs = liftVector abs
+    fromInteger = fromList . return . fromInteger
+
+instance (Eq a, Field a) => Eq (Matrix a) where
+    a == b = rows a == rows b && cols a == cols b && cdat a == cdat b && fdat a == fdat b
+
+instance (Num a, Field a) => Num (Matrix a) where
+    (+) = liftMatrix2 add
+    (-) = liftMatrix2 sub
+    (*) = liftMatrix2 mul
+    signum = liftMatrix signum
+    abs = liftMatrix abs
+    fromInteger = (1><1) . return . fromInteger
+
+---------------------------------------------------
+
+adaptScalar f1 f2 f3 x y
+    | dim x == 1 = f1 (x@>0) y
+    | dim y == 1 = f3 x (y@>0)
+    | otherwise = f2 x y
+
+{-
+subvv = vectorZip 4
+subvc v c = addConstant (-c) v
+subcv c v = addConstant c (scale (-1) v)
+
+mul = vectorZip 1
+
+instance Num (Vector Double) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    (-) = adaptScalar subcv subvv subvc
+    (*) = adaptScalar scale mul (flip scale)
+    abs = vectorMap 3
+    signum = vectorMap 15
+    fromInteger n = fromList [fromInteger n]
+
+----------------------------------------------------
+
+--addConstantC a = gmap (+a)
+--subCvv u v = u `add` scale (-1) v
+subCvv = vectorZipComplex 4 -- faster?
+subCvc v c = addConstantC (-c) v
+subCcv c v = addConstantC c (scale (-1) v)
+
+
+instance Num (Vector (Complex Double)) where
+    (+) = adaptScalar addConstantC add (flip addConstantC)
+    (-) = adaptScalar subCcv subCvv subCvc
+    (*) = adaptScalar scale (vectorZipComplex 1) (flip scale)
+    abs = gmap abs
+    signum = gmap signum
+    fromInteger n = fromList [fromInteger n]
+
+
+-- | adapts a function on two vectors to work on all the elements of two matrices
+liftMatrix2' :: (Vector a -> Vector b -> Vector c) -> Matrix a -> Matrix b -> Matrix c
+liftMatrix2' f m1@(M r1 c1 _) m2@(M r2 c2 _)
+    | sameShape m1 m2 || r1*c1==1 || r2*c2==1
+        = reshape (max c1 c2) $ f (flatten m1) (flatten m2)
+    | otherwise = error "inconsistent matrix dimensions" 
+
+---------------------------------------------------
+
+instance (Eq a, Field a) => Eq (Matrix a) where
+    a == b = rows a == rows b && cdat a == cdat b
+
+instance Num (Matrix Double) where
+    (+) = liftMatrix2' (+)
+    (-) = liftMatrix2' (-)
+    (*) = liftMatrix2' (*)
+    abs = liftMatrix abs
+    signum = liftMatrix signum
+    fromInteger n = fromLists [[fromInteger n]]
+
+----------------------------------------------------
+
+instance Num (Matrix (Complex Double)) where
+    (+) = liftMatrix2' (+)
+    (-) = liftMatrix2' (-)
+    (*) = liftMatrix2' (*)
+    abs = liftMatrix abs
+    signum = liftMatrix signum
+    fromInteger n = fromLists [[fromInteger n]]
+
+------------------------------------------------------
+
+instance Fractional (Vector Double) where
+    fromRational n = fromList [fromRational n]
+    (/) = adaptScalar f (vectorZip 2) g where
+        r `f` v = vectorZip 2 (constant r (dim v)) v
+        v `g` r = scale (recip r) v
+
+-------------------------------------------------------
+
+instance Fractional (Vector (Complex Double)) where
+    fromRational n = fromList [fromRational n]
+    (/) = adaptScalar f (vectorZipComplex 2) g where
+        r `f` v = gmap ((*r).recip) v
+        v `g` r = gmap (/r) v
+
+------------------------------------------------------
+
+instance Fractional (Matrix Double) where
+    fromRational n = fromLists [[fromRational n]]
+    (/) = liftMatrix2' (/)
+
+-------------------------------------------------------
+
+instance Fractional (Matrix (Complex Double)) where
+    fromRational n = fromLists [[fromRational n]]
+    (/) = liftMatrix2' (/)
+
+---------------------------------------------------------
+
+instance Floating (Vector Double) where
+    sin   = vectorMap 0
+    cos   = vectorMap 1
+    tan   = vectorMap 2
+    asin  = vectorMap 4
+    acos  = vectorMap 5
+    atan  = vectorMap 6
+    sinh  = vectorMap 7
+    cosh  = vectorMap 8
+    tanh  = vectorMap 9
+    asinh = vectorMap 10
+    acosh = vectorMap 11
+    atanh = vectorMap 12
+    exp   = vectorMap 13
+    log   = vectorMap 14
+    sqrt  = vectorMap 16
+    (**)  = adaptScalar f (vectorZip 5) g where f s v = constant s (dim v) ** v
+                                                g v s = v ** constant s (dim v)
+    pi    = fromList [pi]
+
+-----------------------------------------------------------
+
+instance Floating (Matrix Double) where
+    sin   = liftMatrix sin
+    cos   = liftMatrix cos
+    tan   = liftMatrix tan
+    asin  = liftMatrix asin
+    acos  = liftMatrix acos
+    atan  = liftMatrix atan
+    sinh  = liftMatrix sinh
+    cosh  = liftMatrix cosh
+    tanh  = liftMatrix tanh
+    asinh = liftMatrix asinh
+    acosh = liftMatrix acosh
+    atanh = liftMatrix atanh
+    exp   = liftMatrix exp
+    log   = liftMatrix log
+    sqrt  = liftMatrix sqrt
+    (**)  = liftMatrix2 (**)
+    pi    = fromLists [[pi]]
+
+-------------------------------------------------------------
+
+instance Floating (Vector (Complex Double)) where
+    sin   = vectorMapComplex 0
+    cos   = vectorMapComplex 1
+    tan   = vectorMapComplex 2
+    asin  = vectorMapComplex 4
+    acos  = vectorMapComplex 5
+    atan  = vectorMapComplex 6
+    sinh  = vectorMapComplex 7
+    cosh  = vectorMapComplex 8
+    tanh  = vectorMapComplex 9
+    asinh = vectorMapComplex 10
+    acosh = vectorMapComplex 11
+    atanh = vectorMapComplex 12
+    exp   = vectorMapComplex 13
+    log   = vectorMapComplex 14
+    sqrt  = vectorMapComplex 16
+    (**)  = adaptScalar f (vectorZipComplex 5) g where f s v = constantC s (dim v) ** v
+                                                       g v s = v ** constantC s (dim v)
+    pi    = fromList [pi]
+
+---------------------------------------------------------------
+
+instance Floating (Matrix (Complex Double)) where
+    sin   = liftMatrix sin
+    cos   = liftMatrix cos
+    tan   = liftMatrix tan
+    asin  = liftMatrix asin
+    acos  = liftMatrix acos
+    atan  = liftMatrix atan
+    sinh  = liftMatrix sinh
+    cosh  = liftMatrix cosh
+    tanh  = liftMatrix tanh
+    asinh = liftMatrix asinh
+    acosh = liftMatrix acosh
+    atanh = liftMatrix atanh
+    exp   = liftMatrix exp
+    log   = liftMatrix log
+    (**)  = liftMatrix2 (**)
+    sqrt  = liftMatrix sqrt
+    pi    = fromLists [[pi]]
+
+---------------------------------------------------------------
+-}
+
+class Contractible a b c | a b -> c where
+ infixl 7 <>
+{- | An overloaded operator for matrix products, matrix-vector and vector-matrix products, dot products and scaling of vectors and matrices. Type consistency is statically checked. Alternatively, you can use the specific functions described below, but using this operator you can automatically combine real and complex objects.
+
+@v  = 'fromList' [1,2,3]    :: Vector Double
+cv = 'fromList' [1+'i',2]
+m  = 'fromLists' [[1,2,3],
+                [4,5,7]] :: Matrix Double
+cm = 'fromLists' [[  1,  2],
+                [3+'i',7*'i'],
+                [  'i',  1]]
+\ 
+\> m \<\> v
+14. 35.
+\ 
+\> cv \<\> m
+9.+1.i  12.+2.i  17.+3.i
+\ 
+\> m \<\> cm
+  7.+5.i   5.+14.i
+19.+12.i  15.+35.i
+\ 
+\> v \<\> 'i'
+1.i  2.i  3.i
+\ 
+\> v \<\> v
+14.0
+\ 
+\> cv \<\> cv
+4.0 :+ 2.0@
+
+-}
+ (<>) :: a -> b -> c
+
+
+instance Contractible Double Double Double where
+ (<>) = (*)
+
+instance Contractible Double (Complex Double) (Complex Double) where
+ a <> b = (a:+0) * b
+
+instance Contractible (Complex Double) Double (Complex Double) where
+ a <> b = a * (b:+0)
+
+instance Contractible (Complex Double) (Complex Double) (Complex Double) where
+ (<>) = (*)
+
+--------------------------------- matrix matrix
+
+instance Contractible (Matrix Double) (Matrix Double) (Matrix Double) where
+ (<>) = mXm
+
+instance Contractible (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where
+ (<>) = mXm
+
+instance Contractible (Matrix (Complex Double)) (Matrix Double) (Matrix (Complex Double)) where
+ c <> r = c <> liftMatrix comp r
+
+instance Contractible (Matrix Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
+ r <> c = liftMatrix comp r <> c
+
+--------------------------------- (Matrix Double) (Vector Double)
+
+instance Contractible (Matrix Double) (Vector Double) (Vector Double) where
+ (<>) = mXv
+
+instance Contractible (Matrix (Complex Double)) (Vector (Complex Double)) (Vector (Complex Double)) where
+ (<>) = mXv
+
+instance Contractible (Matrix (Complex Double)) (Vector Double) (Vector (Complex Double)) where
+ m <> v = m <> comp v
+
+instance Contractible (Matrix Double) (Vector (Complex Double)) (Vector (Complex Double)) where
+ m <> v = liftMatrix comp m <> v
+
+--------------------------------- (Vector Double) (Matrix Double)
+
+instance Contractible (Vector Double) (Matrix Double) (Vector Double) where
+ (<>) = vXm
+
+instance Contractible (Vector (Complex Double)) (Matrix (Complex Double)) (Vector (Complex Double)) where
+ (<>) = vXm
+
+instance Contractible (Vector (Complex Double)) (Matrix Double) (Vector (Complex Double)) where
+ v <> m = v <> liftMatrix comp m
+
+instance Contractible (Vector Double) (Matrix (Complex Double)) (Vector (Complex Double)) where
+ v <> m = comp v <> m
+
+--------------------------------- dot product
+
+instance Contractible (Vector Double) (Vector Double) Double where
+ (<>) = dot
+
+instance Contractible (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where
+ (<>) = dot
+
+instance Contractible (Vector Double) (Vector (Complex Double)) (Complex Double) where
+ a <> b = comp a <> b
+
+instance Contractible (Vector (Complex Double)) (Vector Double) (Complex Double) where
+ (<>) = flip (<>)
+
+--------------------------------- scaling vectors  
+
+instance Contractible Double (Vector Double) (Vector Double) where
+ (<>) = scale
+
+instance Contractible (Vector Double) Double (Vector Double) where
+ (<>) = flip (<>)
+
+instance Contractible (Complex Double) (Vector (Complex Double)) (Vector (Complex Double)) where
+ (<>) = scale
+
+instance Contractible (Vector (Complex Double)) (Complex Double) (Vector (Complex Double)) where
+ (<>) = flip (<>)
+
+instance Contractible Double (Vector (Complex Double)) (Vector (Complex Double)) where
+ a <> v = (a:+0) <> v
+
+instance Contractible (Vector (Complex Double)) Double (Vector (Complex Double)) where
+ (<>) = flip (<>)
+
+instance Contractible (Complex Double) (Vector Double) (Vector (Complex Double)) where
+ a <> v = a <> comp v
+
+instance Contractible (Vector Double) (Complex Double) (Vector (Complex Double)) where
+ (<>) = flip (<>)
+
+--------------------------------- scaling matrices
+
+instance Contractible Double (Matrix Double) (Matrix Double) where
+ (<>) a = liftMatrix (a <>)
+
+instance Contractible (Matrix Double) Double (Matrix Double) where
+ (<>) = flip (<>)
+
+instance Contractible (Complex Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
+ (<>) a = liftMatrix (a <>)
+
+instance Contractible (Matrix (Complex Double)) (Complex Double) (Matrix (Complex Double)) where
+ (<>) = flip (<>)
+
+instance Contractible Double (Matrix (Complex Double)) (Matrix (Complex Double)) where
+ a <> m = (a:+0) <> m
+
+instance Contractible (Matrix (Complex Double)) Double (Matrix (Complex Double)) where
+ (<>) = flip (<>)
+
+instance Contractible (Complex Double) (Matrix Double) (Matrix (Complex Double)) where
+ a <> m = a <> liftMatrix comp m
+
+instance Contractible (Matrix Double) (Complex Double) (Matrix (Complex Double)) where
+ (<>) = flip (<>)
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 9309d1d..dd33943 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -93,6 +93,15 @@ createMatrix order r c = do
     p <- createVector (r*c)
     return (matrixFromVector order c p)
 
+{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns.
+
+@\> reshape 4 ('fromList' [1..12])
+(3><4)
+ [ 1.0,  2.0,  3.0,  4.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 9.0, 10.0, 11.0, 12.0 ]@
+
+-}
 reshape :: (Field t) => Int -> Vector t -> Matrix t
 reshape c v = matrixFromVector RowMajor c v
 
@@ -140,7 +149,6 @@ liftMatrix f m = m { dat = f (dat m), tdat = f (tdat m) } -- check sizes
 
 liftMatrix2 :: (Field t) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
 liftMatrix2 f m1 m2 = reshape (cols m1) (f (cdat m1) (cdat m2)) -- check sizes
-
 ------------------------------------------------------------------
 
 dotL a b = sum (zipWith (*) a b)
@@ -200,6 +208,14 @@ multiplyD order a b
 
 outer' u v = dat (outer u v)
 
+{- | Outer product of two vectors.
+
+@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
+(3><3)
+ [  5.0, 2.0, 3.0
+ , 10.0, 4.0, 6.0
+ , 15.0, 6.0, 9.0 ]@
+-}
 outer :: (Num t, Field t) => Vector t -> Vector t -> Matrix t
 outer u v = multiply RowMajor r c
     where r = matrixFromVector RowMajor 1 u
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 25e848d..f1addf4 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -48,7 +48,7 @@ fromList l = unsafePerformIO $ do
 toList :: Storable a => Vector a -> [a]
 toList v = unsafePerformIO $ peekArray (dim v) (ptr v)
 
-n # l = if length l == n then fromList l else error "# with wrong size"
+n |> l = if length l == n then fromList l else error "|> with wrong size"
 
 at' :: Storable a => Vector a -> Int -> a
 at' v n = unsafePerformIO $ peekElemOff (ptr v) n
@@ -58,7 +58,7 @@ at v n | n >= 0 && n < dim v = at' v n
        | otherwise          = error "vector index out of range"
 
 instance (Show a, Storable a) => (Show (Vector a)) where
-    show v = (show (dim v))++" # " ++ show (toList v)
+    show v = (show (dim v))++" |> " ++ show (toList v)
 
 -- | creates a Vector taking a number of consecutive toList from another Vector
 subVector :: Storable t => Int       -- ^ index of the starting element
@@ -129,3 +129,5 @@ constant x n | isReal id x = scast $ constantR (scast x) n
              | isComp id x = scast $ constantC (scast x) n
              | otherwise   = constantG x n
 
+liftVector  f = fromList . map f . toList
+liftVector2 f u v = fromList $ zipWith f (toList u) (toList v)
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 0f9d998..36bf32e 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -13,11 +13,11 @@
 -----------------------------------------------------------------------------
 
 module Data.Packed.Matrix (
-    Matrix(rows,cols), Field,
+    Matrix(rows,cols),
     fromLists, toLists, (><), (>|<), (@@>),
     trans, conjTrans,
-    reshape, flatten,
-    fromRows, toRows, fromColumns, toColumns,
+    reshape, flatten, asRow, asColumn,
+    fromRows, toRows, fromColumns, toColumns, fromBlocks,
     joinVert, joinHoriz,
     flipud, fliprl,
     liftMatrix, liftMatrix2,
@@ -43,6 +43,22 @@ joinVert ms = case common cols ms of
 joinHoriz :: Field t => [Matrix t] -> Matrix t
 joinHoriz ms = trans. joinVert . map trans $ ms
 
+{- | Creates a matrix from blocks given as a list of lists of matrices:
+
+@\> let a = 'diag' $ 'fromList' [5,7,2]
+\> let b = 'reshape' 4 $ 'constant' (-1) 12
+\> fromBlocks [[a,b],[b,a]]
+(6><7)
+ [  5.0,  0.0,  0.0, -1.0, -1.0, -1.0, -1.0
+ ,  0.0,  7.0,  0.0, -1.0, -1.0, -1.0, -1.0
+ ,  0.0,  0.0,  2.0, -1.0, -1.0, -1.0, -1.0
+ , -1.0, -1.0, -1.0, -1.0,  5.0,  0.0,  0.0
+ , -1.0, -1.0, -1.0, -1.0,  0.0,  7.0,  0.0
+ , -1.0, -1.0, -1.0, -1.0,  0.0,  0.0,  2.0 ]@
+-}
+fromBlocks :: Field t => [[Matrix t]] -> Matrix t
+fromBlocks = joinVert . map joinHoriz 
+
 -- | Reverse rows 
 flipud :: Field t => Matrix t -> Matrix t
 flipud m = fromRows . reverse . toRows $ m
@@ -98,6 +114,11 @@ dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat
 
 ----------------------------------------------------------------
 
+{- | Creates a vector by concatenation of rows
+
+@\> flatten ('ident' 3)
+9 # [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
+-}
 flatten :: Matrix t -> Vector t
 flatten = cdat
 
@@ -106,4 +127,10 @@ fromLists :: Field t => [[t]] -> Matrix t
 fromLists = fromRows . map fromList
 
 conjTrans :: Matrix (Complex Double) -> Matrix (Complex Double)
-conjTrans = trans . liftMatrix conj
\ No newline at end of file
+conjTrans = trans . liftMatrix conj
+
+asRow :: Field a => Vector a -> Matrix a
+asRow v = reshape (dim v) v
+
+asColumn :: Field a => Vector a -> Matrix a
+asColumn v = reshape 1 v
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 9d9d879..94f70be 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -15,7 +15,7 @@
 module Data.Packed.Vector (
     Vector(dim), Field,
     fromList, toList,
-    at,
+    (@>),
     subVector, join,
     constant,
     toComplex, comp,
@@ -26,6 +26,7 @@ module Data.Packed.Vector (
 
 import Data.Packed.Internal
 import Complex
+import GSL.Vector
 
 -- | creates a complex vector from vectors with real and imaginary parts
 toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
@@ -41,10 +42,14 @@ comp v = toComplex (v,constant 0 (dim v))
 
 {- | Creates a real vector containing a range of values:
 
-> > linspace 10 (-2,2)
->-2. -1.556 -1.111 -0.667 -0.222 0.222 0.667 1.111 1.556 2.
-
+@\> linspace 5 (-3,7)
+5 |> [-3.0,-0.5,2.0,4.5,7.0]@
 -}
 linspace :: Int -> (Double, Double) -> Vector Double
 linspace n (a,b) = fromList [a::Double,a+delta .. b]
     where delta = (b-a)/(fromIntegral n -1)
+
+-- | Reads a vector position.
+(@>) :: Field t => Vector t -> Int -> t
+infixl 9 @>
+(@>) = at
diff --git a/lib/GSL.hs b/lib/GSL.hs
new file mode 100644
index 0000000..8e033c3
--- /dev/null
+++ b/lib/GSL.hs
@@ -0,0 +1,45 @@
+{- |
+
+Module      :  GSL
+Copyright   :  (c) Alberto Ruiz 2006-7
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses -fffi and -fglasgow-exts
+
+This module reexports the basic functionality and a collection of utilities (old interface)
+
+-}
+
+module GSL (
+
+module Data.Packed.Vector,
+module Data.Packed.Matrix,
+module Data.Packed.Tensor,
+module Data.Packed.Instances,
+module LinearAlgebra.Algorithms,
+module LAPACK,
+module GSL.Integration,
+module GSL.Differentiation,
+module GSL.Special,
+module GSL.Fourier,
+module GSL.Polynomials,
+module GSL.Minimization,
+module Data.Packed.Plot
+
+) where
+
+import Data.Packed.Vector
+import Data.Packed.Matrix
+import Data.Packed.Tensor
+import Data.Packed.Instances
+import LinearAlgebra.Algorithms
+import LAPACK
+import GSL.Integration
+import GSL.Differentiation
+import GSL.Special
+import GSL.Fourier
+import GSL.Polynomials
+import GSL.Minimization
+import Data.Packed.Plot
diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs
new file mode 100644
index 0000000..919c2d9
--- /dev/null
+++ b/lib/GSL/Matrix.hs
@@ -0,0 +1,311 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  GSL.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- A few linear algebra computations based on the GSL (<http://www.gnu.org/software/gsl>).
+--
+-----------------------------------------------------------------------------
+
+module GSL.Matrix(
+    eigSg, eigHg,
+    svdg,
+    qr,
+    chol,
+    luSolveR, luSolveC,
+    luR, luC,
+    fromFile
+) where
+
+import Data.Packed.Internal
+import Data.Packed.Matrix(fromLists,ident,takeDiag)
+import GSL.Vector
+import Foreign
+import Foreign.C.Types
+import Complex
+import Foreign.C.String
+
+{- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/.
+
+> > let (l,v) = eigS $ 'fromLists' [[1,2],[2,1]]
+> > l
+> 3.000 -1.000
+>
+> > v
+> 0.707 -0.707
+> 0.707  0.707
+>
+> > v <> diag l <> trans v
+> 1.000 2.000
+> 2.000 1.000
+
+-}
+eigSg :: Matrix Double -> (Vector Double, Matrix Double)
+eigSg (m@M {rows = r})
+    | r == 1 = (fromList [cdat m `at` 0], singleton 1)
+    | otherwise = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix RowMajor r r
+        c_eigS // mat cdat m // vec l // mat dat v // check "eigSg" [cdat m]
+        return (l,v)
+foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
+
+------------------------------------------------------------------
+
+
+
+{- | eigendecomposition of a complex hermitian matrix using /gsl_eigen_hermv/
+
+> > let (l,v) = eigH $ 'fromLists' [[1,2+i],[2-i,3]]
+>
+> > l
+> 4.449 -0.449
+>
+> > v
+>         -0.544          0.839
+> (-0.751,0.375) (-0.487,0.243)
+>
+> > v <> diag l <> (conjTrans) v
+>          1.000 (2.000,1.000)
+> (2.000,-1.000)         3.000
+
+-}
+eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double))
+eigHg (m@M {rows = r})
+    | r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1)
+    | otherwise = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix RowMajor r r
+        c_eigH // mat cdat m // vec l // mat dat v // check "eigHg" [cdat m]
+        return (l,v)
+foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
+
+
+{- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/:
+
+@\> let (u,s,v) = svdg $ 'fromLists' [[1,2,3],[-4,1,7]]
+\
+\> u
+0.310 -0.951
+0.951  0.310
+\
+\> s
+8.497 2.792
+\
+\> v
+-0.411 -0.785
+ 0.185 -0.570
+ 0.893 -0.243
+\
+\> u \<\> 'diag' s \<\> 'trans' v
+ 1. 2. 3.
+-4. 1. 7.@
+
+-}
+svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+svdg x@M {rows = r, cols = c} = if r>=c
+    then svd' x
+    else (v, s, u) where (u,s,v) = svd' (trans x)
+
+svd' x@M {rows = r, cols = c} = unsafePerformIO $ do
+    u <- createMatrix RowMajor r c
+    s <- createVector c
+    v <- createMatrix RowMajor c c
+    c_svd // mat cdat x // mat dat u // vec s // mat dat v // check "svdg" [cdat x]
+    return (u,s,v)
+foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
+
+{- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/.
+
+@\> let (q,r) = qr $ 'fromLists' [[1,3,5,7],[2,0,-2,4]]
+\
+\> q
+-0.447 -0.894
+-0.894  0.447
+\
+\> r
+-2.236 -1.342 -0.447 -6.708
+    0. -2.683 -5.367 -4.472
+\
+\> q \<\> r
+1.000 3.000  5.000 7.000
+2.000    0. -2.000 4.000@
+
+-}
+qr :: Matrix Double -> (Matrix Double, Matrix Double)
+qr x@M {rows = r, cols = c} = unsafePerformIO $ do
+    q <- createMatrix RowMajor r r
+    rot <- createMatrix RowMajor r c
+    c_qr // mat cdat x // mat dat q // mat dat rot // check "qr" [cdat x]
+    return (q,rot)
+foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
+
+{- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
+
+@\> let c = chol $ 'fromLists' [[5,4],[4,5]]
+\
+\> c
+2.236    0.
+1.789 1.342
+\
+\> c \<\> 'trans' c
+5.000 4.000
+4.000 5.000@
+
+-}
+chol :: Matrix Double -> Matrix Double
+--chol x@(M r _ p) = createM [p] "chol" r r $ m c_chol x
+chol x@M {rows = r} = unsafePerformIO $ do
+    res <- createMatrix RowMajor r r
+    c_chol // mat cdat x // mat dat res // check "chol" [cdat x]
+    return res
+foreign import ccall "gsl-aux.h chol" c_chol :: TMM
+
+--------------------------------------------------------
+
+{- -| efficient multiplication by the inverse of a matrix (for real matrices)
+-}
+luSolveR :: Matrix Double -> Matrix Double -> Matrix Double
+luSolveR  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
+    | n1==n2 && n1==r = unsafePerformIO $ do
+        s <- createMatrix RowMajor r c
+        c_luSolveR // mat cdat a // mat cdat b // mat dat s // check "luSolveR" [cdat a, cdat b]
+        return s
+    | otherwise = error "luSolveR of nonsquare matrix"
+
+foreign import ccall "gsl-aux.h luSolveR" c_luSolveR ::  TMMM
+
+{- -| efficient multiplication by the inverse of a matrix (for complex matrices). 
+-}
+luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+luSolveC  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
+    | n1==n2 && n1==r = unsafePerformIO $ do
+        s <- createMatrix RowMajor r c
+        c_luSolveC // mat cdat a // mat cdat b // mat dat s // check "luSolveC" [cdat a, cdat b]
+        return s
+    | otherwise = error "luSolveC of nonsquare matrix"
+
+foreign import ccall "gsl-aux.h luSolveC" c_luSolveC ::  TCMCMCM
+
+{- | lu decomposition of real matrix (packed as a vector including l, u, the permutation and sign)
+-}
+luRaux  :: Matrix Double -> Vector Double
+luRaux x@M {rows = r, cols = c} = unsafePerformIO $ do
+    res <- createVector (r*r+r+1)
+    c_luRaux // mat cdat x // vec res // check "luRaux" [cdat x]
+    return res
+foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV
+
+{- | lu decomposition of complex matrix (packed as a vector including l, u, the permutation and sign)
+-}
+luCaux  :: Matrix (Complex Double) -> Vector (Complex Double)
+luCaux x@M {rows = r, cols = c} = unsafePerformIO $ do
+    res <- createVector (r*r+r+1)
+    c_luCaux // mat cdat x // vec res // check "luCaux" [cdat x]
+    return res
+foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV
+
+{- | The LU decomposition of a square matrix. Is based on /gsl_linalg_LU_decomp/ and  /gsl_linalg_complex_LU_decomp/ as described in <http://www.gnu.org/software/gsl/manual/gsl-ref_13.html#SEC223>.
+
+@\> let m = 'fromLists' [[1,2,-3],[2+3*i,-7,0],[1,-i,2*i]]
+\> let (l,u,p,s) = luR m@
+
+L is the lower triangular:
+
+@\> l
+          1.            0.  0.
+0.154-0.231i            1.  0.
+0.154-0.231i  0.624-0.522i  1.@
+
+U is the upper triangular:
+
+@\> u
+2.+3.i           -7.            0.
+    0.  3.077-1.615i           -3.
+    0.            0.  1.873+0.433i@
+
+p is a permutation:
+
+@\> p
+[1,0,2]@
+
+L \* U obtains a permuted version of the original matrix:
+
+@\> 'extractRows' p m
+  2.+3.i   -7.   0.
+      1.    2.  -3.
+      1.  -1.i  2.i
+\ 
+\> l \<\> u
+ 2.+3.i   -7.   0.
+     1.    2.  -3.
+     1.  -1.i  2.i@
+
+s is the sign of the permutation, required to obtain sign of the determinant:
+
+@\> s * product ('toList' $ 'takeDiag' u)
+(-18.0) :+ (-16.000000000000004)
+\> 'LinearAlgebra.Algorithms.det' m
+(-18.0) :+ (-16.000000000000004)@
+
+ -}
+luR :: Matrix Double -> (Matrix Double, Matrix Double, [Int], Double)
+luR m = (l,u,p, fromIntegral s') where
+    r = rows m
+    v = luRaux m
+    lu = reshape r $ subVector 0 (r*r) v
+    s':p = map round . toList . subVector (r*r) (r+1) $ v
+    u = triang r r 0 1`mul` lu
+    l = (triang r r 0 0 `mul` lu) `add` ident r
+    add = liftMatrix2 $ vectorZipR Add
+    mul = liftMatrix2 $ vectorZipR Mul
+
+-- | Complex version of 'luR'.
+luC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double), [Int], Complex Double)
+luC m = (l,u,p, fromIntegral s') where
+    r = rows m
+    v = luCaux m
+    lu = reshape r $ subVector 0 (r*r) v
+    s':p = map (round.realPart) . toList . subVector (r*r) (r+1) $ v
+    u = triang r r 0 1 `mul` lu
+    l = (triang r r 0 0 `mul` lu) `add` ident r
+    add = liftMatrix2 $ vectorZipC Add
+    mul = liftMatrix2 $ vectorZipC Mul
+
+extract l is = [l!!i |i<-is]
+
+{- auxiliary function to get triangular matrices
+-}
+triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]
+    where el i j = if j-i>=h then v else 1 - v
+
+{- | rearranges the rows of a matrix according to the order given in a list of integers. 
+
+> > extractRows [3,3,0,1] (ident 4)
+> 0. 0. 0. 1.
+> 0. 0. 0. 1.
+> 1. 0. 0. 0.
+> 0. 1. 0. 0.
+
+-}
+extractRows :: Field t => [Int] -> Matrix t -> Matrix t
+extractRows l m = fromRows $ extract (toRows $ m) l
+
+--------------------------------------------------------------
+
+-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = do
+    charname <- newCString filename
+    res <- createMatrix RowMajor r c
+    c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
+    --free charname  -- TO DO: free the auxiliary CString
+    return res
+foreign import ccall "gsl-aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
+
+---------------------------------------------------------------------------
diff --git a/lib/GSL/Vector.hs b/lib/GSL/Vector.hs
index a772b34..a074254 100644
--- a/lib/GSL/Vector.hs
+++ b/lib/GSL/Vector.hs
@@ -18,7 +18,7 @@ module GSL.Vector (
     FunCodeV(..), vectorMapR, vectorMapC,
     FunCodeSV(..), vectorMapValR, vectorMapValC,
     FunCodeVV(..), vectorZipR, vectorZipC,
-    scale, addConstant, add, mul,
+    scale, addConstant, add, sub, mul,
 ) where
 
 import Data.Packed.Internal
@@ -84,6 +84,11 @@ add u v | isReal baseOf v = scast $ vectorZipR Add (scast u) (scast v)
                 | isComp baseOf v = scast $ vectorZipC Add (scast u) (scast v)
                 | otherwise   = fromList $ zipWith (+) (toList u) (toList v)
 
+sub :: (Num a, Field a) => Vector a -> Vector a -> Vector a
+sub u v | isReal baseOf v = scast $ vectorZipR Sub (scast u) (scast v)
+                | isComp baseOf v = scast $ vectorZipC Sub (scast u) (scast v)
+                | otherwise   = fromList $ zipWith (-) (toList u) (toList v)
+
 mul :: (Num a, Field a) => Vector a -> Vector a -> Vector a
 mul u v | isReal baseOf v = scast $ vectorZipR Mul (scast u) (scast v)
                 | isComp baseOf v = scast $ vectorZipC Mul (scast u) (scast v)
diff --git a/lib/LinearAlgebra.hs b/lib/LinearAlgebra.hs
index 474f9a6..b0c8b9d 100644
--- a/lib/LinearAlgebra.hs
+++ b/lib/LinearAlgebra.hs
@@ -15,3 +15,4 @@ Some linear algebra algorithms, implemented by means of BLAS, LAPACK or GSL.
 module LinearAlgebra (
 
 ) where
+
diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 680612f..126549a 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -1,3 +1,4 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  LinearAlgebra.Algorithms
@@ -13,5 +14,229 @@ Portability :  uses ffi
 -----------------------------------------------------------------------------
 
 module LinearAlgebra.Algorithms (
-
+    mXm, mXv, vXm,
+    inv,
+    pinv,
+    pinvTol,
+    pinvTolg,
+    Normed(..), NormType(..),
+    det,
+    eps, i
 ) where
+
+
+import Data.Packed.Internal
+import Data.Packed.Matrix
+import GSL.Matrix
+import GSL.Vector
+import LAPACK
+import Complex
+
+{- | Machine precision of a Double.
+
+>> eps
+> 2.22044604925031e-16
+
+(The value used by GNU-Octave)
+
+-}
+eps :: Double
+eps =  2.22044604925031e-16
+
+{- | The imaginary unit
+
+@> 'ident' 3 \<\> i
+1.i   0.   0.
+ 0.  1.i   0.
+ 0.   0.  1.i@
+
+-}
+i :: Complex Double
+i = 0:+1
+
+
+-- | matrix product
+mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t
+mXm = multiply RowMajor
+
+-- | matrix - vector product
+mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t
+mXv m v = flatten $ m `mXm` (asColumn v)
+
+-- | vector - matrix product
+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
+norm2 = toScalarR Norm2
+
+norm1 :: Vector Double -> Double
+norm1 = toScalarR AbsSum
+
+vectorMax :: Vector Double -> Double
+vectorMax = toScalarR Max
+vectorMin :: Vector Double -> Double
+vectorMin = toScalarR Min
+vectorMaxIndex :: Vector Double -> Int
+vectorMaxIndex = round . toScalarR MaxIdx
+vectorMinIndex :: Vector Double -> Int
+vectorMinIndex = round . toScalarR MinIdx
+
+data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int
+
+pnormRV PNorm2 = norm2
+pnormRV PNorm1 = norm1
+pnormRV Infinity = vectorMax . vectorMapR Abs
+--pnormRV _ = error "pnormRV not yet defined"
+
+pnormCV PNorm2 = norm2 . asReal
+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 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 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"
+
+-- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'.
+--pnorm :: (Container t, Field a) => Int -> t a -> Double
+--pnorm = pnormG
+
+class Normed t where
+    pnorm :: NormType -> t -> Double
+    norm :: t -> Double
+    norm = pnorm PNorm2
+
+instance Normed (Vector Double) where
+    pnorm = pnormRV
+
+instance Normed (Vector (Complex Double)) where
+    pnorm = pnormCV
+
+instance Normed (Matrix Double) where
+    pnorm = pnormRM
+
+instance Normed (Matrix (Complex Double)) where
+    pnorm = pnormCM