\section{Densematrix}
============================================================+
================== A MATRIX IMPLEMENTATION ==================
============================================================+
This is an implementation of the types "matrix" and "vector"
and their associated operations. It is most appropriate for
dense matrices and vectors.
The operations involved are:
mmult matvecmult vmmult Multiplication
madd vadd Addition
msub vsub Subtraction
vdot vouter Vector outer and dot products
minverse the inverse of a matrix.
msize vsize Find size of data
{(#rows,#cols) for matrices}
norm Norm of a vector (vdot v v)
swaprow swapcol swaprow a b m ==> swap ath and bth rows
of matrix m
droprow drop the first row of a matrix.
getrow getcol getrow i m ==> get the ith row of a matrix
svmult scalar-vector multiplication
update vupdate update (i,j) x m
vupdate i x v ==> update a singe element
vhd vtl mergevectors treat vectors like [Float]
mkvec mkrvec mkcvec vector creation from [Float].
mkmat mkrmat mkcmat matrix creation from [[Float]]. mkcmat
assumes the arg is column major.
\begin{code}
module Densematrix(Matrix,Vector,
mmult, madd, msub, vouter, vdot, norm ,mneg ,mxpose,mident, msize,
mkmat, mkrmat, mkcmat, mkvec, mkrvec, mkcvec,
vadd, vsub, vsize, vneg,
swaprow, swapcol ,droprow, getrow, getcol,
subscript, vsubscript, vecpart,
update ,vupdate, update2,
vhd, vtl,mergevectors,
matvecmult, vmmult,svmult,
showmatrix,displayvector,
minverse,
veclist, matlist ) where
import List (transpose)
import Utils
type Matrix = [[Float]]
type Vector = [Float]
mmult :: Matrix -> Matrix -> Matrix
madd :: Matrix -> Matrix -> Matrix
msub :: Matrix -> Matrix -> Matrix
vadd :: Vector -> Vector -> Vector
vsub :: Vector -> Vector -> Vector
vsize :: Vector -> Int
vouter :: Vector -> Vector -> Matrix
vdot :: Vector -> Vector -> Float
norm :: Vector -> Float
vneg :: Vector -> Vector
mneg :: Matrix -> Matrix
mxpose :: Matrix -> Matrix
mident :: Int -> Matrix
msize :: Matrix -> (Int,Int)
mkmat :: [[Float]] -> Matrix
mkrmat :: [[Float]] -> Matrix
mkcmat :: [[Float]] -> Matrix
mkvec :: [Float] -> Vector
mkrvec :: [Float] -> Vector
mkcvec :: [Float] -> Vector
veclist :: Vector -> [Float]
matlist :: Matrix -> [[Float]]
swaprow :: Int -> Int -> Matrix -> Matrix
swapcol :: Int -> Int -> Matrix -> Matrix
droprow :: Matrix -> Matrix
getrow :: Int -> Matrix -> Vector
getcol :: Int -> Matrix -> Vector
subscript :: Matrix -> (Int,Int) -> Float
vsubscript :: Vector -> Int -> Float
vecpart :: Int -> Int -> Vector -> Vector
update :: Matrix -> (Int,Int) -> Float -> Matrix
update2 :: Matrix -> (Int,Int,Int) -> (Float,Float) -> Matrix
vupdate :: Vector -> Int -> Float -> Vector
vhd :: Vector -> Float
vtl :: Vector -> Vector
mergevectors :: [Vector] -> Vector
minverse :: Matrix -> Matrix
lu_decomp :: Matrix -> Lu_factor
apply_factor :: Lu_factor -> Matrix -> Matrix
forward :: [[Float]] -> [Int] -> Matrix -> Matrix
backward :: [[Float]] -> [Int] -> Matrix -> Matrix
showmatrix :: Matrix -> [Char]
displayvector :: Vector -> [Char]
matvecmult :: Matrix -> Vector -> Vector
vmmult :: Vector -> Matrix -> Vector
svmult :: Float -> Vector -> Vector
\end{code}
Notes:
vmmult implements an inner product vector matrix
multiplication.
\begin{code}
mmult m n = if compatible then [row i | i <- [0..((length m)-1)]]
else error errmsg
where
compatible = (snd (msize m)) == (fst (msize n))
row i = [item i j | j <- [0..length(head n)-1]]
item i j = vdot (m!!i) (map (!!j) n)
errmsg = "mmult in densematrix: incompatible matrices " ++
(show (msize m)) ++ "*" ++ (show (msize n))
matvecmult m v = map (vdot v) m
vmmult v m = res
where
res' = mmult [v] m
res = concat res'
svmult s v = map (*s) v
vdot u v = sum (map2 (*) u v)
norm v = vdot v v
vneg = map negate
vouter u v = [ map (k*) u | k <- v]
madd a b = map2 vadd a b
vadd u v = map2 (+) u v
msub a b = map2 vsub a b
vsub u v = map2 (-) u v
vsize = (length)
mneg rs = map (map negate) rs
mxpose = transpose
mident n
= [ row k | k <- [1..n] ]
where
row i = (rep (i-1) 0) ++ [1] ++ (rep (n-i) 0)
msize m = (length m , length(head m))
mkmat = id
mkrmat = id
mkcmat = transpose
matlist = id
mkvec = id
mkrvec = id
mkcvec = id
veclist = id
swaprow a b m = swapitems a b m
swapcol a b m = map (swapitems a b) m
swapitems :: Int -> Int -> [a] -> [a]
swapitems a b xs
= if a< b then toa ++ [xs!!b] ++ atob ++ [xs!!a] ++ pastb
else if a > b then swapitems b a xs else xs
where
toa = take a xs
atob = take (b-a-1) (drop (a+1) xs)
pastb = drop (b+1) xs
droprow = tail
getrow i m = m!!i
getcol i m = map (!!i) m
subscript m (i,j) = m !! i !! j
vsubscript v n = v!!n
vecpart start intber
= (take intber) . (drop start)
update m (i,j) val
= (take i m) ++ [(f (m!!i))] ++ (drop (i+1) m)
where
f xs = (take j xs) ++ [val] ++ (drop (j+1) xs)
update2 m (i1,i2,j) (val1,val2)
= (take i1 m) ++ [(f1 (m!!i1))] ++ [(f2 (m!!i2))] ++ (drop (i1+2) m)
where
f1 xs = (take j xs) ++ [val1] ++ (drop (j+1) xs)
f2 xs = (take j xs) ++ [val2] ++ (drop (j+1) xs)
vupdate v i val
= (take i v) ++ [val] ++ (drop (i+1) v)
vhd = head
vtl = tail
mergevectors = concat
showmatrix m = "\n" ++ showmat m ++ "\n"
showmat :: [Vector] -> [Char]
showmat m
= concat (map show_row m)
where
show_row v = (displayvector v) ++ "\n"
displayvector
= concat .
(map (rjustify 13)) .
(map show)
\end{code}
============================================================
============== LU DECOMPOSITION ==============
============== (to do inverses) ==============
============================================================
An lu_factor is a rather contorted representation of an LU
factorization. The first list of vectors represents the L
matrix, the second represents the U part. The third list
contains information about how the matrix was permuted as
it was decomposed.
About the only thing one can do with an Lu_factor is to
use it as an argument to apply_factor.
Notes:
minverse returns the inverse of a matrix.
apply_factor takes one of these lu_factors and applies it
to a matrix. This has the effect of multiplying the inverse
of the factorized matrix by the given matrix.
forward and backward do the forward and backward substitution
on the pieces of an lu_factor.
lu_decomp a
returns ( l-vectors, u-vectors, p-vector )
A baroque representation of an LU factorization.
pivot m pivots the matrix m. It returns the pivoted matrix and
the row on which the pivot was perfomed.
pivoting preserves diagonal elements. two rows and
two columns are swapped.
findpivot m finds the row/column to pivot on in (square) matrix m.
the pivot row/col is the one with the element of greatest
absolulte value on the diagonal.
\begin{code}
data Lu_factor = Lu_fact [[Float]] [[Float]] [Int]
minverse m
= apply_factor (lu_decomp m) (mident (fst(msize m)))
apply_factor factor m
= ((backward uvecs ps) . (forward lvecs ps)) m
where
(Lu_fact lvecs uvecs ps) = factor
forward [] p b = b
forward (l:ls) (p:ps) b
= y1 : (forward ls ps y')
where
y' = msub y (vouter y1 l)
(y1:y) = swaprow 0 p b
backward [[u]] p y = map (map (/u)) y
backward ((u1:u):us) (p:ps) (y:ys)
= x
where
x = swaprow 0 p (xk:nextx)
xk = map (/u1) (vsub y (vmmult u nextx))
nextx = backward us ps ys
lu_decomp a
= if (fst(msize a) == 1) then Lu_fact [] [[only_one]] [1]
else Lu_fact (l:ls) (u:us) (p:ps)
where
Lu_fact ls us ps = lu_decomp a11'
(a',p) = pivot a
u11 = head u
u = getrow 0 a'
m = getcol 0 a'
l = if not(u11 == 0) then tail (map (/ u11) m)
else error ("div by 0 in lu: block is\n "++(show a))
a11 = map tail (tail a')
a11' = msub a11 (vouter (tail u) l)
only_one = subscript a (0,0)
pivot :: Matrix -> (Matrix,Int)
pivot m
= (swaprow 0 p (swapcol 0 p m), p)
where
p = findpivot m
findpivot :: Matrix -> Int
findpivot a
= loc_of_max 0 diag
where
diag = [a `subscript` (i,i) | i <- [0..(fst(msize a))-1] ]
absmax = maxlist (map absfloat diag)
loc_of_max n (x:xs)
= if (absfloat x) == absmax then n
else loc_of_max (n+1) xs
absfloat :: Float -> Float
absfloat n = if n < 0 then -n
else n
maxlist :: [Float] -> Float
maxlist xs = foldl1 max xs
\end{code}
|
