(*:Mathematica Version: 2.0  *)

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

(*:Context: LinearAlgebra`MatrixManipulation`  *)

(*:Title: Basic Matrix Manipulation Functions  *)

(*:Authors:   Kevin McIsaac, Daniel Lichtblau *)

(*:Summary:
This package provides basic matrix manipulation routines for composing a matrix
from block elements and extracting subblocks from a matrix.  Included are
functions for generating special matrices such as Hilbert and Hankel.
Also provided is a routine for extracting a matrix and right-hand-side
vector from a system of linear equations is a given set of variables, and
a routine to find the polar decomposition of a square matrix.

Finally, we have routines to find the norm of a vector, a matrix, or its
inverse. The vector p-norm, for a given vector vec and 1<=p<=Infinity, is
defined to be

Sum[Abs[v[[j]]]^p, {j,Length[p]}]^(1/p)

We implement this for a vector of approximate numbers. For a matrix mat, the
p-norm is defined to be the maximum, over all vectors vec of p-norm one, of the
p-norm of mat.vec. While it is in general not obvious how to compute this, the
cases where p is 1 or Infinity are trivial, and the case where p is 2 is not
hard. For the norms of the inverse matrix one can explicitly compute
Inverse[mat] and work from that, but it is asymptotically about three times
faster to work with the LUDecomposition of mat. Hence we allow either mat or
its LUDecomposition as an argument to InverseMatrixNorm, with the latter case
being significantly faster than the former, which is in turn faster than
MatrixNorm[Inverse[mat]]
*)


(*:Limitations:

The matrix norm is only approximated for the case p=2. The inverse norm is
only approximated for all of the allowable values of p. Hence, while we
use the precision of the input the process of computing the norms, we only
return a machine number, as at best a few digits are correct. This is adequqate
because typically one is interested only in the order of magnitude of the norm.

The techniques used, while quite reasonable, may fail for specialized classes
of input. It should also be noted that any technique used to find the inverse
matrix norm (especially including explicit computation of the inverse matrix)
will fail if the input has inadequate precision. It is easy to recognize this
case, because the log of the condition number will typically be larger than
the precision of the input. In terms of efficiency and correctness, it is
clear from the code that one can get a good value for either the 1-or-infinity
norm of a matrix. If one is concerned that the matrix may be pathological for
the inverse-norm estimators, perhaps the safest ploy is to estimate the 2-norm
of the inverse; this uses the power method, and is fairly accurate. But note
that the 1-and-infinity estimators also use an abbreviated power-method
estimate as a reality check on the main method, so these should not do
anything drastically wrong. Note also that we use a fixed number of iterations
in the power-method computations. One could attain more accuracy (at expense
of run-time) by altering these to use e.g. FixedPoint with a SameTest option
to determine when the desired accuracy has been met.
*)

(*:Sources:   
G. Golub and C. Van Loan, Matrix Computations, The Johns Hopkins
University Press, 1996.

W. Hager. Condition Estimates. SIAM J. Sci. Stat. Comput. 5(2). 311-316. June 19
84

N. Higham. FORTRAN codes for estimating the one-norm of a real or complex matrix
, with applications to condition estimation. ACM Transactions on Mathematical So
ftware 14(4). 381-396. December 1988.
*)

(*:Keywords: matrix, block, submatrix, linear equations,
polar decomposition, matrix norm, matrix condition number *)

(*:Requirements: None.  *)

(*:Package Version: 1.2 *)

(*:History:   Created January 25, 1989 by K. McIsaac.
        Updated February 4, 1989 by John Novak.
	LinearEquationsToMatrices added April 1995 by Daniel Lichtblau
	PolarDecomposition added October 1995 by Daniel Lichtblau
	HankelMatrix correction contributed by Jan Jelinek; minor
	  changes to error checking in LinearEquationsToMatrices and
	  PolarDecomposition, ECM,  Feb. 1997.
	Added VectorNorm, MatrixNorm, MatrixConditionNumber April 1997,
	  Daniel Lichtblau.
   	V1.2: Added LUMatrices, WRI 1998.
*)

