Sept 2015
This library works with two types of numeric arrays: Matrix and Vector. (General multidimensional arrays are provided by other packages like repa or hTensor.)
The Matrix type is a 2-dimensional array with an internal representation suitable for the standard linear algebra BLAS/LAPACK libraries.
The Vector type
is just the Data.Vector.Storable.Vector from the standard vector package.
These arrays are dense, immutable and strict in all elements (unboxed), and are usually manipulated as whole blocks by high-level, indexless functions.
All functions admit double precision real and complex elements (R=Double and C=Complex Double).
Single precision (Float and Complex Float) and machine integers (I=Int32 and Z=Int64) are also supported in arithmetic operations and matrix products. There is experimental support for modular arithmetic in selected functions. Conversion between arrays of different element is explicit.
The instances of the Haskell numeric classes are defined in the element-by-element sense. The matrix product is (<>), the matrix-vector product is (#>), the dot product is (<.>), and the (conjugate) transpose is tr. The general linear system solver is (<\>).
Most operations are "autoconformable": a single row or column is replicated as needed to match the other operand.
The rest of this page is an introduction to the standard interface, which uses run-time dimension checking. See this page for information about the alternative statically checked interface.
To get a quick idea of library usage we first show a couple of simple examples. The following code solves a small overdetermined linear system.
import Numeric.LinearAlgebra -- coefficient matrix of a linear system a = (4><3) [ 1, 2, 3 , 4, 0, 5 , 7, 7, 2 , 3, 3, 1 ] :: Matrix R -- multiple right-hand sides b = (4><2) [ 10, 1 , 20, 2 , 30, 3 , 15, 1 ] :: Matrix R -- least squares solution x = linearSolveLS a b residual = norm_Frob (a <> x - b)
λ> x
(3><2) [ 2.862124248496993, 0.2739078156312625 , 1.0300601202404813, 8.7374749498998e-2 , 1.7134268537074153, 0.1803607214428857 ]
λ> residual
1.7692315735147401
The obtained right-hand sides:
λ> disp 3 (a <> x)
4x2 10.063 0.990 20.016 1.997 30.672 2.890 13.390 1.264
The residual of each column:
λ> map norm_2 $ toColumns (a <> x - b)
[1.7458797760397966,0.28650334786294096]
The above solver is very fast but the coefficient matrix must have full rank. For rank-deficient systems we must use the operator <\> (resembling the backslash of Matlab/Octave), based on the SVD. Examples of different types of linear systems can be found in this page.
The library provides all common matrix decompositions. For example, we can check the SVD on the coefficient matrix of the above system.
λ> let (u,s,v) = thinSVD a
λ> disp 3 u
4x3 -0.251 0.265 -0.930 -0.390 0.852 0.348 -0.812 -0.423 0.114 -0.353 -0.158 0.015
λ> disp 3 (asRow s)
1x3 12.153 4.991 1.842
λ> disp 3 v
3x3 -0.704 0.048 0.708 -0.596 -0.582 -0.553 -0.385 0.812 -0.438
λ> disp 5 (u <> diag s <> tr v)
4x3 1.00000 2.00000 3.00000 4.00000 0.00000 5.00000 7.00000 7.00000 2.00000 3.00000 3.00000 1.00000
The singular values of $A$ are the eigenvalues of $A^TA$:
λ> singularValues a
fromList [12.153316359624688,4.9905623542834805,1.8415180832366043]
λ> sqrt $ eigenvalues (tr a <> a)
fromList [12.153316359624688 :+ 0.0,1.8415180832366078 :+ 0.0,4.9905623542834805 :+ 0.0]
We have a more efficient eigendecomposition specialized for symmetric matrices:
λ> sqrt $ fst $ eigSH (mTm a)
fromList [12.153316359624691,4.990562354283482,1.8415180832366058]
Detailed information about the available matrix decompositions can be found in the Haddock documentation of the package.
Other common linear algebra functions are readily available:
λ> rank a
3
λ> rcond a
0.15152391567411425
λ> nullspace (tr a)
(4><1) [ -3.581291848286757e-2 , -8.95322962071691e-3 , -0.38498887369082685 , 0.9221826509338411 ]
λ> orth a
(4><3) [ -0.2511962379361302, 0.2646741385060122, -0.930357719010583 , -0.39028183686836193, 0.8522734562693702, 0.34818081973118137 , -0.8124981271267512, -0.42266466628276006, 0.11395148277171996 , -0.35274258500061345, -0.15789913242587644, 1.4821977697794364e-2 ]
λ> det (takeRows 3 a)
103.0
λ> pinv a
(3><4) [ -0.34076152304609214, 0.16480961923847692, 8.681362725450899e-2, 2.460921843687378e-2 , 0.26092184368737475, -0.1847695390781563, 5.490981963927855e-2, 3.126252505010018e-2 , 0.2725450901803607, 6.813627254509018e-2, -7.01402805611222e-2, -1.80360721442886e-2 ]
u = vector [1..5] v = vector [10,-3,0,4,5]
λ> u + v
fromList [11.0,-1.0,3.0,8.0,10.0]
λ> u * v
fromList [10.0,-6.0,0.0,16.0,25.0]
λ> u <.> v
45.0
λ> u - 10
fromList [-9.0,-8.0,-7.0,-6.0,-5.0]
λ> let w = vector [1,0,1,0,1,0,1]
λ> norm_2 w
2.0
λ> unitary w
fromList [0.5,0.0,0.5,0.0,0.5,0.0,0.5]
λ> w / scalar (sqrt (w <.> w))
fromList [0.5,0.0,0.5,0.0,0.5,0.0,0.5]
a = (5><3) [1..] :: Matrix R
b = matrix 5
[ 1, 2, 3, 0, -2
, 0, 7, -2, 4, 1 ]
λ> a
(5><3) [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 , 7.0, 8.0, 9.0 , 10.0, 11.0, 12.0 , 13.0, 14.0, 15.0 ]
λ> b
(2><5) [ 1.0, 2.0, 3.0, 0.0, -2.0 , 0.0, 7.0, -2.0, 4.0, 1.0 ]
λ> tr b
(5><2) [ 1.0, 0.0 , 2.0, 7.0 , 3.0, -2.0 , 0.0, 4.0 , -2.0, 1.0 ]
λ> b <> a
(2><3) [ 4.0, 8.0, 12.0 , 67.0, 77.0, 87.0 ]
λ> b #> u
fromList [4.0,29.0]
λ> a * 2
(5><3) [ 2.0, 4.0, 6.0 , 8.0, 10.0, 12.0 , 14.0, 16.0, 18.0 , 20.0, 22.0, 24.0 , 26.0, 28.0, 30.0 ]
λ> a + row [100,200,300]
(5><3) [ 101.0, 202.0, 303.0 , 104.0, 205.0, 306.0 , 107.0, 208.0, 309.0 , 110.0, 211.0, 312.0 , 113.0, 214.0, 315.0 ]
λ> diagl [1..4]
(4><4) [ 1.0, 0.0, 0.0, 0.0 , 0.0, 2.0, 0.0, 0.0 , 0.0, 0.0, 3.0, 0.0 , 0.0, 0.0, 0.0, 4.0 ]
λ> u <.> diagl [1..5] #> v
187.0
λ> vector [1,2,3] `outer` vector [10,20,30,40]
3x4 10 20 30 40 20 40 60 80 30 60 90 120
λ> diagl [1,20] `kronecker` b
4x10 1 2 3 0 -2 0 0 0 0 0 0 7 -2 4 1 0 0 0 0 0 0 0 0 0 0 20 40 60 0 -40 0 0 0 0 0 0 140 -40 80 20
There are several ways to define a vector. The simplest method is fromList, which admits any desired element type.
λ> fromList [1,2,3::R]
fromList [1.0,2.0,3.0]
λ> fromList [2,3+2*iC,-8]
fromList [2.0 :+ 0.0,3.0 :+ 2.0,(-8.0) :+ (-0.0)]
The function vector used above is specialized for double precision real elements (Vector R). It can be useful to avoid a possible type annotation.
λ> vector [1,2,3]
fromList [1.0,2.0,3.0]
The operator |> explicitly takes the length and admits extra elements which are discarded.
λ> 3 |> [0..] :: Vector Z
fromList [0,1,2]
The functions idxs and range create Vector I which can be used to define matrix slices (explained below).
λ> idxs [5,3,2,1]
fromList [5,3,2,1]
λ> range 4
fromList [0,1,2,3]
More examples:
λ> linspace 5 (0,1) :: Vector Float
fromList [0.0,0.25,0.5,0.75,1.0]
λ> build 5 (^2) :: Vector R
fromList [0.0,1.0,4.0,9.0,16.0]
λ> subVector 2 3 (vector [0..9])
fromList [2.0,3.0,4.0]
λ> vjoin[vector[5,3,8], 12, vector[1..3]]
fromList [5.0,3.0,8.0,12.0,1.0,2.0,3.0]
We can also import qualified Data.Storable.Vector as V to use any desired function from the vector package. Many functions provided by hmatrix admit both vectors and matrices.
A matrix can be created from a list of elements in row order, using the operator (><) to specify the dimensions. This method admits any Element type and extra elements are discarded.
λ> (3><4)[0..] :: Matrix Z
(3><4) [ 0, 1, 2, 3 , 4, 5, 6, 7 , 8, 9, 10, 11 ]
There is also a specialized constructor of Matrix R (double precision real elements), which takes only the number of columns of the matrix. (The number of rows is deduced from the length of the list, which in this case cannot be infinite.)
λ> matrix 2 [1..6]
(3><2) [ 1.0, 2.0 , 3.0, 4.0 , 5.0, 6.0 ]
Other examples:
λ> build (5,2) (\i j -> 2*i*iC + 3*j) :: Matrix C
(5><2) [ 0.0 :+ 0.0, 3.0 :+ 0.0 , 0.0 :+ 2.0, 3.0 :+ 2.0 , 0.0 :+ 4.0, 3.0 :+ 4.0 , 0.0 :+ 6.0, 3.0 :+ 6.0 , 0.0 :+ 8.0, 3.0 :+ 8.0 ]
λ> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix R
(5><5) [ 1.0, 0.0, 0.0, 3.0, 0.0 , 0.0, 6.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 1.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 1.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 1.0 ]
λ> repmat (ident 2) 2 3 :: Matrix I
(4><6) [ 1, 0, 1, 0, 1, 0 , 0, 1, 0, 1, 0, 1 , 1, 0, 1, 0, 1, 0 , 0, 1, 0, 1, 0, 1 ]
See Numeric.LinearAlgebra.Data for more information.
λ> fromBlocks [[konst 7 (1,5), 0, ident 3],[0,konst 1 (3,2),8]]
6x10 7 7 7 7 7 0 0 1 0 0 7 7 7 7 7 0 0 0 1 0 7 7 7 7 7 0 0 0 0 1 0 0 0 0 0 1 1 8 8 8 0 0 0 0 0 1 1 8 8 8 0 0 0 0 0 1 1 8 8 8
λ> fromBlocks [[konst 7 (2,5), 0],[0,1]]
3x6 7 7 7 7 7 0 7 7 7 7 7 0 0 0 0 0 0 1
λ> ident 3 ||| 7 :: Matrix I
(3><4) [ 1, 0, 0, 7 , 0, 1, 0, 7 , 0, 0, 1, 7 ]
λ> konst 5 (3,4) ||| row [1,2,3]
3x7 5 5 5 5 1 2 3 5 5 5 5 1 2 3 5 5 5 5 1 2 3
λ> col [10,20,30] ||| konst 5 (3,4)
3x5 10 5 5 5 5 20 5 5 5 5 30 5 5 5 5
λ> row [1,1,1] ||| row [2,2] === row [3,3] ||| row [4,4,4]
2x5 1 1 1 2 2 3 3 4 4 4
λ> diagBlock [konst 7 (2,3), 1]
3x4 7 7 7 0 7 7 7 0 0 0 0 1
λ> let m = (5><10) [0..] :: Matrix I
λ> m
(5><10) [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 , 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 , 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 , 30, 31, 32, 33, 34, 35, 36, 37, 38, 39 , 40, 41, 42, 43, 44, 45, 46, 47, 48, 49 ]
λ> m ?? (All, Take 3)
(5><3) [ 0, 1, 2 , 10, 11, 12 , 20, 21, 22 , 30, 31, 32 , 40, 41, 42 ]
λ> m ?? (DropLast 2, Range 2 1 4)
(3><3) [ 2, 3, 4 , 12, 13, 14 , 22, 23, 24 ]
We can select arbitrary sets of indexes using Pos (idxs creates a Vector I from a list Ints):
λ> m ?? (Pos (idxs [2,3,0]), Pos (4-2*range 3))
(3><3) [ 24, 22, 20 , 34, 32, 30 , 4, 2, 0 ]
Cyclic access is provided by PosCyc:
λ> m ?? (PosCyc (idxs [8,-2]), All)
(2><10) [ 30, 31, 32, 33, 34, 35, 36, 37, 38, 39 , 30, 31, 32, 33, 34, 35, 36, 37, 38, 39 ]
Many data manipulation functions are defined in terms of (??).
λ> m ? [2,1]
(2><10) [ 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 , 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 ]
λ> m ¿ [5,7]
(5><2) [ 5, 7 , 15, 17 , 25, 27 , 35, 37 , 45, 47 ]
Rectangular slices are taken from the original matrix without data copy:
λ> subMatrix (1,3) (3,2) m
(3><2) [ 13, 14 , 23, 24 , 33, 34 ]
(If we only need a small slice of a big matrix it is better to take a fresh array using ?? with the Pos selector, so that the original matrix can be freed.)
Finally, there is also a more general "2D" element extraction:
λ> let r = (2><3) [2,1,2,0,3,3] :: Matrix I
λ> let c = (2><3) [1,1,3,2,0,1] :: Matrix I
λ> r
(2><3) [ 2, 1, 2 , 0, 3, 3 ]
λ> c
(2><3) [ 1, 1, 3 , 2, 0, 1 ]
λ> remap r c m
(2><3) [ 21, 11, 23 , 2, 30, 31 ]
The matrices containing the indexes for remap are autoconformable:
λ> remap (asColumn (range 3)) (asRow (range 4)) m
(3><4) [ 0, 1, 2, 3 , 10, 11, 12, 13 , 20, 21, 22, 23 ]
λ> v <- flatten <$> randn 1 10
λ> v
-0.368 2.427 -2.534 0.703 1.384 1.330 -0.177 -0.452 0.262 0.823
λ> sortVector v
-2.534 -0.452 -0.368 -0.177 0.262 0.703 0.823 1.330 1.384 2.427
λ> sortIndex v
fromList [2,7,0,6,8,3,9,5,4,1]
λ> m <- randn 4 10
λ> disp 2 m
4x10 0.46 0.32 0.40 -0.89 -0.85 1.04 -1.00 1.05 -2.27 0.73 1.19 1.38 -0.13 0.32 0.01 0.33 -0.49 -1.98 -0.18 -0.56 1.33 0.91 0.88 2.37 -0.38 -1.30 1.01 0.91 -0.79 -1.82 1.02 -1.98 0.15 -1.45 -1.04 -0.46 0.49 -0.56 1.00 -1.30
Sort columns by the second row:
λ> disp 2 $ m ?? (All, Pos $ sortIndex (m!1))
4x10 1.05 0.73 -1.00 -2.27 0.40 -0.85 -0.89 1.04 0.46 0.32 -1.98 -0.56 -0.49 -0.18 -0.13 0.01 0.32 0.33 1.19 1.38 0.91 -1.82 1.01 -0.79 0.88 -0.38 2.37 -1.30 1.33 0.91 -0.56 -1.30 0.49 1.00 0.15 -1.04 -1.45 -0.46 1.02 -1.98
λ> b
(2><5) [ 1.0, 2.0, 3.0, 0.0, -2.0 , 0.0, 7.0, -2.0, 4.0, 1.0 ]
λ> disp 4 (b/3)
2x5 0.3333 0.6667 1.0000 0.0000 -0.6667 0.0000 2.3333 -0.6667 1.3333 0.3333
λ> disp 4 b
2x5 1 2 3 0 -2 0 7 -2 4 1
λ> ident 3 :: Matrix C
(3><3) [ 1.0 :+ 0.0, 0.0 :+ 0.0, 0.0 :+ 0.0 , 0.0 :+ 0.0, 1.0 :+ 0.0, 0.0 :+ 0.0 , 0.0 :+ 0.0, 0.0 :+ 0.0, 1.0 :+ 0.0 ]
λ> putStr . dispcf 2 . diag . fromList $ [1,iC,3-2*iC]
3x3 1 0 0 0 i 0 0 0 3-2i
λ> let z = conv2 (konst 1 (3,3)) (diagl [1..10])
λ> z
(12><12) [ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 1.0, 3.0, 3.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 1.0, 3.0, 6.0, 5.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 2.0, 5.0, 9.0, 7.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 3.0, 7.0, 12.0, 9.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 4.0, 9.0, 15.0, 11.0, 6.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 5.0, 11.0, 18.0, 13.0, 7.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 13.0, 21.0, 15.0, 8.0, 0.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.0, 15.0, 24.0, 17.0, 9.0, 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 17.0, 27.0, 19.0, 10.0 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.0, 19.0, 19.0, 10.0 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 10.0 ]
λ> disp 2 (z/2)
12x12 0.50 0.50 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 1.50 1.50 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 1.50 3.00 2.50 1.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2.50 4.50 3.50 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.50 3.50 6.00 4.50 2.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 4.50 7.50 5.50 3.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.50 5.50 9.00 6.50 3.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 6.50 10.50 7.50 4.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.50 7.50 12.00 8.50 4.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.00 8.50 13.50 9.50 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.50 9.50 9.50 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 5.00 5.00
λ> disp 2 z
12x12 1 1 1 0 0 0 0 0 0 0 0 0 1 3 3 2 0 0 0 0 0 0 0 0 1 3 6 5 3 0 0 0 0 0 0 0 0 2 5 9 7 4 0 0 0 0 0 0 0 0 3 7 12 9 5 0 0 0 0 0 0 0 0 4 9 15 11 6 0 0 0 0 0 0 0 0 5 11 18 13 7 0 0 0 0 0 0 0 0 6 13 21 15 8 0 0 0 0 0 0 0 0 7 15 24 17 9 0 0 0 0 0 0 0 0 8 17 27 19 10 0 0 0 0 0 0 0 0 9 19 19 10 0 0 0 0 0 0 0 0 0 10 10 10
λ> dispDots 1 z
1 1 1 . . . . . . . . . 1 3 3 2 . . . . . . . . 1 3 6 5 3 . . . . . . . . 2 5 9 7 4 . . . . . . . . 3 7 12 9 5 . . . . . . . . 4 9 15 11 6 . . . . . . . . 5 11 18 13 7 . . . . . . . . 6 13 21 15 8 . . . . . . . . 7 15 24 17 9 . . . . . . . . 8 17 27 19 10 . . . . . . . . 9 19 19 10 . . . . . . . . . 10 10 10
λ> dispDots 2 $ z/9
0.11 0.11 0.11 . . . . . . . . . 0.11 0.33 0.33 0.22 . . . . . . . . 0.11 0.33 0.67 0.56 0.33 . . . . . . . . 0.22 0.56 1. 0.78 0.44 . . . . . . . . 0.33 0.78 1.33 1. 0.56 . . . . . . . . 0.44 1. 1.67 1.22 0.67 . . . . . . . . 0.56 1.22 2. 1.44 0.78 . . . . . . . . 0.67 1.44 2.33 1.67 0.89 . . . . . . . . 0.78 1.67 2.67 1.89 1. . . . . . . . . 0.89 1.89 3. 2.11 1.11 . . . . . . . . 1. 2.11 2.11 1.11 . . . . . . . . . 1.11 1.11 1.11
λ> dispBlanks 2 $ z/1000
0.00 0.00 0.00
0.00 0.00 0.00 0.00
0.00 0.00 0.01 0.00 0.00
0.00 0.00 0.01 0.01 0.00
0.00 0.01 0.01 0.01 0.00
0.00 0.01 0.02 0.01 0.01
0.00 0.01 0.02 0.01 0.01
0.01 0.01 0.02 0.02 0.01
0.01 0.02 0.02 0.02 0.01
0.01 0.02 0.03 0.02 0.01
0.01 0.02 0.02 0.01
0.01 0.01 0.01
λ> dispDots 2 (diagl [1..15])
1 . . . . . . . . . . . . . . . 2 . . . . . . . . . . . . . . . 3 . . . . . . . . . . . . . . . 4 . . . . . . . . . . . . . . . 5 . . . . . . . . . . . . . . . 6 . . . . . . . . . . . . . . . 7 . . . . . . . . . . . . . . . 8 . . . . . . . . . . . . . . . 9 . . . . . . . . . . . . . . . 10 . . . . . . . . . . . . . . . 11 . . . . . . . . . . . . . . . 12 . . . . . . . . . . . . . . . 13 . . . . . . . . . . . . . . . 14 . . . . . . . . . . . . . . . 15
λ> dispShort 10 10 1 (matrix 30 [1..30*30])
30x30 1.0 2.0 3.0 4.0 5.0 6.0 7.0 .. 29.0 30.0 31.0 32.0 33.0 34.0 35.0 36.0 37.0 .. 59.0 60.0 61.0 62.0 63.0 64.0 65.0 66.0 67.0 .. 89.0 90.0 91.0 92.0 93.0 94.0 95.0 96.0 97.0 .. 119.0 120.0 121.0 122.0 123.0 124.0 125.0 126.0 127.0 .. 149.0 150.0 151.0 152.0 153.0 154.0 155.0 156.0 157.0 .. 179.0 180.0 181.0 182.0 183.0 184.0 185.0 186.0 187.0 .. 209.0 210.0 : : : : : : : :: : : 841.0 842.0 843.0 844.0 845.0 846.0 847.0 .. 869.0 870.0 871.0 872.0 873.0 874.0 875.0 876.0 877.0 .. 899.0 900.0
λ> dispShort 10 10 2 =<< randn 1000 5
1000x5 1.41 -0.72 0.46 1.62 -1.19 -0.45 1.00 0.71 0.12 -0.77 1.31 -1.22 1.41 -0.16 1.04 -0.54 -2.06 1.21 -1.48 -1.63 -0.34 -0.41 1.36 -0.91 -1.91 -0.53 0.33 0.54 0.40 1.08 -0.35 0.38 0.56 -0.71 0.57 : : : : : -0.86 -0.09 0.46 -0.45 0.62 -0.08 1.06 -2.45 -0.43 1.15
λ> let h = build (5,5) (\r c -> r - 2*iC*(c-2))
λ> h
(5><5) [ 0.0 :+ 4.0, 0.0 :+ 2.0, 0.0 :+ 0.0, 0.0 :+ (-2.0), 0.0 :+ (-4.0) , 1.0 :+ 4.0, 1.0 :+ 2.0, 1.0 :+ 0.0, 1.0 :+ (-2.0), 1.0 :+ (-4.0) , 2.0 :+ 4.0, 2.0 :+ 2.0, 2.0 :+ 0.0, 2.0 :+ (-2.0), 2.0 :+ (-4.0) , 3.0 :+ 4.0, 3.0 :+ 2.0, 3.0 :+ 0.0, 3.0 :+ (-2.0), 3.0 :+ (-4.0) , 4.0 :+ 4.0, 4.0 :+ 2.0, 4.0 :+ 0.0, 4.0 :+ (-2.0), 4.0 :+ (-4.0) ]
λ> putStrLn $ dispcf 2 h
5x5 4i 2i 0 -2i -4i 1+4i 1+2i 1 1-2i 1-4i 2+4i 2+2i 2 2-2i 2-4i 3+4i 3+2i 3 3-2i 3-4i 4+4i 4+2i 4 4-2i 4-4i
λ> let latex = putStrLn . latexFormat "bmatrix" . dispcf 2