diff options
-rw-r--r-- | examples/tests.hs | 20 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 2 | ||||
-rw-r--r-- | lib/LAPACK/Internal.hs | 4 |
3 files changed, 25 insertions, 1 deletions
diff --git a/examples/tests.hs b/examples/tests.hs index f167b92..5af33ba 100644 --- a/examples/tests.hs +++ b/examples/tests.hs | |||
@@ -66,6 +66,8 @@ aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length | |||
66 | (|~|) = aprox abs | 66 | (|~|) = aprox abs |
67 | (|~~|) = aprox magnitude | 67 | (|~~|) = aprox magnitude |
68 | 68 | ||
69 | v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2 | ||
70 | |||
69 | eps = 1E-8::Double | 71 | eps = 1E-8::Double |
70 | 72 | ||
71 | asFortran m = (rows m >|< cols m) $ toList (fdat m) | 73 | asFortran m = (rows m >|< cols m) $ toList (fdat m) |
@@ -116,6 +118,15 @@ instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where | |||
116 | --return $ PairM ((a><b) l1) ((b><c) l2) | 118 | --return $ PairM ((a><b) l1) ((b><c) l2) |
117 | coarbitrary = undefined | 119 | coarbitrary = undefined |
118 | 120 | ||
121 | data SqM a = SqM (Matrix a) deriving Show | ||
122 | instance (Field a, Arbitrary a) => Arbitrary (SqM a) where | ||
123 | arbitrary = do | ||
124 | n <- choose (1,10) | ||
125 | l <- vector (n*n) | ||
126 | return $ SqM $ (n><n) l | ||
127 | coarbitrary = undefined | ||
128 | |||
129 | |||
119 | type BaseType = Double | 130 | type BaseType = Double |
120 | 131 | ||
121 | 132 | ||
@@ -133,6 +144,14 @@ svdTestC fun prod m = u <> s' <> (trans v) |~~| m | |||
133 | (<>) = prod | 144 | (<>) = prod |
134 | s' = liftMatrix comp s | 145 | s' = liftMatrix comp s |
135 | 146 | ||
147 | eigTestC fun prod (SqM m) = (m <> v) |~~| (v <> diag s) | ||
148 | && takeDiag ((liftMatrix conj (trans v)) `mulC` v) ~~ constant (rows m) 1 | ||
149 | where (s,v) = fun m | ||
150 | (<>) = prod | ||
151 | |||
152 | takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] | ||
153 | |||
154 | |||
136 | comp v = toComplex (v,constant (dim v) 0) | 155 | comp v = toComplex (v,constant (dim v) 0) |
137 | 156 | ||
138 | main = do | 157 | main = do |
@@ -147,3 +166,4 @@ main = do | |||
147 | quickCheck (svdTestR svdR mulF) | 166 | quickCheck (svdTestR svdR mulF) |
148 | quickCheck (svdTestC svdC mulC) | 167 | quickCheck (svdTestC svdC mulC) |
149 | quickCheck (svdTestC svdC mulF) | 168 | quickCheck (svdTestC svdC mulF) |
169 | quickCheck (eigTestC eigC mulC) | ||
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 0664e9b..b8de245 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -108,6 +108,8 @@ createMatrix order r c = do | |||
108 | 108 | ||
109 | reshape c v = matrixFromVector RowMajor c v | 109 | reshape c v = matrixFromVector RowMajor c v |
110 | 110 | ||
111 | singleton x = reshape 1 (fromList [x]) | ||
112 | |||
111 | transdataG :: Storable a => Int -> Vector a -> Int -> Vector a | 113 | transdataG :: Storable a => Int -> Vector a -> Int -> Vector a |
112 | transdataG c1 d c2 = fromList . concat . transpose . partit c1 . toList $ d | 114 | transdataG c1 d c2 = fromList . concat . transpose . partit c1 . toList $ d |
113 | 115 | ||
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs index e39dd10..e3c9927 100644 --- a/lib/LAPACK/Internal.hs +++ b/lib/LAPACK/Internal.hs | |||
@@ -76,7 +76,9 @@ foreign import ccall "lapack-aux.h eig_l_C" | |||
76 | -- The eigenvectors are the columns of v. | 76 | -- The eigenvectors are the columns of v. |
77 | -- The eigenvalues are not sorted. | 77 | -- The eigenvalues are not sorted. |
78 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) | 78 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) |
79 | eigC (m@M {rows = r}) = unsafePerformIO $ do | 79 | eigC (m@M {rows = r}) |
80 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) | ||
81 | | otherwise = unsafePerformIO $ do | ||
80 | l <- createVector r | 82 | l <- createVector r |
81 | v <- createMatrix ColumnMajor r r | 83 | v <- createMatrix ColumnMajor r r |
82 | dummy <- createMatrix ColumnMajor 1 1 | 84 | dummy <- createMatrix ColumnMajor 1 1 |