BeginPackage["LinearAlgebra`MatrixManipulation`"];

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

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

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

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

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

TakeMatrix::usage=
"TakeMatrix[mat, start, end] gives the submatrix starting at position start
and ending at position end, where each position is expressed as a 
{row, column} pair, e.g., {1, 2}."

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

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

UpperDiagonalMatrix::usage=
"UpperDiagonalMatrix[f, m] gives an upper-diagonal matrix, m x m,
using a function of two arguments f.  The (i,j)-th element of the matrix is
f[i, j] for i >= j and zero otherwise."

LowerDiagonalMatrix::usage=
"LowerDiagonalMatrix[f, m] gives a lower-diagonal matrix, m x m,
using a function of two arguments f.  The (i,j)-th element of the matrix is
f[i, j] for i <= j and zero otherwise."

TridiagonalMatrix::usage=
"TridiagonalMatrix[f, m] gives a tridiagonal matrix, m x m,
using a function of two arguments f.  The (i,j)-th element of the matrix is
f[i,j] for Abs[i - j] <= 1 and zero otherwise."

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

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

HankelMatrix::usage =
"HankelMatrix[column, row] gives 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, 2, ..., n."

LinearEquationsToMatrices::usage =
"LinearEquationsToMatrices[eqns, vars] gives a list of the form {mat, vec}.
The first element is the matrix of coefficients, the second is the vector of
right-hand-sides, such that Thread[mat . vars == vec] gives eqns, possibly
in expanded form.  LinearEquationsToMatrices[eqns, vars, CheckLinearity ->
True] checks that eqns are linear (of total degree 1) in vars.  The check can
slow the conversion significantly, but insures that an incorrect result will
not be returned when eqns are not linear in vars.  The default setting for
CheckLinearity is False."

CheckLinearity::usage =
"CheckLinearity is an option to LinearEquationsToMatrices specifying whether 
or not to check that the equations are linear (of total degree 1) in the
variables."

PolarDecomposition::usage = 
"PolarDecomposition[mat] gives a list {u, s} comprising the polar
decomposition of the matrix mat. Specifically, s is a positive definite matrix 
and u satisfies  u . Transpose[Conjugate[u]] == IdentityMatrix[Length[mat]].
Note that u . s == mat."

VectorNorm::usage =
"VectorNorm[vec] gives the infinity-norm (the largest component) of
the vector vec of approximate numbers. VectorNorm[vec, p] gives the
p-norm of vec."

MatrixNorm::usage =
"MatrixNorm[mat] gives the infinity-norm of the matrix mat of
approximate numbers. MatrixNorm[mat, p] gives the p-norm, where
p must be 1, 2, or Infinity."

InverseMatrixNorm::usage =
"InverseMatrixNorm[mat] gives the infinity-norm of the matrix Inverse[mat],
where mat is a matrix of approximate numbers.  Given an LUDecomposition lu
for mat, InverseMatrixNorm[lu] gives the infinity-norm of Inverse[mat].
InverseMatrixNorm[mat, p] and InverseMatrixNorm[lu, p] compute the p-norms,
where p must be 1, 2, or Infinity."

MatrixConditionNumber::usage =
"MatrixConditionNumber[mat] gives the infinity-norm condition number
of the matrix mat of approximate numbers.  MatrixConditionNumber[mat, p]
gives the condition number in the p-norm, where p must be 1, 2, or Infinity."

LUMatrices::usage =
"LUMatrices[lu] gives a list {lmat, umat} containing the lower and upper
triangular matrices produced by the LU decomposition of a matrix m, where
lu is the combined matrix given by the first element of LUDecomposition[m]."


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, m}, {j, n}]

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

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

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

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

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

