refactoring

20 files changed, 387 insertions, 124 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index 2bd30bc..0dda5d7 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -18,17 +18,23 @@ Extensions:         ForeignFunctionInterface
 --ghc-options:        -Wall
 ghc-options:        -O
 hs-source-dirs:     lib
-Exposed-modules:    Data.Packed.Internal, Data.Packed.Internal.Common,
+Exposed-modules:    Data.Packed.Internal,
+                    Data.Packed.Internal.Common,
                     Data.Packed.Internal.Vector, Data.Packed.Vector,
                     Data.Packed.Internal.Matrix, Data.Packed.Matrix,
-                    Data.Packed.Internal.Tensor,
+                    Data.Packed.Internal.Tensor, Data.Packed.Tensor,
+                    Data.Packed.Plot
+                    -- Data.Packed.Instances
                     LAPACK,
-                    GSL.Vector,
+                    GSL.Vector, 
+                    --GSL.Matrix
                     GSL.Differentiation, GSL.Integration,
                     GSL.Special,
                     GSL.Fourier,
                     GSL.Polynomials,
-                    GSL.Minimization
+                    GSL.Minimization,
+                    LinearAlgebra, LinearAlgebra.Algorithms,
+                    GSL, HSSL
 Other-modules:
 C-sources:          lib/Data/Packed/Internal/aux.c,
                     lib/LAPACK/lapack-aux.c,
diff --git a/examples/pru.hs b/examples/pru.hs
index fb33962..f855bce 100644
--- a/examples/pru.hs
+++ b/examples/pru.hs
@@ -52,12 +52,14 @@ main = do
     print $ foldl part t [("p",1),("r",2),("q",0)]
     print $ foldl part t $ reverse [("p",1),("r",2),("q",0)]
 
-t = T [(4,(Covariant,"p")),(2,(Covariant,"q")),(3,(Contravariant,"r"))] $ fromList [1..24::Double]
+t = T [IdxDesc 4 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 3 Contravariant "r"]
+    $ fromList [1..24::Double]
 
 
-
-t1 = T [(4,(Covariant,"p")),(4,(Contravariant,"q")),(2,(Covariant,"r"))] $ fromList [1..32::Double]
-t2 = T [(4,(Covariant,"p")),(4,(Contravariant,"q"))] $ fromList [1..16::Double]
+t1 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q" ,IdxDesc 2 Covariant "r"]
+    $ fromList [1..32::Double]
+t2 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q"] 
+    $ fromList [1..16::Double]
 
 
 
@@ -72,15 +74,18 @@ e i n = fromList [ delta k i | k <- [1..n]]
 diagl = diag.fromList
 
 scalar x = T [] (fromList [x])
-tensorFromVector idx v = T {dims = [(dim v,idx)], ten = v}
-tensorFromMatrix idxr idxc m = T {dims = [(rows m,idxr),(cols m,idxc)], ten = cdat m}
+tensorFromVector (tp,nm) v = T {dims = [IdxDesc (dim v) tp nm]
+                                       , ten = v}
+tensorFromMatrix (tpr,nmr) (tpc,nmc) m = T {dims = [IdxDesc (rows m) tpr nmr,IdxDesc (cols m) tpc nmc]
+                                           , ten = cdat m}
 
 td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diagl [1..4] :: Tensor Double
 
 tn = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
 tt = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
 
-tq = T [(3,(Covariant,"p")),(2,(Covariant,"q")),(2,(Covariant,"r"))] $ fromList [11 .. 22] :: Tensor Double
+tq = T  [IdxDesc 3 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 2 Covariant "r"]
+    $ fromList [11 .. 22] :: Tensor Double
 
 r1 = contraction tt "j" tq "p"
 r1' = contraction' tt "j" tq "p"
@@ -101,7 +106,7 @@ pru = do
     print $ foldl part t2 [("j'",0),("p",1),("r",1)]
 
 scsig t = scalar (signature (nms t)) `prod` t
-    where nms = map (snd.snd) . dims
+    where nms = map idxName . dims
 
 antisym' t = addT $ map (scsig . flip tridx t) (perms (names t))
 
@@ -115,9 +120,9 @@ antisym' t = addT $ map (scsig . flip tridx t) (perms (names t))
 -}
 
 auxrename (T d v) = T d' v
-    where d' = [(n,(c,show (pos q))) | (n,(c,q)) <- d]
+    where d' = [IdxDesc n c (show (pos q)) | IdxDesc n c q <- d]
           pos n = i where Just i = elemIndex n nms
