summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs9
-rw-r--r--lib/Numeric/LinearAlgebra/Instances.hs276
-rw-r--r--lib/Numeric/LinearAlgebra/Interface.hs118
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK.hs3
-rw-r--r--lib/Numeric/LinearAlgebra/Linear.hs160
5 files changed, 6 insertions, 560 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 8306961..fa6c475 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -77,15 +77,16 @@ module Numeric.LinearAlgebra.Algorithms (
77import Data.Packed.Internal hiding ((//)) 77import Data.Packed.Internal hiding ((//))
78import Data.Packed.Matrix 78import Data.Packed.Matrix
79import Data.Complex 79import Data.Complex
80import Numeric.LinearAlgebra.Linear 80--import Numeric.LinearAlgebra.Linear
81import Numeric.LinearAlgebra.LAPACK as LAPACK 81import Numeric.LinearAlgebra.LAPACK as LAPACK
82import Data.List(foldl1') 82import Data.List(foldl1')
83import Data.Array 83import Data.Array
84import Numeric.Vector 84import Numeric.Container
85import Numeric.Matrix() 85
86constant x = constantD x
86 87
87-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 88-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
88class (Product t, Linear Vector t, Container Vector t, Container Matrix t) => Field t where 89class (Product t, Container Vector t, Container Matrix t) => Field t where
89 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 90 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
90 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 91 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
91 sv' :: Matrix t -> Vector Double 92 sv' :: Matrix t -> Vector Double
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
deleted file mode 100644
index 04a9d88..0000000
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ /dev/null
@@ -1,276 +0,0 @@
1{-# LANGUAGE UndecidableInstances, FlexibleInstances #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.LinearAlgebra.Instances
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : portable
11
12This module exports Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix.
13
14In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
15
16-}
17-----------------------------------------------------------------------------
18
19module Numeric.LinearAlgebra.Instances(
20) where
21
22import Numeric.LinearAlgebra.Linear
23import Numeric.GSL.Vector
24import Data.Packed.Matrix
25import Data.Complex
26import Data.List(transpose,intersperse)
27import Data.Packed.Internal.Vector
28
29#ifndef VECTOR
30import Foreign(Storable)
31#endif
32
33------------------------------------------------------------------
34
35instance (Show a, Element a) => (Show (Matrix a)) where
36 show m = (sizes++) . dsp . map (map show) . toLists $ m
37 where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"
38
39dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
40 where
41 mt = transpose as
42 longs = map (maximum . map length) mt
43 mtp = zipWith (\a b -> map (pad a) b) longs mt
44 pad n str = replicate (n - length str) ' ' ++ str
45 unwords' = concat . intersperse ", "
46
47#ifndef VECTOR
48
49instance (Show a, Storable a) => (Show (Vector a)) where
50 show v = (show (dim v))++" |> " ++ show (toList v)
51
52#endif
53
54------------------------------------------------------------------
55
56instance (Element a, Read a) => Read (Matrix a) where
57 readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
58 where (thing,rest) = breakAt ']' s
59 (dims,listnums) = breakAt ')' thing
60 cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
61 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
62
63#ifdef VECTOR
64
65instance (Element a, Read a) => Read (Vector a) where
66 readsPrec _ s = [(fromList . read $ listnums, rest)]
67 where (thing,trest) = breakAt ']' s
68 (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing)
69 rest = drop 31 trest
70#else
71
72instance (Element a, Read a) => Read (Vector a) where
73 readsPrec _ s = [((d |>) . read $ listnums, rest)]
74 where (thing,rest) = breakAt ']' s
75 (dims,listnums) = breakAt '>' thing
76 d = read . init . fst . breakAt '|' $ dims
77
78#endif
79
80breakAt c l = (a++[c],tail b) where
81 (a,b) = break (==c) l
82
83------------------------------------------------------------------
84
85adaptScalar f1 f2 f3 x y
86 | dim x == 1 = f1 (x@>0) y
87 | dim y == 1 = f3 x (y@>0)
88 | otherwise = f2 x y
89
90#ifndef VECTOR
91
92instance Linear Vector a => Eq (Vector a) where
93 (==) = equal
94
95#endif
96
97instance Num (Vector Float) where
98 (+) = adaptScalar addConstant add (flip addConstant)
99 negate = scale (-1)
100 (*) = adaptScalar scale mul (flip scale)
101 signum = vectorMapF Sign
102 abs = vectorMapF Abs
103 fromInteger = fromList . return . fromInteger
104
105instance Num (Vector Double) where
106 (+) = adaptScalar addConstant add (flip addConstant)
107 negate = scale (-1)
108 (*) = adaptScalar scale mul (flip scale)
109 signum = vectorMapR Sign
110 abs = vectorMapR Abs
111 fromInteger = fromList . return . fromInteger
112
113instance Num (Vector (Complex Double)) where
114 (+) = adaptScalar addConstant add (flip addConstant)
115 negate = scale (-1)
116 (*) = adaptScalar scale mul (flip scale)
117 signum = vectorMapC Sign
118 abs = vectorMapC Abs
119 fromInteger = fromList . return . fromInteger
120
121instance Num (Vector (Complex Float)) where
122 (+) = adaptScalar addConstant add (flip addConstant)
123 negate = scale (-1)
124 (*) = adaptScalar scale mul (flip scale)
125 signum = vectorMapQ Sign
126 abs = vectorMapQ Abs
127 fromInteger = fromList . return . fromInteger
128
129instance Linear Matrix a => Eq (Matrix a) where
130 (==) = equal
131
132instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
133 (+) = liftMatrix2Auto (+)
134 (-) = liftMatrix2Auto (-)
135 negate = liftMatrix negate
136 (*) = liftMatrix2Auto (*)
137 signum = liftMatrix signum
138 abs = liftMatrix abs
139 fromInteger = (1><1) . return . fromInteger
140
141---------------------------------------------------
142
143instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where
144 fromRational n = fromList [fromRational n]
145 (/) = adaptScalar f divide g where
146 r `f` v = scaleRecip r v
147 v `g` r = scale (recip r) v
148
149-------------------------------------------------------
150
151instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
152 fromRational n = (1><1) [fromRational n]
153 (/) = liftMatrix2Auto (/)
154
155---------------------------------------------------------
156
157instance Floating (Vector Float) where
158 sin = vectorMapF Sin
159 cos = vectorMapF Cos
160 tan = vectorMapF Tan
161 asin = vectorMapF ASin
162 acos = vectorMapF ACos
163 atan = vectorMapF ATan
164 sinh = vectorMapF Sinh
165 cosh = vectorMapF Cosh
166 tanh = vectorMapF Tanh
167 asinh = vectorMapF ASinh
168 acosh = vectorMapF ACosh
169 atanh = vectorMapF ATanh
170 exp = vectorMapF Exp
171 log = vectorMapF Log
172 sqrt = vectorMapF Sqrt
173 (**) = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS))
174 pi = fromList [pi]
175
176-------------------------------------------------------------
177
178instance Floating (Vector Double) where
179 sin = vectorMapR Sin
180 cos = vectorMapR Cos
181 tan = vectorMapR Tan
182 asin = vectorMapR ASin
183 acos = vectorMapR ACos
184 atan = vectorMapR ATan
185 sinh = vectorMapR Sinh
186 cosh = vectorMapR Cosh
187 tanh = vectorMapR Tanh
188 asinh = vectorMapR ASinh
189 acosh = vectorMapR ACosh
190 atanh = vectorMapR ATanh
191 exp = vectorMapR Exp
192 log = vectorMapR Log
193 sqrt = vectorMapR Sqrt
194 (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))
195 pi = fromList [pi]
196
197-------------------------------------------------------------
198
199instance Floating (Vector (Complex Double)) where
200 sin = vectorMapC Sin
201 cos = vectorMapC Cos
202 tan = vectorMapC Tan
203 asin = vectorMapC ASin
204 acos = vectorMapC ACos
205 atan = vectorMapC ATan
206 sinh = vectorMapC Sinh
207 cosh = vectorMapC Cosh
208 tanh = vectorMapC Tanh
209 asinh = vectorMapC ASinh
210 acosh = vectorMapC ACosh
211 atanh = vectorMapC ATanh
212 exp = vectorMapC Exp
213 log = vectorMapC Log
214 sqrt = vectorMapC Sqrt
215 (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))
216 pi = fromList [pi]
217
218-----------------------------------------------------------
219
220instance Floating (Vector (Complex Float)) where
221 sin = vectorMapQ Sin
222 cos = vectorMapQ Cos
223 tan = vectorMapQ Tan
224 asin = vectorMapQ ASin
225 acos = vectorMapQ ACos
226 atan = vectorMapQ ATan
227 sinh = vectorMapQ Sinh
228 cosh = vectorMapQ Cosh
229 tanh = vectorMapQ Tanh
230 asinh = vectorMapQ ASinh
231 acosh = vectorMapQ ACosh
232 atanh = vectorMapQ ATanh
233 exp = vectorMapQ Exp
234 log = vectorMapQ Log
235 sqrt = vectorMapQ Sqrt
236 (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
237 pi = fromList [pi]
238
239-----------------------------------------------------------
240
241instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
242 sin = liftMatrix sin
243 cos = liftMatrix cos
244 tan = liftMatrix tan
245 asin = liftMatrix asin
246 acos = liftMatrix acos
247 atan = liftMatrix atan
248 sinh = liftMatrix sinh
249 cosh = liftMatrix cosh
250 tanh = liftMatrix tanh
251 asinh = liftMatrix asinh
252 acosh = liftMatrix acosh
253 atanh = liftMatrix atanh
254 exp = liftMatrix exp
255 log = liftMatrix log
256 (**) = liftMatrix2Auto (**)
257 sqrt = liftMatrix sqrt
258 pi = (1><1) [pi]
259
260---------------------------------------------------------------
261
262-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
263-- mempty = 0 { idim = 0 }
264-- mappend a b = mconcat [a,b]
265-- mconcat = j . filter ((>0).dim)
266-- where j [] = mempty
267-- j l = join l
268
269---------------------------------------------------------------
270
271-- instance (NFData a, Storable a) => NFData (Vector a) where
272-- rnf = rnf . (@>0)
273--
274-- instance (NFData a, Element a) => NFData (Matrix a) where
275-- rnf = rnf . flatten
276
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
deleted file mode 100644
index fa3e209..0000000
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ /dev/null
@@ -1,118 +0,0 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2{-# LANGUAGE UndecidableInstances #-}
3-----------------------------------------------------------------------------
4{- |
5Module : Numeric.LinearAlgebra.Interface
6Copyright : (c) Alberto Ruiz 2007
7License : GPL-style
8
9Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional
11Portability : portable
12
13Some useful operators, and Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix.
14
15In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
16
17
18-}
19-----------------------------------------------------------------------------
20
21module Numeric.LinearAlgebra.Interface(
22 (<>),(<.>),
23 (<\>),
24 (.*),(*/),
25 (<|>),(<->),
26) where
27
28import Numeric.Vector
29import Numeric.Matrix
30import Numeric.LinearAlgebra.Algorithms
31import Numeric.LinearAlgebra.Linear
32
33class Mul a b c | a b -> c where
34 infixl 7 <>
35 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
36 (<>) :: Product t => a t -> b t -> c t
37
38instance Mul Matrix Matrix Matrix where
39 (<>) = mXm
40
41instance Mul Matrix Vector Vector where
42 (<>) m v = flatten $ m <> (asColumn v)
43
44instance Mul Vector Matrix Vector where
45 (<>) v m = flatten $ (asRow v) <> m
46
47---------------------------------------------------
48
49-- | Dot product: @u \<.\> v = dot u v@
50--(<.>) :: (Field t) => Vector t -> Vector t -> t
51(<.>) :: Vectors Vector t => Vector t -> Vector t -> t
52infixl 7 <.>
53(<.>) = dot
54
55----------------------------------------------------
56
57{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
58
59-- -- | @x .* a = scale x a@
60-- (.*) :: (Linear c a) => a -> c a -> c a
61infixl 7 .*
62a .* x = scale a x
63
64----------------------------------------------------
65
66{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
67
68-- -- | @a *\/ x = scale (recip x) a@
69-- (*/) :: (Linear c a) => c a -> a -> c a
70infixl 7 */
71v */ x = scale (recip x) v
72
73-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
74(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
75infixl 7 <\>
76m <\> v = flatten (linearSolveSVD m (reshape 1 v))
77
78------------------------------------------------
79
80{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
81{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
82
83class Joinable a b where
84 joinH :: Element t => a t -> b t -> Matrix t
85 joinV :: Element t => a t -> b t -> Matrix t
86
87instance Joinable Matrix Matrix where
88 joinH m1 m2 = fromBlocks [[m1,m2]]
89 joinV m1 m2 = fromBlocks [[m1],[m2]]
90
91instance Joinable Matrix Vector where
92 joinH m v = joinH m (asColumn v)
93 joinV m v = joinV m (asRow v)
94
95instance Joinable Vector Matrix where
96 joinH v m = joinH (asColumn v) m
97 joinV v m = joinV (asRow v) m
98
99infixl 4 <|>
100infixl 3 <->
101
102{-- - | Horizontal concatenation of matrices and vectors:
103
104@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
105(6><4)
106 [ 1.0, 0.0, 0.0, 1.0
107 , 0.0, 1.0, 0.0, 2.0
108 , 0.0, 0.0, 1.0, 3.0
109 , 3.0, 0.0, 0.0, 4.0
110 , 0.0, 3.0, 0.0, 5.0
111 , 0.0, 0.0, 3.0, 6.0 ]@
112-}
113-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
114a <|> b = joinH a b
115
116-- -- | Vertical concatenation of matrices and vectors.
117-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
118a <-> b = joinV a b
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 5d0154d..288439f 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -44,8 +44,7 @@ module Numeric.LinearAlgebra.LAPACK (
44import Data.Packed.Internal 44import Data.Packed.Internal
45import Data.Packed.Matrix 45import Data.Packed.Matrix
46import Data.Complex 46import Data.Complex
47import Numeric.Vector() 47import Numeric.Conversion
48import Numeric.Container
49import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) 48import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale))
50import Foreign 49import Foreign
51import Foreign.C.Types (CInt) 50import Foreign.C.Types (CInt)
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
deleted file mode 100644
index 775060e..0000000
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ /dev/null
@@ -1,160 +0,0 @@
1{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE TypeFamilies #-}
4-----------------------------------------------------------------------------
5{- |
6Module : Numeric.LinearAlgebra.Linear
7Copyright : (c) Alberto Ruiz 2006-7
8License : GPL-style
9
10Maintainer : Alberto Ruiz (aruiz at um dot es)
11Stability : provisional
12Portability : uses ffi
13
14Basic optimized operations on vectors and matrices.
15
16-}
17-----------------------------------------------------------------------------
18
19module Numeric.LinearAlgebra.Linear (
20 -- * Linear Algebra Typeclasses
21 Vectors(..),
22 -- * Products
23 Product(..),
24 mXm,mXv,vXm,
25 outer, kronecker,
26 -- * Modules
27 --module Numeric.Vector,
28 --module Numeric.Matrix,
29 module Numeric.Container
30) where
31
32import Data.Packed.Internal.Common
33import Data.Packed.Matrix
34import Data.Packed.Vector
35import Data.Complex
36import Numeric.Container
37import Numeric.Vector()
38import Numeric.Matrix()
39import Numeric.GSL.Vector
40import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
41
42-- | Linear algebraic properties of objects
43class Num e => Vectors a e where
44 -- | dot (inner) product
45 dot :: a e -> a e -> e
46 -- | sum of absolute value of elements (differs in complex case from @norm1@
47 absSum :: a e -> e
48 -- | sum of absolute value of elements
49 norm1 :: a e -> e
50 -- | euclidean norm
51 norm2 :: a e -> e
52 -- | element of maximum magnitude
53 normInf :: a e -> e
54
55instance Vectors Vector Float where
56 norm2 = toScalarF Norm2
57 absSum = toScalarF AbsSum
58 dot = dotF
59 norm1 = toScalarF AbsSum
60 normInf = maxElement . vectorMapF Abs
61
62instance Vectors Vector Double where
63 norm2 = toScalarR Norm2
64 absSum = toScalarR AbsSum
65 dot = dotR
66 norm1 = toScalarR AbsSum
67 normInf = maxElement . vectorMapR Abs
68
69instance Vectors Vector (Complex Float) where
70 norm2 = (:+ 0) . toScalarQ Norm2
71 absSum = (:+ 0) . toScalarQ AbsSum
72 dot = dotQ
73 norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs
74 normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs
75
76instance Vectors Vector (Complex Double) where
77 norm2 = (:+ 0) . toScalarC Norm2
78 absSum = (:+ 0) . toScalarC AbsSum
79 dot = dotC
80 norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs
81 normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs
82
83----------------------------------------------------
84
85class Element t => Product t where
86 multiply :: Matrix t -> Matrix t -> Matrix t
87 ctrans :: Matrix t -> Matrix t
88
89instance Product Double where
90 multiply = multiplyR
91 ctrans = trans
92
93instance Product (Complex Double) where
94 multiply = multiplyC
95 ctrans = conj . trans
96
97instance Product Float where
98 multiply = multiplyF
99 ctrans = trans
100
101instance Product (Complex Float) where
102 multiply = multiplyQ
103 ctrans = conj . trans
104
105----------------------------------------------------------
106
107-- synonym for matrix product
108mXm :: Product t => Matrix t -> Matrix t -> Matrix t
109mXm = multiply
110
111-- matrix - vector product
112mXv :: Product t => Matrix t -> Vector t -> Vector t
113mXv m v = flatten $ m `mXm` (asColumn v)
114
115-- vector - matrix product
116vXm :: Product t => Vector t -> Matrix t -> Vector t
117vXm v m = flatten $ (asRow v) `mXm` m
118
119{- | Outer product of two vectors.
120
121@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
122(3><3)
123 [ 5.0, 2.0, 3.0
124 , 10.0, 4.0, 6.0
125 , 15.0, 6.0, 9.0 ]@
126-}
127outer :: (Product t) => Vector t -> Vector t -> Matrix t
128outer u v = asColumn u `multiply` asRow v
129
130{- | Kronecker product of two matrices.
131
132@m1=(2><3)
133 [ 1.0, 2.0, 0.0
134 , 0.0, -1.0, 3.0 ]
135m2=(4><3)
136 [ 1.0, 2.0, 3.0
137 , 4.0, 5.0, 6.0
138 , 7.0, 8.0, 9.0
139 , 10.0, 11.0, 12.0 ]@
140
141@\> kronecker m1 m2
142(8><9)
143 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
144 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
145 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
146 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
147 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
148 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
149 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
150 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
151-}
152kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
153kronecker a b = fromBlocks
154 . splitEvery (cols a)
155 . map (reshape (cols b))
156 . toRows
157 $ flatten a `outer` flatten b
158
159
160-------------------------------------------------------------------