
(*:Mathematica:: V1.2  *)

(*:Context: "LinearAlgebra`MatrixManipulation`"  *)

(*:Title: Basic Matrix Manipulation Functions  *)

(*:Author:   Kevin McIsaac  *)

(*:Summary:  Provides basic matrix manipulation routines for composing matrices 
from block elements and taking subblocks from a matrix.  *)

(*:Keywords: matrix, block, submatrix  *)

(*:Requirements: none.  *)

(*:Sources:   *)

(*:History:   Created January 25, 1989 by K. McIsaac.
                   Updated February 4, 1989 by John Novak.  *)

BeginPackage["LinearAlgebra`MatrixManipulation`"];

AppendColumns::usage=
"AppendColumns[mat1, mat2...] returns a new matrix composed of the
submatrices mat1, mat2, etc., by joining their columns. The submatrices must
all have the same number of columns."

AppendRows::usage=
"AppendRows[mat1, mat2...] returns a new matrix composed of the submatrices
mat1, mat2, etc., by joining their rows. The submatrices must all have the same
number of rows."

SquareMatrixQ::usage=
"SquareMatrixQ[mat] tests whether mat is square." 

TakeRows::usage=
"TakeRows[mat, part] takes the rows of matrix mat give by the partspec,
part. This uses the same notation as Take."

TakeColumns::usage=
"TakeColumns[mat, part] takes the columns of matrix mat give by the
partspec, part. This uses the same notation as Take."

TakeMatrix::usage=
"TakeMatrix[mat, start, end] returns the submatrix starting at start and
ending at end, where the positions are positions in the matrix, e.g., {1,2}."

SubMatrix::usage=
"SubMatrix[mat, start, dim] returns the submatrix of dimension dim starting
at position start in mat."

BlockMatrix::usage=
"BlockMatrix[block] returns a matrix composed of a matrix of matrices."

TakeMatrix::usage=
"TakeMatrix[mat, start, end] returns the submatrix that starts at start and
ends at end, where start and end are expressed as the pair defining a position
in the array (e.g., {1,2})."

UpperDiagonalMatrix::usage=
"UpperDiagonalMatrix[fn, size] creates an upper-diagonal matrix with the i-th,
j-th element being f[i,j] for i>= j and zero otherwise."

LowerDiagonalMatrix::usage=
"LowerDiagonalMatrix[fn, size] creates a tridiagonal matrix with the i-th, j-th
element being f[i,j] for i<= j and zero otherwise."

TridiagonalMatrix::usage=
"TridiagonalMatrix[fn, size] creates a tridiagonal matrix with the i-th, j-th
element being f[i,j] for i= j and zero otherwise."

ZeroMatrix::usage=
"ZeroMatrix[m] returns an m x m block matrix of zeros. ZeroMatrix[m, n]
returns an m x n block matrix of zeros."

HilbertMatrix::usage =
"HilbertMatrix[m,n] returns the m x n Hilbert Matrix, where the
(i,j)-th element is defined as 1/(i + j + 1). HilbertMatrix[m] returns the m x
m Hilbert matrix."  

HankelMatrix::usage =
"HankelMatrix[column,row] returns the Hankel matrix defined by the
given column and row; note that the last element of column must match the
first element of row.  HankelMatrix[col] gives the Hankel matrix whose first
column is col, and where every element under the anti-diagonal is filled with
zeros. HankelMatrix[n] gives the n x n Hankel matrix where the leading column
is the integers 1 - n."

Begin["`Private`"]

AppendColumns[l__?MatrixQ] := Join[l] /; SameColumnSize[{l}]

AppendRows[l__?MatrixQ] := Apply[Join, Transpose[{l}], {1}] /; 
   SameRowSize[{l}]

BlockMatrix[block_] :=
	AppendColumns @@ Apply[AppendRows, block, {1}] 

SquareMatrixQ[mat_?MatrixQ] := SameQ @@ Dimensions[mat]

SameColumnSize[l_List] := (SameQ @@ (Dimensions[#][[2]]& /@ l) )
	SameRowSize[l_List] := (SameQ @@ (Dimensions[#][[1]]& /@ l) )


TakeRows[mat_?MatrixQ, part_] := Take[mat, part]

TakeColumns[mat_?MatrixQ, part_] := Take[#, part]& /@ mat

TakeMatrix[mat_?MatrixQ, start:{startR_Integer, startC_Integer},
end:{endR_Integer, endC_Integer}] :=
	Take[#, {startC, endC}]& /@ Take[mat, {startR, endR}] /;
	And @@ Thread[Dimensions[mat] >= start] && 
	And @@ Thread[Dimensions[mat] >= end]

SubMatrix[mat_List, start:{_Integer, _Integer}, dim:{_Integer,_Integer}] :=
	TakeMatrix[mat, start, start+dim-1] 

TridiagonalMatrix[fn_,size_Integer] :=
	DiagonalMatrix[Map[fn[#,#]&,Range[1,size]]]

(* TridiagonalMatrix[fn_, size_Integer] :=
	Array[If[#1===#2, fn[#1, #2], 0]&, {size, size}] /;
	size >= 0
*)

LowerDiagonalMatrix[fn_, size_Integer] :=
	Array[If[#1>=#2, fn[#1, #2], 0]&, {size, size}] /;
	size >= 0

UpperDiagonalMatrix[fn_, size_Integer] :=
	Array[If[#1<=#2, fn[#1, #2], 0]&, {size, size}] /;
	size >= 0

ZeroMatrix[m_Integer,n_Integer] := 
	Table[0, {m}, {n}] /; m >= 0 && n>=0

ZeroMatrix[m_Integer] := ZeroMatrix[m, m] /; m >= 0 

HilbertMatrix[m_Integer,n_Integer] := 
	Table[1/(i+j+1), {i, 0, m-1}, {j, 0, n-1}]

HilbertMatrix[m_Integer] := HilbertMatrix[m,m]

HankelMatrix[size_Integer] :=
	HankelMatrix[Range[1,size]] /; size >= 0

HankelMatrix[col_List] :=
	Module[{size = Length[col]},
		HankelMatrix[Join[col,Table[0,{size - 1}]],
			{size,size}]
	]

HankelMatrix[col_List,row_List] :=
	HankelMatrix[Join[col,Drop[row,1]],
		{Length[col],Length[row]}]/;
	Last[col] == First[row]

HankelMatrix[vec_List,size:{rows_Integer,cols_Integer}] :=
	Array[(vec[[#1 + #2 - 1]])&,size]/;
		Length[vec] >= (rows + cols - 1)

End[]

EndPackage[]