-          nms = map (snd.snd) d
+          nms = map idxName d
 
 antisym t = T (dims t) (ten (antisym' (auxrename t)))
 
@@ -156,16 +161,16 @@ l1 = vector "p" [0,0,0,1]
 l2 = vector "q" [1,0,0,1]
 l3 = vector "r" [0,1,0,1]
 
-leviCivita n = antisym $ foldl1 prod $ zipWith tcovector (map show [1..]) (toRows (ident n))
+leviCivita n = antisym $ foldl1 prod $ zipWith tcovector (map show [1,2..]) (toRows (ident n))
 
 contractionF t1 t2 = contraction t1 n1 t2 n2
     where n1 = fn t1
           n2 = fn t2
-          fn = snd . snd . head . dims
+          fn = idxName . head . dims
 
 
 dualV vs = foldl' contractionF (leviCivita n) vs
-    where n = fst . head . dims . head $ vs
+    where n = idxDim . head . dims . head $ vs
 
 
 dual1 = foldl' contractionF (leviCivita 3) [u,v]
@@ -184,16 +189,16 @@ x3 = vector "r" [-3,-1,-1]
 x4 = vector "s" [12,0,3]
 
 raise (T d v) = T (map raise' d) v
-    where raise' (n,(Covariant,s)) = (n,(Contravariant,s))
-          raise' (n,(Contravariant,s)) = (n,(Covariant,s))
+    where raise' idx@IdxDesc {idxType = Covariant } = idx {idxType = Contravariant}
+          raise' idx@IdxDesc {idxType = Contravariant } = idx {idxType = Covariant}
 
 dualMV t = prod (foldl' contract1b (lc <*> t) ds) (scalar (recip $ fromIntegral $ fact (length ds)))
     where
         lc = leviCivita n
-        nms1 = map (snd.snd) (dims lc)
-        nms2 = map ((++"'").snd.snd) (dims t)
+        nms1 = map idxName (dims lc)
+        nms2 = map ((++"'").idxName) (dims t)
         ds = zip nms1 nms2
-        n = fst . head . dims $ t
+        n = idxDim . head . dims $ t
 
 -- intersection of two lines :-)
 -- > raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
@@ -211,7 +216,7 @@ y4 = vector "s" [12,0,0,3]
 
 asBase r n = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]
 
-partF t i = part t (name,i) where name = snd . snd . head . dims $ t
+partF t i = part t (name,i) where name = idxName . head . dims $ t
 
 --partL = foldl' partF
 
@@ -219,7 +224,7 @@ niceAS t = filter ((/=0.0).fst) $ zip vals base
     where vals = map ((`at` 0).ten.foldl' partF t) (map (map pred) base)
           base = asBase r n
           r = length (dims t)
-          n = fst . head . dims $ t
+          n = idxDim . head . dims $ t
 
 z1 = vector "p" [0,0,0,1]
 z2 = vector "q" [1,0,0,1]
diff --git a/examples/tests.hs b/examples/tests.hs
index 4adc104..37800cd 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -191,12 +191,12 @@ svdTestC prod m = u <> s' <> (trans v) |~~| m
           s' = liftMatrix comp s
 
 eigTestC prod (SqM m) = (m <> v) |~~| (v <> diag s)
-                        && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant (rows m) 1 --normalized
+                        && takeDiag ((liftMatrix conj (trans 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 (rows m) 1 --normalized ???
+                      -- && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant 1 (rows m) --normalized ???
     where (s,v) = eigR m
           (<>) = prod
 
@@ -274,7 +274,7 @@ integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < eps
 
 arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u)
           `add` vectorMapValR PowVS 2 (vectorMapR Cos u)
-          ~|~ constant (dim u) 1
+          ~|~ constant 1 (dim u)
 
 arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u)
           ~|~ vectorMapR Sin u
diff --git a/lib/Data/Packed/Internal.hs b/lib/Data/Packed/Internal.hs
index a7fca1a..a5a77c5 100644
--- a/lib/Data/Packed/Internal.hs
+++ b/lib/Data/Packed/Internal.hs
@@ -15,9 +15,11 @@
 module Data.Packed.Internal (
     module Data.Packed.Internal.Common,
     module Data.Packed.Internal.Vector,
-    module Data.Packed.Internal.Matrix
+    module Data.Packed.Internal.Matrix,
+    module Data.Packed.Internal.Tensor
 ) where
 
 import Data.Packed.Internal.Common
 import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
+import Data.Packed.Internal.Tensor
\ No newline at end of file
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index bdd7f34..1bfed6d 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -40,6 +40,7 @@ instance (Storable a, RealFloat a) => Storable (Complex a) where    --
     poke p (a :+ b) = pokeArray (castPtr p) [a,b]                   --
 ----------------------------------------------------------------------
 
+on :: (a -> a -> b) -> (t -> a) -> t -> t -> b
 on f g = \x y -> f (g x) (g y)
 
 partit :: Int -> [a] -> [[a]]
@@ -54,12 +55,14 @@ common f = commonval . map f where
     commonval [a] = Just a
     commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
 
+xor :: Bool -> Bool -> Bool
 xor a b = a && not b || b && not a
 
 (//) :: x -> (x -> y) -> y
 infixl 0 //
 (//) = flip ($)
 
+errorCode :: Int -> String
 errorCode 1000 = "bad size"
 errorCode 1001 = "bad function code"
 errorCode 1002 = "memory problem"
@@ -68,6 +71,7 @@ errorCode 1004 = "singular"
 errorCode 1005 = "didn't converge"
 errorCode n    = "code "++show n
 
+check :: String -> [Vector a] -> IO Int -> IO ()
 check msg ls f = do
     err <- f
     when (err/=0) (error (msg++": "++errorCode err))
@@ -77,7 +81,10 @@ check msg ls f = do
 class (Storable a, Typeable a) => Field a
 instance (Storable a, Typeable a) => Field a
 
+isReal :: (Data.Typeable.Typeable a) => (t -> a) -> t -> Bool
 isReal w x   = typeOf (undefined :: Double) == typeOf (w x)
+
+isComp :: (Data.Typeable.Typeable a) => (t -> a) -> t -> Bool
 isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
 
 scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 2925fc0..32dc603 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -194,7 +194,9 @@ multiplyD order a b
 
 ----------------------------------------------------------------------
 
-outer u v = dat (multiply RowMajor r c)
+outer' u v = dat (outer u v)
+
+outer u v = multiply RowMajor r c
     where r = matrixFromVector RowMajor 1 u
           c = matrixFromVector RowMajor (dim v) v
 
@@ -212,8 +214,7 @@ subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do
     r <- createMatrix RowMajor rt ct
     c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat cdat r // check "subMatrixR" [dat r]
     return r
-foreign import ccall "aux.h submatrixR"
-    c_submatrixR :: Int -> Int -> Int -> Int -> TMM
+foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
 
 -- | extraction of a submatrix of a complex matrix
 subMatrixC :: (Int,Int) -- ^ (r0,c0) starting position
diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs
index 27fce6a..c4faf49 100644
--- a/lib/Data/Packed/Internal/Tensor.hs
+++ b/lib/Data/Packed/Internal/Tensor.hs
@@ -14,18 +14,25 @@
 
 module Data.Packed.Internal.Tensor where
 
-import Data.Packed.Internal
+import Data.Packed.Internal.Common
 import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
 import Foreign.Storable
 import Data.List(sort,elemIndex,nub)
 
-data IdxTp = Covariant | Contravariant deriving (Show,Eq)
+data IdxType = Covariant | Contravariant deriving (Show,Eq)
 
-data Tensor t = T { dims   :: [(Int,(IdxTp,String))]
+type IdxName = String
+
+data IdxDesc = IdxDesc { idxDim  :: Int,
+                         idxType :: IdxType,
+                         idxName :: IdxName }
+
+data Tensor t = T { dims   :: [IdxDesc]
                   , ten    :: Vector t
                   }
 
+rank :: Tensor t -> Int
 rank = length . dims
 
 instance (Show a,Storable a) => Show (Tensor a) where
@@ -33,41 +40,49 @@ instance (Show a,Storable a) => Show (Tensor a) where
     show T {dims = ds, ten = t} = "("++shdims ds ++") "++ show (toList t)
 
 
-shdims [(n,(t,name))] = name ++ sym t ++"["++show n++"]"
+shdims :: [IdxDesc] -> String
+shdims [IdxDesc n t name] = name ++ sym t ++"["++show n++"]"
     where sym Covariant     = "_"
           sym Contravariant = "^"
 shdims (d:ds) = shdims [d] ++ "><"++ shdims ds
 
 
-
+findIdx :: (Field t) => IdxName -> Tensor t
+        -> (([IdxDesc], [IdxDesc]), Matrix t)
 findIdx name t = ((d1,d2),m) where
-    (d1,d2) = span (\(_,(_,n)) -> n /=name) (dims t)
-    c = product (map fst d2)
+    (d1,d2) = span (\d -> idxName d /= name) (dims t)
+    c = product (map idxDim d2)
     m = matrixFromVector RowMajor c (ten t)
 
+putFirstIdx :: (Field t) => String -> Tensor t -> ([IdxDesc], Matrix t)
 putFirstIdx name t = (nd,m')
     where ((d1,d2),m) = findIdx name t
           m' = matrixFromVector RowMajor c $ cdat $ trans m
           nd = d2++d1
-          c = dim (ten t) `div` (fst $ head d2)
+          c = dim (ten t) `div` (idxDim $ head d2)
 
+part :: (Field t) => Tensor t -> (IdxName, Int) -> Tensor t
 part t (name,k) = if k<0 || k>=l
-                    then error $ "part "++show (name,k)++" out of range in "++show t
+                    then error $ "part "++show (name,k)++" out of range" -- in "++show t
                     else T {dims = ds, ten = toRows m !! k}
     where (d:ds,m) = putFirstIdx name t
-          (l,_) = d
+          l = idxDim d
 
+parts :: (Field t) => Tensor t -> IdxName -> [Tensor t]
 parts t name = map f (toRows m)
     where (d:ds,m) = putFirstIdx name t
-          (l,_) = d
+          l = idxDim d
           f t = T {dims=ds, ten=t}
 
+concatRename :: [IdxDesc] -> [IdxDesc] -> [IdxDesc]
 concatRename l1 l2 = l1 ++ map ren l2 where
-    ren (n,(t,s)) = if {- s `elem` fs -} True then (n,(t,s++"'")) else (n,(t,s))
-    fs = map (snd.snd) l1
+    ren idx = if {- s `elem` fs -} True then idx {idxName = idxName idx ++ "'"} else idx
+    fs = map idxName l1
 
-prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer v1 v2)
+prod :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t
+prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer' v1 v2)
 
+contraction :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
 contraction t1 n1 t2 n2 =
     if compatIdx t1 n1 t2 n2
         then T (concatRename (tail d1) (tail d2)) (cdat m)
@@ -76,18 +91,22 @@ contraction t1 n1 t2 n2 =
         (d2,m2) = putFirstIdx n2 t2
         m = multiply RowMajor (trans m1) m2
 
+sumT :: (Storable t, Enum t, Num t) => [Tensor t] -> [t]
 sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls)
 
+contract1 :: (Num t, Enum t, Field t) => Tensor t -> IdxName -> IdxName -> Tensor t
 contract1 t name1 name2 = T d $ fromList $ sumT y
     where d = dims (head y)
           x = (map (flip parts name2) (parts t name1))
           y = map head $ zipWith drop [0..] x
 
+contraction' :: (Field t, Enum t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
 contraction' t1 n1 t2 n2 =
     if compatIdx t1 n1 t2 n2
         then contract1 (prod t1 t2) n1 (n2++"'")
         else error "wrong contraction'"
 
+tridx :: (Field t) => [IdxName] -> Tensor t -> Tensor t
 tridx [] t = t
 tridx (name:rest) t = T (d:ds) (join ts) where
     ((_,d:_),_) = findIdx name t
@@ -95,30 +114,38 @@ tridx (name:rest) t = T (d:ds) (join ts) where
     ts = map ten ps
     ds = dims (head ps)
 
-compatIdxAux (n1,(t1,_)) (n2, (t2,_)) = t1 /= t2 && n1 == n2
+compatIdxAux :: IdxDesc -> IdxDesc -> Bool
+compatIdxAux IdxDesc {idxDim = n1, idxType = t1}
+             IdxDesc {idxDim = n2, idxType = t2}
+    = t1 /= t2 && n1 == n2
 
+compatIdx :: (Field t1, Field t) => Tensor t1 -> IdxName -> Tensor t -> IdxName -> Bool
 compatIdx t1 n1 t2 n2 = compatIdxAux d1 d2 where
     d1 = head $ snd $ fst $ findIdx n1 t1
     d2 = head $ snd $ fst $ findIdx n2 t2
 
-names t = sort $ map (snd.snd) (dims t)
+names :: Tensor t -> [IdxName]
+names t = sort $ map idxName (dims t)
 
+normal :: (Field t) => Tensor t -> Tensor t
 normal t = tridx (names t) t
 
+contractions :: (Num t, Field t) => Tensor t -> Tensor t -> [Tensor t]
 contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
 
 -- sent to Haskell-Cafe by Sebastian Sylvan
+perms :: [t] -> [[t]]
 perms [x] = [[x]]
 perms xs = [y:ps | (y,ys) <- selections xs , ps <- perms ys]
 selections []     = []
 selections (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- selections xs]
 
-
+interchanges :: (Ord a) => [a] -> Int
 interchanges ls = sum (map (count ls) ls)
-    where count l p = n
+    where count l p = length $ filter (>p) $ take pel l
               where Just pel = elemIndex p l
-                    n = length $ filter (>p) $ take pel l
 
+signature :: (Num t, Ord a) => [a] -> t
 signature l | length (nub l) < length l =  0
             | even (interchanges l)     =  1
             | otherwise                 = -1
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 8848062..25e848d 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -103,15 +103,15 @@ asComplex :: Vector Double -> Vector (Complex Double)
 asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v), ptr = castPtr (ptr v) }
 
 
-constantG n x = fromList (replicate n x)
+constantG x n = fromList (replicate n x)
 
-constantR :: Int -> Double -> Vector Double
+constantR :: Double -> Int -> Vector Double
 constantR = constantAux cconstantR
 
-constantC :: Int -> Complex Double -> Vector (Complex Double)
+constantC :: Complex Double -> Int -> Vector (Complex Double)
 constantC = constantAux cconstantC
 
-constantAux fun n x = unsafePerformIO $ do
+constantAux fun x n = unsafePerformIO $ do
     v <- createVector n
     px <- newArray [x]
     fun px // vec v // check "constantAux" []
@@ -124,8 +124,8 @@ foreign import ccall safe "aux.h constantR"
 foreign import ccall safe "aux.h constantC"
     cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int
 
-constant :: Field a => Int -> a -> Vector a
-constant n x | isReal id x = scast $ constantR n (scast x)
-             | isComp id x = scast $ constantC n (scast x)
-             | otherwise   = constantG n x
+constant :: Field a => a -> Int -> Vector a
+constant x n | isReal id x = scast $ constantR (scast x) n
+             | isComp id x = scast $ constantC (scast x) n
+             | otherwise   = constantG x n
 
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index ec5744d..c7d5cfa 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -16,12 +16,13 @@ module Data.Packed.Matrix (
     Matrix(rows,cols), Field,
     toLists, (><), (>|<), (@@>),
     trans,
-    reshape,
+    reshape, flatten,
     fromRows, toRows, fromColumns, toColumns,
     joinVert, joinHoriz,
     flipud, fliprl,
     liftMatrix, liftMatrix2,
     multiply,
+    outer,
     subMatrix,
     takeRows, dropRows, takeColumns, dropColumns,
     diag, takeDiag, diagRect, ident
@@ -54,11 +55,11 @@ diagRect s r c
     | r == c    = diag s
     | r < c     = trans $ diagRect s c r
     | r > c     = joinVert  [diag s , zeros (r-c,c)]
-    where zeros (r,c) = reshape c $ constant (r*c) 0
+    where zeros (r,c) = reshape c $ constant 0 (r*c)
 
 takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
 
-ident n = diag (constant n 1)
+ident n = diag (constant 1 n)
 
 r >< c = f where
     f l | dim v == r*c = matrixFromVector RowMajor c v
@@ -88,3 +89,5 @@ dropColumns :: Field t => Int -> Matrix t -> Matrix t
 dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat
 
 ----------------------------------------------------------------
+
+flatten = cdat
diff --git a/lib/Data/Packed/Plot.hs b/lib/Data/Packed/Plot.hs
new file mode 100644
index 0000000..9eddc9f
--- /dev/null
+++ b/lib/Data/Packed/Plot.hs
@@ -0,0 +1,167 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Plot
+-- Copyright   :  (c) Alberto Ruiz 2005
+-- License     :  GPL-style
+-- 
+-- Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+-- Stability   :  provisional
+-- Portability :  uses gnuplot and ImageMagick
+--
+-- Very basic (and provisional) drawing tools.
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Plot(
+
+    gnuplotX, mplot,
+
+    plot, parametricPlot, 
+
+    splot, mesh, mesh', meshdom,
+
+    matrixToPGM, imshow,
+
+) where
+
+import Data.Packed.Vector
+import Data.Packed.Matrix
+import GSL.Vector(FunCodeS(Max,Min),toScalarR)
+import Data.List(intersperse)
+import System
+import Data.IORef
+import System.Exit
+import Foreign hiding (rotate)
+
+
+size = dim
+
+-- | Loads a real matrix from a formatted ASCII text file 
+--fromFile :: FilePath -> IO Matrix
+--fromFile filename = readFile filename >>= return . readMatrix read
+
+-- | Saves a real matrix to a formatted ascii text file
+toFile :: FilePath -> Matrix Double -> IO ()
+toFile filename matrix = writeFile filename (unlines . map unwords. map (map show) . toLists $ matrix)
+
+------------------------------------------------------------------------
+
+
+-- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions.
+meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double)
+meshdom r1 r2 = (outer r1 (constant 1 (size r2)), outer (constant 1 (size r1)) r2)
+
+
+gnuplotX command = do {system cmdstr; return()} where
+    cmdstr = "echo \""++command++"\" | gnuplot -persist"
+
+datafollows = "\\\"-\\\""
+
+prep = (++"e\n\n") . unlines . map (unwords . (map show))
+
+
+{- | Draws a 3D surface representation of a real matrix.
+
+> > mesh (hilb 20)
+
+In certain versions you can interactively rotate the graphic using the mouse.
+
+-}
+mesh :: Matrix Double -> IO ()
+mesh m = gnuplotX (command++dat) where
+    command = "splot "++datafollows++" matrix with lines\n"
+    dat = prep $ toLists $ m
+
+mesh' m = do
+    writeFile "splot-gnu-command" "splot \"splot-tmp.txt\" matrix with lines; pause -1"; 
+    toFile "splot-tmp.txt" m
+    putStr "Press [Return] to close the graphic and continue... "
+    system "gnuplot -persist splot-gnu-command"
+    system "rm splot-tmp.txt splot-gnu-command"
+    return ()
+
+{- | Draws the surface represented by the function f in the desired ranges and number of points, internally using 'mesh'.
+
+> > let f x y = cos (x + y) 
+> > splot f (0,pi) (0,2*pi) 50    
+
+-}
+splot :: (Matrix Double->Matrix Double->Matrix Double) -> (Double,Double) -> (Double,Double) -> Int -> IO () 
+splot f rx ry n = mesh' z where
+    (x,y) = meshdom (linspace n rx) (linspace n ry)
+    z = f x y
+
+{- | plots several vectors against the first one -}
+mplot :: [Vector Double] -> IO ()
+mplot m = gnuplotX (commands++dats) where
+    commands = if length m == 1 then command1 else commandmore
+    command1 = "plot "++datafollows++" with lines\n" ++ dat
+    commandmore = "plot " ++ plots ++ "\n"
+    plots = concat $ intersperse ", " (map cmd [2 .. length m])
+    cmd k = datafollows++" using 1:"++show k++" with lines"
+    dat = prep $ toLists $ fromColumns m
+    dats = concat (replicate (length m-1) dat)
+
+
+
+
+
+
+mplot' m = do
+    writeFile "plot-gnu-command" (commands++endcmd)
+    toFile "plot-tmp.txt" (fromColumns m)
+    putStr "Press [Return] to close the graphic and continue... "
+    system "gnuplot plot-gnu-command"
+    system "rm plot-tmp.txt plot-gnu-command"
+    return ()
+ where
+    commands = if length m == 1 then command1 else commandmore
+    command1 = "plot \"plot-tmp.txt\" with lines\n"
+    commandmore = "plot " ++ plots ++ "\n"
+    plots = concat $ intersperse ", " (map cmd [2 .. length m])
+    cmd k = "\"plot-tmp.txt\" using 1:"++show k++" with lines"
+    endcmd = "pause -1"
+
+-- apply several functions to one object
+mapf fs x = map ($ x) fs
+
+{- | Draws a list of functions over a desired range and with a desired number of points 
+
+> > plot [sin, cos, sin.(3*)] (0,2*pi) 1000
+
+-}
+plot :: [Vector Double->Vector Double] -> (Double,Double) -> Int -> IO ()
+plot fs rx n = mplot (x: mapf fs x)
+    where x = linspace n rx  
+
+{- | Draws a parametric curve. For instance, to draw a spiral we can do something like:
+
+> > parametricPlot (\t->(t * sin t, t * cos t)) (0,10*pi) 1000
+
+-}
+parametricPlot :: (Vector Double->(Vector Double,Vector Double)) -> (Double, Double) -> Int -> IO ()
+parametricPlot f rt n = mplot [fx, fy]
+    where t = linspace n rt
+          (fx,fy) = f t
+
+
+-- | writes a matrix to pgm image file
+matrixToPGM :: Matrix Double -> String
+matrixToPGM m = header ++ unlines (map unwords ll) where
+    c = cols m
+    r = rows m
+    header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"
+    maxgray = 255.0
+    maxval = toScalarR Max $ flatten $ m
+    minval = toScalarR Min $ flatten $ m
+    scale = if (maxval == minval) 
+        then 0.0
+        else maxgray / (maxval - minval)
+    f x = show ( round ( scale *(x - minval) ) :: Int )
+    ll = map (map f) (toLists m)
+
+-- | imshow shows a representation of a matrix as a gray level image using ImageMagick's display.
+imshow :: Matrix Double -> IO ()
+imshow m = do
+    system $ "echo \""++ matrixToPGM m ++"\"| display -antialias -resize 300 - &"
+    return ()
diff --git a/lib/Data/Packed/Tensor.hs b/lib/Data/Packed/Tensor.hs
index 8d1c8b6..75a9288 100644
--- a/lib/Data/Packed/Tensor.hs
+++ b/lib/Data/Packed/Tensor.hs
@@ -1 +1,20 @@
- 
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Tensor
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Tensors
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Tensor (
+
+) where
+
+import Data.Packed.Internal
+import Complex
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 992301a..aa1b489 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -20,7 +20,8 @@ module Data.Packed.Vector (
     constant,
     toComplex, comp,
     conj,
-    dot
+    dot,
+    linspace
 ) where
 
 import Data.Packed.Internal
@@ -35,6 +36,14 @@ conj :: Vector (Complex Double) -> Vector (Complex Double)
 conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1])
     where mulC = multiply RowMajor
 
-comp v = toComplex (v,constant (dim v) 0)
+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 :: Int -> (Double, Double) -> Vector Double
+linspace n (a,b) = fromList [a::Double,a+delta .. b]
+    where delta = (b-a)/(fromIntegral n -1)
diff --git a/lib/GSL/Fourier.hs b/lib/GSL/Fourier.hs
index bf6cd60..9788602 100644
--- a/lib/GSL/Fourier.hs
+++ b/lib/GSL/Fourier.hs
@@ -29,7 +29,7 @@ genfft code v = unsafePerformIO $ do
     c_fft code // vec v // vec r // check "fft" [v]
     return r
 
-foreign import ccall "gsl-aux.h fft" c_fft ::  Int -> TCVCV -- Complex Double :> Complex Double :>  IO Int
+foreign import ccall "gsl-aux.h fft" c_fft ::  Int -> TCVCV
 
 
 {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
diff --git a/lib/GSL/Minimization.hs b/lib/GSL/Minimization.hs
index 32b266f..a59977e 100644
--- a/lib/GSL/Minimization.hs
+++ b/lib/GSL/Minimization.hs
@@ -97,7 +97,7 @@ minimizeNMSimplex f xi sz tol maxit = unsafePerformIO $ do
 
 foreign import ccall "gsl-aux.h minimize"
     c_minimizeNMSimplex:: FunPtr (Int -> Ptr Double -> Double) -> Double -> Int
-                       -> TVVM -- Double :> Double :> Double ::> IO Int
+                       -> TVVM
 
 ----------------------------------------------------------------------------------
 
@@ -165,7 +165,7 @@ foreign import ccall "gsl-aux.h minimizeWithDeriv"
     c_minimizeConjugateGradient :: FunPtr (Int -> Ptr Double -> Double)
                                 -> FunPtr (Int -> Ptr Double -> Ptr Double -> IO ())
                                 -> Double -> Double -> Double -> Int
-                                -> TVM -- Double :> Double ::> IO Int
+                                -> TVM
 
 ---------------------------------------------------------------------
 iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double)
diff --git a/lib/GSL/Polynomials.hs b/lib/GSL/Polynomials.hs
index 365c74e..a87fa56 100644
--- a/lib/GSL/Polynomials.hs
+++ b/lib/GSL/Polynomials.hs
@@ -51,4 +51,4 @@ polySolve' v | dim v > 1 = unsafePerformIO $ do
     return r
              | otherwise = error "polySolve on a polynomial of degree zero"
 
-foreign import ccall "gsl-aux.h polySolve" c_polySolve:: TVCV -- Double :> Complex Double :> IO Int
+foreign import ccall "gsl-aux.h polySolve" c_polySolve:: TVCV
diff --git a/lib/GSL/Vector.hs b/lib/GSL/Vector.hs
index fc80c16..a772b34 100644
--- a/lib/GSL/Vector.hs
+++ b/lib/GSL/Vector.hs
@@ -121,8 +121,7 @@ vectorZipAux fun code u v = unsafePerformIO $ do
 toScalarR :: FunCodeS -> Vector Double -> Double
 toScalarR oper =  toScalarAux c_toScalarR (fromEnum oper)
 
-foreign import ccall safe "gsl-aux.h toScalarR"
-    c_toScalarR :: Int -> TVV -- Double :> Double :> IO Int
+foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: Int -> TVV
 
 ------------------------------------------------------------------
 
@@ -130,15 +129,13 @@ foreign import ccall safe "gsl-aux.h toScalarR"
 vectorMapR :: FunCodeV -> Vector Double -> Vector Double
 vectorMapR = vectorMapAux c_vectorMapR
 
-foreign import ccall safe "gsl-aux.h mapR"
-    c_vectorMapR :: Int -> TVV -- Double :> Double :> IO Int
+foreign import ccall safe "gsl-aux.h mapR" c_vectorMapR :: Int -> TVV
 
 -- | map of complex vectors with given function
 vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double)
 vectorMapC oper = vectorMapAux c_vectorMapC (fromEnum oper)
 
-foreign import ccall safe "gsl-aux.h mapC"
-    c_vectorMapC :: Int -> TCVCV -- Complex Double :> Complex Double :> IO Int
+foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: Int -> TCVCV
 
 -------------------------------------------------------------------
 
@@ -146,15 +143,13 @@ foreign import ccall safe "gsl-aux.h mapC"
 vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double
 vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromEnum oper)
 
