Off[General::spell1];
Clear[Hessenberg,a0];
Hessenberg[a0_] := 
Module[{a, b,n,k,c,beta,y,v,vvt,u}, 
 a = a0;
 n = Length[a];
    norm[V_] :=N[ Sqrt[V.Conjugate[V]/.List->Plus], digit];
    b = Transpose[a];
    qmatrix = IdentityMatrix[n];
 Clear[k];
 For[k=1, k <= n-2, k++,
   c=Table[If[i>=k+1, b[[k,i]],0], {i,1,n}]; (*k-th column*)
   beta= If[b[[k,k]] == 0, 1,   
   N[ -(b[[k,k]]/ Abs[b[[k,k]]]) norm[c],digit] ];
   y =  c - beta Table[If[i != k+1,0,1], {i,1,n}];
   alpha = N[Sqrt[2] / norm[y],digit];
   v = alpha y;
   vvt= Table[v[[i]] Conjugate[v[[j]]],{i,n},{j,n}];
   u = IdentityMatrix[n]-vvt;  
   b = Simplify[u.b.Transpose[u]];
   qmatrix = u.qmatrix;    
 (*Print["R",k,"=" MatrixForm[Chop[Transpose[b]]]];*)
   ]; 
a=Transpose[b];
Return[a]; ];
On[General::spell1];

digit = 8;

Clear[A];
A={{1,2,3,4},{4,5,6,7},{2,1,5,0},{4,2,1,0}}; 
MatrixForm[A]

Clear[H];
H=Hessenberg[A];  
MatrixForm[Chop[H]]

Clear[U];
U=Transpose[qmatrix];  
MatrixForm[U]

MatrixForm[Chop[U.H.Transpose[U]]]

Off[General::spell1]; 
Clear[QRmethod,a0,Steps0];
QRmethod[a0_,Steps0_] := 
Module[{a, n,k,R,Q,M,z}, 
 a = a0;
 M = Steps0;
   a = Hessenberg[a];
   For[j=1, j<= M, j++,
   R = QRfact[a];
   Q = Transpose[q]; Clear[q];
   a = R.Q ;
  (* Print["A",k,"=" MatrixForm[Chop[a]]]; Print["    "]; *)
   ];
 Return[a]; ];
On[General::spell1];

Lambda= QRmethod[A,10];
MatrixForm[N[Chop[Lambda],6]]

Off[General::spell1];
Clear[QRshhift,a0,M0];
QRshift[a0_,M0_] := 
Module[{a, b,n,m,k,R,Q,z},  
 a = a0;
 M = M0;
 n = Length[a];
 lambda = Table[0, {i,1,n}];
   b = Hessenberg[a];
   For[k=1, k <= M, k++,
     m = Length[b]; 
     z= b[[m,m]];
     b = b -z IdentityMatrix[m]; 
     R = QRfact[b];
     Q = Transpose[q];Clear[q];
     b = R.Q + z IdentityMatrix[m];
   If[Abs[b[[m,m-1]]] <= 10^(-digit), 
     lambda[[m]] = b[[m,m]]; 
     b = Table[b[[i,j]], {i,1,m-1},{j,1,m-1}], 
     b = b];
   (Print["A",k,"=" MatrixForm[N[Chop[b],4]]]; 
        Print["    "];)
     ];  
    a=Table[ If[i<= Length[b]&& j <= Length[b], b[[i,j]],
       If[i ==j, lambda[[i]], 0]] , {i,1,n},{j,1,n}  ] ; 
    Return[a]; ];
On[General::spell1];

Clear[Lambda];  
Lambda= QRshift[A,10];
MatrixForm[N[Chop[Lambda],6]]

Clear[A];
A={{190,66,-84,30},{66,303,42,-36},{336,-168,147,-112},
  {30,-36,28,291}};  
MatrixForm[A]

Clear[Lambda];  
Lambda= QRshift[A,10];
MatrixForm[N[Chop[Lambda],6]]

remain=Table[Lambda[[i,j]], {i,1,2},{j,1,2}]; 
Eigenvalues[remain]