#
# Matrix(R)
#
# Creates a Matrix domain where R is a ring.
# This is really just a package of Matrix operations as
# it does not specify matrix dimensions.
#
# Author: MBM 1990 - 1992
#

macro( MatrixInput = `Matrix/Input`, MatrixOutput = `Matrix/Output` );
macro( Construct = `Matrix/Construct` );
macro( MatrixTranspose = `Matrix/Transpose` );
macro( MatrixAdd = `Matrix/Add`, MatrixMultiply = `Matrix/Mul` );
macro( MatrixVectorMultiply = `Matrix/Vector/Multiply` );
macro( VectorMatrixMultiply = `Vector/Matrix/Multiply` );

macro( HNF = readlib(`Matrix/HermiteNormalForm`, ``.gauss.`/MX/HNF.m`) );
macro( GE = readlib(`Matrix/GaussianElimination`, ``.gauss.`/MX/GE.m`) );
macro( GJ = readlib(`Matrix/GaussJordan`, ``.gauss.`/MX/GJ.m`) );
macro( NS = readlib(`Matrix/Nullspace`, ``.gauss.`/MX/NS.m`) );

Matrix := proc() local R,M,env;

    R := args[1];
    if not hasCategory(R,Ring) then ERROR(`1st argument must be a Ring`) fi;
    M := Set();
    M[CoefficientRing] := R;
    env := ['C' = R, 'D' = M];

    defOperation( New, [Integer,Integer,R] &-> M, M );
    defOperation( {Row,Col}, [M,Integer] &-> List(R), M );
    defOperation( {Rows,Cols}, M &-> Integer, M );
    defOperation( Matrix, { List(List(R)) &-> M,
        [Array(R),Integer,Integer] &-> M }, M );
    defOperation( Symmetric, M &-> Boolean, M );
    defOperation( Diagonal, List(R) &-> M, M );
    defOperation( Transpose, M &-> M, M );
    defOperation( ArrayCoeffs, M &-> Array(R), M );
    defOperation( {`+`,`*`}, [M,M] &-> M, M );
    defOperation( `-`, {M &-> M, [M,M] &-> M}, M );
    defOperation( `.`, [R,M] &-> M, M );
    defOperation( `*V`, [M,Vector(R)] &-> Vector(R), M ); 
    defOperation( `V*`, [Vector(R),M] &-> Vector(R), M ); 
    defOperation( RandMatrix, {[Integer,Integer] &-> M,
        [Integer,Integer,[] &-> R] &-> M}, M );

    if hasCategory(R,CommutativeRing) then
        defOperation( Rank, M &-> Integer, M );
    fi;
    if hasCategory(R,EuclideanDomain) then
        defOperation( {Hermite,Smith}, M &-> M, M );
        M[Hermite] := subs(env, proc(a) HNF(a,D) end);
    fi;
    if hasCategory(R,Field) then
        defOperations( {GaussianElimination, GaussJordan},
            [M,Name,Name] &-> M, M );
        defOperation( Nullspace, [M,Name] &-> List(Vector(R)), M );
        M[GaussianElimination] := subs(env, proc() GE(D,args) end);
        M[Rank] := subs(env, proc(a) local r; GE(D,a,'r'); r end);
        M[GaussJordan] := subs(env, proc() GJ(D,args) end);
        M[Nullspace] := subs(env, proc() NS(D,args) end);
    fi;

    if hasProperty(R,UniquelyRepresented) then M[`=`] := <evalb(x=y)> fi;
    M[Input] := subs(env, proc(a) MatrixInput(D,C,a) end);
    M[Output] := subs(env, proc(a) MatrixOutput(D,C,a) end);
    M[Type] := subs(env, proc(a) local r,x;
	if not type(a,list(list)) then RETURN(false) fi;
	for r in a do for x in r do
	    if not C[Type](x) then RETURN(false) fi;
	od od;
	true
	end);
    M[CheckDimensions] := proc() local n;
        for n in [args] do
            if not type(n,integer) or n < 1 then
                ERROR(`matrix dimensions must be a positive`,n) fi
        od;
        end;
    M[New] := subs(env,proc(m,n,x) D[CheckDimensions](m,n); [[x$n]$m] end);
    M[RandMatrix] := subs(env, proc(m,n) local i,j,f;
        if nargs = 3 then f := args[3] else f := C[Random] fi;
        [ seq([seq(f(), i=1..n)], j=1..m) ]
        end);
    M[Random] := subs(env, proc()
        D[RandMatrix](2+irem(rand(),3),2+irem(rand(),3))
        end);
    M[Row] := proc(x,r) x[r] end;
    M[Rows] := proc(x) nops(x) end;
    M[Col] := proc(x,c) local i,j,a;
        for i to nops(x) do a[i] := x[i][c] od;
        [ seq(a[j], j=1..nops(x)) ]
        end;
    M[Cols] := proc(x) nops(x[1]) end;
    M[Matrix] := proc() Construct(args) end;
    M[Diagonal] := subs(env, proc(d) local i,j,a,n;
        n := nops(d);
        a := array(1..n,1..n);
        for i to n do
            for j to n do
                if i = j then a[i,i] := d[i] else a[i,j] := C[0] fi
            od;
        od;
        D[Matrix](a,n,n)
        end);
    M[Symmetric] := subs( env, proc(x) local i,j,m,n;
        m := nops(x);
        n := nops(x[1]);
        if m <> n then RETURN( false ) fi;
        for i to m do for j to n do
            if C[`<>`](x[i][j],x[j][i]) then RETURN( false ) fi;
        od od;
        true
        end);

    M[ArrayCoeffs] := proc(x) local a,i,j,m,n,r;
        m := nops(x); n := nops(x[1]); a := array(1..m,1..n);
        for i to m do r := x[i]; for j to n do a[i,j] := r[j] od od;
        op(a)
        end;
    M[Error] := proc() ERROR(`matrix dimensions incompatible`) end;
    M[`+`] := subs(env, proc(x,y) MatrixAdd(D,C,x,y) end);
    M[`-`] := subs(env, proc(x,y)
        if nargs = 1 then D[`R*`](-1,x) else D[`+`](x,D[`-`](y)) fi
        end);
    M[`*`] := subs(env, proc(x,y) MatrixMultiply(D,C,args) end);
    M[`.`] := subs(env, proc(x,a) local m,n,i,j,c;
        m := D[Rows](a); n := D[Cols](a); c := array(1..m,1..n);
        for i to m do for j to n do c[i,j] := C[`*`](x,a[i][j]) od od;
        D[Matrix](c,m,n)
        end);
    M[Transpose] := subs(env, proc(a) MatrixTranspose(D,a) end);
    M[`*V`] := subs(env, proc(x,v) MatrixVectorMultiply(D,C,x,v) end);
    M[`V*`] := subs(env, proc(x,v) VectorMatrixMultiply(D,C,x,v) end);
    op(M)