-foreign import ccall safe "gsl-aux.h mapValR"
-    c_vectorMapValR :: Int -> Ptr Double -> TVV -- Double :> Double :> IO Int
+foreign import ccall safe "gsl-aux.h mapValR" c_vectorMapValR :: Int -> Ptr Double -> TVV
 
 -- | map of complex vectors with given function
 vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
 vectorMapValC = vectorMapValAux c_vectorMapValC
 
-foreign import ccall safe "gsl-aux.h mapValC"
-    c_vectorMapValC :: Int -> Ptr (Complex Double) -> TCVCV -- Complex Double :> Complex Double :> IO Int
+foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: Int -> Ptr (Complex Double) -> TCVCV
 
 -------------------------------------------------------------------
 
@@ -162,14 +157,10 @@ foreign import ccall safe "gsl-aux.h mapValC"
 vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double
 vectorZipR = vectorZipAux c_vectorZipR
 
-foreign import ccall safe "gsl-aux.h zipR"
-    c_vectorZipR :: Int -> TVVV -- Double :> Double :> Double :> IO Int
+foreign import ccall safe "gsl-aux.h zipR" c_vectorZipR :: Int -> TVVV
 
 -- | elementwise operation on complex vectors
 vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
 vectorZipC = vectorZipAux c_vectorZipC
 
