#
# Gaussian elimination and related routines for computing
# the determinant and inverse of a matrix over a field F
# and solving a linear system over F
# 
# Author: MBM 1989
#

ForwardElimination := proc(C,A,m,n) local i,j,k,s,t;
	
    s := 1;
    for k to m do 

        for i from k to m while A[i,k] = C[0] do od;
        if i > m then RETURN(FAIL) fi;
        for j from i+1 to m do
	    if A[j,k] = C[0] then next fi;
	    if length(A[j,k]) < length(A[i,k]) then i := j fi
	od;
        if i <> k then
	    s := -s;
	    for j from k to n do t := A[k,j]; A[k,j] := A[i,j]; A[i,j] := t od
        fi;

	userinfo(2,GaussianElimination,`forward elimination at row`,k);
        for i from k+1 to m do
	    if A[i,k] = C[0] then next fi;
	    t := C[`/`](A[i,k],A[k,k]);
	    for j from k+1 to n do A[i,j] := C[`-`](A[i,j],C[`*`](t,A[k,j])) od;
	    A[i,k] := 0
        od

    od;
    s

end:

BackElimination := proc(C,A,m,n) local i,j,k,t;
	
    for k from m by -1 to 1 do
	userinfo(2,GaussianElimination,`back elimination at row`,k);
	t := C[Inv](A[k,k]);
	for j from m+1 to n do A[k,j] := C[`*`](t,A[k,j]) od;
	for i to k-1 do
	    for j from m+1 to n do
		A[i,j] := C[`-`](A[i,j],C[`*`](A[i,k],A[k,j]))
	    od;
	    A[i,k] := C[0];
	od;
	A[k,k] := C[1]
    od;
    NULL

end:

MatrixDeterminant := proc(C,D,a) local d,i,n,A;

    n := nops(a);
    A := D[ArrayCoeffs](a);
    if n = 2 then RETURN(
	C[`-`](C[`*`](A[1,1],A[2,2]),C[`*`](A[1,2],A[2,1])) ) fi;
    d := ForwardElimination(C,A,n,n);
    if d = FAIL then RETURN(C[0]) fi;
    for i to n do d := C[`*`](d,A[i,i]) od;
    d

end:

MatrixInverse := proc(C,D,a) local A,i,j,n,d,I;

    n := nops(a);
    A := array(1..n,1..2*n);
    for i to n do
        for j to n do A[i,j] := a[i][j] od;
        for j from n+1 to 2*n do A[i,j] := C[0] od;
        A[i,n+i] := C[1]
    od;
    d := ForwardElimination(C,A,n,2*n);
    if d = FAIL then RETURN(FAIL) fi;
    d := BackElimination(C,A,n,2*n);
    I := array(1..n,1..n);
    for i to n do for j to n do I[i,j] := A[i,n+j] od od;
    D[Matrix](I,n,n)

end:

MatrixSolve := proc(C,D,a,b) local A,i,j,k,n,d;

    n := nops(a);
    for i to n do for j to n do A[i,j] := a[i][j] od; A[i,n+1] := b[i] od;
    d := ForwardElimination(C,A,n,n+1);
    if d = FAIL then RETURN(FAIL) fi;
    d := BackElimination(C,A,n,n+1);
    [ seq(A[k,n+1], k=1..n) ]

end:
	

save `GE.m`;
quit
