(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



QFun[v_,i_]:=Module[
{m,w,Q},
m=Length[v];
w=v;
w[[i]]+=Norm[v];
w=w/Norm[v];
Q=IdentityMatrix[m]-2/w.Conjugate[w] Outer[Times,w,Conjugate[w]]
]


ZapEntries[v_,i_] := Module[
{w,j},
w=v;
For[j=1, j<i, j++,
	w[[j]]=0.0;
];
w
]


MakeMatrixUpperTriangle[A_, m_]:=Module[
{x, Q, R, Qinside, i},
(*Initialize the variables that will end up being output*)
R=A;
Q=IdentityMatrix[m];
(*Go through each column and make everything below the diagonal 0*)
For[i=1, i<m,i++,
x=R[[All,i]];
(*ZapEntries kills the top of the vector x*)
x=ZapEntries[x, i];
Qinside= QFun[x, i];
(*Accumulate the effects of the operations*)
R=Qinside.R;
Q=Q.Transpose[Qinside];
];
{Q,R}
]


BackSolve[R_, b_]:=Module[
{stuff, x, i,j,m,n},
(* *)
{m,n}=Dimensions[R];
If[ n!= m, Return[Print["Non square matrix!"]]];
If[ n!= Length[b], Return[Print["RHS wrong size!"]]];
m = Length[b];
x=Table[0,{m}];
(*Loop through the rows (i) and then the columns (j) to solve for each x*)
For[i=m, i>=1,i--,
stuff=0.0;
For[j=m, j>i,j--,
stuff =stuff +(R[[i,j]]x[[j]]);
];
x[[i]]=(b[[i]]-stuff)/R[[i,i]];
];
(*Return x*)
x
]


DoBackSolve[R_, b_]:=Module[
{stuff, x, i,j,m,n},
(* *)
{m,n}=Dimensions[R];
If[ n!= m, Return[Print["Non square matrix!"]]];
If[ n!= Length[b], Return[Print["RHS wrong size!"]]];
m = Length[b];
x=Table[0,{m}];

(*Loop through the rows (i) and then the columns (j) to solve for each x*)
Do[stuff=0.0;
Do[stuff =stuff +(R[[i,j]]x[[j]]),
{j,m,i+1,-1}
];
x[[i]]=(b[[i]]-stuff)/R[[i,i]],
{i,m,1,-1}
];
(*Return x*)
x
]



