lapack
Numerical Linear Algebra using LAPACK
https://hub.darcs.net/thielema/lapack/
Version on this page:  0.5.1 
LTS Haskell 22.33:  0.5.2 
Stackage Nightly 20240909:  0.5.2 
Latest on Hackage:  0.5.2 
lapack0.5.1@sha256:26f71569a13435f9700cfc6498a3841d1505f152382e0d722c3ee110a472f2ed,8349
Module documentation for 0.5.1
 Numeric
 Numeric.LAPACK
 Numeric.LAPACK.Example
 Numeric.LAPACK.Format
 Numeric.LAPACK.Linear
 Numeric.LAPACK.Matrix
 Numeric.LAPACK.Matrix.Array
 Numeric.LAPACK.Matrix.Banded
 Numeric.LAPACK.Matrix.BandedHermitian
 Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite
 Numeric.LAPACK.Matrix.Block
 Numeric.LAPACK.Matrix.Diagonal
 Numeric.LAPACK.Matrix.Extent
 Numeric.LAPACK.Matrix.Full
 Numeric.LAPACK.Matrix.Function
 Numeric.LAPACK.Matrix.Hermitian
 Numeric.LAPACK.Matrix.HermitianPositiveDefinite
 Numeric.LAPACK.Matrix.Layout
 Numeric.LAPACK.Matrix.Permutation
 Numeric.LAPACK.Matrix.Quadratic
 Numeric.LAPACK.Matrix.Shape
 Numeric.LAPACK.Matrix.Special
 Numeric.LAPACK.Matrix.Square
 Numeric.LAPACK.Matrix.Superscript
 Numeric.LAPACK.Matrix.Symmetric
 Numeric.LAPACK.Matrix.Triangular
 Numeric.LAPACK.Matrix.Type
 Numeric.LAPACK.Orthogonal
 Numeric.LAPACK.Output
 Numeric.LAPACK.Permutation
 Numeric.LAPACK.Scalar
 Numeric.LAPACK.Shape
 Numeric.LAPACK.Singular
 Numeric.LAPACK.Vector
 Numeric.LAPACK