-foreign import ccall safe "gsl-aux.h zipC"
-    c_vectorZipC :: Int -> TCVCVCV -- Complex Double :> Complex Double :> Complex Double :> IO Int
-
-
+foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: Int -> TCVCVCV
diff --git a/lib/HSSL.hs b/lib/HSSL.hs
new file mode 100644
index 0000000..4506386
--- /dev/null
+++ b/lib/HSSL.hs
@@ -0,0 +1,17 @@
+{- |
+
+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.
+
+-}
+
+module HSSL (
+
+) where
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index dc9eda1..602f5df 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -26,13 +26,12 @@ import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
 import Data.Packed.Vector
 import Data.Packed.Matrix
+import GSL.Vector(scale)
 import Complex
 import Foreign
 
 -----------------------------------------------------------------------------
--- dgesvd
-foreign import ccall "LAPACK/lapack-aux.h svd_l_R"
-    dgesvd :: TMMVM -- Double ::> Double ::> (Double :> Double ::> IO Int)
+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.
 --
@@ -48,9 +47,7 @@ svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
     return (u,s,trans v)
 
 -----------------------------------------------------------------------------
--- dgesdd
-foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd"
-    dgesdd :: TMMVM --Double ::> Double ::> (Double :> Double ::> IO Int)
+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.
 --
