\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}