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)