(* :Name: LinearAlgebra`GJSF` *)

(* :Title: Elementary Row Operations and Gauss Jordan Standard Form *)

(* :Author: Gary S. Stoudt *)

(* :Summary:

	This package is intended for educational purposes.  The package 

	performs the elementary row or column operations on a matrix.

	In addition, it will find elementary matrices that correspond

	to the row operations.  The Gauss Jordan Standard Form of a

	matrix can be found (reduced row-echelon form), and the matrix

	needed to get to that form is also produced.

  *)

(* :Context: LinearAlgebra`GJSF` *) 

(* :Package Version: Version 1 *)


(* :Copyright: Copyright 1994 Gary S. Stoudt *)

(* :History: Created Summer, 1994 *)

(* :Keywords: Elementary Matrix, Elementary Row Operation,

	      Gauss Jordan Elimination

  *)


(* :Source: Any elementary linear algebra textbook. *)

(* :Mathematica Version: 2.2 *)

(* :Limitation: *)



BeginPackage["LinearAlgebra`GJSF`"]



Type1::usage = 

"Type1[A_,i_,j_] performs a Type I row operation on the mxn matrix A;

that is, rows i and j are interchanged."



Type2::usage = 

"Type2[A_,i_,x_] performs a Type II row operation on the mxn matrix A;

that is, row i is replaced by x times row i."



Type3::usage = 

"Type3[A_,i_,j_,x_] performs a Type III row operation on the mxn matrix A;

that is, row j is replaced by row j plus x times row i."



E1::usage =

"E1[i_,j_,size_] yields a 'size' by 'size' square matrix that when 

left multiplied by a matrix performs a Type 1 row operation, and when

right multiplied performs a Type 1 column operation."



E2::usage =

"E2[i_,x_,size_] yields a 'size' by 'size' square matrix that when 

left multiplied by a matrix performs a Type 2 row operation, and when

right multiplied performs a Type 2 column operation."



E3::usage =

"E3[i_,j_,x_,size_] yields a 'size' by 'size' square matrix that when 

left multiplied by a matrix performs a Type 3 row operation, and when

right multiplied performs a Type 3 column operation."


slu::usage =

"slu[a_, l_, u_] yields a lower triangular matrix l and an upper triangular matrix u

such that a=lu.  This is LU factorization without pivoting."


Begin["LinearAlgebra`GJSF`Private`"]



Type1[A_,i_,j_]:=

     Module[{t=A,m=A},m[[i]]=t[[j]];m[[j]]=t[[i]];m];



E1[i_,j_,size_]:=Type1[IdentityMatrix[size],i,j]



Type2[A_,i_,x_]:=Module[{t=A,m=A},m[[i]]=x t[[i]];m];



E2[i_,x_,size_]:=Type2[IdentityMatrix[size],i,x]



Type3[A_,i_,j_,x_]:=

	Module[{t=A,m=A},m[[j]]=x t[[i]]+t[[j]];m];



E3[i_,j_,x_,size_]:=Type3[IdentityMatrix[size],i,j,x]

slu[a_,l_,u_]:=
  Module[{n,tol,i,j,k,aa},
    n=Dimensions[a][[1]]; tol=1 10^-6;
    low=Table[0,{i,n},{j,n}];up=Table[0,{i,n},{j,n}];    
    aa=a;
    For[k=1,k<=n, k++,
      low[[k,k]]=1;
	  For[i=k+1,i<=n, i++,
 	   low[[i,k]]= aa[[i,k]]/aa[[k,k]];
     	 For[j=k+1,j<=n, j++,
           aa[[i,j]] = aa[[i,j]] - low[[i,k]]aa[[k,j]]
   	     ];
	  ];
	  For[j=k,j<=n, j++,
      	     up[[k,j]] = aa[[k,j]]
  	  ];
    ];
      l=low;
      u=up;
   ];    
      


End[ ]

End[ ]