diff options
author | Alberto Ruiz <aruiz@um.es> | 2015-06-13 19:18:16 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2015-06-13 19:18:16 +0200 |
commit | 717c680a4b65a2226b0dd6fc13f7c63e7bc0431d (patch) | |
tree | 1775c3c363a0b61f5f6a6ec1f22fe9b7d5864dc4 /packages/base/src/Internal/Modular.hs | |
parent | 4b3e29097aa272d429f8005fe17b459cf0c049c8 (diff) |
setRect, general luPacked' based on luST
Diffstat (limited to 'packages/base/src/Internal/Modular.hs')
-rw-r--r-- | packages/base/src/Internal/Modular.hs | 75 |
1 files changed, 61 insertions, 14 deletions
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 824fc57..3b27310 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs | |||
@@ -33,12 +33,15 @@ import Internal.Element | |||
33 | import Internal.Container | 33 | import Internal.Container |
34 | import Internal.Vectorized (prodI,sumI,prodL,sumL) | 34 | import Internal.Vectorized (prodI,sumI,prodL,sumL) |
35 | import Internal.LAPACK (multiplyI, multiplyL) | 35 | import Internal.LAPACK (multiplyI, multiplyL) |
36 | import Internal.Util(Indexable(..),gaussElim) | 36 | import Internal.Algorithms(luFact) |
37 | import Internal.Util(Normed(..),Indexable(..),gaussElim, gaussElim_1, gaussElim_2,luST, magnit) | ||
38 | import Internal.ST(mutable) | ||
37 | import GHC.TypeLits | 39 | import GHC.TypeLits |
38 | import Data.Proxy(Proxy) | 40 | import Data.Proxy(Proxy) |
39 | import Foreign.ForeignPtr(castForeignPtr) | 41 | import Foreign.ForeignPtr(castForeignPtr) |
40 | import Foreign.Storable | 42 | import Foreign.Storable |
41 | import Data.Ratio | 43 | import Data.Ratio |
44 | import Data.Complex | ||
42 | 45 | ||
43 | 46 | ||
44 | 47 | ||
@@ -116,6 +119,7 @@ instance KnownNat m => Element (Mod m I) | |||
116 | transdata n v m = i2f (transdata n (f2i v) m) | 119 | transdata n v m = i2f (transdata n (f2i v) m) |
117 | constantD x n = i2f (constantD (unMod x) n) | 120 | constantD x n = i2f (constantD (unMod x) n) |
118 | extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js | 121 | extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js |
122 | setRect i j m x = setRect i j (f2iM m) (f2iM x) | ||
119 | sortI = sortI . f2i | 123 | sortI = sortI . f2i |
120 | sortV = i2f . sortV . f2i | 124 | sortV = i2f . sortV . f2i |
121 | compareV u v = compareV (f2i u) (f2i v) | 125 | compareV u v = compareV (f2i u) (f2i v) |
@@ -130,6 +134,7 @@ instance KnownNat m => Element (Mod m Z) | |||
130 | transdata n v m = i2f (transdata n (f2i v) m) | 134 | transdata n v m = i2f (transdata n (f2i v) m) |
131 | constantD x n = i2f (constantD (unMod x) n) | 135 | constantD x n = i2f (constantD (unMod x) n) |
132 | extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js | 136 | extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js |
137 | setRect i j m x = setRect i j (f2iM m) (f2iM x) | ||
133 | sortI = sortI . f2i | 138 | sortI = sortI . f2i |
134 | sortV = i2f . sortV . f2i | 139 | sortV = i2f . sortV . f2i |
135 | compareV u v = compareV (f2i u) (f2i v) | 140 | compareV u v = compareV (f2i u) (f2i v) |
@@ -139,18 +144,6 @@ instance KnownNat m => Element (Mod m Z) | |||
139 | where | 144 | where |
140 | m' = fromIntegral . natVal $ (undefined :: Proxy m) | 145 | m' = fromIntegral . natVal $ (undefined :: Proxy m) |
141 | 146 | ||
142 | {- | ||
143 | instance (Ord t, Element t) => Element (Mod m t) | ||
144 | where | ||
145 | transdata n v m = i2f (transdata n (f2i v) m) | ||
146 | constantD x n = i2f (constantD (unMod x) n) | ||
147 | extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js | ||
148 | sortI = sortI . f2i | ||
149 | sortV = i2f . sortV . f2i | ||
150 | compareV u v = compareV (f2i u) (f2i v) | ||
151 | selectV c l e g = i2f (selectV c (f2i l) (f2i e) (f2i g)) | ||
152 | remapM i j m = i2fM (remap i j (f2iM m)) | ||
153 | -} | ||
154 | 147 | ||
155 | instance forall m . KnownNat m => Container Vector (Mod m I) | 148 | instance forall m . KnownNat m => Container Vector (Mod m I) |
156 | where | 149 | where |
@@ -258,6 +251,20 @@ instance KnownNat m => Product (Mod m Z) where | |||
258 | where | 251 | where |
259 | m' = fromIntegral . natVal $ (undefined :: Proxy m) | 252 | m' = fromIntegral . natVal $ (undefined :: Proxy m) |
260 | 253 | ||
254 | instance KnownNat m => Normed (Vector (Mod m I)) | ||
255 | where | ||
256 | norm_0 = norm_0 . toInt | ||
257 | norm_1 = norm_1 . toInt | ||
258 | norm_2 = norm_2 . toInt | ||
259 | norm_Inf = norm_Inf . toInt | ||
260 | |||
261 | instance KnownNat m => Normed (Vector (Mod m Z)) | ||
262 | where | ||
263 | norm_0 = norm_0 . toZ | ||
264 | norm_1 = norm_1 . toZ | ||
265 | norm_2 = norm_2 . toZ | ||
266 | norm_Inf = norm_Inf . toZ | ||
267 | |||
261 | 268 | ||
262 | instance KnownNat m => Numeric (Mod m I) | 269 | instance KnownNat m => Numeric (Mod m I) |
263 | instance KnownNat m => Numeric (Mod m Z) | 270 | instance KnownNat m => Numeric (Mod m Z) |
@@ -334,6 +341,15 @@ test = (ok, info) | |||
334 | lg = (3><3) (repeat (3*10^(9::Int))) :: Matrix Z | 341 | lg = (3><3) (repeat (3*10^(9::Int))) :: Matrix Z |
335 | lgm = fromZ lg :: Matrix (Mod 10000000000 Z) | 342 | lgm = fromZ lg :: Matrix (Mod 10000000000 Z) |
336 | 343 | ||
344 | gen n = diagRect 1 (konst 5 n) n n :: Numeric t => Matrix t | ||
345 | |||
346 | checkGen x = norm_Inf $ flatten $ invg x <> x - ident (rows x) | ||
347 | |||
348 | invg t = gaussElim t (ident (rows t)) | ||
349 | |||
350 | checkLU okf t = norm_Inf $ flatten (l <> u <> p - t) | ||
351 | where | ||
352 | (l,u,p,_ :: Int) = luFact $ mutable (luST okf) t | ||
337 | 353 | ||
338 | info = do | 354 | info = do |
339 | print v | 355 | print v |
@@ -356,11 +372,42 @@ test = (ok, info) | |||
356 | print $ lg <> lg | 372 | print $ lg <> lg |
357 | print lgm | 373 | print lgm |
358 | print $ lgm <> lgm | 374 | print $ lgm <> lgm |
375 | |||
376 | print (checkGen (gen 5 :: Matrix R)) | ||
377 | print (checkGen (gen 5 :: Matrix C)) | ||
378 | print (checkGen (gen 5 :: Matrix Float)) | ||
379 | print (checkGen (gen 5 :: Matrix (Complex Float))) | ||
380 | print (invg (gen 5) :: Matrix (Mod 7 I)) | ||
381 | print (invg (gen 5) :: Matrix (Mod 7 Z)) | ||
382 | |||
383 | print $ mutable (luST (const True)) (gen 5 :: Matrix R) | ||
384 | print $ mutable (luST (const True)) (gen 5 :: Matrix (Mod 11 Z)) | ||
385 | |||
386 | print $ checkLU (magnit 0) (gen 5 :: Matrix R) | ||
387 | print $ checkLU (magnit 0) (gen 5 :: Matrix Float) | ||
388 | print $ checkLU (magnit 0) (gen 5 :: Matrix C) | ||
389 | print $ checkLU (magnit 0) (gen 5 :: Matrix (Complex Float)) | ||
390 | print $ checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 I)) | ||
391 | print $ checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 Z)) | ||
359 | 392 | ||
360 | 393 | ||
361 | ok = and | 394 | ok = and |
362 | [ toInt (m #> v) == cmod 11 (toInt m #> toInt v ) | 395 | [ toInt (m #> v) == cmod 11 (toInt m #> toInt v ) |
363 | , am <> gaussElim am bm == bm | 396 | , am <> gaussElim_1 am bm == bm |
397 | , am <> gaussElim_2 am bm == bm | ||
398 | , am <> gaussElim am bm == bm | ||
399 | , (checkGen (gen 5 :: Matrix R)) < 1E-15 | ||
400 | , (checkGen (gen 5 :: Matrix Float)) < 1E-7 | ||
401 | , (checkGen (gen 5 :: Matrix C)) < 1E-15 | ||
402 | , (checkGen (gen 5 :: Matrix (Complex Float))) < 1E-7 | ||
403 | , (checkGen (gen 5 :: Matrix (Mod 7 I))) == 0 | ||
404 | , (checkGen (gen 5 :: Matrix (Mod 7 Z))) == 0 | ||
405 | , (checkLU (magnit 1E-10) (gen 5 :: Matrix R)) < 1E-15 | ||
406 | , (checkLU (magnit 1E-5) (gen 5 :: Matrix Float)) < 1E-6 | ||
407 | , (checkLU (magnit 1E-10) (gen 5 :: Matrix C)) < 1E-15 | ||
408 | , (checkLU (magnit 1E-5) (gen 5 :: Matrix (Complex Float))) < 1E-6 | ||
409 | , (checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 I))) == 0 | ||
410 | , (checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 Z))) == 0 | ||
364 | , prodElements (konst (9:: Mod 10 I) (12::Int)) == product (replicate 12 (9:: Mod 10 I)) | 411 | , prodElements (konst (9:: Mod 10 I) (12::Int)) == product (replicate 12 (9:: Mod 10 I)) |
365 | , gm <> gm == konst 0 (3,3) | 412 | , gm <> gm == konst 0 (3,3) |
366 | , lgm <> lgm == konst 0 (3,3) | 413 | , lgm <> lgm == konst 0 (3,3) |