@@ -66,9 +63,7 @@ svdRdd' x@M {rows = r, cols = c} = unsafePerformIO $ do
     return (u,s,trans v)
 
 -----------------------------------------------------------------------------
--- zgesvd
-foreign import ccall "LAPACK/lapack-aux.h svd_l_C"
-    zgesvd :: TCMCMVCM -- (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
+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.
 --
@@ -85,9 +80,7 @@ svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
     return (u,s,trans v)
 
 -----------------------------------------------------------------------------
--- zgeev
-foreign import ccall "LAPACK/lapack-aux.h eig_l_C"
-    zgeev :: TCMCMCVCM -- (Complex Double) ::> (Complex Double) ::> ((Complex Double) :> (Complex Double) ::> IO Int)
+foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
 
 -- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
 --
@@ -106,9 +99,7 @@ eigC (m@M {rows = r})
         return (l,v)
 
 -----------------------------------------------------------------------------
--- dgeev
-foreign import ccall "LAPACK/lapack-aux.h eig_l_R"
-    dgeev :: TMMCVM -- Double ::> Double ::> ((Complex Double) :> Double ::> IO Int)
+foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
 
 -- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
 --
@@ -139,12 +130,10 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
     | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs
     | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs)
 
-scale r v = fromList [r] `outer` v
+-- scale r v = fromList [r] `outer` v
 
 -----------------------------------------------------------------------------