end:

`Matrix/Input` := proc(M,R,a) local i,j,m,n,r;
    if type(a,matrix) then
        RETURN( `Matrix/Input`(M,R,convert(a,listlist)) )
    fi;
    if type(a,list) then m := nops(a) else RETURN(FAIL) fi;
    if m>0 and type(a[1],list) then n := nops(a[1]) else RETURN(FAIL) fi; 
    for i to m do
        if not type(a[i],list) or nops(a[i]) <> n then RETURN(FAIL) fi;
        r[i] := map(R[Input],a[i]);
        if member(FAIL,r[i]) then RETURN(FAIL) fi
    od;
    M[Matrix]([seq(r[j] , j=1..m)])
end:

`Matrix/Output` := proc(D,C,a) local i,j,m,n,r;
    m := D[Rows](a); n := D[Cols](a);
    for i to m do
        for j to n do
        r[i,j] := C[Output](a[i][j]);
        od
    od;
    [seq([seq(r[i,j], j=1..n)], i=1..m)];
end:

`Matrix/Construct` := proc(x,m,n) local i,j;
    if type(x,list) then x
    elif type(x,array) and nargs = 3 then
        if type(m,integer) and m > 0 and type(n,integer) and n > 0
        then [ seq([ seq(x[i,j], j=1..n) ], i=1..m) ]
            else ERROR(`invalid dimensions`)
        fi
    else ERROR(`1st argument must be a list of lists or an array`)
    fi
end:

`Matrix/Add` := proc(M,R,a,b) local i,j,m,n,c;

    m := M[Rows](a); n := M[Cols](a); c := array(1..m,1..n);
    if m <> M[Rows](b) or n <> M[Cols](b) then
        ERROR(`incompatible dimensions`) fi;
    for i to m do for j to n do c[i,j] := R[`+`](a[i][j],b[i][j]) od od;
    M[Matrix](c,m,n)
end:

`Matrix/Mul` := proc(M,R,x,y) local a,b,c,i,j,k,l,m,n,add,mul;

    if nargs = 3 then RETURN(x) fi;
    if nargs > 4 then
        RETURN( procname(M,R,procname(M,R,x,y),args[5..nargs]) ) fi;
    if type(x,integer) then ERROR(`not implemented`) fi;

    l := M[Cols](x); m := M[Rows](x); n := M[Rows](y);
    if m <> M[Cols](y) then ERROR(`incompatible dimensions`) fi;
    a := M[ArrayCoeffs](x); b := M[ArrayCoeffs](y); c := array(1..l,1..n);
    add := eval(R[`+`]):
    mul := eval(R[`*`]):
    for i to l do
        for j to n do
        #c[i,j] := R[`*`](a[i,1],b[1,j]);
        c[i,j] := mul(a[i,1],b[1,j]);
        #c[i,j] := a[i,1]*b[1,j];
        for k from 2 to m do
            #c[i,j] := R[`+`]( R[`*`](a[i,k],b[k,j]), c[i,j] )
            c[i,j] := add( mul(a[i,k],b[k,j]), c[i,j] )
            #c[i,j] := a[i,k]*b[k,j]+c[i,j]
        od
        od
    od;
    M[Matrix](c,l,n)
end:

`Matrix/Vector/Multiply` := proc(M,R,x,v) local a,i,j,k,m,n,r;
	
	m := M[Rows](x); n := M[Cols](x);
	if n <> nops(v) then ERROR(`incompatible dimensions`) fi;
	a := array(1..m);
	for i to m do
	    r := M[Row](x,i);
	    a[i] := R[`+`]( seq(R[`*`](r[j],v[j]), j=1..n) );
	od;
	[seq(a[k], k=1..m)]
end:

`Vector/Matrix/Multiply` := proc(M,R,v,x) local a,i,j,k,m,n,r;
	
	m := M[Rows](x); n := M[Cols](x);
	if m <> nops(v) then ERROR(`incompatible dimensions`) fi;
	a := array(1..n);
	for i to n do
	    r := M[Col](x,i);
	    a[i] := R[`+`]( seq(R[`*`](r[j],v[j]), j=1..m) );
	od;
	[seq(a[k], k=1..n)]
end:

`Matrix/Transpose` := proc(M,a) local i,j,m,n,t;
	m := M[Rows](a); n := M[Cols](a); t := array(1..n,1..m);
	for i to m do for j to n do t[j,i] := a[i][j] od od;
	M[Matrix](t,n,m)
end:

save `Matrix.m`;
quit
