Gauss Elimination
This is -*- Literate Haskell -*-
\begin{code}
module MathObj.Matrix.Gauss where
import MathObj.Matrix.MatArr
import Data.Array ((!), (//), bounds, Ix, listArray)
import MathObj.Matrix.Residue
import MathObj.Generic.VulgarFraction (Fraction, maxFarey)
#define __EDEN__ 1
#ifdef __EDEN__
import Parallel.Eden
#endif
import MathObj.Generic.VulgarFraction ((/:/))
import Debug.Trace (trace)
\end{code}
Straightforward Gauss:
n=length(A)
for i=1:n
for j=i+1:n
y=A(j, i)/A(i, i);
A(j, i)=y;
A(j, i+1:n)=A(j, i+1:n)-y*A(i, i+1:n);
end
end
\begin{code}
gaussIterator (i, bound, matrix) = let ((ln, lm), (n, m)) = bound
zs = [ ((j, k), makeEl i j k)
| j<-[i+1..n], k<-[i..n] ]
makeY (i, j) = matrix!(j, i)/matrix!(i, i)
makeEl i j k = matrix!(j, k)makeY(i, j)*matrix!(i,k)
result = (i+1, bound, matrix//zs)
in result
gauss :: (Num x, Fractional x) => MatArr Int x -> MatArr Int x
gauss matrix = let ((ln, lm), (n, m)) = bounds matrix
(_, _, result) = last $ take (nln+1)
$ iterate gaussIterator (ln, ((ln, lm), (n, m)), matrix)
in result
\end{code}
\begin{code}
diag matrix = let ((ln, lm), (n, m)) = bounds matrix
in [ matrix!(i,i) | i<-[ln..n] ]
\end{code}
\begin{seqcode}
det :: (Num x, Fractional x) => MatArr Int x -> x
det = product . diag . gauss
\end{seqcode}
Now for something completly different!
\begin{seqcode}
gaussResidue = lift1pm gauss
detResidue ps x = let y = gaussResidue ps x
in product $ map unmaybe $ diag $ y
\end{seqcode}
Parallele Version.
\begin{code}
unmaybe (Just x) = x
unmaybe Nothing = error "Nothing an der falschen Stelle!"
gaussResidueG myMap = lift1gpm myMap gauss
gaussResidueGC c myMap = lift1gpmc c myMap gauss
gaussResidueGCD c myMap = lift1tss c myMap gauss
gaussResidueGCL c myMap = lift1tsl c myMap gauss
gaussResidueGCD2 c myMap = lift1tss2 c myMap gauss
detResidueG myMap ps x = let y = trace ("Calling Gauss...") (gaussResidueG myMap ps x)
dy = diag y
(a, b) = countJust dy 0 0
countJust [] a b = (a, b)
countJust ((Just _):xs) a b = countJust xs (a+1) b
countJust (Nothing:xs) a b = countJust xs a (b+1)
d = trace ("Diag: " ++ (show a) ++ " Justs, "
++ (show b) ++ " Nothings") $
d1
d1 = map unmaybe $ dy
res = rnf d `seq` product d
bound = maximum $ map maxFarey d
in trace ("Computed bound: " ++ (show bound))
$ res
\end{code}
\begin{type}
type Map a = (a -> a) -> [a] -> [a]
detResidueGC :: (Trans a, Trans x, Num x) =>
Int -> Map a -> [x] -> MatArr Int x -> x
\end{type}
\begin{code}
detResidueGC c myMap ps x
= let dy = trace "Calling chunking Gauss..." $
gaussResidueGC c myMap ps x
d = map unmaybe $ diag dy
res = product d
in res
\end{code}
\begin{codex}
detResidueGCD c myMap ps x
= let dy = trace "Calling diagonal Gauss..." $
gaussResidueGCD c myMap ps x -- returns diagonals
d1 = map unmaybe $ dy
res = rnf d `seq` product d
in res
\end{codex}
\begin{type}
detResidueGCD :: (Trans a, Trans b, Trans n, Trans x, Num n, Num x) =>
Int -> Map a b -> [n] -> MatArr Int x -> x
\end{type}
\begin{code}
detResidueGCD c myMap ps x
= let ys = trace "Calling Gauss..." $
gaussResidueGCD c myMap ps x
ds = map unmaybe ys
res = rnf ds `seq`
product ds
in res
detResidueGCL c myMap ps x
= let ys = trace "Calling Gauss..." $
gaussResidueGCL c myMap ps x
ds = map unmaybe ys
res = rnf ds `seq`
product ds
in res
detResidueGCD2 c myMap ps x
= let ys = trace "Calling Gauss..." $
gaussResidueGCD2 c myMap ps x
ds = map unmaybe ys
res = rnf ds `seq`
product ds
in res
\end{code}