polytotaldegree[poly_, vars_List] :=
    Module[{j, len=Length[poly], nvars=Length[vars], hed=Head[poly],
      base, expon, deg=0, d, out},
	If [nvars == 0, Return[0]];
	For [j=1, j<=nvars, j++, If [SameQ[poly, vars[[j]]], Return[1]]];
	If [hed===Power,
	    If [len=!=2,
	        Message[LinearEquationsToMatrices::notpoly, poly, vars];
		Return[False]];
	    expon = poly[[2]];
	    j = polytotaldegree[expon, vars];
	    If [j=!=0,
	        Message[LinearEquationsToMatrices::notpoly, poly, vars];
		Return[False]];
	    base = poly[[1]];
	    d = polytotaldegree[base, vars];
	    If [d===False,
	        (error = True; Message[ToMatrix::notpoly, poly, vars]; Break)];
	    If [d===0, (deg = 0; Break)];
	    If [!(IntegerQ[expon] && expon>=0),
	        Message[LinearEquationsToMatrices::notpoly, poly, vars];
		Return[False]];
	    deg = d*expon;
	    ];
	If [hed===Plus,
	    For [j=1, j<=len, j++,
	        d = polytotaldegree[poly[[j]], vars];
		If [d===False,
	            Message[LinearEquationsToMatrices::notpoly, poly, vars];
		    Return[False]];
		deg = Max[deg, d];
		]
	    ];
	If [hed===Times,
	    For [j=1, j<=len, j++,
	        d = polytotaldegree[poly[[j]], vars];
		If [d===False,
	            Message[LinearEquationsToMatrices::notpoly, poly, vars];
		    Return[False]];
		deg += d;
		]
	    ];
    deg
    ]

LinearEquationsToMatrices::noteqn = "`1` is not an equation."
LinearEquationsToMatrices::notpoly = "`1` is not a polynomial in `2`."
LinearEquationsToMatrices::notlin = "`1` is not linear in `2`."

Options[LinearEquationsToMatrices] = {CheckLinearity -> False}

LinearEquationsToMatrices[oeqns_, vars_List, opts___?OptionQ] :=
	(
	{{},{}}
	) /; Length[oeqns]==0

