Code

I would like to state that because I didn't have a working Q from last week, I am starting with the 'Q code' that you posted online.

Householder and QR stuff
Week 1 and 2

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]]
]

This will make the first i-1 entries in v equal to 0

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 Stuff
Week 3

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
]

Testing

Testing the QRFactorization code

m=3;
A=RandomReal[{-2.0,2.0},{m,m}];
{Q,R}=MakeMatrixUpperTriangle[A, m];
MatrixForm[R]
MatrixForm[Q.R-A]
MatrixForm[Q.Transpose[Q]]

R comes out upper triangular, Q.r is the same as A and Q is Orth.

Testing the BackSub code

I am using the code from earlier to amke upper triangular matrices and comparing my backsolve to the output from LinearSolve. 

Here I am solving Rx = b

m=140;
A=RandomReal[{-2.0,2.0},{m,m}];
{Q,R}=MakeMatrixUpperTriangle[A, m];
b=RandomReal[{-2.,2.}, {m}];
Jenx=BackSolve[R,b];
Mmax=LinearSolve[R,b];
ListPlot[{Jenx- Mmax}]

Here I am solving Ax = b  by doing A = QR and the Rx=Q*b.

m=23;
A=RandomReal[{-2.0,2.0},{m,m}];
b=RandomReal[{-2.,2.}, {m}];
{Q,R}=MakeMatrixUpperTriangle[A, m];
Jenx=BackSolve[R,Transpose[Q].b];
Mmax=LinearSolve[A,b];
ListPlot[{Jenx, Mmax}]
Norm[Jenx- Mmax]

m=140;
A=RandomReal[{-2.0,2.0},{m,m}];
{Q,R}=MakeMatrixUpperTriangle[A, m];
b=RandomReal[{-2.,2.}, {m}];
Jenx=DoBackSolve[R,b];
Mmax=LinearSolve[R,b];
ListPlot[{Jenx- Mmax}]