{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Copyright : (c) Mikael Johansson 2006
Maintainer : mik@math.uni-jena.de
Stability : provisional
Portability : requires multi-parameter type classes
Routines and abstractions for Matrices and
basic linear algebra over fields or rings.
-}
module MathObj.Matrix where
import qualified Algebra.Module as Module
import qualified Algebra.Vector as Vector
import qualified Algebra.Ring as Ring
import qualified Algebra.Additive as Additive
import Algebra.Module((*>), )
import Algebra.Ring((*), fromInteger, scalarProduct, )
import Algebra.Additive((+), (-), zero, subtract, )
import Data.Array (Array, listArray, elems, bounds, (!), ixmap, range, )
import qualified Data.List as List
import Control.Monad (liftM2, )
import Control.Exception (assert, )
import Data.Tuple.HT (swap, )
import Data.List.HT (outerProduct, )
import NumericPrelude (Integer, )
import PreludeBase hiding (zipWith, )
{- |
A matrix is a twodimensional array of ring elements, indexed by integers.
-}
data {-(Ring.C a) =>-}
T a = Cons (Array (Integer, Integer) a) deriving (Eq,Ord,Read)
{- |
Transposition of matrices is just transposition in the sense of
Data.List.
-}
transpose :: T a -> T a
transpose (Cons m) =
let (lower,upper) = bounds m
in Cons (ixmap (swap lower, swap upper) swap m)
rows :: T a -> [[a]]
rows (Cons m) =
let ((lr,lc), (ur,uc)) = bounds m
in outerProduct (curry(m!)) (range (lr,ur)) (range (lc,uc))
columns :: T a -> [[a]]
columns (Cons m) =
let ((lr,lc), (ur,uc)) = bounds m
in outerProduct (curry(m!)) (range (lc,uc)) (range (lr,ur))
fromList :: Integer -> Integer -> [a] -> T a
fromList m n xs = Cons (listArray ((1,1),(m,n)) xs)
instance (Ring.C a, Show a) => Show (T a) where
show m = List.unlines $ map (\r -> "(" ++ r ++ ")")
$ map (unwords . map show) $ rows m
dimension :: T a -> (Integer,Integer)
dimension (Cons m) = uncurry subtract (bounds m) + (1,1)
numRows :: T a -> Integer
numRows = fst . dimension
numColumns :: T a -> Integer
numColumns = snd . dimension
-- These implementations may benefit from a better exception than
-- just assertions to validate dimensionalities
instance (Additive.C a) => Additive.C (T a) where
(+) = zipWith (+)
(-) = zipWith (-)
zero = zeroMatrix 1 1
zipWith :: (a -> b -> c) -> T a -> T b -> T c
zipWith op mM@(Cons m) nM@(Cons n) =
let d = dimension mM
em = elems m
en = elems n
in assert (d == dimension nM) $
uncurry fromList d (List.zipWith op em en)
zeroMatrix :: (Additive.C a) => Integer -> Integer -> T a
zeroMatrix m n = fromList m n $
List.repeat zero
-- List.replicate (fromInteger (m*n)) zero
instance (Ring.C a) => Ring.C (T a) where
mM * nM = assert (numRows mM == numColumns nM) $
fromList (numColumns mM) (numRows nM)
(liftM2 scalarProduct (rows mM) (columns nM))
fromInteger n = fromList 1 1 [fromInteger n]
instance Functor T where
fmap f (Cons m) = Cons (fmap f m)
instance Vector.C T where
zero = zero
(<+>) = (+)
(*>) = Vector.functorScale
instance Module.C a b => Module.C a (T b) where
x *> m = fmap (x*>) m
{- |
What more do we need from our matrix class? We have addition,
subtraction and multiplication, and thus composition of generic
free-module-maps. We're going to want to solve linear equations with
or without fields underneath, so we're going to want an implementation
of the Gaussian algorithm as well as most probably Smith normal
form. Determinants are cool, and these are to be calculated either
with the Gaussian algorithm or some other goodish method.
-}
{- |
We'll want generic linear equation solving, returning one solution,
any solution really, or nothing. Basically, this is asking for the
preimage of a given vector over the given map, so
a_11 x_1 + .. + a_1n x_n = y_1
...
a_m1 x_1 + .. + a_mn a_n = y_m
has really x_1,...,x_n as a preimage of the vector y_1,..,y_m under
the map (a_ij), since obviously y_1,..,y_m = (a_ij) x_1,..,x_n
So, generic linear equation solving boils down to the function
-}
preimage :: (Ring.C a) => T a -> T a -> Maybe (T a)
preimage a y = assert
(numRows a == numRows y && -- they match
numColumns y == 1) -- and y is a column vector
Nothing
{-
Cf. /usr/lib/hugs/demos/Matrix.hs
-}