diff options
-rw-r--r-- | examples/tests.hs | 6 | ||||
-rw-r--r-- | lib/Data/Packed.hs | 2 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 36 | ||||
-rw-r--r-- | lib/Data/Packed/Matrix.hs | 44 | ||||
-rw-r--r-- | lib/GSLHaskell.hs | 327 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 46 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Instances.hs | 8 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Interface.hs | 14 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Linear.hs | 6 |
9 files changed, 81 insertions, 408 deletions
diff --git a/examples/tests.hs b/examples/tests.hs index 3fea81e..ad70266 100644 --- a/examples/tests.hs +++ b/examples/tests.hs | |||
@@ -39,7 +39,7 @@ instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where | |||
39 | return (r:+i) | 39 | return (r:+i) |
40 | coarbitrary = undefined | 40 | coarbitrary = undefined |
41 | 41 | ||
42 | instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where | 42 | instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where |
43 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | 43 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) |
44 | m <- choose (1,maxdim) | 44 | m <- choose (1,maxdim) |
45 | n <- choose (1,maxdim) | 45 | n <- choose (1,maxdim) |
@@ -52,7 +52,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where | |||
52 | coarbitrary = undefined | 52 | coarbitrary = undefined |
53 | 53 | ||
54 | data PairM a = PairM (Matrix a) (Matrix a) deriving Show | 54 | data PairM a = PairM (Matrix a) (Matrix a) deriving Show |
55 | instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where | 55 | instance (Num a, Element a, Arbitrary a) => Arbitrary (PairM a) where |
56 | arbitrary = do | 56 | arbitrary = do |
57 | a <- choose (1,maxdim) | 57 | a <- choose (1,maxdim) |
58 | b <- choose (1,maxdim) | 58 | b <- choose (1,maxdim) |
@@ -65,7 +65,7 @@ instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where | |||
65 | 65 | ||
66 | data SqM a = SqM (Matrix a) deriving Show | 66 | data SqM a = SqM (Matrix a) deriving Show |
67 | sqm (SqM a) = a | 67 | sqm (SqM a) = a |
68 | instance (Field a, Arbitrary a) => Arbitrary (SqM a) where | 68 | instance (Element a, Arbitrary a) => Arbitrary (SqM a) where |
69 | arbitrary = do | 69 | arbitrary = do |
70 | n <- choose (1,maxdim) | 70 | n <- choose (1,maxdim) |
71 | l <- vector (n*n) | 71 | l <- vector (n*n) |
diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs index 668d2f7..53aced9 100644 --- a/lib/Data/Packed.hs +++ b/lib/Data/Packed.hs | |||
@@ -27,7 +27,7 @@ import Data.Complex | |||
27 | import Data.Packed.Internal | 27 | import Data.Packed.Internal |
28 | 28 | ||
29 | -- | conversion utilities | 29 | -- | conversion utilities |
30 | class (Field e) => Container c e where | 30 | class (Element e) => Container c e where |
31 | toComplex :: RealFloat e => (c e, c e) -> c (Complex e) | 31 | toComplex :: RealFloat e => (c e, c e) -> c (Complex e) |
32 | fromComplex :: RealFloat e => c (Complex e) -> (c e, c e) | 32 | fromComplex :: RealFloat e => c (Complex e) -> (c e, c e) |
33 | comp :: RealFloat e => c e -> c (Complex e) | 33 | comp :: RealFloat e => c e -> c (Complex e) |
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index f63ee52..fbab33c 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -84,17 +84,17 @@ type Mt t s = Int -> Int -> Ptr t -> s | |||
84 | -- type t ::> s = Mt t s | 84 | -- type t ::> s = Mt t s |
85 | 85 | ||
86 | -- | the inverse of 'Data.Packed.Matrix.fromLists' | 86 | -- | the inverse of 'Data.Packed.Matrix.fromLists' |
87 | toLists :: (Field t) => Matrix t -> [[t]] | 87 | toLists :: (Element t) => Matrix t -> [[t]] |
88 | toLists m = partit (cols m) . toList . cdat $ m | 88 | toLists m = partit (cols m) . toList . cdat $ m |
89 | 89 | ||
90 | -- | creates a Matrix from a list of vectors | 90 | -- | creates a Matrix from a list of vectors |
91 | fromRows :: Field t => [Vector t] -> Matrix t | 91 | fromRows :: Element t => [Vector t] -> Matrix t |
92 | fromRows vs = case common dim vs of | 92 | fromRows vs = case common dim vs of |
93 | Nothing -> error "fromRows applied to [] or to vectors with different sizes" | 93 | Nothing -> error "fromRows applied to [] or to vectors with different sizes" |
94 | Just c -> reshape c (join vs) | 94 | Just c -> reshape c (join vs) |
95 | 95 | ||
96 | -- | extracts the rows of a matrix as a list of vectors | 96 | -- | extracts the rows of a matrix as a list of vectors |
97 | toRows :: Field t => Matrix t -> [Vector t] | 97 | toRows :: Element t => Matrix t -> [Vector t] |
98 | toRows m = toRows' 0 where | 98 | toRows m = toRows' 0 where |
99 | v = cdat m | 99 | v = cdat m |
100 | r = rows m | 100 | r = rows m |
@@ -103,11 +103,11 @@ toRows m = toRows' 0 where | |||
103 | | otherwise = subVector k c v : toRows' (k+c) | 103 | | otherwise = subVector k c v : toRows' (k+c) |
104 | 104 | ||
105 | -- | Creates a matrix from a list of vectors, as columns | 105 | -- | Creates a matrix from a list of vectors, as columns |
106 | fromColumns :: Field t => [Vector t] -> Matrix t | 106 | fromColumns :: Element t => [Vector t] -> Matrix t |
107 | fromColumns m = trans . fromRows $ m | 107 | fromColumns m = trans . fromRows $ m |
108 | 108 | ||
109 | -- | Creates a list of vectors from the columns of a matrix | 109 | -- | Creates a list of vectors from the columns of a matrix |
110 | toColumns :: Field t => Matrix t -> [Vector t] | 110 | toColumns :: Element t => Matrix t -> [Vector t] |
111 | toColumns m = toRows . trans $ m | 111 | toColumns m = toRows . trans $ m |
112 | 112 | ||
113 | 113 | ||
@@ -152,18 +152,18 @@ where r is the desired number of rows.) | |||
152 | , 9.0, 10.0, 11.0, 12.0 ]@ | 152 | , 9.0, 10.0, 11.0, 12.0 ]@ |
153 | 153 | ||
154 | -} | 154 | -} |
155 | reshape :: Field t => Int -> Vector t -> Matrix t | 155 | reshape :: Element t => Int -> Vector t -> Matrix t |
156 | reshape c v = matrixFromVector RowMajor c v | 156 | reshape c v = matrixFromVector RowMajor c v |
157 | 157 | ||
158 | singleton x = reshape 1 (fromList [x]) | 158 | singleton x = reshape 1 (fromList [x]) |
159 | 159 | ||
160 | -- | application of a vector function on the flattened matrix elements | 160 | -- | application of a vector function on the flattened matrix elements |
161 | liftMatrix :: (Field a, Field b) => (Vector a -> Vector b) -> Matrix a -> Matrix b | 161 | liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b |
162 | liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor c (f d) | 162 | liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor c (f d) |
163 | liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) | 163 | liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) |
164 | 164 | ||
165 | -- | application of a vector function on the flattened matrices elements | 165 | -- | application of a vector function on the flattened matrices elements |
166 | liftMatrix2 :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t | 166 | liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t |
167 | liftMatrix2 f m1 m2 | 167 | liftMatrix2 f m1 m2 |
168 | | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2" | 168 | | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2" |
169 | | otherwise = case m1 of | 169 | | otherwise = case m1 of |
@@ -176,8 +176,8 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 | |||
176 | 176 | ||
177 | ---------------------------------------------------------------- | 177 | ---------------------------------------------------------------- |
178 | 178 | ||
179 | -- | Optimized matrix computations are provided for elements in the Field class. | 179 | -- | Optimized matrix computations are provided for elements in the Element class. |
180 | class (Storable a, Floating a) => Field a where | 180 | class (Storable a, Floating a) => Element a where |
181 | constantD :: a -> Int -> Vector a | 181 | constantD :: a -> Int -> Vector a |
182 | transdata :: Int -> Vector a -> Int -> Vector a | 182 | transdata :: Int -> Vector a -> Int -> Vector a |
183 | multiplyD :: Matrix a -> Matrix a -> Matrix a | 183 | multiplyD :: Matrix a -> Matrix a -> Matrix a |
@@ -186,14 +186,14 @@ class (Storable a, Floating a) => Field a where | |||
186 | -> Matrix a -> Matrix a | 186 | -> Matrix a -> Matrix a |
187 | diagD :: Vector a -> Matrix a | 187 | diagD :: Vector a -> Matrix a |
188 | 188 | ||
189 | instance Field Double where | 189 | instance Element Double where |
190 | constantD = constantR | 190 | constantD = constantR |
191 | transdata = transdataR | 191 | transdata = transdataR |
192 | multiplyD = multiplyR | 192 | multiplyD = multiplyR |
193 | subMatrixD = subMatrixR | 193 | subMatrixD = subMatrixR |
194 | diagD = diagR | 194 | diagD = diagR |
195 | 195 | ||
196 | instance Field (Complex Double) where | 196 | instance Element (Complex Double) where |
197 | constantD = constantC | 197 | constantD = constantC |
198 | transdata = transdataC | 198 | transdata = transdataC |
199 | multiplyD = multiplyC | 199 | multiplyD = multiplyC |
@@ -202,7 +202,7 @@ instance Field (Complex Double) where | |||
202 | 202 | ||
203 | ------------------------------------------------------------------ | 203 | ------------------------------------------------------------------ |
204 | 204 | ||
205 | (>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a | 205 | (>|<) :: (Element a) => Int -> Int -> [a] -> Matrix a |
206 | r >|< c = f where | 206 | r >|< c = f where |
207 | f l | dim v == r*c = matrixFromVector ColumnMajor c v | 207 | f l | dim v == r*c = matrixFromVector ColumnMajor c v |
208 | | otherwise = error $ "inconsistent list size = " | 208 | | otherwise = error $ "inconsistent list size = " |
@@ -260,13 +260,13 @@ foreign import ccall safe "auxi.h multiplyC" | |||
260 | -> Int -> Int -> Ptr (Complex Double) | 260 | -> Int -> Int -> Ptr (Complex Double) |
261 | -> IO Int | 261 | -> IO Int |
262 | 262 | ||
263 | multiply' :: (Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a | 263 | multiply' :: (Element a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a |
264 | multiply' RowMajor a b = multiplyD a b | 264 | multiply' RowMajor a b = multiplyD a b |
265 | multiply' ColumnMajor a b = trans $ multiplyD (trans b) (trans a) | 265 | multiply' ColumnMajor a b = trans $ multiplyD (trans b) (trans a) |
266 | 266 | ||
267 | 267 | ||
268 | -- | matrix product | 268 | -- | matrix product |
269 | multiply :: (Field a) => Matrix a -> Matrix a -> Matrix a | 269 | multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a |
270 | multiply = multiplyD | 270 | multiply = multiplyD |
271 | 271 | ||
272 | ---------------------------------------------------------------------- | 272 | ---------------------------------------------------------------------- |
@@ -287,7 +287,7 @@ subMatrixC (r0,c0) (rt,ct) x = | |||
287 | reshape (2*cols x) . asReal . cdat $ x | 287 | reshape (2*cols x) . asReal . cdat $ x |
288 | 288 | ||
289 | -- | Extracts a submatrix from a matrix. | 289 | -- | Extracts a submatrix from a matrix. |
290 | subMatrix :: Field a | 290 | subMatrix :: Element a |
291 | => (Int,Int) -- ^ (r0,c0) starting position | 291 | => (Int,Int) -- ^ (r0,c0) starting position |
292 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | 292 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix |
293 | -> Matrix a -- ^ input matrix | 293 | -> Matrix a -- ^ input matrix |
@@ -313,7 +313,7 @@ diagC = diagAux c_diagC "diagC" | |||
313 | foreign import ccall "auxi.h diagC" c_diagC :: TCVCM | 313 | foreign import ccall "auxi.h diagC" c_diagC :: TCVCM |
314 | 314 | ||
315 | -- | creates a square matrix with the given diagonal | 315 | -- | creates a square matrix with the given diagonal |
316 | diag :: Field a => Vector a -> Matrix a | 316 | diag :: Element a => Vector a -> Matrix a |
317 | diag = diagD | 317 | diag = diagD |
318 | 318 | ||
319 | ------------------------------------------------------------------------ | 319 | ------------------------------------------------------------------------ |
@@ -340,7 +340,7 @@ foreign import ccall safe "auxi.h constantC" | |||
340 | @> constant 2 7 | 340 | @> constant 2 7 |
341 | 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ | 341 | 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ |
342 | -} | 342 | -} |
343 | constant :: Field a => a -> Int -> Vector a | 343 | constant :: Element a => a -> Int -> Vector a |
344 | constant = constantD | 344 | constant = constantD |
345 | 345 | ||
346 | -------------------------------------------------------------------------- | 346 | -------------------------------------------------------------------------- |
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs index a705975..e96500f 100644 --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs | |||
@@ -14,7 +14,7 @@ | |||
14 | ----------------------------------------------------------------------------- | 14 | ----------------------------------------------------------------------------- |
15 | 15 | ||
16 | module Data.Packed.Matrix ( | 16 | module Data.Packed.Matrix ( |
17 | Field, | 17 | Element, |
18 | Matrix,rows,cols, | 18 | Matrix,rows,cols, |
19 | (><), | 19 | (><), |
20 | trans, | 20 | trans, |
@@ -41,13 +41,13 @@ import Data.List(transpose,intersperse) | |||
41 | import Data.Array | 41 | import Data.Array |
42 | 42 | ||
43 | -- | creates a matrix from a vertical list of matrices | 43 | -- | creates a matrix from a vertical list of matrices |
44 | joinVert :: Field t => [Matrix t] -> Matrix t | 44 | joinVert :: Element t => [Matrix t] -> Matrix t |
45 | joinVert ms = case common cols ms of | 45 | joinVert ms = case common cols ms of |
46 | Nothing -> error "joinVert on matrices with different number of columns" | 46 | Nothing -> error "joinVert on matrices with different number of columns" |
47 | Just c -> reshape c $ join (map cdat ms) | 47 | Just c -> reshape c $ join (map cdat ms) |
48 | 48 | ||
49 | -- | creates a matrix from a horizontal list of matrices | 49 | -- | creates a matrix from a horizontal list of matrices |
50 | joinHoriz :: Field t => [Matrix t] -> Matrix t | 50 | joinHoriz :: Element t => [Matrix t] -> Matrix t |
51 | joinHoriz ms = trans. joinVert . map trans $ ms | 51 | joinHoriz ms = trans. joinVert . map trans $ ms |
52 | 52 | ||
53 | {- | Creates a matrix from blocks given as a list of lists of matrices: | 53 | {- | Creates a matrix from blocks given as a list of lists of matrices: |
@@ -63,15 +63,15 @@ joinHoriz ms = trans. joinVert . map trans $ ms | |||
63 | , -1.0, -1.0, -1.0, -1.0, 0.0, 7.0, 0.0 | 63 | , -1.0, -1.0, -1.0, -1.0, 0.0, 7.0, 0.0 |
64 | , -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 2.0 ]@ | 64 | , -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 2.0 ]@ |
65 | -} | 65 | -} |
66 | fromBlocks :: Field t => [[Matrix t]] -> Matrix t | 66 | fromBlocks :: Element t => [[Matrix t]] -> Matrix t |
67 | fromBlocks = joinVert . map joinHoriz | 67 | fromBlocks = joinVert . map joinHoriz |
68 | 68 | ||
69 | -- | Reverse rows | 69 | -- | Reverse rows |
70 | flipud :: Field t => Matrix t -> Matrix t | 70 | flipud :: Element t => Matrix t -> Matrix t |
71 | flipud m = fromRows . reverse . toRows $ m | 71 | flipud m = fromRows . reverse . toRows $ m |
72 | 72 | ||
73 | -- | Reverse columns | 73 | -- | Reverse columns |
74 | fliprl :: Field t => Matrix t -> Matrix t | 74 | fliprl :: Element t => Matrix t -> Matrix t |
75 | fliprl m = fromColumns . reverse . toColumns $ m | 75 | fliprl m = fromColumns . reverse . toColumns $ m |
76 | 76 | ||
77 | ------------------------------------------------------------ | 77 | ------------------------------------------------------------ |
@@ -84,7 +84,7 @@ fliprl m = fromColumns . reverse . toColumns $ m | |||
84 | , 0.0, 5.0, 0.0, 0.0 | 84 | , 0.0, 5.0, 0.0, 0.0 |
85 | , 0.0, 0.0, 5.0, 0.0 ]@ | 85 | , 0.0, 0.0, 5.0, 0.0 ]@ |
86 | -} | 86 | -} |
87 | diagRect :: (Field t, Num t) => Vector t -> Int -> Int -> Matrix t | 87 | diagRect :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t |
88 | diagRect s r c | 88 | diagRect s r c |
89 | | dim s < min r c = error "diagRect" | 89 | | dim s < min r c = error "diagRect" |
90 | | r == c = diag s | 90 | | r == c = diag s |
@@ -93,11 +93,11 @@ diagRect s r c | |||
93 | where zeros (r,c) = reshape c $ constantD 0 (r*c) | 93 | where zeros (r,c) = reshape c $ constantD 0 (r*c) |
94 | 94 | ||
95 | -- | extracts the diagonal from a rectangular matrix | 95 | -- | extracts the diagonal from a rectangular matrix |
96 | takeDiag :: (Field t) => Matrix t -> Vector t | 96 | takeDiag :: (Element t) => Matrix t -> Vector t |
97 | takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] | 97 | takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] |
98 | 98 | ||
99 | -- | creates the identity matrix of given dimension | 99 | -- | creates the identity matrix of given dimension |
100 | ident :: Field a => Int -> Matrix a | 100 | ident :: Element a => Int -> Matrix a |
101 | ident n = diag (constant 1 n) | 101 | ident n = diag (constant 1 n) |
102 | 102 | ||
103 | ------------------------------------------------------------ | 103 | ------------------------------------------------------------ |
@@ -112,7 +112,7 @@ ident n = diag (constant 1 n) | |||
112 | This is the format produced by the instances of Show (Matrix a), which | 112 | This is the format produced by the instances of Show (Matrix a), which |
113 | can also be used for input. | 113 | can also be used for input. |
114 | -} | 114 | -} |
115 | (><) :: (Field a) => Int -> Int -> [a] -> Matrix a | 115 | (><) :: (Element a) => Int -> Int -> [a] -> Matrix a |
116 | r >< c = f where | 116 | r >< c = f where |
117 | f l | dim v == r*c = matrixFromVector RowMajor c v | 117 | f l | dim v == r*c = matrixFromVector RowMajor c v |
118 | | otherwise = error $ "inconsistent list size = " | 118 | | otherwise = error $ "inconsistent list size = " |
@@ -122,16 +122,16 @@ r >< c = f where | |||
122 | ---------------------------------------------------------------- | 122 | ---------------------------------------------------------------- |
123 | 123 | ||
124 | -- | Creates a matrix with the first n rows of another matrix | 124 | -- | Creates a matrix with the first n rows of another matrix |
125 | takeRows :: Field t => Int -> Matrix t -> Matrix t | 125 | takeRows :: Element t => Int -> Matrix t -> Matrix t |
126 | takeRows n mat = subMatrix (0,0) (n, cols mat) mat | 126 | takeRows n mat = subMatrix (0,0) (n, cols mat) mat |
127 | -- | Creates a copy of a matrix without the first n rows | 127 | -- | Creates a copy of a matrix without the first n rows |
128 | dropRows :: Field t => Int -> Matrix t -> Matrix t | 128 | dropRows :: Element t => Int -> Matrix t -> Matrix t |
129 | dropRows n mat = subMatrix (n,0) (rows mat - n, cols mat) mat | 129 | dropRows n mat = subMatrix (n,0) (rows mat - n, cols mat) mat |
130 | -- |Creates a matrix with the first n columns of another matrix | 130 | -- |Creates a matrix with the first n columns of another matrix |
131 | takeColumns :: Field t => Int -> Matrix t -> Matrix t | 131 | takeColumns :: Element t => Int -> Matrix t -> Matrix t |
132 | takeColumns n mat = subMatrix (0,0) (rows mat, n) mat | 132 | takeColumns n mat = subMatrix (0,0) (rows mat, n) mat |
133 | -- | Creates a copy of a matrix without the first n columns | 133 | -- | Creates a copy of a matrix without the first n columns |
134 | dropColumns :: Field t => Int -> Matrix t -> Matrix t | 134 | dropColumns :: Element t => Int -> Matrix t -> Matrix t |
135 | dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat | 135 | dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat |
136 | 136 | ||
137 | ---------------------------------------------------------------- | 137 | ---------------------------------------------------------------- |
@@ -141,7 +141,7 @@ dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat | |||
141 | @\> flatten ('ident' 3) | 141 | @\> flatten ('ident' 3) |
142 | 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ | 142 | 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ |
143 | -} | 143 | -} |
144 | flatten :: Field t => Matrix t -> Vector t | 144 | flatten :: Element t => Matrix t -> Vector t |
145 | flatten = cdat | 145 | flatten = cdat |
146 | 146 | ||
147 | {- | Creates a 'Matrix' from a list of lists (considered as rows). | 147 | {- | Creates a 'Matrix' from a list of lists (considered as rows). |
@@ -152,20 +152,20 @@ flatten = cdat | |||
152 | , 3.0, 4.0 | 152 | , 3.0, 4.0 |
153 | , 5.0, 6.0 ]@ | 153 | , 5.0, 6.0 ]@ |
154 | -} | 154 | -} |
155 | fromLists :: Field t => [[t]] -> Matrix t | 155 | fromLists :: Element t => [[t]] -> Matrix t |
156 | fromLists = fromRows . map fromList | 156 | fromLists = fromRows . map fromList |
157 | 157 | ||
158 | -- | creates a 1-row matrix from a vector | 158 | -- | creates a 1-row matrix from a vector |
159 | asRow :: Field a => Vector a -> Matrix a | 159 | asRow :: Element a => Vector a -> Matrix a |
160 | asRow v = reshape (dim v) v | 160 | asRow v = reshape (dim v) v |
161 | 161 | ||
162 | -- | creates a 1-column matrix from a vector | 162 | -- | creates a 1-column matrix from a vector |
163 | asColumn :: Field a => Vector a -> Matrix a | 163 | asColumn :: Element a => Vector a -> Matrix a |
164 | asColumn v = reshape 1 v | 164 | asColumn v = reshape 1 v |
165 | 165 | ||
166 | ----------------------------------------------------- | 166 | ----------------------------------------------------- |
167 | 167 | ||
168 | fromArray2D :: (Field e) => Array (Int, Int) e -> Matrix e | 168 | fromArray2D :: (Element e) => Array (Int, Int) e -> Matrix e |
169 | fromArray2D m = (r><c) (elems m) | 169 | fromArray2D m = (r><c) (elems m) |
170 | where ((r0,c0),(r1,c1)) = bounds m | 170 | where ((r0,c0),(r1,c1)) = bounds m |
171 | r = r1-r0+1 | 171 | r = r1-r0+1 |
@@ -201,7 +201,7 @@ this function the user can easily define any desired display function: | |||
201 | @disp = putStrLn . format \" \" (printf \"%.2f\")@ | 201 | @disp = putStrLn . format \" \" (printf \"%.2f\")@ |
202 | 202 | ||
203 | -} | 203 | -} |
204 | format :: (Field t) => String -> (t -> String) -> Matrix t -> String | 204 | format :: (Element t) => String -> (t -> String) -> Matrix t -> String |
205 | format sep f m = dsp' sep . map (map f) . toLists $ m | 205 | format sep f m = dsp' sep . map (map f) . toLists $ m |
206 | 206 | ||
207 | disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m | 207 | disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m |
@@ -217,7 +217,7 @@ readMatrix :: String -> Matrix Double | |||
217 | readMatrix = fromLists . map (map read). map words . filter (not.null) . lines | 217 | readMatrix = fromLists . map (map read). map words . filter (not.null) . lines |
218 | 218 | ||
219 | -- | rearranges the rows of a matrix according to the order given in a list of integers. | 219 | -- | rearranges the rows of a matrix according to the order given in a list of integers. |
220 | extractRows :: Field t => [Int] -> Matrix t -> Matrix t | 220 | extractRows :: Element t => [Int] -> Matrix t -> Matrix t |
221 | extractRows l m = fromRows $ extract (toRows $ m) l | 221 | extractRows l m = fromRows $ extract (toRows $ m) l |
222 | where extract l is = [l!!i |i<-is] | 222 | where extract l is = [l!!i |i<-is] |
223 | 223 | ||
@@ -231,5 +231,5 @@ extractRows l m = fromRows $ extract (toRows $ m) l | |||
231 | , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]@ | 231 | , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]@ |
232 | 232 | ||
233 | -} | 233 | -} |
234 | repmat :: (Field t) => Matrix t -> Int -> Int -> Matrix t | 234 | repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t |
235 | repmat m r c = fromBlocks $ partit c $ replicate (r*c) m \ No newline at end of file | 235 | repmat m r c = fromBlocks $ partit c $ replicate (r*c) m \ No newline at end of file |
diff --git a/lib/GSLHaskell.hs b/lib/GSLHaskell.hs deleted file mode 100644 index 254a957..0000000 --- a/lib/GSLHaskell.hs +++ /dev/null | |||
@@ -1,327 +0,0 @@ | |||
1 | {-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} | ||
2 | ----------------------------------------------------------------------------- | ||
3 | {- | | ||
4 | Module : GSLHaskell | ||
5 | Copyright : (c) Alberto Ruiz 2006 | ||
6 | License : GPL-style | ||
7 | |||
8 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
9 | Stability : provisional | ||
10 | Portability : uses -fffi and -fglasgow-exts | ||
11 | |||
12 | Old GSLHaskell interface. | ||
13 | |||
14 | -} | ||
15 | ----------------------------------------------------------------------------- | ||
16 | |||
17 | module GSLHaskell( | ||
18 | module Data.Packed.Vector, | ||
19 | module Data.Packed.Matrix, | ||
20 | module LinearAlgebra.Algorithms, | ||
21 | module LinearAlgebra.Linear, | ||
22 | module LAPACK, | ||
23 | module GSL.Integration, | ||
24 | module GSL.Differentiation, | ||
25 | module GSL.Fourier, | ||
26 | module GSL.Polynomials, | ||
27 | module GSL.Minimization, | ||
28 | module GSL.Matrix, | ||
29 | module GSL.Special, | ||
30 | module Graphics.Plot, | ||
31 | module Complex, | ||
32 | Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->), | ||
33 | fromArray2D, GSLHaskell.pnorm, | ||
34 | ) where | ||
35 | |||
36 | import GSL.Integration | ||
37 | import GSL.Differentiation | ||
38 | import GSL.Fourier | ||
39 | import GSL.Polynomials | ||
40 | import GSL.Minimization | ||
41 | import Graphics.Plot | ||
42 | import Complex | ||
43 | import GSL.Special(setErrorHandlerOff, | ||
44 | erf, | ||
45 | erf_Z, | ||
46 | bessel_J0_e, | ||
47 | exp_e10_e, | ||
48 | gamma) | ||
49 | import Data.Packed.Vector | ||
50 | import Data.Packed.Matrix | ||
51 | import Data.Packed.Matrix hiding ((><)) | ||
52 | import GSL.Vector | ||
53 | import qualified LinearAlgebra.Algorithms | ||
54 | import LAPACK | ||
55 | import GSL.Matrix | ||
56 | import LinearAlgebra.Algorithms hiding (pnorm) | ||
57 | import LinearAlgebra.Linear hiding (Mul,(<>)) | ||
58 | import Data.Packed.Internal.Matrix(multiply) | ||
59 | import Complex | ||
60 | import Numeric(showGFloat) | ||
61 | import Data.List(transpose,intersperse) | ||
62 | import Foreign(Storable) | ||
63 | import Data.Array | ||
64 | import LinearAlgebra.Instances | ||
65 | |||
66 | |||
67 | class Mul a b c | a b -> c where | ||
68 | infixl 7 <> | ||
69 | {- | 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. | ||
70 | |||
71 | @v = 'fromList' [1,2,3] :: Vector Double | ||
72 | cv = 'fromList' [1+'i',2] | ||
73 | m = 'fromLists' [[1,2,3], | ||
74 | [4,5,7]] :: Matrix Double | ||
75 | cm = 'fromLists' [[ 1, 2], | ||
76 | [3+'i',7*'i'], | ||
77 | [ 'i', 1]] | ||
78 | \ | ||
79 | \> m \<\> v | ||
80 | 14. 35. | ||
81 | \ | ||
82 | \> cv \<\> m | ||
83 | 9.+1.i 12.+2.i 17.+3.i | ||
84 | \ | ||
85 | \> m \<\> cm | ||
86 | 7.+5.i 5.+14.i | ||
87 | 19.+12.i 15.+35.i | ||
88 | \ | ||
89 | \> v \<\> 'i' | ||
90 | 1.i 2.i 3.i | ||
91 | \ | ||
92 | \> v \<\> v | ||
93 | 14.0 | ||
94 | \ | ||
95 | \> cv \<\> cv | ||
96 | 4.0 :+ 2.0@ | ||
97 | |||
98 | -} | ||
99 | (<>) :: a -> b -> c | ||
100 | |||
101 | |||
102 | instance Mul Double Double Double where | ||
103 | (<>) = (*) | ||
104 | |||
105 | instance Mul Double (Complex Double) (Complex Double) where | ||
106 | a <> b = (a:+0) * b | ||
107 | |||
108 | instance Mul (Complex Double) Double (Complex Double) where | ||
109 | a <> b = a * (b:+0) | ||
110 | |||
111 | instance Mul (Complex Double) (Complex Double) (Complex Double) where | ||
112 | (<>) = (*) | ||
113 | |||
114 | --------------------------------- matrix matrix | ||
115 | |||
116 | instance Mul (Matrix Double) (Matrix Double) (Matrix Double) where | ||
117 | (<>) = multiply | ||
118 | |||
119 | instance Mul (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where | ||
120 | (<>) = multiply | ||
121 | |||
122 | instance Mul (Matrix (Complex Double)) (Matrix Double) (Matrix (Complex Double)) where | ||
123 | c <> r = c <> liftMatrix comp r | ||
124 | |||
125 | instance Mul (Matrix Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where | ||
126 | r <> c = liftMatrix comp r <> c | ||
127 | |||
128 | --------------------------------- (Matrix Double) (Vector Double) | ||
129 | |||
130 | instance Mul (Matrix Double) (Vector Double) (Vector Double) where | ||
131 | (<>) = mXv | ||
132 | |||
133 | instance Mul (Matrix (Complex Double)) (Vector (Complex Double)) (Vector (Complex Double)) where | ||
134 | (<>) = mXv | ||
135 | |||
136 | instance Mul (Matrix (Complex Double)) (Vector Double) (Vector (Complex Double)) where | ||
137 | m <> v = m <> comp v | ||
138 | |||
139 | instance Mul (Matrix Double) (Vector (Complex Double)) (Vector (Complex Double)) where | ||
140 | m <> v = liftMatrix comp m <> v | ||
141 | |||
142 | --------------------------------- (Vector Double) (Matrix Double) | ||
143 | |||
144 | instance Mul (Vector Double) (Matrix Double) (Vector Double) where | ||
145 | (<>) = vXm | ||
146 | |||
147 | instance Mul (Vector (Complex Double)) (Matrix (Complex Double)) (Vector (Complex Double)) where | ||
148 | (<>) = vXm | ||
149 | |||
150 | instance Mul (Vector (Complex Double)) (Matrix Double) (Vector (Complex Double)) where | ||
151 | v <> m = v <> liftMatrix comp m | ||
152 | |||
153 | instance Mul (Vector Double) (Matrix (Complex Double)) (Vector (Complex Double)) where | ||
154 | v <> m = comp v <> m | ||
155 | |||
156 | --------------------------------- dot product | ||
157 | |||
158 | instance Mul (Vector Double) (Vector Double) Double where | ||
159 | (<>) = dot | ||
160 | |||
161 | instance Mul (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where | ||
162 | (<>) = dot | ||
163 | |||
164 | instance Mul (Vector Double) (Vector (Complex Double)) (Complex Double) where | ||
165 | a <> b = comp a <> b | ||
166 | |||
167 | instance Mul (Vector (Complex Double)) (Vector Double) (Complex Double) where | ||
168 | (<>) = flip (<>) | ||
169 | |||
170 | --------------------------------- scaling vectors | ||
171 | |||
172 | instance Mul Double (Vector Double) (Vector Double) where | ||
173 | (<>) = scale | ||
174 | |||
175 | instance Mul (Vector Double) Double (Vector Double) where | ||
176 | (<>) = flip (<>) | ||
177 | |||
178 | instance Mul (Complex Double) (Vector (Complex Double)) (Vector (Complex Double)) where | ||
179 | (<>) = scale | ||
180 | |||
181 | instance Mul (Vector (Complex Double)) (Complex Double) (Vector (Complex Double)) where | ||
182 | (<>) = flip (<>) | ||
183 | |||
184 | instance Mul Double (Vector (Complex Double)) (Vector (Complex Double)) where | ||
185 | a <> v = (a:+0) <> v | ||
186 | |||
187 | instance Mul (Vector (Complex Double)) Double (Vector (Complex Double)) where | ||
188 | (<>) = flip (<>) | ||
189 | |||
190 | instance Mul (Complex Double) (Vector Double) (Vector (Complex Double)) where | ||
191 | a <> v = a <> comp v | ||
192 | |||
193 | instance Mul (Vector Double) (Complex Double) (Vector (Complex Double)) where | ||
194 | (<>) = flip (<>) | ||
195 | |||
196 | --------------------------------- scaling matrices | ||
197 | |||
198 | instance Mul Double (Matrix Double) (Matrix Double) where | ||
199 | (<>) a = liftMatrix (a <>) | ||
200 | |||
201 | instance Mul (Matrix Double) Double (Matrix Double) where | ||
202 | (<>) = flip (<>) | ||
203 | |||
204 | instance Mul (Complex Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where | ||
205 | (<>) a = liftMatrix (a <>) | ||
206 | |||
207 | instance Mul (Matrix (Complex Double)) (Complex Double) (Matrix (Complex Double)) where | ||
208 | (<>) = flip (<>) | ||
209 | |||
210 | instance Mul Double (Matrix (Complex Double)) (Matrix (Complex Double)) where | ||
211 | a <> m = (a:+0) <> m | ||
212 | |||
213 | instance Mul (Matrix (Complex Double)) Double (Matrix (Complex Double)) where | ||
214 | (<>) = flip (<>) | ||
215 | |||
216 | instance Mul (Complex Double) (Matrix Double) (Matrix (Complex Double)) where | ||
217 | a <> m = a <> liftMatrix comp m | ||
218 | |||
219 | instance Mul (Matrix Double) (Complex Double) (Matrix (Complex Double)) where | ||
220 | (<>) = flip (<>) | ||
221 | |||
222 | ----------------------------------------------------------------------------------- | ||
223 | |||
224 | size :: Vector a -> Int | ||
225 | size = dim | ||
226 | |||
227 | gmap :: (Storable a, Storable b) => (a->b) -> Vector a -> Vector b | ||
228 | gmap f v = liftVector f v | ||
229 | |||
230 | -- shows a Double with n digits after the decimal point | ||
231 | shf :: (RealFloat a) => Int -> a -> String | ||
232 | shf dec n | abs n < 1e-10 = "0." | ||
233 | | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"." | ||
234 | | otherwise = showGFloat (Just dec) n "" | ||
235 | -- shows a Complex Double as a pair, with n digits after the decimal point | ||
236 | shfc n z@ (a:+b) | ||
237 | | magnitude z <1e-10 = "0." | ||
238 | | abs b < 1e-10 = shf n a | ||
239 | | abs a < 1e-10 = shf n b ++"i" | ||
240 | | b > 0 = shf n a ++"+"++shf n b ++"i" | ||
241 | | otherwise = shf n a ++shf n b ++"i" | ||
242 | |||
243 | dsp :: String -> [[String]] -> String | ||
244 | dsp sep as = unlines . map unwords' $ transpose mtp where | ||
245 | mt = transpose as | ||
246 | longs = map (maximum . map length) mt | ||
247 | mtp = zipWith (\a b -> map (pad a) b) longs mt | ||
248 | pad n str = replicate (n - length str) ' ' ++ str | ||
249 | unwords' = concat . intersperse sep | ||
250 | |||
251 | format :: (Field t) => String -> (t -> String) -> Matrix t -> String | ||
252 | format sep f m = dsp sep . map (map f) . toLists $ m | ||
253 | |||
254 | disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m | ||
255 | |||
256 | dispR :: Int -> Matrix Double -> IO () | ||
257 | dispR d m = disp m (shf d) | ||
258 | |||
259 | dispC :: Int -> Matrix (Complex Double) -> IO () | ||
260 | dispC d m = disp m (shfc d) | ||
261 | |||
262 | -- | creates a matrix from a table of numbers. | ||
263 | readMatrix :: String -> Matrix Double | ||
264 | readMatrix = fromLists . map (map read). map words . filter (not.null) . lines | ||
265 | |||
266 | ------------------------------------------------------------- | ||
267 | |||
268 | class Joinable a b c | a b -> c where | ||
269 | joinH :: a -> b -> c | ||
270 | joinV :: a -> b -> c | ||
271 | |||
272 | instance Joinable (Matrix Double) (Vector Double) (Matrix Double) where | ||
273 | joinH m v = fromBlocks [[m,reshape 1 v]] | ||
274 | joinV m v = fromBlocks [[m],[reshape (size v) v]] | ||
275 | |||
276 | instance Joinable (Vector Double) (Matrix Double) (Matrix Double) where | ||
277 | joinH v m = fromBlocks [[reshape 1 v,m]] | ||
278 | joinV v m = fromBlocks [[reshape (size v) v],[m]] | ||
279 | |||
280 | instance Joinable (Matrix Double) (Matrix Double) (Matrix Double) where | ||
281 | joinH m1 m2 = fromBlocks [[m1,m2]] | ||
282 | joinV m1 m2 = fromBlocks [[m1],[m2]] | ||
283 | |||
284 | instance Joinable (Matrix (Complex Double)) (Vector (Complex Double)) (Matrix (Complex Double)) where | ||
285 | joinH m v = fromBlocks [[m,reshape 1 v]] | ||
286 | joinV m v = fromBlocks [[m],[reshape (size v) v]] | ||
287 | |||
288 | instance Joinable (Vector (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where | ||
289 | joinH v m = fromBlocks [[reshape 1 v,m]] | ||
290 | joinV v m = fromBlocks [[reshape (size v) v],[m]] | ||
291 | |||
292 | instance Joinable (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where | ||
293 | joinH m1 m2 = fromBlocks [[m1,m2]] | ||
294 | joinV m1 m2 = fromBlocks [[m1],[m2]] | ||
295 | |||
296 | infixl 3 <|>, <-> | ||
297 | |||
298 | {- | Horizontal concatenation of matrices and vectors: | ||
299 | |||
300 | @\> 'ident' 3 \<-\> i\<\>'ident' 3 \<|\> 'fromList' [1..6] | ||
301 | 1. 0. 0. 1. | ||
302 | 0. 1. 0. 2. | ||
303 | 0. 0. 1. 3. | ||
304 | 1.i 0. 0. 4. | ||
305 | 0. 1.i 0. 5. | ||
306 | 0. 0. 1.i 6.@ | ||
307 | -} | ||
308 | (<|>) :: (Joinable a b c) => a -> b -> c | ||
309 | a <|> b = joinH a b | ||
310 | |||
311 | -- | Vertical concatenation of matrices and vectors. | ||
312 | (<->) :: (Joinable a b c) => a -> b -> c | ||
313 | a <-> b = joinV a b | ||
314 | |||
315 | ---------------------------------------------------------- | ||
316 | |||
317 | fromArray2D :: (Field e) => Array (Int, Int) e -> Matrix e | ||
318 | fromArray2D m = (r><c) (elems m) | ||
319 | where ((r0,c0),(r1,c1)) = bounds m | ||
320 | r = r1-r0+1 | ||
321 | c = c1-c0+1 | ||
322 | |||
323 | |||
324 | pnorm :: (Normed t1, Num t) => t -> t1 -> Double | ||
325 | pnorm 0 = LinearAlgebra.Algorithms.pnorm Infinity | ||
326 | pnorm 1 = LinearAlgebra.Algorithms.pnorm PNorm1 | ||
327 | pnorm 2 = LinearAlgebra.Algorithms.pnorm PNorm2 \ No newline at end of file | ||
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 794ef69..b7e208a 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -51,7 +51,7 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
51 | -- * Util | 51 | -- * Util |
52 | haussholder, | 52 | haussholder, |
53 | unpackQR, unpackHess, | 53 | unpackQR, unpackHess, |
54 | GenMat(linearSolveSVD,lu,eigSH',cholSH) | 54 | Field(linearSolveSVD,lu,eigSH',cholSH) |
55 | ) where | 55 | ) where |
56 | 56 | ||
57 | 57 | ||
@@ -65,7 +65,7 @@ import Numeric.LinearAlgebra.Linear | |||
65 | import Data.List(foldl1') | 65 | import Data.List(foldl1') |
66 | 66 | ||
67 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. | 67 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. |
68 | class (Normed (Matrix t), Linear Matrix t) => GenMat t where | 68 | class (Normed (Matrix t), Linear Matrix t) => Field t where |
69 | -- | Singular value decomposition using lapack's dgesvd or zgesvd. | 69 | -- | Singular value decomposition using lapack's dgesvd or zgesvd. |
70 | svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 70 | svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
71 | lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) | 71 | lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) |
@@ -103,7 +103,7 @@ class (Normed (Matrix t), Linear Matrix t) => GenMat t where | |||
103 | ctrans :: Matrix t -> Matrix t | 103 | ctrans :: Matrix t -> Matrix t |
104 | 104 | ||
105 | 105 | ||
106 | instance GenMat Double where | 106 | instance Field Double where |
107 | svd = svdR | 107 | svd = svdR |
108 | lu = GSL.luR | 108 | lu = GSL.luR |
109 | linearSolve = linearSolveR | 109 | linearSolve = linearSolveR |
@@ -116,7 +116,7 @@ instance GenMat Double where | |||
116 | hess = unpackHess hessR | 116 | hess = unpackHess hessR |
117 | schur = schurR | 117 | schur = schurR |
118 | 118 | ||
119 | instance GenMat (Complex Double) where | 119 | instance Field (Complex Double) where |
120 | svd = svdC | 120 | svd = svdC |
121 | lu = GSL.luC | 121 | lu = GSL.luC |
122 | linearSolve = linearSolveC | 122 | linearSolve = linearSolveC |
@@ -132,37 +132,37 @@ instance GenMat (Complex Double) where | |||
132 | -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. | 132 | -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. |
133 | -- | 133 | -- |
134 | -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ | 134 | -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ |
135 | eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) | 135 | eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) |
136 | eigSH m | m `equal` ctrans m = eigSH' m | 136 | eigSH m | m `equal` ctrans m = eigSH' m |
137 | | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" | 137 | | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" |
138 | 138 | ||
139 | -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. | 139 | -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. |
140 | -- | 140 | -- |
141 | -- If @c = chol m@ then @m == c \<> ctrans c@. | 141 | -- If @c = chol m@ then @m == c \<> ctrans c@. |
142 | chol :: GenMat t => Matrix t -> Matrix t | 142 | chol :: Field t => Matrix t -> Matrix t |
143 | chol m | m `equal` ctrans m = cholSH m | 143 | chol m | m `equal` ctrans m = cholSH m |
144 | | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" | 144 | | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" |
145 | 145 | ||
146 | square m = rows m == cols m | 146 | square m = rows m == cols m |
147 | 147 | ||
148 | det :: GenMat t => Matrix t -> t | 148 | det :: Field t => Matrix t -> t |
149 | det m | square m = s * (product $ toList $ takeDiag $ u) | 149 | det m | square m = s * (product $ toList $ takeDiag $ u) |
150 | | otherwise = error "det of nonsquare matrix" | 150 | | otherwise = error "det of nonsquare matrix" |
151 | where (_,u,_,s) = lu m | 151 | where (_,u,_,s) = lu m |
152 | 152 | ||
153 | -- | Inverse of a square matrix using lapacks' dgesv and zgesv. | 153 | -- | Inverse of a square matrix using lapacks' dgesv and zgesv. |
154 | inv :: GenMat t => Matrix t -> Matrix t | 154 | inv :: Field t => Matrix t -> Matrix t |
155 | inv m | square m = m `linearSolve` ident (rows m) | 155 | inv m | square m = m `linearSolve` ident (rows m) |
156 | | otherwise = error "inv of nonsquare matrix" | 156 | | otherwise = error "inv of nonsquare matrix" |
157 | 157 | ||
158 | -- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss. | 158 | -- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss. |
159 | pinv :: GenMat t => Matrix t -> Matrix t | 159 | pinv :: Field t => Matrix t -> Matrix t |
160 | pinv m = linearSolveSVD m (ident (rows m)) | 160 | pinv m = linearSolveSVD m (ident (rows m)) |
161 | 161 | ||
162 | -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. | 162 | -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. |
163 | -- | 163 | -- |
164 | -- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. | 164 | -- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. |
165 | full :: Field t | 165 | full :: Element t |
166 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) | 166 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) |
167 | full svd m = (u, d ,v) where | 167 | full svd m = (u, d ,v) where |
168 | (u,s,v) = svd m | 168 | (u,s,v) = svd m |
@@ -173,7 +173,7 @@ full svd m = (u, d ,v) where | |||
173 | -- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. | 173 | -- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. |
174 | -- | 174 | -- |
175 | -- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. | 175 | -- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. |
176 | economy :: Field t | 176 | economy :: Element t |
177 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) | 177 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) |
178 | economy svd m = (u', subVector 0 d s, v') where | 178 | economy svd m = (u', subVector 0 d s, v') where |
179 | (u,s,v) = svd m | 179 | (u,s,v) = svd m |
@@ -198,15 +198,15 @@ i = 0:+1 | |||
198 | 198 | ||
199 | 199 | ||
200 | -- matrix product | 200 | -- matrix product |
201 | mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t | 201 | mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t |
202 | mXm = multiply | 202 | mXm = multiply |
203 | 203 | ||
204 | -- matrix - vector product | 204 | -- matrix - vector product |
205 | mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t | 205 | mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t |
206 | mXv m v = flatten $ m `mXm` (asColumn v) | 206 | mXv m v = flatten $ m `mXm` (asColumn v) |
207 | 207 | ||
208 | -- vector - matrix product | 208 | -- vector - matrix product |
209 | vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t | 209 | vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t |
210 | vXm v m = flatten $ (asRow v) `mXm` m | 210 | vXm v m = flatten $ (asRow v) `mXm` m |
211 | 211 | ||
212 | 212 | ||
@@ -264,7 +264,7 @@ instance Normed (Matrix (Complex Double)) where | |||
264 | ----------------------------------------------------------------------- | 264 | ----------------------------------------------------------------------- |
265 | 265 | ||
266 | -- | The nullspace of a matrix from its SVD decomposition. | 266 | -- | The nullspace of a matrix from its SVD decomposition. |
267 | nullspacePrec :: GenMat t | 267 | nullspacePrec :: Field t |
268 | => Double -- ^ relative tolerance in 'eps' units | 268 | => Double -- ^ relative tolerance in 'eps' units |
269 | -> Matrix t -- ^ input matrix | 269 | -> Matrix t -- ^ input matrix |
270 | -> [Vector t] -- ^ list of unitary vectors spanning the nullspace | 270 | -> [Vector t] -- ^ list of unitary vectors spanning the nullspace |
@@ -276,7 +276,7 @@ nullspacePrec t m = ns where | |||
276 | ns = drop rank $ toRows $ ctrans v | 276 | ns = drop rank $ toRows $ ctrans v |
277 | 277 | ||
278 | -- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@). | 278 | -- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@). |
279 | nullVector :: GenMat t => Matrix t -> Vector t | 279 | nullVector :: Field t => Matrix t -> Vector t |
280 | nullVector = last . nullspacePrec 1 | 280 | nullVector = last . nullspacePrec 1 |
281 | 281 | ||
282 | ------------------------------------------------------------------------ | 282 | ------------------------------------------------------------------------ |
@@ -316,7 +316,7 @@ pinvTol t m = v' `mXm` diag s' `mXm` trans u' where | |||
316 | 316 | ||
317 | -- many thanks, quickcheck! | 317 | -- many thanks, quickcheck! |
318 | 318 | ||
319 | haussholder :: (GenMat a) => a -> Vector a -> Matrix a | 319 | haussholder :: (Field a) => a -> Vector a -> Matrix a |
320 | haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) | 320 | haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) |
321 | where w = asColumn v | 321 | where w = asColumn v |
322 | 322 | ||
@@ -328,7 +328,7 @@ zt 0 v = v | |||
328 | zt k v = join [subVector 0 (dim v - k) v, constant 0 k] | 328 | zt k v = join [subVector 0 (dim v - k) v, constant 0 k] |
329 | 329 | ||
330 | 330 | ||
331 | unpackQR :: (GenMat t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) | 331 | unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) |
332 | unpackQR (pq, tau) = (q,r) | 332 | unpackQR (pq, tau) = (q,r) |
333 | where cs = toColumns pq | 333 | where cs = toColumns pq |
334 | m = rows pq | 334 | m = rows pq |
@@ -339,7 +339,7 @@ unpackQR (pq, tau) = (q,r) | |||
339 | hs = zipWith haussholder (toList tau) vs | 339 | hs = zipWith haussholder (toList tau) vs |
340 | q = foldl1' mXm hs | 340 | q = foldl1' mXm hs |
341 | 341 | ||
342 | unpackHess :: (GenMat t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) | 342 | unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) |
343 | unpackHess hf m | 343 | unpackHess hf m |
344 | | rows m == 1 = ((1><1)[1],m) | 344 | | rows m == 1 = ((1><1)[1],m) |
345 | | otherwise = (uH . hf) m | 345 | | otherwise = (uH . hf) m |
@@ -357,13 +357,13 @@ uH (pq, tau) = (p,h) | |||
357 | -------------------------------------------------------------------------- | 357 | -------------------------------------------------------------------------- |
358 | 358 | ||
359 | -- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. | 359 | -- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. |
360 | rcond :: GenMat t => Matrix t -> Double | 360 | rcond :: Field t => Matrix t -> Double |
361 | rcond m = last s / head s | 361 | rcond m = last s / head s |
362 | where (_,s',_) = svd m | 362 | where (_,s',_) = svd m |
363 | s = toList s' | 363 | s = toList s' |
364 | 364 | ||
365 | -- | Number of linearly independent rows or columns. | 365 | -- | Number of linearly independent rows or columns. |
366 | rank :: GenMat t => Matrix t -> Int | 366 | rank :: Field t => Matrix t -> Int |
367 | rank m | pnorm PNorm1 m < eps = 0 | 367 | rank m | pnorm PNorm1 m < eps = 0 |
368 | | otherwise = dim s where (_,s,_) = economy svd m | 368 | | otherwise = dim s where (_,s,_) = economy svd m |
369 | 369 | ||
@@ -381,7 +381,7 @@ diagonalize m = if rank v == n | |||
381 | else eig m | 381 | else eig m |
382 | 382 | ||
383 | -- | Generic matrix functions for diagonalizable matrices. | 383 | -- | Generic matrix functions for diagonalizable matrices. |
384 | matFunc :: GenMat t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) | 384 | matFunc :: Field t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) |
385 | matFunc f m = case diagonalize (complex m) of | 385 | matFunc f m = case diagonalize (complex m) of |
386 | Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v | 386 | Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v |
387 | Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix" | 387 | Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix" |
@@ -420,5 +420,5 @@ expGolub m = iterate msq f !! j | |||
420 | {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, | 420 | {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, |
421 | based on a scaled Pade approximation. | 421 | based on a scaled Pade approximation. |
422 | -} | 422 | -} |
423 | expm :: GenMat t => Matrix t -> Matrix t | 423 | expm :: Field t => Matrix t -> Matrix t |
424 | expm = expGolub | 424 | expm = expGolub |
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs index 4ee576f..3f295bf 100644 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Instances.hs | |||
@@ -29,7 +29,7 @@ import Foreign(Storable) | |||
29 | 29 | ||
30 | ------------------------------------------------------------------ | 30 | ------------------------------------------------------------------ |
31 | 31 | ||
32 | instance (Show a, Field a) => (Show (Matrix a)) where | 32 | instance (Show a, Element a) => (Show (Matrix a)) where |
33 | show m = (sizes++) . dsp . map (map show) . toLists $ m | 33 | show m = (sizes++) . dsp . map (map show) . toLists $ m |
34 | where sizes = "("++show (rows m)++"><"++show (cols m)++")\n" | 34 | where sizes = "("++show (rows m)++"><"++show (cols m)++")\n" |
35 | 35 | ||
@@ -51,7 +51,7 @@ adaptScalar f1 f2 f3 x y | |||
51 | | dim y == 1 = f3 x (y@>0) | 51 | | dim y == 1 = f3 x (y@>0) |
52 | | otherwise = f2 x y | 52 | | otherwise = f2 x y |
53 | 53 | ||
54 | liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t | 54 | liftMatrix2' :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t |
55 | liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2)) | 55 | liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2)) |
56 | | otherwise = error "nonconformant matrices in liftMatrix2'" | 56 | | otherwise = error "nonconformant matrices in liftMatrix2'" |
57 | 57 | ||
@@ -60,7 +60,7 @@ compat' m1 m2 = rows m1 == 1 && cols m1 == 1 | |||
60 | || rows m2 == 1 && cols m2 == 1 | 60 | || rows m2 == 1 && cols m2 == 1 |
61 | || rows m1 == rows m2 && cols m1 == cols m2 | 61 | || rows m1 == rows m2 && cols m1 == cols m2 |
62 | 62 | ||
63 | instance (Eq a, Field a) => Eq (Vector a) where | 63 | instance (Eq a, Element a) => Eq (Vector a) where |
64 | a == b = dim a == dim b && toList a == toList b | 64 | a == b = dim a == dim b && toList a == toList b |
65 | 65 | ||
66 | instance (Linear Vector a) => Num (Vector a) where | 66 | instance (Linear Vector a) => Num (Vector a) where |
@@ -71,7 +71,7 @@ instance (Linear Vector a) => Num (Vector a) where | |||
71 | abs = liftVector abs | 71 | abs = liftVector abs |
72 | fromInteger = fromList . return . fromInteger | 72 | fromInteger = fromList . return . fromInteger |
73 | 73 | ||
74 | instance (Eq a, Field a) => Eq (Matrix a) where | 74 | instance (Eq a, Element a) => Eq (Matrix a) where |
75 | a == b = cols a == cols b && flatten a == flatten b | 75 | a == b = cols a == cols b && flatten a == flatten b |
76 | 76 | ||
77 | instance (Linear Vector a) => Num (Matrix a) where | 77 | instance (Linear Vector a) => Num (Matrix a) where |
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs index fd076ec..4a9b309 100644 --- a/lib/Numeric/LinearAlgebra/Interface.hs +++ b/lib/Numeric/LinearAlgebra/Interface.hs | |||
@@ -29,7 +29,7 @@ import Numeric.LinearAlgebra.Algorithms | |||
29 | class Mul a b c | a b -> c where | 29 | class Mul a b c | a b -> c where |
30 | infixl 7 <> | 30 | infixl 7 <> |
31 | -- | matrix product | 31 | -- | matrix product |
32 | (<>) :: Field t => a t -> b t -> c t | 32 | (<>) :: Element t => a t -> b t -> c t |
33 | 33 | ||
34 | instance Mul Matrix Matrix Matrix where | 34 | instance Mul Matrix Matrix Matrix where |
35 | (<>) = multiply | 35 | (<>) = multiply |
@@ -43,7 +43,7 @@ instance Mul Vector Matrix Vector where | |||
43 | --------------------------------------------------- | 43 | --------------------------------------------------- |
44 | 44 | ||
45 | -- | @u \<.\> v = dot u v@ | 45 | -- | @u \<.\> v = dot u v@ |
46 | (<.>) :: (Field t) => Vector t -> Vector t -> t | 46 | (<.>) :: (Element t) => Vector t -> Vector t -> t |
47 | infixl 7 <.> | 47 | infixl 7 <.> |
48 | (<.>) = dot | 48 | (<.>) = dot |
49 | 49 | ||
@@ -62,15 +62,15 @@ infixl 7 */ | |||
62 | v */ x = scale (recip x) v | 62 | v */ x = scale (recip x) v |
63 | 63 | ||
64 | -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). | 64 | -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). |
65 | (<\>) :: (GenMat a) => Matrix a -> Vector a -> Vector a | 65 | (<\>) :: (Field a) => Matrix a -> Vector a -> Vector a |
66 | infixl 7 <\> | 66 | infixl 7 <\> |
67 | m <\> v = flatten (linearSolveSVD m (reshape 1 v)) | 67 | m <\> v = flatten (linearSolveSVD m (reshape 1 v)) |
68 | 68 | ||
69 | ------------------------------------------------ | 69 | ------------------------------------------------ |
70 | 70 | ||
71 | class Joinable a b where | 71 | class Joinable a b where |
72 | joinH :: Field t => a t -> b t -> Matrix t | 72 | joinH :: Element t => a t -> b t -> Matrix t |
73 | joinV :: Field t => a t -> b t -> Matrix t | 73 | joinV :: Element t => a t -> b t -> Matrix t |
74 | 74 | ||
75 | instance Joinable Matrix Matrix where | 75 | instance Joinable Matrix Matrix where |
76 | joinH m1 m2 = fromBlocks [[m1,m2]] | 76 | joinH m1 m2 = fromBlocks [[m1,m2]] |
@@ -98,9 +98,9 @@ infixl 3 <-> | |||
98 | , 0.0, 3.0, 0.0, 5.0 | 98 | , 0.0, 3.0, 0.0, 5.0 |
99 | , 0.0, 0.0, 3.0, 6.0 ]@ | 99 | , 0.0, 0.0, 3.0, 6.0 ]@ |
100 | -} | 100 | -} |
101 | (<|>) :: (Field t, Joinable a b) => a t -> b t -> Matrix t | 101 | (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t |
102 | a <|> b = joinH a b | 102 | a <|> b = joinH a b |
103 | 103 | ||
104 | -- | Vertical concatenation of matrices and vectors. | 104 | -- | Vertical concatenation of matrices and vectors. |
105 | (<->) :: (Field t, Joinable a b) => a t -> b t -> Matrix t | 105 | (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t |
106 | a <-> b = joinV a b | 106 | a <-> b = joinV a b |
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs index 94f6958..13d69ab 100644 --- a/lib/Numeric/LinearAlgebra/Linear.hs +++ b/lib/Numeric/LinearAlgebra/Linear.hs | |||
@@ -84,7 +84,7 @@ instance Linear Matrix (Complex Double) where | |||
84 | -------------------------------------------------- | 84 | -------------------------------------------------- |
85 | 85 | ||
86 | -- | euclidean inner product | 86 | -- | euclidean inner product |
87 | dot :: (Field t) => Vector t -> Vector t -> t | 87 | dot :: (Element t) => Vector t -> Vector t -> t |
88 | dot u v = dat (multiply r c) `at` 0 | 88 | dot u v = dat (multiply r c) `at` 0 |
89 | where r = asRow u | 89 | where r = asRow u |
90 | c = asColumn v | 90 | c = asColumn v |
@@ -98,7 +98,7 @@ dot u v = dat (multiply r c) `at` 0 | |||
98 | , 10.0, 4.0, 6.0 | 98 | , 10.0, 4.0, 6.0 |
99 | , 15.0, 6.0, 9.0 ]@ | 99 | , 15.0, 6.0, 9.0 ]@ |
100 | -} | 100 | -} |
101 | outer :: (Field t) => Vector t -> Vector t -> Matrix t | 101 | outer :: (Element t) => Vector t -> Vector t -> Matrix t |
102 | outer u v = asColumn u `multiply` asRow v | 102 | outer u v = asColumn u `multiply` asRow v |
103 | 103 | ||
104 | {- | Kronecker product of two matrices. | 104 | {- | Kronecker product of two matrices. |
@@ -123,7 +123,7 @@ m2=(4><3) | |||
123 | , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 | 123 | , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 |
124 | , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ | 124 | , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ |
125 | -} | 125 | -} |
126 | kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t | 126 | kronecker :: (Element t) => Matrix t -> Matrix t -> Matrix t |
127 | kronecker a b = fromBlocks | 127 | kronecker a b = fromBlocks |
128 | . partit (cols a) | 128 | . partit (cols a) |
129 | . map (reshape (cols b)) | 129 | . map (reshape (cols b)) |