LinearEquationsToMatrices[oeqns_, vars_List, opts___?OptionQ] :=
    Module[{eqns=oeqns, len=Length[oeqns], polys, mat, rules, j, poly, eqn,
		chk, giveResult},
      (
	    
	rules = Thread[vars -> 0];
	polys = Map[(#[[1]]-#[[2]])&, eqns];
	vec = polys /. rules;
	mat = Outer[D, polys, vars];
        {mat, -vec}
      ) /; ( chk = TrueQ[CheckLinearity /. {opts} /.
                         Options[LinearEquationsToMatrices]];
             If [Head[eqns]===And, eqns = Apply[List, eqns]];
             If [Head[eqns]=!=List, eqns = Thread[eqns]];
             If [Head[eqns]=!=List, eqns = {eqns}];
	     giveResult = True;
             If [chk,
                 For[j=1, j<=len, j++,
                     eqn = eqns[[j]];
                     If [Head[eqn]=!=Equal || Length[eqn]=!=2,
                         Message[LinearEquationsToMatrices::noteqn, eqn];
			 giveResult = False,
                         poly = eqn[[1]] - eqn[[2]];
                         If [giveResult && !PolynomialQ[poly, vars],
                	     Message[LinearEquationsToMatrices::notpoly,
				 rhs, vars];
                             giveResult = False,
                             If [giveResult && polytotaldegree[poly, vars]=!=1,
                                 Message[LinearEquationsToMatrices::notlin,
					 poly, vars];
			         giveResult = False
			     ] (* end If *)
			 ] (* end If *)
		     ] (* end If *)	
                 ] (* end For *)
             ]; (* end If *)
	     giveResult 
	   )

    ]

PolarDecomposition::nexist = "No polar decomposition exists."
PolarDecomposition::fail = "Polar decomposition failed."

PolarDecomposition[mat_?MatrixQ] :=
  Module[{giveResult, matsquare, vals, vecs, tmat, diag, inv, U, S},
    (
    U = mat . Inverse[S];
    {U, S}
    ) /; (
	  giveResult = True;
	  matsquare = Transpose[Conjugate[mat]] . mat;
          {vals, vecs} = Eigensystem[matsquare];
          tmat = Transpose[vecs];
          diag = DiagonalMatrix[Sqrt[vals]];
          If[Cases[vals, _?(#==0&)] =!= {},
	     Message[PolarDecomposition::nexist];
	     giveResult = False,
             inv = Inverse[tmat];
             If[Head[inv]===Inverse,
		Message[PolarDecomposition::fail];
 		giveResult = False,
                S = tmat . diag . inv;
                If[TrueQ[Det[S] == 0],
		   Message[PolarDecomposition::fail];
		   giveResult = False
		]
	     ]
	  ];
	  giveResult	
	 )
    ]


matnorminf[m_] := N[Max[Apply[Plus, Abs[m], 1]]]

matnorm1[m_] := matnorminf[Transpose[m]]

matnorm2[mat_?MatrixQ, prec_]:= Module[
	{symmat, len=Length[mat], x, y},
	symmat = mat.Transpose[Conjugate[mat]];
	y = Table[Random[Real, {10, 11}, prec+10], {len}];
	For [j=1, j<=5, j++,
		x = y;
		y = SetPrecision[symmat.x, prec + 10]
		];
	N[(y.Conjugate[y] / x.Conjugate[x])^(1/4)]
	]

LUPrecision[lu_] := Precision[lu[[1]]]

LUOkay[lu_] := Length[lu]==3 && lu[[3]]!=-1 && lu[[3]]!=Infinity && 
  Precision[lu[[3]]]!=Infinity 

LUSingularQ[lu_] := Length[lu]==3 && (lu[[3]]===Infinity || lu[[3]]===-1)

LUDim[lu_] := Length[lu[[2]]]

matinvnorm2[lu_] := Module[
	{len=LUDim[lu], prec=LUPrecision[lu]+10, x, y},
	x = Table[Random[Real, {10, 11}, prec], {len}];
	For [j=1, j<=5, j++,
		y = x;
		x = SetPrecision[
		LUBackSubstitution[lu, LUBackSubstitution[lu, y],
		  Transpose->True, Conjugate->True], prec]
		];
	N[(x.Conjugate[x] / y.Conjugate[y])^(1/4)]
	]

vecnorminf[vec_] := N[Max[Abs[vec]]]

vecnorminfPosition[vec_] := Module[
	{absvec=Abs[vec], len=Length[vec], pos, j},
	pos = 1;
	For [j=2, j<=len, j++, If [absvec[[j]] > absvec[[pos]], pos = j]];
	pos
	]

vecnorm1[vec_] := N[Plus @@ Abs[vec]]

matinvnorm1inf[lu_, p_:1] := Module[
	{len=LUDim[lu], x, y, z, w, oldw, j, jj, norm, oldnorm},
	x = Table[1, {jj, len}] / len;
	If [p==1, y = LUBackSubstitution[lu, x],
		y = Experimental`ExtendedLUBackSubstitution[lu, x, Transpose]];
	norm = vecnorm1[y];
	If [len==1, Return[norm]];
	w = Table[sgn = y[[jj]]; If[sgn!=0, sgn, 1], {jj, len}];
	z = Re[Experimental`ExtendedLUBackSubstitution[lu, w,
	  And[Conjugate, If[p==1, Transpose, True]]]];
	Do[	j = vecnorminfPosition[z];
		x = Table[If[k==j, 1, 0], {k, len}];
		If [p==1, y = LUBackSubstitution[lu, x],
			y = Experimental`ExtendedLUBackSubstitution[lu, x, Transpose]];
		oldnorm = norm;
		norm = vecnorm1[y];
		If [norm<=oldnorm, norm=oldnorm; Break[]];
		oldw = w;
		w = Table[sgn = y[[jj]]; If[sgn!=0, sgn, 1], {jj, len}];
		If [oldw==w, Break[]];
		z = Re[Experimental`ExtendedLUBackSubstitution[lu, w,
		  And[Conjugate, If[p==1, Transpose, True]]]];
		If [vecnorminf[z]==Abs[z[[j]]], Break[]];
	  , {3}];
	y = Table[Random[Real, {10, 11}, LUPrecision[lu]+10], {len}];
	If [p==1, y = LUBackSubstitution[lu, y],
		y = Experimental`ExtendedLUBackSubstitution[lu, y, Transpose]];
	y = y / If [p==1,vecnorm1[y],vecnorminf[y]]; (* make y have norm 1 *)
	If [p==1, x = LUBackSubstitution[lu, y],
		x = Experimental`ExtendedLUBackSubstitution[lu, y, Transpose]];
	N[Max[norm, 2*If[p==1,vecnorm1[w],vecnorminf[w]]/(3*len),
	  If[p==1,vecnorm1[x],vecnorminf[x]]]]
	]

VectorNorm[vec_?VectorQ, p_:Infinity] :=
  Module[{prec = Precision[vec], bad=Select[vec,(!NumberQ[#])&]},
   (	
	If [p==1, Return[vecnorm1[vec]]];
	If [p==Infinity, Return[vecnorminf[vec]]];
	(Apply[Plus, Abs[vec]^p])^(1/p)
   ) /; (bad==={} &&
	 If[prec === Infinity,
		Message[VectorNorm::prec];  False, True])
  ] /; (p===Infinity || (NumberQ[p] && p>=1))
	
VectorNorm::prec =
"VectorNorm has received a vector with infinite precision."
 
MatrixNorm[mat_?MatrixQ, p_:Infinity] /; (p==1 || p==2 || p==Infinity) :=
  Module[{prec=Precision[mat], bad=Select[Flatten[mat],(!NumberQ[#])&]},
   (
	If[p==1, matnorm1[mat],
           If[p==Infinity, matnorminf[mat],
               matnorm2[mat, prec]
	   ]
        ]
   ) /; (bad==={} &&
	 If[prec === Infinity,
	    Message[MatrixNorm::prec];  False, True])
  ]

MatrixNorm::prec =
"MatrixNorm has received a matrix with infinite precision."
 

InverseMatrixNorm[lu_, p_:Infinity] /; (p==1 || p==2 || p==Infinity) :=
  Module[{},
	If[p==1, matinvnorm1inf[lu, 1],
           If[p==Infinity, matinvnorm1inf[lu, Infinity],
              matinvnorm2[lu]
	   ]
        ]
  ] /; LUOkay[lu]

InverseMatrixNorm[lu_, p_:Infinity] :=
	Infinity /; ((p==1 || p==2 || p==Infinity) && LUSingularQ[lu])

InverseMatrixNorm[mat_?MatrixQ, p_:Infinity] /;
  (p==1 || p==2 || p==Infinity) :=
  Module[{bad=Select[Flatten[mat],(!NumberQ[#])&], prec=Precision[mat]},
	InverseMatrixNorm[LUDecomposition[mat],p] /;
	 (bad==={} &&
	  If[prec === Infinity,
		Message[InverseMatrixNorm::prec];  False, True])
  ]

InverseMatrixNorm::prec =
"InverseMatrixNorm has received a matrix with infinite precision."

MatrixConditionNumber[mat_?MatrixQ, lu_, p_:Infinity] /; LUOkay[lu] :=
	Module[{prec = Precision[mat]},
         (
	  MatrixNorm[mat,p] * InverseMatrixNorm[lu,p]
	 ) /; If[prec === Infinity,
		Message[MatrixConditionNumber::prec]; False, True]
	]

MatrixConditionNumber[mat_?MatrixQ,p_:Infinity] :=
	Module[{prec = Precision[mat]},
	 (
	  MatrixNorm[mat,p] * InverseMatrixNorm[mat,p]
	 ) /; If[prec === Infinity,
                Message[MatrixConditionNumber::prec]; False, True]
	]

MatrixConditionNumber::prec = 
"MatrixConditionNumber has received a matrix with infinite precision."

LUMatrices[lu_?MatrixQ]:=
  Module[{l, u, k = Length[lu]},
        u = lu Table[If[i<=j,1,0],{i,k},{j,k}];
        l = lu-u+IdentityMatrix[k];
        {l, u}
  ]


End[]

EndPackage[]

