diff options
Diffstat (limited to 'packages')
4 files changed, 40 insertions, 1 deletions
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs index b834b20..c9f5d0f 100644 --- a/packages/base/src/Internal/LAPACK.hs +++ b/packages/base/src/Internal/LAPACK.hs | |||
@@ -302,6 +302,14 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | |||
302 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) | 302 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) |
303 | fixeig _ _ = error "fixeig with impossible inputs" | 303 | fixeig _ _ = error "fixeig with impossible inputs" |
304 | 304 | ||
305 | -- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0. | ||
306 | -- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1) | ||
307 | -- can be different. Therefore old 'fixeig' would fail for 'eigG'. | ||
308 | fixeigG [] _ = [] | ||
309 | fixeigG [_] [v] = [comp' v] | ||
310 | fixeigG ((ar1:+ai1) : an : as) (v1:v2:vs) | ||
311 | | abs ai1 > 1e-13 = toComplex' (v1, v2) : toComplex' (v1, mapVector negate v2) : fixeigG as vs | ||
312 | | otherwise = comp' v1 : fixeigG (an:as) (v2:vs) | ||
305 | 313 | ||
306 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. | 314 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. |
307 | -- The eigenvalues are not sorted. | 315 | -- The eigenvalues are not sorted. |
@@ -316,7 +324,7 @@ eigG a b = (alpha', beta, v'') | |||
316 | (alpha, beta, v) = eigGaux dggev a b "eigG" | 324 | (alpha, beta, v) = eigGaux dggev a b "eigG" |
317 | alpha' = fixeig1 alpha | 325 | alpha' = fixeig1 alpha |
318 | v' = toRows $ trans v | 326 | v' = toRows $ trans v |
319 | v'' = fromColumns $ fixeig (toList alpha') v' | 327 | v'' = fromColumns $ fixeigG (toList alpha') v' |
320 | 328 | ||
321 | eigGaux f ma mb st = unsafePerformIO $ do | 329 | eigGaux f ma mb st = unsafePerformIO $ do |
322 | a <- copy ColumnMajor ma | 330 | a <- copy ColumnMajor ma |
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index 4ed1462..c0c151a 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs | |||
@@ -756,6 +756,9 @@ runTests n = do | |||
756 | test (eigSHProp2 . cHer) | 756 | test (eigSHProp2 . cHer) |
757 | test (eigProp2 . rSq) | 757 | test (eigProp2 . rSq) |
758 | test (eigProp2 . cSq) | 758 | test (eigProp2 . cSq) |
759 | putStrLn "------ geig" | ||
760 | test (uncurry geigProp . rSq2WC) | ||
761 | test (uncurry geigProp . cSq2WC) | ||
759 | putStrLn "------ nullSpace" | 762 | putStrLn "------ nullSpace" |
760 | test (nullspaceProp . rM) | 763 | test (nullspaceProp . rM) |
761 | test (nullspaceProp . cM) | 764 | test (nullspaceProp . cM) |
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs index 59230e0..97cfd01 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs | |||
@@ -17,6 +17,7 @@ Arbitrary instances for vectors, matrices. | |||
17 | 17 | ||
18 | module Numeric.LinearAlgebra.Tests.Instances( | 18 | module Numeric.LinearAlgebra.Tests.Instances( |
19 | Sq(..), rSq,cSq, | 19 | Sq(..), rSq,cSq, |
20 | Sq2WC(..), rSq2WC,cSq2WC, | ||
20 | Rot(..), rRot,cRot, | 21 | Rot(..), rRot,cRot, |
21 | rHer,cHer, | 22 | rHer,cHer, |
22 | WC(..), rWC,cWC, | 23 | WC(..), rWC,cWC, |
@@ -105,6 +106,23 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where | |||
105 | 106 | ||
106 | shrink (Sq a) = [ Sq b | b <- shrink a ] | 107 | shrink (Sq a) = [ Sq b | b <- shrink a ] |
107 | 108 | ||
109 | -- a pair of square matrices | ||
110 | newtype (Sq2WC a) = Sq2WC (Matrix a, Matrix a) deriving Show | ||
111 | instance (ArbitraryField a, Numeric a) => Arbitrary (Sq2WC a) where | ||
112 | arbitrary = do | ||
113 | n <- chooseDim | ||
114 | l <- vector (n*n) | ||
115 | r <- vector (n*n) | ||
116 | l' <- makeWC $ (n><n) l | ||
117 | r' <- makeWC $ (n><n) r | ||
118 | return $ Sq2WC (l', r') | ||
119 | where | ||
120 | makeWC m = do | ||
121 | let (u,_,v) = svd m | ||
122 | n = rows m | ||
123 | sv' <- replicateM n (choose (1,100)) | ||
124 | let s = diag (fromList sv') | ||
125 | return $ u <> real s <> tr v | ||
108 | 126 | ||
109 | -- a unitary matrix | 127 | -- a unitary matrix |
110 | newtype (Rot a) = Rot (Matrix a) deriving Show | 128 | newtype (Rot a) = Rot (Matrix a) deriving Show |
@@ -204,6 +222,9 @@ cRot (Rot m) = m :: CM | |||
204 | rSq (Sq m) = m :: RM | 222 | rSq (Sq m) = m :: RM |
205 | cSq (Sq m) = m :: CM | 223 | cSq (Sq m) = m :: CM |
206 | 224 | ||
225 | rSq2WC (Sq2WC (a, b)) = (a, b) :: (RM, RM) | ||
226 | cSq2WC (Sq2WC (a, b)) = (a, b) :: (CM, CM) | ||
227 | |||
207 | rWC (WC m) = m :: RM | 228 | rWC (WC m) = m :: RM |
208 | cWC (WC m) = m :: CM | 229 | cWC (WC m) = m :: CM |
209 | 230 | ||
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs index 6cd3a9e..38aa977 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs | |||
@@ -36,6 +36,7 @@ module Numeric.LinearAlgebra.Tests.Properties ( | |||
36 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, | 36 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, |
37 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, | 37 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, |
38 | eigProp, eigSHProp, eigProp2, eigSHProp2, | 38 | eigProp, eigSHProp, eigProp2, eigSHProp2, |
39 | geigProp, | ||
39 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, | 40 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, |
40 | hessProp, | 41 | hessProp, |
41 | schurProp1, schurProp2, | 42 | schurProp1, schurProp2, |
@@ -237,6 +238,12 @@ eigProp2 m = fst (eig m) |~| eigenvalues m | |||
237 | 238 | ||
238 | eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m | 239 | eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m |
239 | 240 | ||
241 | geigProp a b = a' <> v <> diag betas' |~| b' <> v <> diag alphas | ||
242 | where (alphas, betas, v) = geig a b | ||
243 | betas' = complex betas | ||
244 | a' = complex a | ||
245 | b' = complex b | ||
246 | |||
240 | ------------------------------------------------------------------ | 247 | ------------------------------------------------------------------ |
241 | 248 | ||
242 | qrProp m = q <> r |~| m && unitary q && upperTriang r | 249 | qrProp m = q <> r |~| m && unitary q && upperTriang r |