I would like to put two rather prosaic things into Haskell 1.3. They almost
fall into the "syntactic sugar" class, but they would make my life easier.
The first is that I would like to see arrays be a class instead of whatever
they are. I wanted to construct a subclass of arrays that were constrained
to have lower bounds equal to one, but after fooling around for some time I
just gave up. Maybe it's easy, and I just don't know the right way to hold
my mouth. I would also like to be able to construct a sub-class of one-
dimensional array that is a vector, and a sub-class of two-dimensional
array that is a matrix, and overload "*" to mean "inner product".
The second thing I would like is an array section notation. In many operations
of linear algebra, one needs to view a matrix sometimes as an array of row
vectors, and sometimes as an array of column vectors. This arose in development
of a function that implements Crout's method to factor a matrix. (Crout's
method is especially attractive for functional languages because each element
of the factor is written exactly once. That is not the case with Gauss-like
methods.) I ended up writing three functions, one that computes the inner
product of two vectors, another that computes the inner product of a row and
column of a single matrix, and another that computes the inner product of a
row of one matrix with a column of another. Others would need functions that
compute the inner product of a vector with a row or column of a matrix. It
would be easier to write one function that computes the inner product of two
vectors, and create vectors out of pieces of a matrix by using a section
notation. For example, I might write a Crout reduction with no pivoting as:
lu = array b
([(i,1) := a!(i,1) | i <- [1..m]] ++
[(1,j) := a!(1,j)/a!(1,1) | j <- [2..n]] ++
[(i,j) := (a!(i,j) - dot lu!(i,1..j-1) lu!(1..j-1,j))
| i <- [2..m], j <- [2..i]] ++
[(i,j) := (a!(i,j) - dot lu!(i,1..i-1) lu!(1..i-1,j)) / lu!(i,i)
| i <- [2..m], j <- [i+1,m]])
where ((_,_),(m,n)) = b = bounds a.
BTW, I have developed a Crout reduction that uses pivoting, but I _think_
it's hitting something that's a little too strict -- the run-time system
insists there's a black hole, but if I run the code "by hand" I'm always able
to find an order such that data are available -- there aren't any circular
dependencies on un-computed data. Maybe somebody can tell me where I've
gone wrong, or recommend a change in Haskell 1.3 to cope with the problem if
it's real. The Crout reduction with pivoting follows. If anybody wants to try
it through the compiler. you'll need a test harness, which I'll be happy to
send, but I don't think I ought to waste net bandwidth to post it.
It's also unfortunate that array bounds were defined to be ((array of low
bounds),(array of high bounds)) [I know the arrays are really tuples] instead
of array of tuples (low,high). The latter could be used with the inRange
function from the Prelude, while the former cannot. But it'd probably be
_really_ hard on a lot of people to change this now.
Best regards,
Van Snyder
------------------------ Crout reduction with pivoting ----------------------
module Croutp (croutp) where
-- Crout method for LU factorization
import Linalg (ipamaxc,rcdot2)
-- ipamaxc a j p m n returns the index x of the element of p(m..n) such that
-- a!(p!(x),j) is the element of column j of a having the largest
-- absolute value.
-- rcdot2 a b i j m n computes the inner product of a(i,m..n) and b(m..n,j)
-- croutp takes matrix a and returns l and u factors in one matrix lu.
-- performs pivoting.
-- calculates values of lu from values of a and lu.
croutp :: (RealFrac v) => Array (Int,Int) v -> (Array (Int,Int) v,
Array Int (Array Int Int), Array Int Int, Array Int Int)
croutp a = if k==1 && l==1 && m<=n then (lu,p,mx,mk)
else error "crout: lower bounds not 1 or #rows > #columns"
where
b = bounds a
((k,l),(m,n)) = b
--t :: (RealFrac v) => Array (Int,Int) v
t = array ((1,1),(m,m))
([let k = p!1!i in (k,1) := a!(k,1) | i <- [1..m]] ++
[let k = p!s!i
in (k,s) := a!(k,s) - rcdot2 t lu k s 1 (s-1)
| s <- [2..m], i <- [s..m]])
--p :: Array Int (Array Int Int)
p = array (1,m)
([1 := array (1,m) [i := i | i <- [1..m]]]++
[s := let u = s-1
k = mk!u in
if u == k then p!u
else p!u // [u := mx!u, k := (p!u)!u] | s <- [2..m]])
--mk :: Array Int Int
--With the first definition of mk active, run-time insists there's a black hole.
--With the second, things work, but the function does no pivoting.
mk = array (1,m) [s := ipamaxc t s (p!s) s m | s <- [1..m]]
--mk = array (1,m) [s := s | s <- [1..m]]
mx :: Array Int Int
mx = array (1,m) [s := (p!s)!(mk!s) | s <- [1..m]]
--lu :: (RealFrac v) => Array (Int,Int) v
lu = array b
([(s,j) := t!(mx!s,j) | s <- [1..m], j <- [1..s]]++
[(s,j) := let k = mx!s in
(a!(k,j) - rcdot2 t lu k j 1 (s-1)) / lu!(s,s)
| s <- [1..m], j <- [s+1..n]])
module Linalg (iamax,iamaxc,iamaxr,ipamax,ipamaxc,ipamaxr,rcdot,rcdot2) where
-- iamax a returns the index of the element of a of maximum absolute value.
iamax :: (Real v) => Array Int v -> Int
iamax a = kamax l
where (l,h) = bounds a
kamax k = if k >= h then h else
let u = kamax (k+1) in if abs(a!u)>abs(a!k) then u else k
-- ipamax a p m n returns the index x of the element of p(m..n) such that
-- a!(p!(x)) is the element of a having largest absolute value.
ipamax :: (Real v) => Array Int v -> Array Int Int -> Int -> Int -> Int
ipamax a p m n = kamax m
where kamax k = if k >= n then n else
let u = kamax (k+1)
in if abs(a!(p!u)) > abs(a!(p!k)) then u else k
-- iamaxc a j returns the index of the element of
-- column j of a of maximum absolute value.
iamaxc :: (Real v) => Array (Int,Int) v -> Int -> Int
iamaxc a j = if inRange (l2,h2) j then kamax l1
else error "iamaxc: column number out of range"
where ((l1,l2),(h1,h2)) = bounds a
kamax k = if k >= h1 then h1 else
let u = kamax (k+1)
in if abs(a!(u,j)) > abs(a!(k,j)) then u else k
-- ipamaxc a j p m n returns the index x of the element of p(m..n) such that
-- a!(p!(x),j) is the element of column j of a having the largest
-- absolute value.
ipamaxc :: (Real v) =>
Array (Int,Int) v -> Int -> Array Int Int -> Int -> Int -> Int
ipamaxc a j p m n = kamax m
where kamax k = if k >= n then n else
let u = kamax (k+1)
in if abs(a!(p!(u),j)) > abs(a!(p!(k),j)) then u else k
-- iamaxr a i returns the index of the element of
-- row i of a of maximum absolute value.
iamaxr :: (Real v) => Array (Int,Int) v -> Int -> Int
iamaxr a i = if inRange (l1,h1) i then kamax l2
else error "iamaxr: row number out of range"
where ((l1,l2),(h1,h2)) = bounds a
kamax k = if k >= h2 then h2 else
let u = kamax (k+1)
in if abs(a!(i,u)) > abs(a!(i,k)) then u else k
-- ipamaxr a j p m n returns the index x of the element of p(m..n) such that
-- a!(j,p!(x)) is the element of row j of a having the largest
-- absolute value.
ipamaxr :: (Real v) =>
Array (Int,Int) v -> Int -> Array Int Int -> Int -> Int -> Int
ipamaxr a j p m n = kamax m
where kamax k = if k >= n then n else
let u = kamax (k+1)
in if abs(a!(j,p!(u))) > abs(a!(j,p!(k))) then u else k
-- rcdot a i j m n computes the inner product of a(i,m..n) and a(m..n,j)
rcdot :: (RealFrac v) => Array (Int,Int) v -> Int -> Int -> Int -> Int -> v
rcdot a i j m n = let ((l1,l2),(h1,h2)) = bounds a
d1 = (l1,h1); d2 = (l2,h2) in
if m > n then fromInteger 0 else
if not (inRange d1 i && inRange d1 m && inRange d1 n
&& inRange d2 j && inRange d2 m && inRange d2 n)
then error ("rcdot: index out of range\n"++
"d1 d2 i j m n = "++(show d1)++" "++(show d2)++" "++
(show i)++" "++(show j)++" "++(show m)++" "++(show n)++"\n")
else
if m == n then a!(i,m) * a!(m,j)
else let k = (m+n) `div` 2 in
(rcdot a i j m k) + (rcdot a i j (k+1) n)
-- rcdot2 a b i j m n computes the inner product of a(i,m..n) and b(m..n,j)
rcdot2 :: (RealFrac v) =>
Array (Int,Int) v -> Array (Int,Int) v -> Int -> Int -> Int -> Int -> v
rcdot2 a b i j m n = let ((la1,la2),(ha1,ha2)) = bounds a
((lb1,lb2),(hb1,hb2)) = bounds b
da1 = (la1,ha1); da2 = (la2,ha2)
db1 = (lb1,hb1); db2 = (lb2,hb2)
in
if m > n then fromInteger 0 else
if not (inRange da1 i && inRange db1 m && inRange db1 n
&& inRange db2 j && inRange da2 m && inRange da2 n)
then error "rcdot2: index out of range" else
if m == n then a!(i,m) * b!(m,j)
else let k = (m+n) `div` 2 in
(rcdot2 a b i j m k) + (rcdot2 a b i j (k+1) n)
