diff options
Diffstat (limited to 'packages/base/src/Numeric/LinearAlgebra')
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Data.hs | 11 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Real.hs | 337 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util.hs | 20 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util/CG.hs | 86 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util/Static.hs | 70 |
5 files changed, 419 insertions, 105 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index 3128a24..3417a5e 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs | |||
@@ -49,13 +49,9 @@ module Numeric.LinearAlgebra.Data( | |||
49 | find, maxIndex, minIndex, maxElement, minElement, atIndex, | 49 | find, maxIndex, minIndex, maxElement, minElement, atIndex, |
50 | 50 | ||
51 | -- * Sparse | 51 | -- * Sparse |
52 | SMatrix, AssocMatrix, mkCSR, toDense, | 52 | GMatrix, AssocMatrix, mkSparse, toDense, |
53 | mkDiag, | 53 | mkDiagR, dense, |
54 | |||
55 | -- * Static dimensions | ||
56 | 54 | ||
57 | Static, ddata, R, vect0, sScalar, vect2, vect3, (&), | ||
58 | |||
59 | -- * IO | 55 | -- * IO |
60 | disp, | 56 | disp, |
61 | loadMatrix, saveMatrix, | 57 | loadMatrix, saveMatrix, |
@@ -79,9 +75,8 @@ module Numeric.LinearAlgebra.Data( | |||
79 | import Data.Packed.Vector | 75 | import Data.Packed.Vector |
80 | import Data.Packed.Matrix | 76 | import Data.Packed.Matrix |
81 | import Data.Packed.Numeric | 77 | import Data.Packed.Numeric |
82 | import Numeric.LinearAlgebra.Util hiding ((&)) | 78 | import Numeric.LinearAlgebra.Util hiding ((&),(#)) |
83 | import Data.Complex | 79 | import Data.Complex |
84 | import Numeric.Sparse | 80 | import Numeric.Sparse |
85 | import Numeric.LinearAlgebra.Util.Static | ||
86 | 81 | ||
87 | 82 | ||
diff --git a/packages/base/src/Numeric/LinearAlgebra/Real.hs b/packages/base/src/Numeric/LinearAlgebra/Real.hs new file mode 100644 index 0000000..db15705 --- /dev/null +++ b/packages/base/src/Numeric/LinearAlgebra/Real.hs | |||
@@ -0,0 +1,337 @@ | |||
1 | {-# LANGUAGE DataKinds #-} | ||
2 | {-# LANGUAGE KindSignatures #-} | ||
3 | {-# LANGUAGE GeneralizedNewtypeDeriving #-} | ||
4 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
5 | {-# LANGUAGE FunctionalDependencies #-} | ||
6 | {-# LANGUAGE FlexibleContexts #-} | ||
7 | {-# LANGUAGE ScopedTypeVariables #-} | ||
8 | {-# LANGUAGE EmptyDataDecls #-} | ||
9 | {-# LANGUAGE Rank2Types #-} | ||
10 | {-# LANGUAGE FlexibleInstances #-} | ||
11 | {-# LANGUAGE TypeOperators #-} | ||
12 | {-# LANGUAGE ViewPatterns #-} | ||
13 | {-# LANGUAGE GADTs #-} | ||
14 | |||
15 | |||
16 | {- | | ||
17 | Module : Numeric.LinearAlgebra.Real | ||
18 | Copyright : (c) Alberto Ruiz 2006-14 | ||
19 | License : BSD3 | ||
20 | Stability : provisional | ||
21 | |||
22 | Experimental interface for real arrays with statically checked dimensions. | ||
23 | |||
24 | -} | ||
25 | |||
26 | module Numeric.LinearAlgebra.Real( | ||
27 | -- * Vector | ||
28 | R, | ||
29 | vec2, vec3, vec4, ๐ง, (&), | ||
30 | -- * Matrix | ||
31 | L, Sq, | ||
32 | ๐, | ||
33 | (#),(ยฆ),(โโ), | ||
34 | Konst(..), | ||
35 | eye, | ||
36 | diagR, diag, | ||
37 | blockAt, | ||
38 | -- * Products | ||
39 | (<>),(#>),(<ยท>), | ||
40 | -- * Pretty printing | ||
41 | Disp(..), | ||
42 | -- * Misc | ||
43 | Dim, unDim, | ||
44 | module Numeric.HMatrix | ||
45 | ) where | ||
46 | |||
47 | |||
48 | import GHC.TypeLits | ||
49 | import Numeric.HMatrix hiding ((<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โโ)) | ||
50 | import qualified Numeric.HMatrix as LA | ||
51 | import Data.Packed.ST | ||
52 | |||
53 | newtype Dim (n :: Nat) t = Dim t | ||
54 | deriving Show | ||
55 | |||
56 | unDim :: Dim n t -> t | ||
57 | unDim (Dim x) = x | ||
58 | |||
59 | data Proxy :: Nat -> * | ||
60 | |||
61 | |||
62 | lift1F | ||
63 | :: (c t -> c t) | ||
64 | -> Dim n (c t) -> Dim n (c t) | ||
65 | lift1F f (Dim v) = Dim (f v) | ||
66 | |||
67 | lift2F | ||
68 | :: (c t -> c t -> c t) | ||
69 | -> Dim n (c t) -> Dim n (c t) -> Dim n (c t) | ||
70 | lift2F f (Dim u) (Dim v) = Dim (f u v) | ||
71 | |||
72 | |||
73 | |||
74 | type R n = Dim n (Vector โ) | ||
75 | |||
76 | type L m n = Dim m (Dim n (Matrix โ)) | ||
77 | |||
78 | |||
79 | infixl 4 & | ||
80 | (&) :: forall n . KnownNat n | ||
81 | => R n -> โ -> R (n+1) | ||
82 | Dim v & x = Dim (vjoin [v', scalar x]) | ||
83 | where | ||
84 | d = fromIntegral . natVal $ (undefined :: Proxy n) | ||
85 | v' | d > 1 && size v == 1 = LA.konst (v!0) d | ||
86 | | otherwise = v | ||
87 | |||
88 | |||
89 | -- vect0 :: R 0 | ||
90 | -- vect0 = Dim (fromList[]) | ||
91 | |||
92 | ๐ง :: โ -> R 1 | ||
93 | ๐ง = Dim . scalar | ||
94 | |||
95 | |||
96 | vec2 :: โ -> โ -> R 2 | ||
97 | vec2 a b = Dim $ runSTVector $ do | ||
98 | v <- newUndefinedVector 2 | ||
99 | writeVector v 0 a | ||
100 | writeVector v 1 b | ||
101 | return v | ||
102 | |||
103 | vec3 :: โ -> โ -> โ -> R 3 | ||
104 | vec3 a b c = Dim $ runSTVector $ do | ||
105 | v <- newUndefinedVector 3 | ||
106 | writeVector v 0 a | ||
107 | writeVector v 1 b | ||
108 | writeVector v 2 c | ||
109 | return v | ||
110 | |||
111 | |||
112 | vec4 :: โ -> โ -> โ -> โ -> R 4 | ||
113 | vec4 a b c d = Dim $ runSTVector $ do | ||
114 | v <- newUndefinedVector 4 | ||
115 | writeVector v 0 a | ||
116 | writeVector v 1 b | ||
117 | writeVector v 2 c | ||
118 | writeVector v 3 d | ||
119 | return v | ||
120 | |||
121 | |||
122 | |||
123 | |||
124 | instance forall n t . (Num (Vector t), Numeric t )=> Num (Dim n (Vector t)) | ||
125 | where | ||
126 | (+) = lift2F (+) | ||
127 | (*) = lift2F (*) | ||
128 | (-) = lift2F (-) | ||
129 | abs = lift1F abs | ||
130 | signum = lift1F signum | ||
131 | negate = lift1F negate | ||
132 | fromInteger x = Dim (fromInteger x) | ||
133 | |||
134 | instance (Num (Matrix t), Numeric t) => Num (Dim m (Dim n (Matrix t))) | ||
135 | where | ||
136 | (+) = (lift2F . lift2F) (+) | ||
137 | (*) = (lift2F . lift2F) (*) | ||
138 | (-) = (lift2F . lift2F) (-) | ||
139 | abs = (lift1F . lift1F) abs | ||
140 | signum = (lift1F . lift1F) signum | ||
141 | negate = (lift1F . lift1F) negate | ||
142 | fromInteger x = Dim (Dim (fromInteger x)) | ||
143 | |||
144 | -------------------------------------------------------------------------------- | ||
145 | |||
146 | class Konst t | ||
147 | where | ||
148 | konst :: โ -> t | ||
149 | |||
150 | instance forall n. KnownNat n => Konst (R n) | ||
151 | where | ||
152 | konst x = Dim (LA.konst x d) | ||
153 | where | ||
154 | d = fromIntegral . natVal $ (undefined :: Proxy n) | ||
155 | |||
156 | instance forall m n . (KnownNat m, KnownNat n) => Konst (L m n) | ||
157 | where | ||
158 | konst x = Dim (Dim (LA.konst x (m',n'))) | ||
159 | where | ||
160 | m' = fromIntegral . natVal $ (undefined :: Proxy m) | ||
161 | n' = fromIntegral . natVal $ (undefined :: Proxy n) | ||
162 | |||
163 | -------------------------------------------------------------------------------- | ||
164 | |||
165 | diagR :: forall m n k . (KnownNat m, KnownNat n) => โ -> R k -> L m n | ||
166 | diagR x v = Dim (Dim (diagRect x (unDim v) m' n')) | ||
167 | where | ||
168 | m' = fromIntegral . natVal $ (undefined :: Proxy m) | ||
169 | n' = fromIntegral . natVal $ (undefined :: Proxy n) | ||
170 | |||
171 | diag :: KnownNat n => R n -> Sq n | ||
172 | diag = diagR 0 | ||
173 | |||
174 | -------------------------------------------------------------------------------- | ||
175 | |||
176 | blockAt :: forall m n . (KnownNat m, KnownNat n) => โ -> Int -> Int -> Matrix Double -> L m n | ||
177 | blockAt x r c a = Dim (Dim res) | ||
178 | where | ||
179 | z = scalar x | ||
180 | z1 = LA.konst x (r,c) | ||
181 | z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c))) | ||
182 | ra = min (rows a) . max 0 $ m'-r | ||
183 | ca = min (cols a) . max 0 $ n'-c | ||
184 | sa = subMatrix (0,0) (ra, ca) a | ||
185 | m' = fromIntegral . natVal $ (undefined :: Proxy m) | ||
186 | n' = fromIntegral . natVal $ (undefined :: Proxy n) | ||
187 | res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]] | ||
188 | |||
189 | {- | ||
190 | matrix :: (KnownNat m, KnownNat n) => Matrix Double -> L n m | ||
191 | matrix = blockAt 0 0 0 | ||
192 | -} | ||
193 | |||
194 | -------------------------------------------------------------------------------- | ||
195 | |||
196 | class Disp t | ||
197 | where | ||
198 | disp :: Int -> t -> IO () | ||
199 | |||
200 | instance Disp (L n m) | ||
201 | where | ||
202 | disp n (d2 -> a) = do | ||
203 | if rows a == 1 && cols a == 1 | ||
204 | then putStrLn $ "Const " ++ (last . words . LA.dispf n $ a) | ||
205 | else putStr "Dim " >> LA.disp n a | ||
206 | |||
207 | instance Disp (R n) | ||
208 | where | ||
209 | disp n (unDim -> v) = do | ||
210 | let su = LA.dispf n (asRow v) | ||
211 | if LA.size v == 1 | ||
212 | then putStrLn $ "Const " ++ (last . words $ su ) | ||
213 | else putStr "Dim " >> putStr (tail . dropWhile (/='x') $ su) | ||
214 | |||
215 | -------------------------------------------------------------------------------- | ||
216 | |||
217 | infixl 3 # | ||
218 | (#) :: L r c -> R c -> L (r+1) c | ||
219 | Dim (Dim m) # Dim v = Dim (Dim (m LA.โโ asRow v)) | ||
220 | |||
221 | |||
222 | ๐ :: forall n . KnownNat n => L 0 n | ||
223 | ๐ = Dim (Dim (LA.konst 0 (0,d))) | ||
224 | where | ||
225 | d = fromIntegral . natVal $ (undefined :: Proxy n) | ||
226 | |||
227 | infixl 3 ยฆ | ||
228 | (ยฆ) :: L r c1 -> L r c2 -> L r (c1+c2) | ||
229 | Dim (Dim a) ยฆ Dim (Dim b) = Dim (Dim (a LA.ยฆ b)) | ||
230 | |||
231 | infixl 2 โโ | ||
232 | (โโ) :: L r1 c -> L r2 c -> L (r1+r2) c | ||
233 | Dim (Dim a) โโ Dim (Dim b) = Dim (Dim (a LA.โโ b)) | ||
234 | |||
235 | |||
236 | {- | ||
237 | |||
238 | -} | ||
239 | |||
240 | type Sq n = L n n | ||
241 | |||
242 | type GL = (KnownNat n, KnownNat m) => L m n | ||
243 | type GSq = KnownNat n => Sq n | ||
244 | |||
245 | infixr 8 <> | ||
246 | (<>) :: L m k -> L k n -> L m n | ||
247 | (d2 -> a) <> (d2 -> b) = Dim (Dim (a LA.<> b)) | ||
248 | |||
249 | infixr 8 #> | ||
250 | (#>) :: L m n -> R n -> R m | ||
251 | (d2 -> m) #> (unDim -> v) = Dim (m LA.#> v) | ||
252 | |||
253 | infixr 8 <ยท> | ||
254 | (<ยท>) :: R n -> R n -> โ | ||
255 | (unDim -> u) <ยท> (unDim -> v) = udot u v | ||
256 | |||
257 | |||
258 | d2 :: forall c (n :: Nat) (n1 :: Nat). Dim n1 (Dim n c) -> c | ||
259 | d2 = unDim . unDim | ||
260 | |||
261 | |||
262 | instance Transposable (L m n) (L n m) | ||
263 | where | ||
264 | tr (Dim (Dim a)) = Dim (Dim (tr a)) | ||
265 | |||
266 | |||
267 | eye :: forall n . KnownNat n => Sq n | ||
268 | eye = Dim (Dim (ident d)) | ||
269 | where | ||
270 | d = fromIntegral . natVal $ (undefined :: Proxy n) | ||
271 | |||
272 | |||
273 | -------------------------------------------------------------------------------- | ||
274 | |||
275 | test :: (Bool, IO ()) | ||
276 | test = (ok,info) | ||
277 | where | ||
278 | ok = d2 (eye :: Sq 5) == ident 5 | ||
279 | && d2 (mTm sm :: Sq 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..] | ||
280 | && d2 (tm :: L 3 5) == mat 5 [1..15] | ||
281 | && thingS == thingD | ||
282 | && precS == precD | ||
283 | |||
284 | info = do | ||
285 | print $ u | ||
286 | print $ v | ||
287 | print (eye :: Sq 3) | ||
288 | print $ ((u & 5) + 1) <ยท> v | ||
289 | print (tm :: L 2 5) | ||
290 | print (tm <> sm :: L 2 3) | ||
291 | print thingS | ||
292 | print thingD | ||
293 | print precS | ||
294 | print precD | ||
295 | |||
296 | u = vec2 3 5 | ||
297 | |||
298 | v = ๐ง 2 & 4 & 7 | ||
299 | |||
300 | mTm :: L n m -> Sq m | ||
301 | mTm a = tr a <> a | ||
302 | |||
303 | tm :: GL | ||
304 | tm = lmat 0 [1..] | ||
305 | |||
306 | lmat :: forall m n . (KnownNat m, KnownNat n) => โ -> [โ] -> L m n | ||
307 | lmat z xs = Dim . Dim . reshape n' . fromList . take (m'*n') $ xs ++ repeat z | ||
308 | where | ||
309 | m' = fromIntegral . natVal $ (undefined :: Proxy m) | ||
310 | n' = fromIntegral . natVal $ (undefined :: Proxy n) | ||
311 | |||
312 | sm :: GSq | ||
313 | sm = lmat 0 [1..] | ||
314 | |||
315 | thingS = (u & 1) <ยท> tr q #> q #> v | ||
316 | where | ||
317 | q = tm :: L 10 3 | ||
318 | |||
319 | thingD = vjoin [unDim u, 1] LA.<ยท> tr m LA.#> m LA.#> unDim v | ||
320 | where | ||
321 | m = mat 3 [1..30] | ||
322 | |||
323 | precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <ยท> konst 2 #> v | ||
324 | precD = 1 + 2 * vjoin[unDim u, 6] LA.<ยท> LA.konst 2 (size (unDim u) +1, size (unDim v)) LA.#> unDim v | ||
325 | |||
326 | |||
327 | instance (KnownNat n', KnownNat m') => Testable (L n' m') | ||
328 | where | ||
329 | checkT _ = test | ||
330 | |||
331 | {- | ||
332 | do (snd test) | ||
333 | fst test | ||
334 | -} | ||
335 | |||
336 | |||
337 | |||
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs index a319785..47b1090 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs | |||
@@ -32,7 +32,7 @@ module Numeric.LinearAlgebra.Util( | |||
32 | rand, randn, | 32 | rand, randn, |
33 | cross, | 33 | cross, |
34 | norm, | 34 | norm, |
35 | โ,โค,โ,โ,โn,โn,๐,i_C, --โ | 35 | โ,โค,โ,โ,๐,i_C, --โ |
36 | norm_1, norm_2, norm_0, norm_Inf, norm_Frob, norm_nuclear, | 36 | norm_1, norm_2, norm_0, norm_Inf, norm_Frob, norm_nuclear, |
37 | mnorm_1, mnorm_2, mnorm_0, mnorm_Inf, | 37 | mnorm_1, mnorm_2, mnorm_0, mnorm_Inf, |
38 | unitary, | 38 | unitary, |
@@ -70,8 +70,8 @@ type โ = Double | |||
70 | type โ = Int | 70 | type โ = Int |
71 | type โค = Int | 71 | type โค = Int |
72 | type โ = Complex Double | 72 | type โ = Complex Double |
73 | type โn = Vector โ | 73 | --type โn = Vector โ |
74 | type โn = Vector โ | 74 | --type โn = Vector โ |
75 | --newtype โ m = H m | 75 | --newtype โ m = H m |
76 | 76 | ||
77 | i_C, ๐ :: โ | 77 | i_C, ๐ :: โ |
@@ -84,7 +84,7 @@ i_C = ๐ | |||
84 | fromList [1.0,2.0,3.0,4.0,5.0] | 84 | fromList [1.0,2.0,3.0,4.0,5.0] |
85 | 85 | ||
86 | -} | 86 | -} |
87 | vect :: [โ] -> โn | 87 | vect :: [โ] -> Vector โ |
88 | vect = fromList | 88 | vect = fromList |
89 | 89 | ||
90 | {- | create a real matrix | 90 | {- | create a real matrix |
@@ -103,18 +103,6 @@ mat | |||
103 | mat c = reshape c . fromList | 103 | mat c = reshape c . fromList |
104 | 104 | ||
105 | 105 | ||
106 | |||
107 | class ( Container Vector t | ||
108 | , Container Matrix t | ||
109 | , Konst t Int Vector | ||
110 | , Konst t (Int,Int) Matrix | ||
111 | ) => Numeric t | ||
112 | |||
113 | instance Numeric Double | ||
114 | instance Numeric (Complex Double) | ||
115 | |||
116 | |||
117 | |||
118 | {- | print a real matrix with given number of digits after the decimal point | 106 | {- | print a real matrix with given number of digits after the decimal point |
119 | 107 | ||
120 | >>> disp 5 $ ident 2 / 3 | 108 | >>> disp 5 $ ident 2 / 3 |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs index 5e2ea84..50372f1 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs | |||
@@ -3,11 +3,14 @@ | |||
3 | 3 | ||
4 | module Numeric.LinearAlgebra.Util.CG( | 4 | module Numeric.LinearAlgebra.Util.CG( |
5 | cgSolve, cgSolve', | 5 | cgSolve, cgSolve', |
6 | CGMat, CGState(..), R, V | 6 | CGState(..), R, V |
7 | ) where | 7 | ) where |
8 | 8 | ||
9 | import Data.Packed.Numeric | 9 | import Data.Packed.Numeric |
10 | import Numeric.Sparse | ||
10 | import Numeric.Vector() | 11 | import Numeric.Vector() |
12 | import Numeric.LinearAlgebra.Algorithms(linearSolveLS, relativeError, NormType(..)) | ||
13 | import Control.Arrow((***)) | ||
11 | 14 | ||
12 | {- | 15 | {- |
13 | import Util.Misc(debug, debugMat) | 16 | import Util.Misc(debug, debugMat) |
@@ -51,7 +54,7 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx | |||
51 | rdx = norm2 dx / max 1 (norm2 x) | 54 | rdx = norm2 dx / max 1 (norm2 x) |
52 | 55 | ||
53 | conjugrad | 56 | conjugrad |
54 | :: (Transposable m, Contraction m V V) | 57 | :: (Transposable m mt, Contraction m V V, Contraction mt V V) |
55 | => Bool -> m -> V -> V -> R -> R -> [CGState] | 58 | => Bool -> m -> V -> V -> R -> R -> [CGState] |
56 | conjugrad sym a b = solveG (tr a โ) (a โ) (cg sym) b | 59 | conjugrad sym a b = solveG (tr a โ) (a โ) (cg sym) b |
57 | 60 | ||
@@ -82,27 +85,88 @@ takeUntil q xs = a++ take 1 b | |||
82 | where | 85 | where |
83 | (a,b) = break q xs | 86 | (a,b) = break q xs |
84 | 87 | ||
85 | class (Transposable m, Contraction m V V) => CGMat m | ||
86 | |||
87 | cgSolve | 88 | cgSolve |
88 | :: CGMat m | 89 | :: Bool -- ^ is symmetric |
89 | => Bool -- ^ is symmetric | 90 | -> GMatrix -- ^ coefficient matrix |
90 | -> m -- ^ coefficient matrix | ||
91 | -> Vector Double -- ^ right-hand side | 91 | -> Vector Double -- ^ right-hand side |
92 | -> Vector Double -- ^ solution | 92 | -> Vector Double -- ^ solution |
93 | cgSolve sym a b = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0 | 93 | cgSolve sym a b = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0 |
94 | where | 94 | where |
95 | n = max 10 (round $ sqrt (fromIntegral (dim b) :: Double)) | 95 | n = max 10 (round $ sqrt (fromIntegral (dim b) :: Double)) |
96 | 96 | ||
97 | cgSolve' | 97 | cgSolve' |
98 | :: CGMat m | 98 | :: Bool -- ^ symmetric |
99 | => Bool -- ^ symmetric | ||
100 | -> R -- ^ relative tolerance for the residual (e.g. 1E-4) | 99 | -> R -- ^ relative tolerance for the residual (e.g. 1E-4) |
101 | -> R -- ^ relative tolerance for ฮดx (e.g. 1E-3) | 100 | -> R -- ^ relative tolerance for ฮดx (e.g. 1E-3) |
102 | -> Int -- ^ maximum number of iterations | 101 | -> Int -- ^ maximum number of iterations |
103 | -> m -- ^ coefficient matrix | 102 | -> GMatrix -- ^ coefficient matrix |
104 | -> V -- ^ initial solution | 103 | -> V -- ^ initial solution |
105 | -> V -- ^ right-hand side | 104 | -> V -- ^ right-hand side |
106 | -> [CGState] -- ^ solution | 105 | -> [CGState] -- ^ solution |
107 | cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es | 106 | cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es |
108 | 107 | ||
108 | |||
109 | -------------------------------------------------------------------------------- | ||
110 | |||
111 | instance Testable GMatrix | ||
112 | where | ||
113 | checkT _ = (ok,info) | ||
114 | where | ||
115 | sma = convo2 20 3 | ||
116 | x1 = vect [1..20] | ||
117 | x2 = vect [1..40] | ||
118 | sm = mkSparse sma | ||
119 | dm = toDense sma | ||
120 | |||
121 | s1 = sm !#> x1 | ||
122 | d1 = dm #> x1 | ||
123 | |||
124 | s2 = tr sm !#> x2 | ||
125 | d2 = tr dm #> x2 | ||
126 | |||
127 | sdia = mkDiagR 40 20 (vect [1..10]) | ||
128 | s3 = sdia !#> x1 | ||
129 | s4 = tr sdia !#> x2 | ||
130 | ddia = diagRect 0 (vect [1..10]) 40 20 | ||
131 | d3 = ddia #> x1 | ||
132 | d4 = tr ddia #> x2 | ||
133 | |||
134 | v = testb 40 | ||
135 | s5 = cgSolve False sm v | ||
136 | d5 = denseSolve dm v | ||
137 | |||
138 | info = do | ||
139 | print sm | ||
140 | disp (toDense sma) | ||
141 | print s1; print d1 | ||
142 | print s2; print d2 | ||
143 | print s3; print d3 | ||
144 | print s4; print d4 | ||
145 | print s5; print d5 | ||
146 | print $ relativeError Infinity s5 d5 | ||
147 | |||
148 | ok = s1==d1 | ||
149 | && s2==d2 | ||
150 | && s3==d3 | ||
151 | && s4==d4 | ||
152 | && relativeError Infinity s5 d5 < 1E-10 | ||
153 | |||
154 | disp = putStr . dispf 2 | ||
155 | |||
156 | vect = fromList :: [Double] -> Vector Double | ||
157 | |||
158 | convomat :: Int -> Int -> AssocMatrix | ||
159 | convomat n k = [ ((i,j `mod` n),1) | i<-[0..n-1], j <- [i..i+k-1]] | ||
160 | |||
161 | convo2 :: Int -> Int -> AssocMatrix | ||
162 | convo2 n k = m1 ++ m2 | ||
163 | where | ||
164 | m1 = convomat n k | ||
165 | m2 = map (((+n) *** id) *** id) m1 | ||
166 | |||
167 | testb n = vect $ take n $ cycle ([0..10]++[9,8..1]) | ||
168 | |||
169 | denseSolve a = flatten . linearSolveLS a . asColumn | ||
170 | |||
171 | -- mkDiag v = mkDiagR (dim v) (dim v) v | ||
172 | |||
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Static.hs deleted file mode 100644 index a3f8eb0..0000000 --- a/packages/base/src/Numeric/LinearAlgebra/Util/Static.hs +++ /dev/null | |||
@@ -1,70 +0,0 @@ | |||
1 | {-# LANGUAGE DataKinds #-} | ||
2 | {-# LANGUAGE KindSignatures #-} | ||
3 | {-# LANGUAGE GeneralizedNewtypeDeriving #-} | ||
4 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
5 | {-# LANGUAGE FlexibleContexts #-} | ||
6 | {-# LANGUAGE ScopedTypeVariables #-} | ||
7 | {-# LANGUAGE EmptyDataDecls #-} | ||
8 | {-# LANGUAGE Rank2Types #-} | ||
9 | {-# LANGUAGE FlexibleInstances #-} | ||
10 | {-# LANGUAGE TypeOperators #-} | ||
11 | |||
12 | module Numeric.LinearAlgebra.Util.Static( | ||
13 | Static (ddata), | ||
14 | R, | ||
15 | vect0, sScalar, vect2, vect3, (&) | ||
16 | ) where | ||
17 | |||
18 | |||
19 | import GHC.TypeLits | ||
20 | import Data.Packed.Numeric | ||
21 | import Numeric.Vector() | ||
22 | import Numeric.LinearAlgebra.Util(Numeric,โ) | ||
23 | |||
24 | lift1F :: (Vector t -> Vector t) -> Static n (Vector t) -> Static n (Vector t) | ||
25 | lift1F f (Static v) = Static (f v) | ||
26 | |||
27 | lift2F :: (Vector t -> Vector t -> Vector t) -> Static n (Vector t) -> Static n (Vector t) -> Static n (Vector t) | ||
28 | lift2F f (Static u) (Static v) = Static (f u v) | ||
29 | |||
30 | newtype Static (n :: Nat) t = Static { ddata :: t } deriving Show | ||
31 | |||
32 | type R n = Static n (Vector โ) | ||
33 | |||
34 | |||
35 | infixl 4 & | ||
36 | (&) :: R n -> โ -> R (n+1) | ||
37 | Static v & x = Static (vjoin [v, scalar x]) | ||
38 | |||
39 | vect0 :: R 0 | ||
40 | vect0 = Static (fromList[]) | ||
41 | |||
42 | sScalar :: โ -> R 1 | ||
43 | sScalar = Static . scalar | ||
44 | |||
45 | |||
46 | vect2 :: โ -> โ -> R 2 | ||
47 | vect2 x1 x2 = Static (fromList [x1,x2]) | ||
48 | |||
49 | vect3 :: โ -> โ -> โ -> R 3 | ||
50 | vect3 x1 x2 x3 = Static (fromList [x1,x2,x3]) | ||
51 | |||
52 | |||
53 | |||
54 | |||
55 | |||
56 | |||
57 | instance forall n t . (KnownNat n, Num (Vector t), Numeric t )=> Num (Static n (Vector t)) | ||
58 | where | ||
59 | (+) = lift2F add | ||
60 | (*) = lift2F mul | ||
61 | (-) = lift2F sub | ||
62 | abs = lift1F abs | ||
63 | signum = lift1F signum | ||
64 | negate = lift1F (scale (-1)) | ||
65 | fromInteger x = Static (konst (fromInteger x) d) | ||
66 | where | ||
67 | d = fromIntegral . natVal $ (undefined :: Proxy n) | ||
68 | |||
69 | data Proxy :: Nat -> * | ||
70 | |||