--- dsyev
-foreign import ccall "LAPACK/lapack-aux.h eig_l_S"
-    dsyev :: TMVM -- Double ::> (Double :> Double ::> IO Int)
+foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
 
 -- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
 --
@@ -166,9 +155,7 @@ eigS' (m@M {rows = r})
         return (l,v)
 
 -----------------------------------------------------------------------------
--- zheev
-foreign import ccall "LAPACK/lapack-aux.h eig_l_H"
-    zheev :: TCMVCM -- (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
+foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
 
 -- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix:
 --
@@ -190,9 +177,7 @@ eigH' (m@M {rows = r})
         return (l,v)
 
 -----------------------------------------------------------------------------
--- dgesv
-foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l"
-    dgesv :: TMMM -- Double ::> Double ::> Double ::> IO Int
+foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
 
 -- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition.
 linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
@@ -204,9 +189,7 @@ linearSolveR  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
     | otherwise = error "linearSolveR of nonsquare matrix"
 
 -----------------------------------------------------------------------------
--- zgesv
-foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l"
-    zgesv :: TCMCMCM -- (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
+foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
 
 -- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition.
 linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
@@ -218,9 +201,7 @@ linearSolveC  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
     | otherwise = error "linearSolveC of nonsquare matrix"
 
 -----------------------------------------------------------------------------------
--- dgels
-foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l"
-    dgels :: TMMM -- Double ::> Double ::> Double ::> IO Int
+foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM
 
 -- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'.
 linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
@@ -232,9 +213,7 @@ linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformI
     return r
 
 -----------------------------------------------------------------------------------
--- zgels
-foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l"
-    zgels :: TCMCMCM -- (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
+foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM
 
 -- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'.
 linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
@@ -246,9 +225,7 @@ linearSolveLSC_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformI
     return r
 
 -----------------------------------------------------------------------------------
--- dgelss
-foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l"
-    dgelss :: Double -> TMMM -- Double ::> Double ::> Double ::> IO Int
+foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM
 
 -- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
 linearSolveSVDR :: Maybe Double   -- ^ rcond
@@ -264,9 +241,7 @@ linearSolveSVDR_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafeP
     return r
 
 -----------------------------------------------------------------------------------
--- zgelss
-foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l"
-    zgelss :: Double -> TCMCMCM -- (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
+foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
 
 -- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
 linearSolveSVDC :: Maybe Double            -- ^ rcond
diff --git a/lib/LinearAlgebra.hs b/lib/LinearAlgebra.hs
new file mode 100644
index 0000000..474f9a6
--- /dev/null
+++ b/lib/LinearAlgebra.hs
@@ -0,0 +1,17 @@
+-----------------------------------------------------------------------------
+{- |
+Module      :  LinearAlgebra
+Copyright   :  (c) Alberto Ruiz 2006-7
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+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
new file mode 100644
index 0000000..680612f
--- /dev/null
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -0,0 +1,17 @@
+-----------------------------------------------------------------------------
+{- |
+Module      :  LinearAlgebra.Algorithms
+Copyright   :  (c) Alberto Ruiz 2006-7
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+
+-}
+-----------------------------------------------------------------------------
+
+module LinearAlgebra.Algorithms (
+
+) where