(* :Title: 	NCMatMult.m // Mathematica 1.2 and 2.0 *)

(* :Author: 	Unknown. *)

(* :Context: 	NCMatMult` *)

(* :Summary:
*)

(* :Alias:
*)

(* :Warnings: 
*)

(* :History: 
   :4/23/93: Added a few usage statements.(mstankus)
   :11/8/94: Added NCMToMatMultSub, BlockDecompose and
             BlockDiagonal.(mstankus)
    07/18/97  Bill added SchurComplement and coMat
*)

BeginPackage["NCMatMult`","NonCommutativeMultiply`"];

Clear[MatMult];

MatMult::usage=
     "MatMult[x,y] (where x and y are matrices) gives the matrix \
      multiplication of x and y using NonCommutativeMultiply \
      rather than Times as Dot[] does.";

Clear[NCMToMatMult];

NCMToMatMult::usage = 
     "The definition of NCMToMatMult is \
      NCMToMatMult[expr_] := expr/.NonCommutativeMultiply->MatMult";

Clear[tpMat];

tpMat::usage =
     "Usage note for tpMat not written yet.";

Clear[ajMat];

ajMat::usage =
     "Usage note for ajMat not written yet.";


Clear[coMat];

coMat::usage =
     "coMat take complex conjugate \"co\" of each entry.";


Clear[GaussElimination];

GaussElimination::usage = 
     "Usage note for GaussElimination not written yet.";

Clear[SchurComplementTop];

SchurComplementTop ::usage = 
     "Usage note for SchurComplementTop not written yet.";

Clear[SchurComplementBtm];

SchurComplementBtm ::usage = 
     "Usage note for SchurComplementBtm not written yet.";


Clear[invMat2];

invMat2::usage = 
     "Usage note for invMat2 not written yet.";

Clear[DilationHalmos];

DilationHalmos::usage = 
     "Usage note for DilationHalmos not written yet.";

Clear[Cascade];

Cascade::usage = 
     "Usage note for Cascade not written yet.";

Clear[Chain];

Chain::usage = 
     "Usage note for Chain not written yet.";

Clear[Redheffer];

Redheffer::usage = 
     "Usage note for Redheffer not written yet.";

Clear[NCMToMatMultSub];

NCMToMatMultSub::usage = 
     "Usage note for NCMToMatMultSub not written yet.";

Clear[BlockDecompose];

BlockDecompose::usage =
     "BlockDecompose";

Clear[BlockDiagonal];

BlockDiagonal::usage = 
     "BlockDiagonal";

Clear[DeepMatMult];

DeepMatMult::usage = 
     "DeepMatMult[x,y] (where x and y are matrices) gives the matrix \
      multiplication of x and y using NonCommutativeMultiply \
      rather than Times as Dot[] does. It attempts to take \
      the product correctly if the entries are matrices so \
      that one has a block matrix.";

Begin["`Private`"];
(* -------------------------------------------------------------- *)
(*  This defines block matrix multiplication and their transpose  *)
(* -------------------------------------------------------------- *)
(*  Long comment ------------
The experienced matrix analyst should always remember that the 
Mathematia convention for handling vectors is tricky. 

         v={{1,2,4}} -- is a 1x3 matrix or a row vector

         v={{1},{2},{4}} ---is a 3x1 matrix or a row vector

         v={1,2,4}  ----is a vector but NOT A MATRIX. Indeed whether it 
                        is a row or column vector depends on the context. 
                        DON'T USE IT. DON'T USE IT. Always remember to 
                        put TWO curly brackets on your vectors or 
                        there will probably be trouble. 
End of Long comment ----- *)


MatMult[x_List,y_List]:=Inner[NonCommutativeMultiply,x,y,Plus];

MatMult[x_,y_,z__]:=MatMult[MatMult[x,y],z];

(* ------------------------------------------------------------------ *)
(*  Often it is convenient to manipulate matrices with NCM's and      *)
(*  then change to matrix multiply at the last moment.                *)
(* ------------------------------------------------------------------ *)

NCMToMatMult[expr_]:=expr/.NonCommutativeMultiply->MatMult;


(* --------------------------------------------------------------------- *)
(*  Changed line. This allows one to do  x = tpMat[u];u={{0,2},{3,4}};   *)
(* --------------------------------------------------------------------- *)
tpMat[u_]:=Transpose[Map[tp,u,{2}]] /; Length[Dimensions[u]] >=2

ajMat[u_]:=Transpose[Map[aj,u,{2}]] /; Length[Dimensions[u]] >=2

coMat[u_]:= Map[co,u,{2}] /; Length[Dimensions[u]] >=2


(* ---------------------------------------------------------------- *)
(*  This is the formula for 2 by 2 block gauss elimination          *)
(*  It assumes that a the top diagonal entry is invertible. That    *)
(*  is the default. Also we make one for the bottom                 *)
(* ---------------------------------------------------------------- *)

GaussElimination[{{a_,b_},{c_,d_}}]:= GaussElimination[{{a,b},{c,d}},top];

(* ------------------------------------------------------------------ *)
(*     left pivot, diag, right pivot                                  *)
(* ------------------------------------------------------------------ *)
GaussElimination[{{a_,b_},{c_,d_}},top]:=
         {{{Id,0},{c**inv[a],Id}},
          {{a,0},{0,d-c**inv[a]**b}},
          {{Id,inv[a]**b},{0,Id}}};

GaussElimination[{{a_,b_},{c_,d_}},btm]:=
         {{{Id,b**inv[d]},{0,Id}},
          {{a-b**inv[d]**c,0},{0,d}},
          {{Id,0},{d**c,Id}}};

(* ---------------------------------------------------------------- *)
(*   This is the formula for the 2 by 2 block matrix SchurComplement*)
(*  One  assumes that a the top diagonal entry is invertible.       *)
(*                  Also we make one for the bottom                 *)
(* ---------------------------------------------------------------- *)



SchurComplementTop[{{a_,b_},{c_,d_}}]:=    d-c**inv[a]**b;

SchurComplementBtm[{{a_,b_},{c_,d_}}]:=    a-b**inv[d]**c;








(* --------------------------------------------------------------- *)
(* This is the formula for the inverse of a 2x2 block matrix.      *)
(* --------------------------------------------------------------- *)

invMat2[mat_]:=
Block[{GE1,GE2,GE3,GE12,GE32,GE21,GE22,GE},

   GE=GaussElimination[mat];
   GE1=GE[[1]];
   GE2=GE[[2]];
   GE3=GE[[3]];


   GE12=-GE1[[2,1]];
   GE1[[2,1]]=GE12;

   GE2[[1,1]]=inv[GE2[[1,1]]];
   GE2[[2,2]]=inv[GE2[[2,2]]];

   GE32=-GE3[[1,2]];
   GE3[[1,2]]=GE32;

   invGE=MatMult[GE3,MatMult[GE2,GE1]]; 
   Return[invGE]
];
 
(* --------------------------------------------------------------- *)
(* This is the HALMOS DILATION of a. It is designed so that        *)
(* DilationHalmos[a] is a unitary operator whenever 1 - a**tp[a]   *)
(* is positive.                                                    *) 
(* --------------------------------------------------------------- *)

DilationHalmos[a_]:={{a,rt[Id-a**tp[a]]},
                     {rt[Id-tp[a]**a],-tp[a]}};


Cascade[P_,K_]:=P[[1,1]]+P[[1,2]]**K**inv[Id-P[[2,2]]**K]**P[[2,1]];

(* --------------------------------------------------------------- *)
(* frequency response functions grow from this.                    *)
(* --------------------------------------------------------------- *)

(* --------------------------------------------------------------- *)
(* The chain and Redheffer matrix for a given 2x2 block matrix P   *) 
(* --------------------------------------------------------------- *)
 
Chain[P_]:=
Block[{temp},
  temp={{0,0},{0,0}};
  temp[[1,1]]=P[[1,2]]-P[[1,1]]**inv[P[[2,1]]]**P[[2,2]];
  temp[[1,2]]=P[[1,1]]**inv[P[[2,1]]];
  temp[[2,1]]=-inv[P[[2,1]]]**P[[2,2]];
  temp[[2,2]]=inv[P[[2,1]]];
  Return[temp]
];

Redheffer[P_]:=
Block[{temp},
  temp = {{0,0},{0,0}}; 
  temp[[1,1]] = P[[1,2]]**inv[P[[2,1]]];
  temp[[1,2]] = P[[1,1]]-P[[1,2]]**inv[P[[2,2]]**P[[2,1]]];
  temp[[2,1]] = inv[P[[2,2]]];
  temp[[2,2]] = -inv[P[[2,2]]**P[[2,1]]];
  Return[temp]
];


NCMToMatMultSub[x_,rule_Rule] := 
       NCMToMatMultSub[x,{rule}];

subRule = 
{
MatrixHold[y_?NumberQ] + MatrixHold[mat_List] :>
                       MatrixHold[IdentityMatrix[Length[mat]]] +
                       MatrixHold[mat]
};

NCMToMatMultSub[x_List,rules:{___Rule}] := 
   Map[NCMToMatMultSub[#,rules]&,x];

NCMToMatMultSub[x_,rules:{___Rule}] := 
Block[{result,dummy,y}, 
   result = NCMToDummy[x,MatMultHold];
   result = result//.rules;
   result = result//.subRule;
   result = result//.MatMultHold->MyMatMult;
   Return[result];
];

MyMatMult[x_,y__] := MatMult[x,y];

MyMatMult[x_] := x;

MyMatMult[x___] := Abort[];

NCMToDummy[x_+y_,d_] := NCMToDummy[x,d] + NCMToDummy[y,d];

NCMToDummy[c_?NumberQ,d_] := d[c];

NCMToDummy[a_ b_,d_] := NCMToDummy[a,d] NCMToDummy[b,d];

NCMToDummy[a_NonCommutativeMultiply,d_] := Apply[d,a];

NCMToDummy[a_,d_] := d[a];

NCMToDummy[x___] := Abort[];

BlockDecompose[rows_,cols_] := BlockDecompose[rows,cols,"BD"];

BlockDecompose[rows_,cols_,str_] :=
Module[{i,j,result},
   result = Table[Unique[str],{i,1,rows},{j,1,cols}];
   Apply[SetNonCommutative,Flatten[result]];
   Return[result];
];
 
BlockDiagonal[rows_,rows_] := BlockDiagonal[rows,rows,"BD"];

BlockDiagonal[rows_,rows_,str_] :=
Module[{i,j,fun,result},
  fun[i_,i_] := Unique[str];
  fun[i_,j_] := 0;
  result = Table[fun[i,j],{i,1,rows},{j,1,rows}];
  Apply[SetNonCommutative,Flatten[result]];
  Return[result];
];
  

(*
ans = NCMToMatMultSub[1+c+a**b,
{
a->BlockDecompose[2,2,"a"],
b->BlockDecompose[2,2,"b"],
c->BlockDecompose[2,2,"c"]
}
];
*)

DeepMatMult[x:{{___}...},y:{{___}...}] :=
Module[{dim1,dim2,j,w,k,numberSummands,result},
   dim1 = Dimensions[x];
   dim1 = Take[dim1,{1,2}];
   dim2 = Dimensions[y];
   dim2 = Take[dim2,{1,2}];
   numberSummands = dim1[[2]];
   result = Table[Sum[DeepMatMult[x[[j,w]],y[[w,k]]],{w,1,numberSummands}]
                      ,{j,1,dim1[[1]]},{k,1,dim2[[2]]}];
   Return[result];
];

DeepMatMult[x_,y_,z__]:=DeepMatMult[DeepMatMult[x,y],z];

DeepMatMult[x_?NumberQ,y_] := x y;

DeepMatMult[x,y_?NumberQ] := x y;

DeepMatMult[x_,y_] := x**y;

DeepMatMult[x___] := BadCall["DeepMatMult",x];

(* 
Example

SetNonCommutative[a11,a12,a21,a22];
mat = {{,},{,}};
mat[[1,1]] = {{a11,a12},{a21,a22}};
mat[[1,2]] = {{0},{0}};
mat[[2,1]] = {{0,0}};
(* 
Could type
mat[[2,2]] = 1; 
OR
*)
mat[[2,2]] = {{1}};
ans = DeepMatMult[mat,mat];

Print["Type MatrixForm[mat]"];
Print["Type ans[[1,1]]"];
Print["Type ans[[1,2]]"];
Print["Type ans[[2,1]]"];
Print["Type ans[[2,2]]"];
*)

End[];
EndPackage[];
