(* :Name: LinearAlgebra`GaussianElimination` *)

(* :Title: Gaussian Elimination and Back Substitution *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
This package implements LU decomposition with partial pivoting.
It is not intended for computational purposes, but rather for
educational purposes.  This package can be used in conjunction
with NumericalMath`ComputerArithmetic`.
*)

(* :Context: LinearAlgebra`GaussianElimination` *) 

(* :Package Version: Mathematica 2.0 *)

(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc.
*)

(* :History:
	Original version by Jerry B. Keiper, May 1990.
	Revised by Jerry B. Keiper, December 1990.
*)

(* :Keywords:
	Gaussian elimination, back substitution, partial pivoting,
	LU decomposition
*)

(* :Source: Any elementary numerical analysis textbook. *)

(* :Mathematica Version: 2.0 *)

(* :Limitation: *)

BeginPackage["LinearAlgebra`GaussianElimination`"]

LUFactor::usage = 
"LUFactor[mat] gives the LU decomposition along with a pivot list of the
matrix mat.  The calculation is done using Gaussian elimination with
partial pivoting.  The result is returned as a data object with a head of LU.
LUFactor may be used in the context of ComputerArithmetic.m or
IntervalArithmetic.m."

LUSolve::usage =
"LUSolve[lu, b] solves the linear system represented by lu (the matrix) and
right-hand side b.  The \"matrix\" lu must be the data object with a head of
LU resulting from LUFactor[mat] and b must be an ordinary list of the
appropriate size.  LUSolve may be used in the context of
ComputerArithmetic.m or IntervalArithmetic.m."

LU::usage =
"LU[a, pivots] is the data object returned by LUFactor and is to be given
as the first argument to LUSolve.  The matrix a represents the lower and
upper triangular matrices in the LU decomposition, and pivots is a record
of the pivots used in the decomposition."

Unprotect[LUFactor, LUSolve, LU];

Begin["LinearAlgebra`GaussianElimination`Private`"]

LUFactor[aa_?MatrixQ] :=
    Module[{ans}, ans /; Head[ans = lufactor[aa]] == LU] /;
	(Length[aa] == Length[aa[[1]]])

LUSolve[ap_LU, bb_List] :=
    Module[{ans}, ans /; VectorQ[ans = lusolve[ap, bb]]] /;
	(Length[ap[[1]]] == Length[bb])

divide[a_, b_] :=
    If[FreeQ[{a, b}, NumericalMath`ComputerArithmetic`ComputerNumber],
	a/b, NumericalMath`ComputerArithmetic`IdealDivide[a, b]];

lufactor[aa_] :=
    Module[{a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, n=Length[aa], tmp},
	pivot = Table[i, {i, n}];
	For[ii=1, ii<=n, ii++,	(* for each row do ... *)
	    (* find a pivot *)
	    mpiv = Abs[a[[pivot[[ii]], ii]]];
	    k = ii;
	    For[i=ii+1, i<=n, i++,
		tmp = Abs[a[[pivot[[i]], ii]]];
		If[tmp > mpiv, mpiv = tmp; k = i]
		];
	    tmp = pivot[[ii]];
	    pivot[[ii]] = iip = pivot[[k]];
	    pivot[[k]] = tmp;
	    mpiv = a[[iip, ii]];
	    (* calculate and store the multipliers and reduce the other elements. *)
	    For[i=ii+1, i<=n, i++,
		ip = pivot[[i]];
		m = a[[ip, ii]] = divide[a[[ip, ii]], mpiv];
		For[j=ii+1, j<=n, j++, a[[ip, j]] -= m*a[[iip, j]]]
		];
	    ];
	LU[a, pivot]
	]
			
lusolve[ap_, bb_] :=
    Module[{a, b = bb, pivot, ii, i, m, n, tmp},
	a = ap[[1]];
	pivot = ap[[2]];
	n = Length[a];
	(* forward elimination *)
	For[ii=1, ii<=n, ii++,  (* for each row do ... *)
	    For[i=ii+1, i<=n, i++,
		b[[pivot[[i]]]] -= a[[pivot[[i]], ii]]*b[[pivot[[ii]]]]
		];
	    ];
	(* back substitution *)
	Clear[i];
	For[ii=n, ii>0, ii--,
	    tmp = Sum[a[[pivot[[ii]], i]]*b[[pivot[[i]]]], {i, ii+1, n}];
	    tmp = b[[pivot[[ii]]]] - tmp;
	    b[[pivot[[ii]]]] = divide[tmp, a[[pivot[[ii]], ii]]]
	    ];
	tmp = Table[Null,{n}];
	For[ii=1, ii<=n, ii++, tmp[[ii]] = b[[pivot[[ii]]]]];
	tmp
	]

End[ ] (* "LinearAlgebra`GaussianElimination`Private`" *)

Protect[LUFactor, LUSolve, LU];

EndPackage[ ] (* "LinearAlgebra`GaussianElimination`" *)

(*
test3 := Module[{a, x, b, r, i, j}, 
	a = Table[N[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]; 
	b = Table[N[Random[Integer, {-10, 10}]], {i, 3}]; 
	x = LUSolve[LUFactor[a], b];
	r = a . x - b;
	{a, b, x, r}
	]

test5 := Module[{a, b, i, j}, 
	a = Table[N[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]; 
	b = Table[N[Random[Integer, {-10, 10}]], {i, 5}]; 
	a . LUSolve[LUFactor[a], b] - b
	]

cntest3 := Module[{a, x, b, r, i, j}, 
	a = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]; 
	b = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 3}]; 
	x = LUSolve[LUFactor[a], b];
	r = a . x - b;
	{a, b, x, r}
	]

cntest5 := Module[{a, b, i, j}, 
	a = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]; 
	b = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 5}]; 
	a . LUSolve[LUFactor[a], b] - b
	]

inttest3 := Module[{a, x, b, r, i, j}, 
	a = Table[Interval[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]; 
	b = Table[Interval[Random[Integer, {-10, 10}]], {i, 3}]; 
	x = LUSolve[LUFactor[a], b];
	r = a . x - b;
	{a, b, x, r}
	]

inttest5 := Module[{a, b, i, j}, 
	a = Table[Interval[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]; 
	b = Table[Interval[Random[Integer, {-10, 10}]], {i, 5}]; 
	a . LUSolve[LUFactor[a], b] - b
	]

ninttest3 := Module[{a, x, b, r, i, j}, 
	a = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]];
	b = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 3}]];
	x = LUSolve[LUFactor[a], b];
	r = a . x - b;
	{a, b, x, r}
	]

ninttest5 := Module[{a, b, i, j}, 
	a = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]];
	b = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 5}]];
	a . LUSolve[LUFactor[a], b] - b
	]
*)
