summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2007-06-22 17:33:17 +0000
committerAlberto Ruiz <aruiz@um.es>2007-06-22 17:33:17 +0000
commit978e6d038239af50d70bae2c303f4e45b1879b7a (patch)
tree571b2060f388d0693820f808b40089acb100a5d9 /lib
parent989bdf7e88c13500bd1986dcde36f6cc4f467efb (diff)
refactoring
Diffstat (limited to 'lib')
-rw-r--r--lib/Data/Packed/Internal.hs4
-rw-r--r--lib/Data/Packed/Internal/Common.hs7
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs7
-rw-r--r--lib/Data/Packed/Internal/Tensor.hs65
-rw-r--r--lib/Data/Packed/Internal/Vector.hs16
-rw-r--r--lib/Data/Packed/Matrix.hs9
-rw-r--r--lib/Data/Packed/Plot.hs167
-rw-r--r--lib/Data/Packed/Tensor.hs21
-rw-r--r--lib/Data/Packed/Vector.hs13
-rw-r--r--lib/GSL/Fourier.hs2
-rw-r--r--lib/GSL/Minimization.hs4
-rw-r--r--lib/GSL/Polynomials.hs2
-rw-r--r--lib/GSL/Vector.hs23
-rw-r--r--lib/HSSL.hs17
-rw-r--r--lib/LAPACK.hs55
-rw-r--r--lib/LinearAlgebra.hs17
-rw-r--r--lib/LinearAlgebra/Algorithms.hs17
17 files changed, 349 insertions, 97 deletions
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 @@
15module Data.Packed.Internal ( 15module Data.Packed.Internal (
16 module Data.Packed.Internal.Common, 16 module Data.Packed.Internal.Common,
17 module Data.Packed.Internal.Vector, 17 module Data.Packed.Internal.Vector,
18 module Data.Packed.Internal.Matrix 18 module Data.Packed.Internal.Matrix,
19 module Data.Packed.Internal.Tensor
19) where 20) where
20 21
21import Data.Packed.Internal.Common 22import Data.Packed.Internal.Common
22import Data.Packed.Internal.Vector 23import Data.Packed.Internal.Vector
23import Data.Packed.Internal.Matrix 24import Data.Packed.Internal.Matrix
25import 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 --
40 poke p (a :+ b) = pokeArray (castPtr p) [a,b] -- 40 poke p (a :+ b) = pokeArray (castPtr p) [a,b] --
41---------------------------------------------------------------------- 41----------------------------------------------------------------------
42 42
43on :: (a -> a -> b) -> (t -> a) -> t -> t -> b
43on f g = \x y -> f (g x) (g y) 44on f g = \x y -> f (g x) (g y)
44 45
45partit :: Int -> [a] -> [[a]] 46partit :: Int -> [a] -> [[a]]
@@ -54,12 +55,14 @@ common f = commonval . map f where
54 commonval [a] = Just a 55 commonval [a] = Just a
55 commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing 56 commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
56 57
58xor :: Bool -> Bool -> Bool
57xor a b = a && not b || b && not a 59xor a b = a && not b || b && not a
58 60
59(//) :: x -> (x -> y) -> y 61(//) :: x -> (x -> y) -> y
60infixl 0 // 62infixl 0 //
61(//) = flip ($) 63(//) = flip ($)
62 64
65errorCode :: Int -> String
63errorCode 1000 = "bad size" 66errorCode 1000 = "bad size"
64errorCode 1001 = "bad function code" 67errorCode 1001 = "bad function code"
65errorCode 1002 = "memory problem" 68errorCode 1002 = "memory problem"
@@ -68,6 +71,7 @@ errorCode 1004 = "singular"
68errorCode 1005 = "didn't converge" 71errorCode 1005 = "didn't converge"
69errorCode n = "code "++show n 72errorCode n = "code "++show n
70 73
74check :: String -> [Vector a] -> IO Int -> IO ()
71check msg ls f = do 75check msg ls f = do
72 err <- f 76 err <- f
73 when (err/=0) (error (msg++": "++errorCode err)) 77 when (err/=0) (error (msg++": "++errorCode err))
@@ -77,7 +81,10 @@ check msg ls f = do
77class (Storable a, Typeable a) => Field a 81class (Storable a, Typeable a) => Field a
78instance (Storable a, Typeable a) => Field a 82instance (Storable a, Typeable a) => Field a
79 83
84isReal :: (Data.Typeable.Typeable a) => (t -> a) -> t -> Bool
80isReal w x = typeOf (undefined :: Double) == typeOf (w x) 85isReal w x = typeOf (undefined :: Double) == typeOf (w x)
86
87isComp :: (Data.Typeable.Typeable a) => (t -> a) -> t -> Bool
81isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x) 88isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
82 89
83scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b 90scast :: 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
194 194
195---------------------------------------------------------------------- 195----------------------------------------------------------------------
196 196
197outer u v = dat (multiply RowMajor r c) 197outer' u v = dat (outer u v)
198
199outer u v = multiply RowMajor r c
198 where r = matrixFromVector RowMajor 1 u 200 where r = matrixFromVector RowMajor 1 u
199 c = matrixFromVector RowMajor (dim v) v 201 c = matrixFromVector RowMajor (dim v) v
200 202
@@ -212,8 +214,7 @@ subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do
212 r <- createMatrix RowMajor rt ct 214 r <- createMatrix RowMajor rt ct
213 c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat cdat r // check "subMatrixR" [dat r] 215 c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat cdat r // check "subMatrixR" [dat r]
214 return r 216 return r
215foreign import ccall "aux.h submatrixR" 217foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
216 c_submatrixR :: Int -> Int -> Int -> Int -> TMM
217 218
218-- | extraction of a submatrix of a complex matrix 219-- | extraction of a submatrix of a complex matrix
219subMatrixC :: (Int,Int) -- ^ (r0,c0) starting position 220subMatrixC :: (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 @@
14 14
15module Data.Packed.Internal.Tensor where 15module Data.Packed.Internal.Tensor where
16 16
17import Data.Packed.Internal 17import Data.Packed.Internal.Common
18import Data.Packed.Internal.Vector 18import Data.Packed.Internal.Vector
19import Data.Packed.Internal.Matrix 19import Data.Packed.Internal.Matrix
20import Foreign.Storable 20import Foreign.Storable
21import Data.List(sort,elemIndex,nub) 21import Data.List(sort,elemIndex,nub)
22 22
23data IdxTp = Covariant | Contravariant deriving (Show,Eq) 23data IdxType = Covariant | Contravariant deriving (Show,Eq)
24 24
25data Tensor t = T { dims :: [(Int,(IdxTp,String))] 25type IdxName = String
26
27data IdxDesc = IdxDesc { idxDim :: Int,
28 idxType :: IdxType,
29 idxName :: IdxName }
30
31data Tensor t = T { dims :: [IdxDesc]
26 , ten :: Vector t 32 , ten :: Vector t
27 } 33 }
28 34
35rank :: Tensor t -> Int
29rank = length . dims 36rank = length . dims
30 37
31instance (Show a,Storable a) => Show (Tensor a) where 38instance (Show a,Storable a) => Show (Tensor a) where
@@ -33,41 +40,49 @@ instance (Show a,Storable a) => Show (Tensor a) where
33 show T {dims = ds, ten = t} = "("++shdims ds ++") "++ show (toList t) 40 show T {dims = ds, ten = t} = "("++shdims ds ++") "++ show (toList t)
34 41
35 42
36shdims [(n,(t,name))] = name ++ sym t ++"["++show n++"]" 43shdims :: [IdxDesc] -> String
44shdims [IdxDesc n t name] = name ++ sym t ++"["++show n++"]"
37 where sym Covariant = "_" 45 where sym Covariant = "_"
38 sym Contravariant = "^" 46 sym Contravariant = "^"
39shdims (d:ds) = shdims [d] ++ "><"++ shdims ds 47shdims (d:ds) = shdims [d] ++ "><"++ shdims ds
40 48
41 49
42 50findIdx :: (Field t) => IdxName -> Tensor t
51 -> (([IdxDesc], [IdxDesc]), Matrix t)
43findIdx name t = ((d1,d2),m) where 52findIdx name t = ((d1,d2),m) where
44 (d1,d2) = span (\(_,(_,n)) -> n /=name) (dims t) 53 (d1,d2) = span (\d -> idxName d /= name) (dims t)
45 c = product (map fst d2) 54 c = product (map idxDim d2)
46 m = matrixFromVector RowMajor c (ten t) 55 m = matrixFromVector RowMajor c (ten t)
47 56
57putFirstIdx :: (Field t) => String -> Tensor t -> ([IdxDesc], Matrix t)
48putFirstIdx name t = (nd,m') 58putFirstIdx name t = (nd,m')
49 where ((d1,d2),m) = findIdx name t 59 where ((d1,d2),m) = findIdx name t
50 m' = matrixFromVector RowMajor c $ cdat $ trans m 60 m' = matrixFromVector RowMajor c $ cdat $ trans m
51 nd = d2++d1 61 nd = d2++d1
52 c = dim (ten t) `div` (fst $ head d2) 62 c = dim (ten t) `div` (idxDim $ head d2)
53 63
64part :: (Field t) => Tensor t -> (IdxName, Int) -> Tensor t
54part t (name,k) = if k<0 || k>=l 65part t (name,k) = if k<0 || k>=l
55 then error $ "part "++show (name,k)++" out of range in "++show t 66 then error $ "part "++show (name,k)++" out of range" -- in "++show t
56 else T {dims = ds, ten = toRows m !! k} 67 else T {dims = ds, ten = toRows m !! k}
57 where (d:ds,m) = putFirstIdx name t 68 where (d:ds,m) = putFirstIdx name t
58 (l,_) = d 69 l = idxDim d
59 70
71parts :: (Field t) => Tensor t -> IdxName -> [Tensor t]
60parts t name = map f (toRows m) 72parts t name = map f (toRows m)
61 where (d:ds,m) = putFirstIdx name t 73 where (d:ds,m) = putFirstIdx name t
62 (l,_) = d 74 l = idxDim d
63 f t = T {dims=ds, ten=t} 75 f t = T {dims=ds, ten=t}
64 76
77concatRename :: [IdxDesc] -> [IdxDesc] -> [IdxDesc]
65concatRename l1 l2 = l1 ++ map ren l2 where 78concatRename l1 l2 = l1 ++ map ren l2 where
66 ren (n,(t,s)) = if {- s `elem` fs -} True then (n,(t,s++"'")) else (n,(t,s)) 79 ren idx = if {- s `elem` fs -} True then idx {idxName = idxName idx ++ "'"} else idx
67 fs = map (snd.snd) l1 80 fs = map idxName l1
68 81
69prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer v1 v2) 82prod :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t
83prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer' v1 v2)
70 84
85contraction :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
71contraction t1 n1 t2 n2 = 86contraction t1 n1 t2 n2 =
72 if compatIdx t1 n1 t2 n2 87 if compatIdx t1 n1 t2 n2
73 then T (concatRename (tail d1) (tail d2)) (cdat m) 88 then T (concatRename (tail d1) (tail d2)) (cdat m)
@@ -76,18 +91,22 @@ contraction t1 n1 t2 n2 =
76 (d2,m2) = putFirstIdx n2 t2 91 (d2,m2) = putFirstIdx n2 t2
77 m = multiply RowMajor (trans m1) m2 92 m = multiply RowMajor (trans m1) m2
78 93
94sumT :: (Storable t, Enum t, Num t) => [Tensor t] -> [t]
79sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls) 95sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls)
80 96
97contract1 :: (Num t, Enum t, Field t) => Tensor t -> IdxName -> IdxName -> Tensor t
81contract1 t name1 name2 = T d $ fromList $ sumT y 98contract1 t name1 name2 = T d $ fromList $ sumT y
82 where d = dims (head y) 99 where d = dims (head y)
83 x = (map (flip parts name2) (parts t name1)) 100 x = (map (flip parts name2) (parts t name1))
84 y = map head $ zipWith drop [0..] x 101 y = map head $ zipWith drop [0..] x
85 102
103contraction' :: (Field t, Enum t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
86contraction' t1 n1 t2 n2 = 104contraction' t1 n1 t2 n2 =
87 if compatIdx t1 n1 t2 n2 105 if compatIdx t1 n1 t2 n2
88 then contract1 (prod t1 t2) n1 (n2++"'") 106 then contract1 (prod t1 t2) n1 (n2++"'")
89 else error "wrong contraction'" 107 else error "wrong contraction'"
90 108
109tridx :: (Field t) => [IdxName] -> Tensor t -> Tensor t
91tridx [] t = t 110tridx [] t = t
92tridx (name:rest) t = T (d:ds) (join ts) where 111tridx (name:rest) t = T (d:ds) (join ts) where
93 ((_,d:_),_) = findIdx name t 112 ((_,d:_),_) = findIdx name t
@@ -95,30 +114,38 @@ tridx (name:rest) t = T (d:ds) (join ts) where
95 ts = map ten ps 114 ts = map ten ps
96 ds = dims (head ps) 115 ds = dims (head ps)
97 116
98compatIdxAux (n1,(t1,_)) (n2, (t2,_)) = t1 /= t2 && n1 == n2 117compatIdxAux :: IdxDesc -> IdxDesc -> Bool
118compatIdxAux IdxDesc {idxDim = n1, idxType = t1}
119 IdxDesc {idxDim = n2, idxType = t2}
120 = t1 /= t2 && n1 == n2
99 121
122compatIdx :: (Field t1, Field t) => Tensor t1 -> IdxName -> Tensor t -> IdxName -> Bool
100compatIdx t1 n1 t2 n2 = compatIdxAux d1 d2 where 123compatIdx t1 n1 t2 n2 = compatIdxAux d1 d2 where
101 d1 = head $ snd $ fst $ findIdx n1 t1 124 d1 = head $ snd $ fst $ findIdx n1 t1
102 d2 = head $ snd $ fst $ findIdx n2 t2 125 d2 = head $ snd $ fst $ findIdx n2 t2
103 126
104names t = sort $ map (snd.snd) (dims t) 127names :: Tensor t -> [IdxName]
128names t = sort $ map idxName (dims t)
105 129
130normal :: (Field t) => Tensor t -> Tensor t
106normal t = tridx (names t) t 131normal t = tridx (names t) t
107 132
133contractions :: (Num t, Field t) => Tensor t -> Tensor t -> [Tensor t]
108contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ] 134contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
109 135
110-- sent to Haskell-Cafe by Sebastian Sylvan 136-- sent to Haskell-Cafe by Sebastian Sylvan
137perms :: [t] -> [[t]]
111perms [x] = [[x]] 138perms [x] = [[x]]
112perms xs = [y:ps | (y,ys) <- selections xs , ps <- perms ys] 139perms xs = [y:ps | (y,ys) <- selections xs , ps <- perms ys]
113selections [] = [] 140selections [] = []
114selections (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- selections xs] 141selections (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- selections xs]
115 142
116 143interchanges :: (Ord a) => [a] -> Int
117interchanges ls = sum (map (count ls) ls) 144interchanges ls = sum (map (count ls) ls)
118 where count l p = n 145 where count l p = length $ filter (>p) $ take pel l
119 where Just pel = elemIndex p l 146 where Just pel = elemIndex p l
120 n = length $ filter (>p) $ take pel l
121 147
148signature :: (Num t, Ord a) => [a] -> t
122signature l | length (nub l) < length l = 0 149signature l | length (nub l) < length l = 0
123 | even (interchanges l) = 1 150 | even (interchanges l) = 1
124 | otherwise = -1 151 | 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)
103asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v), ptr = castPtr (ptr v) } 103asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v), ptr = castPtr (ptr v) }
104 104
105 105
106constantG n x = fromList (replicate n x) 106constantG x n = fromList (replicate n x)
107 107
108constantR :: Int -> Double -> Vector Double 108constantR :: Double -> Int -> Vector Double
109constantR = constantAux cconstantR 109constantR = constantAux cconstantR
110 110
111constantC :: Int -> Complex Double -> Vector (Complex Double) 111constantC :: Complex Double -> Int -> Vector (Complex Double)
112constantC = constantAux cconstantC 112constantC = constantAux cconstantC
113 113
114constantAux fun n x = unsafePerformIO $ do 114constantAux fun x n = unsafePerformIO $ do
115 v <- createVector n 115 v <- createVector n
116 px <- newArray [x] 116 px <- newArray [x]
117 fun px // vec v // check "constantAux" [] 117 fun px // vec v // check "constantAux" []
@@ -124,8 +124,8 @@ foreign import ccall safe "aux.h constantR"
124foreign import ccall safe "aux.h constantC" 124foreign import ccall safe "aux.h constantC"
125 cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int 125 cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int
126 126
127constant :: Field a => Int -> a -> Vector a 127constant :: Field a => a -> Int -> Vector a
128constant n x | isReal id x = scast $ constantR n (scast x) 128constant x n | isReal id x = scast $ constantR (scast x) n
129 | isComp id x = scast $ constantC n (scast x) 129 | isComp id x = scast $ constantC (scast x) n
130 | otherwise = constantG n x 130 | otherwise = constantG x n
131 131
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 (
16 Matrix(rows,cols), Field, 16 Matrix(rows,cols), Field,
17 toLists, (><), (>|<), (@@>), 17 toLists, (><), (>|<), (@@>),
18 trans, 18 trans,
19 reshape, 19 reshape, flatten,
20 fromRows, toRows, fromColumns, toColumns, 20 fromRows, toRows, fromColumns, toColumns,
21 joinVert, joinHoriz, 21 joinVert, joinHoriz,
22 flipud, fliprl, 22 flipud, fliprl,
23 liftMatrix, liftMatrix2, 23 liftMatrix, liftMatrix2,
24 multiply, 24 multiply,
25 outer,
25 subMatrix, 26 subMatrix,
26 takeRows, dropRows, takeColumns, dropColumns, 27 takeRows, dropRows, takeColumns, dropColumns,
27 diag, takeDiag, diagRect, ident 28 diag, takeDiag, diagRect, ident
@@ -54,11 +55,11 @@ diagRect s r c
54 | r == c = diag s 55 | r == c = diag s
55 | r < c = trans $ diagRect s c r 56 | r < c = trans $ diagRect s c r
56 | r > c = joinVert [diag s , zeros (r-c,c)] 57 | r > c = joinVert [diag s , zeros (r-c,c)]
57 where zeros (r,c) = reshape c $ constant (r*c) 0 58 where zeros (r,c) = reshape c $ constant 0 (r*c)
58 59
59takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] 60takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
60 61
61ident n = diag (constant n 1) 62ident n = diag (constant 1 n)
62 63
63r >< c = f where 64r >< c = f where
64 f l | dim v == r*c = matrixFromVector RowMajor c v 65 f l | dim v == r*c = matrixFromVector RowMajor c v
@@ -88,3 +89,5 @@ dropColumns :: Field t => Int -> Matrix t -> Matrix t
88dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat 89dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat
89 90
90---------------------------------------------------------------- 91----------------------------------------------------------------
92
93flatten = 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 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Plot
4-- Copyright : (c) Alberto Ruiz 2005
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz (aruiz at um dot es)
8-- Stability : provisional
9-- Portability : uses gnuplot and ImageMagick
10--
11-- Very basic (and provisional) drawing tools.
12--
13-----------------------------------------------------------------------------
14
15module Data.Packed.Plot(
16
17 gnuplotX, mplot,
18
19 plot, parametricPlot,
20
21 splot, mesh, mesh', meshdom,
22
23 matrixToPGM, imshow,
24
25) where
26
27import Data.Packed.Vector
28import Data.Packed.Matrix
29import GSL.Vector(FunCodeS(Max,Min),toScalarR)
30import Data.List(intersperse)
31import System
32import Data.IORef
33import System.Exit
34import Foreign hiding (rotate)
35
36
37size = dim
38
39-- | Loads a real matrix from a formatted ASCII text file
40--fromFile :: FilePath -> IO Matrix
41--fromFile filename = readFile filename >>= return . readMatrix read
42
43-- | Saves a real matrix to a formatted ascii text file
44toFile :: FilePath -> Matrix Double -> IO ()
45toFile filename matrix = writeFile filename (unlines . map unwords. map (map show) . toLists $ matrix)
46
47------------------------------------------------------------------------
48
49
50-- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions.
51meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double)
52meshdom r1 r2 = (outer r1 (constant 1 (size r2)), outer (constant 1 (size r1)) r2)
53
54
55gnuplotX command = do {system cmdstr; return()} where
56 cmdstr = "echo \""++command++"\" | gnuplot -persist"
57
58datafollows = "\\\"-\\\""
59
60prep = (++"e\n\n") . unlines . map (unwords . (map show))
61
62
63{- | Draws a 3D surface representation of a real matrix.
64
65> > mesh (hilb 20)
66
67In certain versions you can interactively rotate the graphic using the mouse.
68
69-}
70mesh :: Matrix Double -> IO ()
71mesh m = gnuplotX (command++dat) where
72 command = "splot "++datafollows++" matrix with lines\n"
73 dat = prep $ toLists $ m
74
75mesh' m = do
76 writeFile "splot-gnu-command" "splot \"splot-tmp.txt\" matrix with lines; pause -1";
77 toFile "splot-tmp.txt" m
78 putStr "Press [Return] to close the graphic and continue... "
79 system "gnuplot -persist splot-gnu-command"
80 system "rm splot-tmp.txt splot-gnu-command"
81 return ()
82
83{- | Draws the surface represented by the function f in the desired ranges and number of points, internally using 'mesh'.
84
85> > let f x y = cos (x + y)
86> > splot f (0,pi) (0,2*pi) 50
87
88-}
89splot :: (Matrix Double->Matrix Double->Matrix Double) -> (Double,Double) -> (Double,Double) -> Int -> IO ()
90splot f rx ry n = mesh' z where
91 (x,y) = meshdom (linspace n rx) (linspace n ry)
92 z = f x y
93
94{- | plots several vectors against the first one -}
95mplot :: [Vector Double] -> IO ()
96mplot m = gnuplotX (commands++dats) where
97 commands = if length m == 1 then command1 else commandmore
98 command1 = "plot "++datafollows++" with lines\n" ++ dat
99 commandmore = "plot " ++ plots ++ "\n"
100 plots = concat $ intersperse ", " (map cmd [2 .. length m])
101 cmd k = datafollows++" using 1:"++show k++" with lines"
102 dat = prep $ toLists $ fromColumns m
103 dats = concat (replicate (length m-1) dat)
104
105
106
107
108
109
110mplot' m = do
111 writeFile "plot-gnu-command" (commands++endcmd)
112 toFile "plot-tmp.txt" (fromColumns m)
113 putStr "Press [Return] to close the graphic and continue... "
114 system "gnuplot plot-gnu-command"
115 system "rm plot-tmp.txt plot-gnu-command"
116 return ()
117 where
118 commands = if length m == 1 then command1 else commandmore
119 command1 = "plot \"plot-tmp.txt\" with lines\n"
120 commandmore = "plot " ++ plots ++ "\n"
121 plots = concat $ intersperse ", " (map cmd [2 .. length m])
122 cmd k = "\"plot-tmp.txt\" using 1:"++show k++" with lines"
123 endcmd = "pause -1"
124
125-- apply several functions to one object
126mapf fs x = map ($ x) fs
127
128{- | Draws a list of functions over a desired range and with a desired number of points
129
130> > plot [sin, cos, sin.(3*)] (0,2*pi) 1000
131
132-}
133plot :: [Vector Double->Vector Double] -> (Double,Double) -> Int -> IO ()
134plot fs rx n = mplot (x: mapf fs x)
135 where x = linspace n rx
136
137{- | Draws a parametric curve. For instance, to draw a spiral we can do something like:
138
139> > parametricPlot (\t->(t * sin t, t * cos t)) (0,10*pi) 1000
140
141-}
142parametricPlot :: (Vector Double->(Vector Double,Vector Double)) -> (Double, Double) -> Int -> IO ()
143parametricPlot f rt n = mplot [fx, fy]
144 where t = linspace n rt
145 (fx,fy) = f t
146
147
148-- | writes a matrix to pgm image file
149matrixToPGM :: Matrix Double -> String
150matrixToPGM m = header ++ unlines (map unwords ll) where
151 c = cols m
152 r = rows m
153 header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"
154 maxgray = 255.0
155 maxval = toScalarR Max $ flatten $ m
156 minval = toScalarR Min $ flatten $ m
157 scale = if (maxval == minval)
158 then 0.0
159 else maxgray / (maxval - minval)
160 f x = show ( round ( scale *(x - minval) ) :: Int )
161 ll = map (map f) (toLists m)
162
163-- | imshow shows a representation of a matrix as a gray level image using ImageMagick's display.
164imshow :: Matrix Double -> IO ()
165imshow m = do
166 system $ "echo \""++ matrixToPGM m ++"\"| display -antialias -resize 300 - &"
167 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 @@
1 1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Tensor
4-- Copyright : (c) Alberto Ruiz 2007
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional
9-- Portability : portable
10--
11-- Tensors
12--
13-----------------------------------------------------------------------------
14
15module Data.Packed.Tensor (
16
17) where
18
19import Data.Packed.Internal
20import 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 (
20 constant, 20 constant,
21 toComplex, comp, 21 toComplex, comp,
22 conj, 22 conj,
23 dot 23 dot,
24 linspace
24) where 25) where
25 26
26import Data.Packed.Internal 27import Data.Packed.Internal
@@ -35,6 +36,14 @@ conj :: Vector (Complex Double) -> Vector (Complex Double)
35conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1]) 36conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1])
36 where mulC = multiply RowMajor 37 where mulC = multiply RowMajor
37 38
38comp v = toComplex (v,constant (dim v) 0) 39comp v = toComplex (v,constant 0 (dim v))
39 40
41{- | Creates a real vector containing a range of values:
40 42
43> > linspace 10 (-2,2)
44>-2. -1.556 -1.111 -0.667 -0.222 0.222 0.667 1.111 1.556 2.
45
46-}
47linspace :: Int -> (Double, Double) -> Vector Double
48linspace n (a,b) = fromList [a::Double,a+delta .. b]
49 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
29 c_fft code // vec v // vec r // check "fft" [v] 29 c_fft code // vec v // vec r // check "fft" [v]
30 return r 30 return r
31 31
32foreign import ccall "gsl-aux.h fft" c_fft :: Int -> TCVCV -- Complex Double :> Complex Double :> IO Int 32foreign import ccall "gsl-aux.h fft" c_fft :: Int -> TCVCV
33 33
34 34
35{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave. 35{- | 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
97 97
98foreign import ccall "gsl-aux.h minimize" 98foreign import ccall "gsl-aux.h minimize"
99 c_minimizeNMSimplex:: FunPtr (Int -> Ptr Double -> Double) -> Double -> Int 99 c_minimizeNMSimplex:: FunPtr (Int -> Ptr Double -> Double) -> Double -> Int
100 -> TVVM -- Double :> Double :> Double ::> IO Int 100 -> TVVM
101 101
102---------------------------------------------------------------------------------- 102----------------------------------------------------------------------------------
103 103
@@ -165,7 +165,7 @@ foreign import ccall "gsl-aux.h minimizeWithDeriv"
165 c_minimizeConjugateGradient :: FunPtr (Int -> Ptr Double -> Double) 165 c_minimizeConjugateGradient :: FunPtr (Int -> Ptr Double -> Double)
166 -> FunPtr (Int -> Ptr Double -> Ptr Double -> IO ()) 166 -> FunPtr (Int -> Ptr Double -> Ptr Double -> IO ())
167 -> Double -> Double -> Double -> Int 167 -> Double -> Double -> Double -> Int
168 -> TVM -- Double :> Double ::> IO Int 168 -> TVM
169 169
170--------------------------------------------------------------------- 170---------------------------------------------------------------------
171iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double) 171iv :: (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
51 return r 51 return r
52 | otherwise = error "polySolve on a polynomial of degree zero" 52 | otherwise = error "polySolve on a polynomial of degree zero"
53 53
54foreign import ccall "gsl-aux.h polySolve" c_polySolve:: TVCV -- Double :> Complex Double :> IO Int 54foreign 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
121toScalarR :: FunCodeS -> Vector Double -> Double 121toScalarR :: FunCodeS -> Vector Double -> Double
122toScalarR oper = toScalarAux c_toScalarR (fromEnum oper) 122toScalarR oper = toScalarAux c_toScalarR (fromEnum oper)
123 123
124foreign import ccall safe "gsl-aux.h toScalarR" 124foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: Int -> TVV
125 c_toScalarR :: Int -> TVV -- Double :> Double :> IO Int
126 125
127------------------------------------------------------------------ 126------------------------------------------------------------------
128 127
@@ -130,15 +129,13 @@ foreign import ccall safe "gsl-aux.h toScalarR"
130vectorMapR :: FunCodeV -> Vector Double -> Vector Double 129vectorMapR :: FunCodeV -> Vector Double -> Vector Double
131vectorMapR = vectorMapAux c_vectorMapR 130vectorMapR = vectorMapAux c_vectorMapR
132 131
133foreign import ccall safe "gsl-aux.h mapR" 132foreign import ccall safe "gsl-aux.h mapR" c_vectorMapR :: Int -> TVV
134 c_vectorMapR :: Int -> TVV -- Double :> Double :> IO Int
135 133
136-- | map of complex vectors with given function 134-- | map of complex vectors with given function
137vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double) 135vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double)
138vectorMapC oper = vectorMapAux c_vectorMapC (fromEnum oper) 136vectorMapC oper = vectorMapAux c_vectorMapC (fromEnum oper)
139 137
140foreign import ccall safe "gsl-aux.h mapC" 138foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: Int -> TCVCV
141 c_vectorMapC :: Int -> TCVCV -- Complex Double :> Complex Double :> IO Int
142 139
143------------------------------------------------------------------- 140-------------------------------------------------------------------
144 141
@@ -146,15 +143,13 @@ foreign import ccall safe "gsl-aux.h mapC"
146vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double 143vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double
147vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromEnum oper) 144vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromEnum oper)
148 145
149foreign import ccall safe "gsl-aux.h mapValR" 146foreign import ccall safe "gsl-aux.h mapValR" c_vectorMapValR :: Int -> Ptr Double -> TVV
150 c_vectorMapValR :: Int -> Ptr Double -> TVV -- Double :> Double :> IO Int
151 147
152-- | map of complex vectors with given function 148-- | map of complex vectors with given function
153vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double) 149vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
154vectorMapValC = vectorMapValAux c_vectorMapValC 150vectorMapValC = vectorMapValAux c_vectorMapValC
155 151
156foreign import ccall safe "gsl-aux.h mapValC" 152foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: Int -> Ptr (Complex Double) -> TCVCV
157 c_vectorMapValC :: Int -> Ptr (Complex Double) -> TCVCV -- Complex Double :> Complex Double :> IO Int
158 153
159------------------------------------------------------------------- 154-------------------------------------------------------------------
160 155
@@ -162,14 +157,10 @@ foreign import ccall safe "gsl-aux.h mapValC"
162vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double 157vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double
163vectorZipR = vectorZipAux c_vectorZipR 158vectorZipR = vectorZipAux c_vectorZipR
164 159
165foreign import ccall safe "gsl-aux.h zipR" 160foreign import ccall safe "gsl-aux.h zipR" c_vectorZipR :: Int -> TVVV
166 c_vectorZipR :: Int -> TVVV -- Double :> Double :> Double :> IO Int
167 161
168-- | elementwise operation on complex vectors 162-- | elementwise operation on complex vectors
169vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) 163vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
170vectorZipC = vectorZipAux c_vectorZipC 164vectorZipC = vectorZipAux c_vectorZipC
171 165
172foreign import ccall safe "gsl-aux.h zipC" 166foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: Int -> TCVCVCV
173 c_vectorZipC :: Int -> TCVCVCV -- Complex Double :> Complex Double :> Complex Double :> IO Int
174
175
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 @@
1{- |
2
3Module : GSL
4Copyright : (c) Alberto Ruiz 2006-7
5License : GPL-style
6
7Maintainer : Alberto Ruiz (aruiz at um dot es)
8Stability : provisional
9Portability : uses -fffi and -fglasgow-exts
10
11This module reexports the basic functionality and a collection of utilities.
12
13-}
14
15module HSSL (
16
17) 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
26import Data.Packed.Internal.Matrix 26import Data.Packed.Internal.Matrix
27import Data.Packed.Vector 27import Data.Packed.Vector
28import Data.Packed.Matrix 28import Data.Packed.Matrix
29import GSL.Vector(scale)
29import Complex 30import Complex
30import Foreign 31import Foreign
31 32
32----------------------------------------------------------------------------- 33-----------------------------------------------------------------------------
33-- dgesvd 34foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM
34foreign import ccall "LAPACK/lapack-aux.h svd_l_R"
35 dgesvd :: TMMVM -- Double ::> Double ::> (Double :> Double ::> IO Int)
36 35
37-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. 36-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
38-- 37--
@@ -48,9 +47,7 @@ svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
48 return (u,s,trans v) 47 return (u,s,trans v)
49 48
50----------------------------------------------------------------------------- 49-----------------------------------------------------------------------------
51-- dgesdd 50foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM
52foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd"
53 dgesdd :: TMMVM --Double ::> Double ::> (Double :> Double ::> IO Int)
54 51
55-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix. 52-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
56-- 53--
@@ -66,9 +63,7 @@ svdRdd' x@M {rows = r, cols = c} = unsafePerformIO $ do
66 return (u,s,trans v) 63 return (u,s,trans v)
67 64
68----------------------------------------------------------------------------- 65-----------------------------------------------------------------------------
69-- zgesvd 66foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
70foreign import ccall "LAPACK/lapack-aux.h svd_l_C"
71 zgesvd :: TCMCMVCM -- (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
72 67
73-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix. 68-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
74-- 69--
@@ -85,9 +80,7 @@ svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
85 return (u,s,trans v) 80 return (u,s,trans v)
86 81
87----------------------------------------------------------------------------- 82-----------------------------------------------------------------------------
88-- zgeev 83foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
89foreign import ccall "LAPACK/lapack-aux.h eig_l_C"
90 zgeev :: TCMCMCVCM -- (Complex Double) ::> (Complex Double) ::> ((Complex Double) :> (Complex Double) ::> IO Int)
91 84
92-- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix: 85-- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
93-- 86--
@@ -106,9 +99,7 @@ eigC (m@M {rows = r})
106 return (l,v) 99 return (l,v)
107 100
108----------------------------------------------------------------------------- 101-----------------------------------------------------------------------------
109-- dgeev 102foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
110foreign import ccall "LAPACK/lapack-aux.h eig_l_R"
111 dgeev :: TMMCVM -- Double ::> Double ::> ((Complex Double) :> Double ::> IO Int)
112 103
113-- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix: 104-- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
114-- 105--
@@ -139,12 +130,10 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
139 | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs 130 | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs
140 | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs) 131 | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs)
141 132
142scale r v = fromList [r] `outer` v 133-- scale r v = fromList [r] `outer` v
143 134
144----------------------------------------------------------------------------- 135-----------------------------------------------------------------------------
145-- dsyev 136foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
146foreign import ccall "LAPACK/lapack-aux.h eig_l_S"
147 dsyev :: TMVM -- Double ::> (Double :> Double ::> IO Int)
148 137
149-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix: 138-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
150-- 139--
@@ -166,9 +155,7 @@ eigS' (m@M {rows = r})
166 return (l,v) 155 return (l,v)
167 156
168----------------------------------------------------------------------------- 157-----------------------------------------------------------------------------
169-- zheev 158foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
170foreign import ccall "LAPACK/lapack-aux.h eig_l_H"
171 zheev :: TCMVCM -- (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
172 159
173-- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix: 160-- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix:
174-- 161--
@@ -190,9 +177,7 @@ eigH' (m@M {rows = r})
190 return (l,v) 177 return (l,v)
191 178
192----------------------------------------------------------------------------- 179-----------------------------------------------------------------------------
193-- dgesv 180foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
194foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l"
195 dgesv :: TMMM -- Double ::> Double ::> Double ::> IO Int
196 181
197-- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition. 182-- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition.
198linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double 183linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
@@ -204,9 +189,7 @@ linearSolveR a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
204 | otherwise = error "linearSolveR of nonsquare matrix" 189 | otherwise = error "linearSolveR of nonsquare matrix"
205 190
206----------------------------------------------------------------------------- 191-----------------------------------------------------------------------------
207-- zgesv 192foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
208foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l"
209 zgesv :: TCMCMCM -- (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
210 193
211-- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition. 194-- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition.
212linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 195linearSolveC :: 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})
218 | otherwise = error "linearSolveC of nonsquare matrix" 201 | otherwise = error "linearSolveC of nonsquare matrix"
219 202
220----------------------------------------------------------------------------------- 203-----------------------------------------------------------------------------------
221-- dgels 204foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM
222foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l"
223 dgels :: TMMM -- Double ::> Double ::> Double ::> IO Int
224 205
225-- | 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'. 206-- | 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'.
226linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double 207linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
@@ -232,9 +213,7 @@ linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformI
232 return r 213 return r
233 214
234----------------------------------------------------------------------------------- 215-----------------------------------------------------------------------------------
235-- zgels 216foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM
236foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l"
237 zgels :: TCMCMCM -- (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
238 217
239-- | 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'. 218-- | 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'.
240linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 219linearSolveLSC :: 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
246 return r 225 return r
247 226
248----------------------------------------------------------------------------------- 227-----------------------------------------------------------------------------------
249-- dgelss 228foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM
250foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l"
251 dgelss :: Double -> TMMM -- Double ::> Double ::> Double ::> IO Int
252 229
253-- | 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. 230-- | 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.
254linearSolveSVDR :: Maybe Double -- ^ rcond 231linearSolveSVDR :: Maybe Double -- ^ rcond
@@ -264,9 +241,7 @@ linearSolveSVDR_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafeP
264 return r 241 return r
265 242
266----------------------------------------------------------------------------------- 243-----------------------------------------------------------------------------------
267-- zgelss 244foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
268foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l"
269 zgelss :: Double -> TCMCMCM -- (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
270 245
271-- | 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. 246-- | 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.
272linearSolveSVDC :: Maybe Double -- ^ rcond 247linearSolveSVDC :: 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 @@
1-----------------------------------------------------------------------------
2{- |
3Module : LinearAlgebra
4Copyright : (c) Alberto Ruiz 2006-7
5License : GPL-style
6
7Maintainer : Alberto Ruiz (aruiz at um dot es)
8Stability : provisional
9Portability : uses ffi
10
11Some linear algebra algorithms, implemented by means of BLAS, LAPACK or GSL.
12
13-}
14-----------------------------------------------------------------------------
15module LinearAlgebra (
16
17) 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 @@
1-----------------------------------------------------------------------------
2{- |
3Module : LinearAlgebra.Algorithms
4Copyright : (c) Alberto Ruiz 2006-7
5License : GPL-style
6
7Maintainer : Alberto Ruiz (aruiz at um dot es)
8Stability : provisional
9Portability : uses ffi
10
11
12-}
13-----------------------------------------------------------------------------
14
15module LinearAlgebra.Algorithms (
16
17) where