(*
	Matrix operations: row and columns operations,
	the exponential of a matrix
	
	Bruce Ikenaga
	Department of Mathematics
	Case Western Reserve University
	Cleveland, Ohio 44106
	bxi@po.cwru.edu
	
	Last revision: 4-9-91
*)

BeginPackage["MatrixOperations`"]

RowSwap::usage = "RowSwap[m, r1, r2] swaps rows
r1 and r2 of matrix m. Note: RowSwap[] returns the
resulting matrix; the original matrix m is not
affected.";

RowMult::usage = "RowMult[m, r, x] multiplies row
r of matrix m by x. Note: RowMult[] returns the
resulting matrix; the original matrix m is not
affected.";

RowComb::usage = "RowComb[m, r1, r2, x] replaces row
r1 of matrix m by r1 plus x times row r2. Note:
RowComb[] returns the resulting matrix; the original
matrix m is not affected.";

ColSwap::usage = "ColSwap[m, c1, c2] swaps columns c1
and c2 of matrix m. Note: ColSwap[] returns the
resulting matrix; the original matrix m is not
affected.";

ColMult::usage = "ColMult[m, c, x] multiplies column
c of matrix m by x. Note: ColMult[] returns the
resulting matrix; the original matrix m is not
affected.";

ColComb::usage = "ColComb[m, c1, c2, x] replaces column
c1 of matrix m by c1 plus x times column c2. Note:
ColComb[] returns the resulting matrix; the original
matrix m is not affected.";

Convolution::usage = "Convolution[f, g, t]
computes the convolution (f * g)(t).";

MatrixExponential::usage = "MatrixExponential[A, t]
computes the exponential E^(A t), where A is a
square matrix. MatrixExponential calls Eigenvalues
to find the eigenvalues of A.";

Begin["`Private`"]

RowSwap[ matrix_, row1_, row2_ ] :=
	Block[ {temp = matrix},
		temp[[row1]] = matrix[[row2]];
		temp[[row2]] = matrix[[row1]];
		temp
	];
	
RowMult[ matrix_, row_, factor_ ] :=
	Block[ {temp = matrix},
		temp[[row]] = factor matrix[[row]];
		temp
	];
	
RowComb[ matrix_, row1_, row2_, factor_ ] :=
	Block[ {temp = matrix},
		temp[[row1]] =
			matrix[[row1]] + factor matrix[[row2]];
		temp
	];
	
ColSwap[ matrix_, col1_, col2_ ] :=
	Block[ {temp = Transpose[matrix]},
		temp[[col1]] = Transpose[matrix][[col2]];
		temp[[col2]] = Transpose[matrix][[col1]];
		Transpose[temp]
	];
ColMult[ matrix_, col_, factor_ ] :=
	Block[ {temp = Transpose[matrix]},
		temp[[col]] =
			factor Transpose[matrix][[col]];
		Transpose[temp]
	];
ColComb[ matrix_, col1_, col2_, factor_ ] :=
	Block[ {temp = Transpose[matrix]},
		temp[[col1]] =
			Transpose[matrix][[col1]] +
				factor Transpose[matrix][[col2]];
		Transpose[temp]
	];

Convolution[ f_, g_, var_ ] :=
	Block[ {u, v},
		Integrate[ f[u - v] g[v], {v, 0, u} ]
			/. u->var ];

MatrixExponential[ arr_, var_ ] :=
	Block[
		{ charlist = Eigenvalues[ arr ],
			id = IdentityMatrix[ Length[ arr ]],
			w, temp = 0 * id, tempa, tempb },
		tempa = (E^(charlist[[1]] #) &);
		tempb = id;
		temp = temp + tempa[w] * tempb;
		Do[	tempa = Function[t,
			Release[Convolution[ E^(charlist[[i]] #)&,
				tempa, t]]];
			tempb = tempb . (arr - charlist[[i-1]] id);
			temp = temp + tempa[w] * tempb,
			{i, 2, Length[arr]} ];
		temp /. w->var
	];

End[]

EndPackage[]