Jul 24, 2015
In principle, any matrix computation can be defined using the available tools for manipulation of immutable data; many algorithms can be expressed as recursive procedures working on lists of row or column vectors or matrix blocks. There is some allocation overhead and we may need to force evaluation of some intermediate results to avoid space leaks, but this is usually not a big problem.
However, in certain situations more traditional destructive operations without memory allocation can be very advantageous. This library also provides the imperative building blocks typically used in the literature to describe matrix algorithms.
The development module exports tools for inplace computations in the ST monad, which allow to define pure functions using imperative style. The most recent version includes new high level operations (matrix products, elementary row operations, etc.) for all Element types.
mutable (\(r,c) x -> setMatrix x 2 1 (col [10,20,30])) (ident 5)
((5><5) [ 1.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 1.0, 0.0, 0.0, 0.0 , 0.0, 10.0, 1.0, 0.0, 0.0 , 0.0, 20.0, 0.0, 1.0, 0.0 , 0.0, 30.0, 0.0, 0.0, 1.0 ],())
mutable (\(r,c) x -> rowOper (SWAP 0 4 AllCols) x) (ident 5 :: Matrix Z)
((5><5) [ 0, 0, 0, 0, 1 , 0, 1, 0, 0, 0 , 0, 0, 1, 0, 0 , 0, 0, 0, 1, 0 , 1, 0, 0, 0, 0 ],())
mutable (\(r,c) x -> rowOper (AXPY 3 0 1 (FromCol 2)) x) (konst 1 (5,7) :: Matrix I)
((5><7) [ 1, 1, 1, 1, 1, 1, 1 , 1, 1, 4, 4, 4, 4, 4 , 1, 1, 1, 1, 1, 1, 1 , 1, 1, 1, 1, 1, 1, 1 , 1, 1, 1, 1, 1, 1, 1 ],())
The mutable helper encapsulates the thawing and freezing auxiliary operations and allows
a possible additional result (not used here).
These tools can be used to define linear algebra computations with element types not supported by LAPACK. For example, consider the forward substitution step of a solver based on the LU decomposition. (This is just Gaussian elimination working with a unit lower triangular matrix.) The standard imperative algorithm destructively replaces the right-hand side by the solution, and can be directly implemented using inplace matrix products of rectangular slices:
forwSubs lup rhs = fst $ mutable f rhs
where
f (r,c) x = do
l <- thawMatrix lup
let go k = gemmm 1 (Slice x k 0 1 c) (-1) (Slice l k 0 1 k) (Slice x 0 0 k c)
mapM_ go [0..r-1]
We test this code with a small triangular system $AX=B$ with real elements:
a
(3><3) [ 1.0, 0.0, 0.0 , 2.0, 1.0, 0.0 , 3.0, 4.0, 1.0 ]
b
(3><2) [ 1.0, 2.0 , 5.0, 8.0 , 20.0, 28.0 ]
forwSubs a b
(3><2) [ 1.0, 2.0 , 3.0, 4.0 , 5.0, 6.0 ]
While the standard solvers (linearSolve, <\>, etc.) only admit Double or Complex Double elements, forwSubs works with any Element type. Division is not needed to solve these simple unit triangular systems, so we can even accept integer elements. For example, we can solve the above test case with elements in Z=Int64 (note that toZ performs truncation):
let a' = toZ a; b' = toZ b
forwSubs a' b'
(3><2) [ 1, 2 , 3, 4 , 5, 6 ]
The function forwSubs is actually used in a complete experimental solver written with this kind of inplace operations. It is of course slower than the optimized LAPACK routines with Double elements, but it can solve modular systems of dimension around 1000 in a few seconds.
This page shows some examples of modular linear systems.