Installation
Before installing the Haskell bindings you need to install the BLAS and LAPACK packages. Please note, that additionally to the reference implementation in FORTRAN 77 there are alternative optimized implementations like OpenBLAS, ATLAS, Intel MKL.
Debian, Ubuntu
sudo aptget install libblasdev liblapackdev
MacOS
You may install pkgconfig and LAPACK via MacPorts:
sudo port install pkgconfig lapack
However, the pkgconfig files for LAPACK seem to be installed in a nonstandard location. You must make them visible to pkgconfig by
export PKG_CONFIG_PATH=/opt/local/lib/lapack/pkgconfig:$PKG_CONFIG_PATH
You may set the search PATH permanently by adding
the above command line to your $HOME/.profile
file.
Alternatively, a solution for all users of your machine would be to set symbolic links:
sudo ln s /opt/local/lib/lapack/pkgconfig/blas.pc /opt/local/lib/pkgconfig/blas.pc
sudo ln s /opt/local/lib/lapack/pkgconfig/lapack.pc /opt/local/lib/pkgconfig/lapack.pc
Introduction
Here is a small example for constructing and formatting matrices:
Prelude> import qualified Numeric.LAPACK.Matrix as Matrix
Prelude Matrix> import Numeric.LAPACK.Format as Fmt ((##))
Prelude Matrix Fmt> let a = Matrix.fromList (Matrix.shapeInt 3) (Matrix.shapeInt 4) [(0::Float)..]
Prelude Matrix Fmt> a ## "%.4f"
0.0000 1.0000 2.0000 3.0000
4.0000 5.0000 6.0000 7.0000
8.0000 9.0000 10.0000 11.0000
Prelude Matrix Fmt> import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
Prelude Matrix Fmt MatrixShape> import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
Prelude Matrix Fmt MatrixShape Triangular> let u = Triangular.upperFromList MatrixShape.RowMajor (Matrix.shapeInt 4) [(0::Float)..]
Prelude Matrix Fmt MatrixShape Triangular> (u, Triangular.transpose u) ## "%.4f"
0.0000 1.0000 2.0000 3.0000
4.0000 5.0000 6.0000
7.0000 8.0000
9.0000
0.0000
1.0000 4.0000
2.0000 5.0000 7.0000
3.0000 6.0000 8.0000 9.0000
You may find a more complex introductory example at: http://code.henningthielemann.de/bob2019/main.pdf
Formatting
We do not try to do fancy formatting for the Show
instance.
The Show
instances of matrices emit plain valid Haskell code in one line
where all numbers are printed in full precision.
If matrices are part of larger Haskell data structures
this kind of formatting works best.
For humanfriendly formatting in GHCi you need to add something like ## "%.4f"
after a matrix or vector expression.
It means: Print all numbers in fixed point representation
using four digits for the fractional part.
You can use the formatting placeholders provided by printf
.
The matrices have Hyper
instances for easy usage in
HyperHaskell.
Formatting is based on the Output
type class
that currently supports output as
Text boxes
for GHCi and
HTML
for HyperHaskell.
You may tell GHCi to use Format
class instead of Show
:
Fmt> let lapackPrint x = x##"%.3f"
Fmt> :set interactiveprint lapackPrint
You may permanently configure this one in your local .ghci
file.
If you want to display values via Show
class,
you can always fall back by:
Fmt> print "Hello"
Matrix vs. Vector
Vectors are Storable.Array
s from the
comfortarray package.
An array can have a fancy shape
like a shape defined by an enumeration type,
the shape of two appended arrays,
a shape that is compatible to a Haskell container type,
a rectangular or triangular shape.
All operations check dynamically
whether corresponding shapes match structurally.
E.g. ab === cd
composes a matrix from four quadrants a
, b
, c
, d
.
It is not enough that ab
and cd
have the same width,
but the widths of a
and c
as well as b
and d
must match.
The type variables for shapes show which dimensions must be compatible.
We recommend to use type variables for the shapes as long as possible
because this increases type safety even
if you eventually use only ShapeInt
.
If you use statically sized shapes you get static size checks.
A matrix can have any of these shapes as height and as width.
E.g. it is no problem to define a matrix
that maps a triangular shaped array to a rectangular shaped one.
There are actually practical applications to such matrices.
A matrix can be treated as vector, but there are limitations.
E.g. if you scale a Hermitian matrix by a complex factor
it will in general be no longer Hermitian.
Another problem: Two equally sized rectangular matrices
may differ in the element order (row major vs. column major).
You cannot simply add them by adding the flattened arrays elementwise.
Thus if you want to perform vector operations on a matrix
the package requires you to “unpack” a matrix to a vector
using Matrix.Array.toVector
.
This conversion is almost a noop and preserves most of the shape information.
The reverse operation is Matrix.Array.fromVector
.
There are more matrix types that are not based on a single array. E.g. we provide a symbolic inverse, a scaling matrix, a permutation matrix. We also support arrays that represent factors of a matrix factorization. You obtain these by LU and QR decompositions. You can extract the matrix factors of it, but you can also multiply the factors to other matrices or use the decompositions for solving simultaneous linear equations.
Matrix type parameters
LAPACK supports a variety of special matrix types,
e.g. dense, banded, triangular, symmetric, Hermitian matrices
and our Haskell interface supports them, too.
There are two layers:
The low level layer addresses how matrices are stored for LAPACK.
Matrices and vectors are stored in the Array
type
from comfortarray:Data.Array.Comfort.Storable
using shape types from Numeric.LAPACK.Matrix.Layout
.
The high level layer provides a matrix type
in Numeric.LAPACK.Matrix.Array
with mathematically relevant type parameters.
The matrix type is:
ArrayMatrix pack property lower upper meas vert horiz height width a
The type parameters are from right to left:

a
 the element type 
height
andwidth
are the vertical and horizontal shapes of the matrix 
meas vert horiz
form a group with following possible assignments:meas vert horiz name condition Shape
Small
Small
Square matrix height == width Size
Small
Small
Liberal square size height == size width Size
Big
Small
Tall matrix size height >= size width Size
Small
Big
Wide matrix size height <= size width Size
Big
Big
General matrix arbitrary height and width Think of
meas
as the measurement that goes into the relation of dimensions. You can either compare shapes (meas ~ Shape
) or their sizes (meas ~ Size
). Forvert
andhoriz
the possible values mean:
Small
: The corresponding dimension is equal to the minimum ofheight
andwidth
. 
Big
: The corresponding dimension has no further restrictions, but it is of course at least the minimum ofheight
andwidth
.
The names
Small
andBig
fit best to tall and wide matrices. The remaining combinationsSmall Small
forSquare
andBig Big
forGeneral
appear to be arbitrary, but they help to e.g. treat square and tall matrices the same way, where sensible. TurningShape
intoSize
orSmall
intoBig
relaxes a dimension relation. 

lower upper
count the numbers of nonzero offdiagonals.Off course, stored offdiagonals can consist entirely of zeros. Thus more precisely we have to say, that
lower
andupper
tell that all values outside the numbered bands are zero.lower
andupper
can be:
Filled
 no restriction on the number of offdiagonals. 
Bands n
, wheren
is a natural number unarily encoded in types.Empty
is a synonym forBands U0
.


property
can be
Arbitrary
 this type does not make any further promises about the matrix elements 
Unit
 matrix is triangular with a unit diagonalIt can be used for matrices that always have a unit diagonal by construction. This property is preserved by some operations and enables optimizations by LAPACK. Solving with a unit triangular matrix does not require division and thus cannot fail due to division by zero.

Symmetric
 matrix is symmetric 
Hermitian
 matrix is Hermitian (also supported for real elements)The internal
Hermitian
property also has three type tagsneg zero pos
to restrict the range of values of bilinearforms. This way you can denote positive definiteness and semidefiniteness.


pack
canPacked
orUnpacked
.Unpacked
means that the full matrix bounded byheight
andwidth
is stored.Packed
format is supported for triangular, symmetric, Hermitian and banded matrices.For banded matrices you should always prefer the packed format. For triangular, symmetric and Hermitian matrices LAPACK does not always support packed format natively and our Haskell interface converts forth and back silently. I also think that unpacked triangular formats enjoy better support by vectorized block algorithms. Thus, the unpacked triangular format may be better for performance although it requires twice as much space as the packed format. The packed triangular format however is still the default format for conversion to and from lists, because this prevents the user from declaring nonzero values in the empty area of a triangular matrix.
Let us examine some examples:

Diagonal matrix:
ArrayMatrix Packed Arbitrary Empty Empty Shape Small Small sh sh a

Packed upper triangular matrix:
ArrayMatrix Packed Arbitrary Empty Filled Shape Small Small sh sh a

Unpacked unit lower triangular matrix:
ArrayMatrix Unpacked Unit Filled Empty Shape Small Small sh sh a

Complexvalued symmetric matrix:
ArrayMatrix Packed Symmetric Filled Filled Shape Small Small sh sh (Complex a)

Tall banded matrix:
ArrayMatrix Packed Arbitrary (Bands sub) (Bands super) Size Big Small height width a
The type tags have a mathematical meaning and this pays off for operations on matrices. E.g. matrix multiplication adds offdiagonals. Matrix inversion fills nonzero triangular matrix parts. The five supported relations for matrix dimensions are transitive, and thus matrix multiplication maintains a dimension relation, e.g. tall times tall is tall.
Please note, that not all type parameter combinations are supported.
Some restrictions are dictated by mathematics,
e.g. Hermitian matrices must always be square,
matrices with unit diagonal must always be triangular and so on.
Some combinations are simply not supported, because they do not add value.
E.g. a (square) diagonal matrix is always symmetric
but we allow Symmetric
only together with Filled
.
Forbidden combinations are often not prevented at the type level,
but you will not be able to construct a matrix of a forbidden type.
Infix operators
The package provides fancy infix operators like #*
and \*#
.
They symbolize both operands and operations.
E.g. in #*
the hash means Matrix, the star means Multiplication
and the bar means Column Vector.
Possible operations are:

a_*_b
a
multiplied byb

a_/_b
a
multiplied byinverse b

a_\_b
inverse a
multiplied byb
Possible operands are:

#
 a matrix that is generalized through a type class 
##
 a full matrix 
\
 a diagonal matrix represented by aVector


 a row vector 

 a column vector 
.
 a scalar
For multiplication of equally shaped matrices
we also provide instances of Semigroup.<>
.
Precedence of the operators is chosen analogously to plain *
and /
.
Associativity is chosen such that the same operator can be applied
multiple times without parentheses.
But sometimes this may mean that you have to mix
left and right associative operators,
and thus you may still need parentheses.
Type errors
You might encounter cryptic type errors that refer to the encoding of particular matrix types via matrix type parameters.
E.g. the error
Couldn't match type `Numeric.LAPACK.Matrix.Extent.Big`
with `Numeric.LAPACK.Matrix.Extent.Small`
may mean that you passed Square
where General
or Tall
was expected.
You may solve the problem with a function
like Square.toFull
or Square.fromFull
.
The error
Couldn't match type `Type.Data.Bool.False`
with `Type.Data.Bool.True`
most likely refers to nonmatching definiteness warranties
in a Hermitian
matrix.
You may try a function like Hermitian.assureFullRank
or Hermitian.relaxIndefinite
to fix the issue.
Changes
Change log for the lapack
package
0.5

Matrix.Block
for Block matrices. Add*Extra
constraint families to many type classes in order to handle the data stored in the extra type parameters ofMatrix
. 
Format.format
now uses custom typeConfig
instead of a plain formatString
.
0.4

Unified
Matrix
type that provides the same type parameters across all special types. This reduces the use of type functions and improves type inference. 
Unified
transpose
andadjoint
functions enabled by the newMatrix
type. 
Unpacked
format: We now support data type and according functions for unpacked triangular, symmetric and Hermitian matrices. Enables declaration e.g. of Hessenberg matrices. 
There are now two types of square matrices:

Square
: height and width shapes match exactly 
LiberalSquare
: only the sizes of height and width match


Hermitian
: Definiteness properties in the type 
eigensystem
,Householder.fromMatrix
,LowerUpper.fromMatrix
etc.: We use the new classShape.Permutable
for shapes where permutation of indices seems to make sense. We tried using liberal squares matrix factors, but this would require extra parameters and consistency checks for the shapes of the factor matrices. 
Square.fromGeneral
>fromFull

Orthogonal.affineKernelFromSpan
>affineFiberFromFrame
,Orthogonal.affineSpanFromKernel
>affineFrameFromFiber

Matrix.Function
: New module providing generalized algebraic and transcendent functions likesqrt
,exp
,log
. 
Matrix.Superscript
: Experimental module for eyecandy notationa#^T
for transposition anda#^Inv
for inverse.
0.3.2

Orthogonal
:project
,affineKernelFromSpan
,affineSpanFromKernel
,leastSquaresConstraint
,gaussMarkovLinearModel

Symmetric.fromHermitian
,Hermitian.fromSymmetric

instance Monoid Matrix
, especiallymempty
for matrices with static shapes. 
Extent.Dimensions
: turn from type family to data family 
Start using
doctestextract
for simple tests
0.3.1
Matrix.Symmetric
: You can now import many functions for symmetric matrices from this module. This is more natural than importing them fromTriangular
.
0.3

Matrix data family

Matrix
:ZeroInt
>ShapeInt
,zeroInt
>shapeInt

Hermitian
,BandedHermitian
:covariance
>gramian

Square.eigensystem
: Return left eigenvectors as rows of the last matrix. This is adjoint with respect to the definition inlapack0.2
but it is consistent with the other eigenvalue and singular value decompositions.