function[]=compqr(kstart, maxit)
% This Matlab routine performs maxit tests with the
% structured QR eigenvalue algorithm applied to an initial
% companion matrix A of size m=2^(1+j), j=1, .., kstart.
% At  each QR  iteration the  current iterate  is represented
% by using approximately 9 m parameters.
errt=zeros(1, kstart);
condt=zeros(1, kstart);
numit=zeros(1,kstart);
timev=zeros(1,kstart);
for itglb=1:maxit
    disp('&&&&&&&&&&&&&&&&&&&&&&&&&');
    disp('test');
    %itglb
    disp('&&&&&&&&&&&&&&&&&&&&&&&&&');
    for nmdim=1:kstart
        disp('&&&&&&&&&&&&&&&&&&&&&&&&&');
        disp('size');
        m=2^(1+nmdim);
        m
        msize=m;
        err=zeros(m, 1);
        err1=zeros(m, 1);
        dapp=zeros(1,m);
        da=zeros(1, m);
        %dh=zeros(1,m);
        subda=ones(1,m-1);
        q=zeros(2,m-1);
        q(1,1)=-1.;
        t=zeros(2,m-1);
        t(1,m-1)=1.;
        b=zeros(2, 2*(m-2));
        for k=1:m-2
            b(1, 2*k-1)=1.;
        end
        u=zeros(m,1);
        w=zeros(m,1);
        w(m)=1;
        %%%%%%%MODIFY THIS PART TO CHANGE THE INPUT POLYNOMIAL%%%%%
        %%%%%CURRENTLY THE PROGROM COMPUTES THE ROOTS OF UNITY %%%%
        for k=1:m
            u(k)=cos(2*pi*k/m) + sqrt(-1)*sin(2*pi*k/m);
        end
        %plot(real(u), imag(u));
        eigtrue=u;
        %uu=poly(u);
        uu(1)=1.;
        for k=2:m
            uu(k)=0.;
        end
        uu(m+1)=-1.;
        u=-uu(m+1:-1:2);
        da(m)=u(m);
        u(1)=u(1)+1;
        %explicit representation of A
        A=zeros(m);
        for k=1:m-1
            A(k+1,k) =subda(k);
        end
        for k=1:m
            A(k,k)=da(k);
        end
        A(1,m)=u(1)-1.;
        for k=2:m-1
            A(k,m)=u(k);
        end
        %Atrue=A;
        % disp('condition number');
        ccond=max(condeig(A));
        % disp('norm of the matrix');
        nnorm=norm(A);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        condt(nmdim)=ccond*nnorm+condt(nmdim);
        n=m;
        fail=0;
        %timeqr=0;
        %timeqrs=0;
        % auxiliary data structures
        d=zeros(1, m); qtil=zeros(2,m-1);
        %qhat=zeros(2,1);
        util=zeros(1,m);
        g=zeros(1, m-1);
        %qaux=zeros(2,1);
        qhatn=zeros(3, m-1); f=zeros(3, 3*(m-2));
        %subdn=zeros(1, m-1);
        wtil=zeros(1,m); giv=zeros(2,m-1); ztil=zeros(3,m-1);
        hhtil=zeros(1,m); ggtil=zeros(1,m-1); giv3=zeros(2, m-3); giv2=zeros(2,m-2);
        ptilf=zeros(4,m); zhatf=zeros(4,m); phatf=zeros(4,m-1); qf=zeros(1,m-1);
        rf=zeros(1,m-2); an=zeros(1,m);  bn=zeros(1,m-1); recs=zeros(1, m-1);
        recss=zeros(1,m-2);  faipr=zeros(1, m-2);
        subdn=subda;
        tttime=cputime;
        while (n>1 && fail==0)
            iter=1;
            shstart=0;
            rshcount=0.;
            stage=1;
            alastnew= da(n);
            shift= 1. + sqrt(-1)* 1.;
            shift=shift/abs(shift);
            while ((n>=2 && abs(subdn(n-1))>eps * abs(alastnew)))
                %                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %                 %%%%%%%%%%%%%%%%%%%SHIFT STRATEGY%%%%%%%%
                %                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                if (iter>1)
                    if (abs(alastold-alastnew)<0.3 *abs(alastold))
                        shstart=1;
                    end
                end
                if (shstart==1)
                    sh=alastnew;
                    rshcount=rshcount+1;
                else
                    sh=shift;
                end
                if (rshcount>15)
                    if (stage==1)
                        sh=1.0+sqrt(-1)*2.0;
                        %sh=sh/abs(sh);
                        shstart=0.;
                        rshcount=0.;
                        stage=stage+1;
                    else
                        fail=1;
                        disp('failure');
                        break
                    end
                end
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%% FAST QR ITERATION %%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%----------------------------------%%
                %%%%% REDUCTION INTO A TRIANGULAR FORM%%%%
                %                 for k=1:n
                %                     d(k)=da(k)-sh;
                %                 end
                d(1:n)=da(1:n)-sh* ones(1,n);
                subda(1:n-1)=subda(1:n-1); 
                qhat=q(:,1);
                uhat=u(1);
                Gf=zeros(2); 
                dddf=zeros(2); 
                if (n>2)
                    for k=1:n-2
                        %[r,ae,be]=giv1(d(k),subda(k));
                        %[Gf,Yf]=planerot([d(k);  subda(k)]);
                        xf=[d(k);  subda(k)];
                        if xf(2) ~= 0
                            Yf1 = norm(xf);
                            %Gf = [xf'; -xf(2) xf(1)]/Yf1;
                            %Yf=[r; 0];
                            %dddf=diag([xf(2), -conj(xf(2))])/abs(xf(2));
                            %Gf=dddf*Gf;
                            %Gf
                            %r=Yf1; 
                            %Yf1=dddf(1,1)*Yf1
                             ccvv=(xf(2)/abs(xf(2)));
                             ae=-(conj(xf(1))/Yf1)* ccvv; 
                             be=real((conj(xf(2))/Yf1)*ccvv);  
                             Yf1=ccvv*Yf1; 
                            %x = [r; 0];
                        else
                            Gf = eye(2,class(xf));
                            Yf1=xf(1);
                            ae=-1.; 
                            be=0.; 
                        end
                        %dddf=diag([conj(Gf(, exp(sqrt(-1)*angle(Gf(1,2)))]);
                        %Gf=dddf*Gf;
                        %ae=-Gf(1,1);
                        %be=real(Gf(1,2));
                        %Yf=dddf*Yf;
                        r=Yf1;
                        giv(1,k)=ae;
                        giv(2,k)=be;
                        d(k)=r;
                        g1=transpose(qhat)*t(:, k) + uhat* conj(w(k+1));
                        g(k)=-ae*g1+be*d(k+1);
                        d(k+1)=be*g1+conj(ae)*d(k+1);
                        qaux=transpose(b(:, [2*k-1,2*k]))*qhat;
                        qtil(:, k)=-ae*qaux + be*q(:,k+1);
                        qhat=be*qaux + conj(ae)*q(:,k+1);
                        util(k)=-ae*uhat + be*u(k+1);
                        uhat =be*uhat +conj(ae)*u(k+1);
                    end
                end
                %[r,ae,be]=giv1(d(n-1),subda(n-1));
                [Gf,Yf]=planerot([d(n-1); subda(n-1)]);
                dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                Gf=dddf*Gf;
                ae=-Gf(1,1);
                be=real(Gf(1,2));
                Yf=dddf*Yf;
                r=Yf(1);
                giv(1,n-1)=ae;
                giv(2,n-1)=be;
                d(n-1)=r;
                g1=transpose(qhat)*t(:, n-1) + uhat* conj(w(n));
                g(n-1)=-ae*g1+be*d(n);
                d(n)=be*g1 + conj(ae)*d(n);
                util(n-1)=-ae*uhat + be*u(n);
                util(n) =be*uhat +conj(ae)*u(n);
                %    VERIFY THE TRIANGULARIZATION
                %                 qm=eye(n);
                %                 aprova=A([1:n], [1:n])-sh * eye(n);
                %                 for k=1:n-1
                %                     givm=[-giv(1,k), giv(2,k); giv(2,k), conj(giv(1,k))];
                %                     qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                %                             zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                %                     aprova=qma*aprova;
                %                     qm=qma*qm;
                %                 end
                %                 aprova
                %                 eig(aprova*transpose(conj(qm)) +sh *eye(n))
                %                 aprova*transpose(conj(qm)) +sh *eye(n)
                %                 norm(qm*qm'-eye(n))
                %                 rprova=zeros(n);
                %                 for k=1:n
                %                     rprova(k,k)=d(k);
                %                 end
                %                 for k=1:n-1
                %                     rprova(k, k+1)=g(k);
                %                 end
                %                 for k=1:n-2
                %                     E=eye(2);
                %                     for j=k+2:n
                %                         rprova(k,j)=transpose(qtil(:,k))*E*t(:,j-1) + util(k)*conj(w(j));
                %                         if(j<m)
                %                             E=E*b(:,[2*j-3,2*j-2]);
                %                         end
                %                     end
                %                 end
                %                 rprova
                %                 eig(rprova*transpose(conj(qm)) +sh *eye(n))
                %
                %%----------------------------------%%
                %%%%% COMPUTING A PARAMETRIZATION OF R*Q %%%%
                g1=d(1)*giv(2,1)+ g(1)*giv(1,1);
                d(1)=-conj(giv(1,1))*d(1)+ g(1)*giv(2,1);
                subdn(1)=giv(2,1)*d(2);
                d(2)=giv(1,1)*d(2);
                wtil(1)=-conj(giv(1,1))*conj(w(1)) + giv(2,1)*conj(w(2));
                what=giv(2,1)* conj(w(1))  + giv(1,1)*conj(w(2));
%                 f(:, :)=f(:,:); 
%                 qhatn(:,:)=qhatn(:,:); 
%                 qtil(:,:)=qtil(:,:); 
%                 util(:)=util(:); 
%                 giv(:,:)=giv(:,:); 
%                 d(:)=d(:); 
%                 subdn(:)=subdn(:); 
%                 w(:)=w(:); 
                if(n>3)
                    for k=2:n-2
                        g2=g1-util(k-1)*what;
                        qhatn(:, k-1)=[qtil(:,k-1); g2];
                        ztil(:,k-1)=[giv(2,k)*t(:,k); -conj(giv(1,k))];
                        f(:, 3*k-5)=[b(:,2*k-1); 0];
                        f(:, 3*k-4)=[b(:,2*k); 0];
                        f(:, 3*k-3)=[giv(1,k)*t(:,k); giv(2,k)];
                        g1=d(k)*giv(2,k)+ g(k)*giv(1,k);
                        d(k)=-conj(giv(1,k))*d(k)+ g(k)*giv(2,k);
                        subdn(k)=giv(2,k)*d(k+1);
                        d(k+1)=giv(1,k)*d(k+1);
                        wtil(k)=-conj(giv(1,k))*what+ giv(2,k)*conj(w(k+1));
                        what=giv(2,k)* what  + giv(1,k)*conj(w(k+1));
                    end
                end
                if (n>=3)
                    g2=g1-util(n-2)*what;
                    qhatn(:, n-2)=[qtil(:,n-2); g2];
                    ztil(:,n-2)=[giv(2,n-1)*t(:,n-1); -conj(giv(1,n-1))];
                    f(:, 3*n-8)=[1; 0; 0];
                    f(:, 3*n-7)=[0; 1; 0];
                    f(:, 3*n-6)=[giv(1,n-1)*t(:,n-1); giv(2,n-1)];
                    wtil(n-1)=-conj(giv(1,n-1))*what+ giv(2,n-1)*conj(w(n));
                    g1=d(n-1)*giv(2,n-1)+ g(n-1)*giv(1,n-1);
                    d(n-1)=-conj(giv(1,n-1))*d(n-1)+ g(n-1)*giv(2,n-1);
                    subdn(n-1)=giv(2,n-1)*d(n);
                    d(n)=giv(1,n-1)*d(n);
                end
                ztil(:,n-1)= [0; 0; 1];
                if(n>2)
                    wtil(n)=giv(2,n-1)* what  + giv(1,n-1)*conj(w(n));
                else
                    wtil(n)=what;
                end
                g2=g1-util(n-1)*wtil(n);
                qhatn(:, n-1)=[0; 0; g2];
                %%%% VERIFY THE PARAMETRIZATION OF R*Q
                %                 rprova*transpose(conj(qm))
                %                 eig(rprova*transpose(conj(qm)) + sh *eye(n))
                %                 s=zeros(n);
                %                 for k=1:n
                %                     s(k,k)=d(k);
                %                 end
                %                 for k=1:n-1
                %                     s(k+1,k)=subdn(k);
                %                 end
                %                 for k=1:n-1
                %                     E=eye(3);
                %                     for j=k+1:n
                %                         s(k,j)=transpose(qhatn(:,k))*E*ztil(:,j-1) + util(k)*wtil(j);
                %                         if(j<m)
                %                             E=E*f(:,[3*j-5,3*j-4, 3*j-3]);
                %                         end
                %                     end
                %                 end
                %                 s
                %%%%COMPUTE A  PARAMETRIZATION OF H%%%%%%
                %                 for k=1:n
                %                     hhtil(k)=d(k)+sh-util(k)*wtil(k);
                %                 end
                hhtil(1:n)=d(1:n)+sh*ones(1,n) -util(1:n).*wtil(1:n);
                %                 for k=1:n-1
                %                     ggtil(k)=subdn(k)-util(k+1)*wtil(k);
                %                 end
                ggtil(1:n-1)=subdn(1:n-1)-util(2:n).*wtil(1:n-1);
                %%%% VERIFY THAT H IS UNITARY%%%%%%%%%%%%%%%%%%%%%
                %                 if (n<m)
                %                     hh=zeros(n+1,n);
                %                     for k=1:n
                %                         hh(k,k)=hhtil(k);
                %                     end
                %                     for k=1:n-1
                %                         hh(k+1, k)=ggtil(k);
                %                     end
                %                     if (n>=3)
                %                         for k=3:n
                %                             for j=1:k-2
                %                                 hh(k,j)=-util(k)*wtil(j);
                %                             end
                %                         end
                %                     end
                %                     for k=1:n-1
                %                         E=eye(3);
                %                         for j=k+1:n
                %                             hh(k,j)=transpose(qhatn(:,k))*E*ztil(:,j-1);
                %                             if(j<m)
                %                                 E=E*f(:,[3*j-5,3*j-4, 3*j-3]);
                %
                %                             end
                %                         end
                %                     end
                %                     for k=1:n
                %                         hh(n+1,k)=-ucum*wtil(k);
                %                     end
                %                     disp('try the orthogonality')
                %                     hhpr=hh'*hh;
                %                     hhpr=hhpr-eye(n);
                %                     norm (hhpr)
                %                     [ypr,ipr]=max(abs(hhpr));
                %                     [yypr, iipr]=max(ypr);
                %                     yypr
                %                     ipr(iipr)
                %                     iipr
                %                 end
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %% COMPUTE THE QR FACTORIZATION OF H %%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %--------------------------------------%
                %%% REDUCTION INTO A BANDED FORM %%%%%%%
                %                 for k=1:n
                %                     u(k)=util(k);
                %                     w(k)=conj(wtil(k));
                %                 end
                u(1:n)=util(1:n);
                w(1:n)=conj(wtil(1:n));
                if (n==m)
                    uf=u(n);
                    gf=ggtil(n-1);
                    hf=hhtil(n);
                    if(n>3)
                        for k=n-1:-1:3
                            %[r,ae,be]=giv1(u(k),uf);
                            %[Gf,Yf]=planerot([u(k);uf]);
                            xf=[u(k); uf];
                            if xf(2) ~= 0
                                Yf1 = norm(xf);
%                                 Gf = [xf'; -xf(2) xf(1)]/Yf1;
%                                 %Yf=[r; 0];
%                                 dddf=diag([xf(2), -conj(xf(2))])/abs(xf(2));
%                                 Gf=dddf*Gf;
%                                 %Gf
%                                 Yf1=dddf(1,1)*Yf1;
                                %x = [r; 0];
                              ccvv=(xf(2)/abs(xf(2)));
                             ae=-(conj(xf(1))/Yf1)* ccvv; 
                             be=real((conj(xf(2))/Yf1)*ccvv);  
                             Yf1=ccvv*Yf1; 
                            else
                                Gf = eye(2,class(xf));
                                Yf1=xf(1);
                                ae=-1.; 
                                be=0.; 
                            end
                            %dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                            %  Gf=dddf*Gf;
                            %ae=-Gf(1,1);
                            %be=real(Gf(1,2));
                            % Yf=dddf*Yf;
                            r=Yf1;
                            giv3(1,k-2)=ae;
                            giv3(2,k-2)=be;
                            qf(k)=be*hhtil(k) +conj(ae)*gf;
                            uf1=uf*wtil(k-1);
                            faipr(k-1)=-uf1;
                            rf(k-1)=be*ggtil(k-1)-conj(ae)*uf1;
                            zhatf(:, k+1)=[ztil(:,k); hf];
                            hf=-ae*hhtil(k) +be *gf;
                            gf=-ae*ggtil(k-1) -be*uf1;
                            uf=r;
                            ptilf(:, k+1)=[be*qhatn(:, k); conj(ae)];
                            phatf(:,k)=[-ae*qhatn(:, k); be];
                        end
                    end
                    if(n>=3)
                        qf(2)=gf;
                        uf1=uf*wtil(1);
                        faipr(1)=-uf1;
                        rf(1)=-uf1;
                        zhatf(:, 3)=[ztil(:,2); hf];
                        ptilf(:,3)=[zeros(3,1); 1];
                    end
                    hf=-hhtil(2);
                    gf=-ggtil(1);
                    phatf(:,2)=[-qhatn(:, 2); 0];
                    qf(1)=-gf;
                    zhatf(:, 2)=[ztil(:,1); hf];
                    ptilf(:,2)=[zeros(3,1); -1];
                    phatf(:,1)=[qhatn(:, 1); 0];
                    zhatf(:, 1)=[zeros(3,1); hhtil(1)];
                    ptilf(:,1)=[zeros(3,1); 1];
                    %%%%%VERIFY THE REDUCTION%%%%%%%%%%
                    %                     pp=zeros(n);
                    %                     if(n>=3)
                    %                         for k=1:n-2
                    %                             pp(k+2,k)=rf(k);
                    %                         end
                    %                     end
                    %                     for k=1:n-1
                    %                         pp(k+1, k)=qf(k);
                    %                     end
                    %                     for k=1:n
                    %                         E=eye(4);
                    %                         for j=k:n
                    %                             pp(k,j)=transpose(ptilf(:,k))*E*zhatf(:,j);
                    %                             if(j<m)
                    %                                 if(j==1)
                    %                                     E=E*[eye(3), zeros(3,1);transpose(phatf(:,1))];
                    %                                 else
                    %                                     E=E*[f(:, [3*j-5, 3*j-4, 3*j-3]), zeros(3,1);transpose(phatf(:,j))];
                    %                                 end
                    %                             end
                    %                         end
                    %                     end
                    %                     pp
                    %                     qm=eye(n);
                    %                     aprova=hh;
                    %                     if(n>3)
                    %                         for k=n-1:-1:3
                    %                             givm=[-giv3(1,k-2), giv3(2,k-2); giv3(2,k-2), conj(giv3(1,k-2))];
                    %                             qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                    %                                     zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                    %                             aprova=qma*aprova;
                    %                             qm=qma*qm;
                    %                         end
                    %                     end
                    %                     qm1=qm;
                    %                     aprova
                    %                     aprova -pp
                    %%%%REDUCTION INTO A HESSENBERG FORM%%%%
                    an(1)=transpose(ptilf(:,1))*zhatf(:,1);
                    zz=transpose( ptilf(:,1))* [eye(3), zeros(3,1);transpose(phatf(:,1))];
                    ptilf(:,1)=transpose(zz);
                    if(n>=3)
                        for k=2:n-1
                            %[r,ae,be]=giv1(qf(k-1),rf(k-1));
                            %[Gf,Yf]=planerot([qf(k-1); rf(k-1)]);
                            xf=[qf(k-1); rf(k-1)];
                            if xf(2) ~= 0
                                Yf1 = norm(xf);
%                                 Gf = [xf'; -xf(2) xf(1)]/Yf1;
%                                 %Yf=[r; 0];
%                                 dddf=diag([xf(2), -conj(xf(2))])/abs(xf(2));
%                                 Gf=dddf*Gf;
%                                 %Gf
%                                 Yf1=dddf(1,1)*Yf1;
                                %x = [r; 0];
                             ccvv=(xf(2)/abs(xf(2)));
                             ae=-(conj(xf(1))/Yf1)* ccvv; 
                             be=real((conj(xf(2))/Yf1)*ccvv);  
                             Yf1=ccvv*Yf1; 
                            else
                                Gf = eye(2,class(xf));
                                Yf1=xf(1);
                                ae=-1.; 
                                be=0.; 
                            end
                            %  dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                            %  Gf=dddf*Gf;
                            %ae=-Gf(1,1);
                            %be=real(Gf(1,2));
                            %    Yf=dddf*Yf;
                            r=Yf1;
                            giv2(1, k-1)=ae;
                            giv2(2,k-1)=be;
                            bn(k-1)=r;
                            uf1=transpose(ptilf(:,k))*zhatf(:,k);
                            an(k)=-ae*uf1+be*qf(k);
                            qf(k)=be*uf1+conj(ae)*qf(k);
                            zz=transpose( ptilf(:,k))*[f(:, [3*k-5, 3*k-4, 3*k-3]), zeros(3,1);transpose(phatf(:,k))];
                            ptilf(:,k)=-ae*transpose(zz)+ be*ptilf(:,k+1);
                            ptilf(:,k+1)=be*transpose(zz)+ conj(ae)*ptilf(:,k+1);
                        end
                    end
                    an(n)=transpose(ptilf(:,n))*zhatf(:,n); bn(n-1)=qf(n-1);
                    %%%% VERIFY THE REDUCTION %%%
                    %                     pp=zeros(n);
                    %                     for k=1:n-1
                    %                         pp(k+1,k)=bn(k);
                    %                     end
                    %                     for k=1:n
                    %                         pp(k, k)=an(k);
                    %                     end
                    %                     for k=1:n-1
                    %                         E=eye(4);
                    %                         for j=k+1:n
                    %                             pp(k,j)=transpose(ptilf(:,k))*E*zhatf(:,j);
                    %                             if(j<m)
                    %                                 E=E*[f(:, [3*j-5, 3*j-4, 3*j-3]), zeros(3,1);transpose(phatf(:,j))];
                    %                             end
                    %                         end
                    %                     end
                    %                     pp
                    %                     %aprova
                    %                     qm=eye(n);
                    %                     if(n>=3)
                    %                         for k=2:n-1
                    %                             givm=[-giv2(1,k-1), giv2(2,k-1); giv2(2,k-1), conj(giv2(1,k-1))];
                    %                             qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                    %                                     zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                    %                             aprova=qma*aprova;
                    %                             %aprova
                    %                             qm=qma*qm;
                    %                         end
                    %                     end
                    %                     qm2=qm;
                    %                     aprova
                    %%%% REDUCTION INTO A DIAGONAL FORM %%%
                    if(n>=3)
                        for k=1:n-2
                            uf1=sqrt(abs(an(k))^2 + abs(bn(k))^2);
                            giv(1,k)=-conj(an(k))/uf1;
                            giv(2,k)=conj(bn(k))/uf1;
                            zz=transpose(ptilf(:,k))*zhatf(:,k+1);
                            an(k+1)=conj(giv(2,k))*zz+conj(giv(1,k))*an(k+1);
                            ptilf(:,k)=transpose(transpose(ptilf(:,k))* [f(:, [3*k-2, 3*k-1, 3*k]), zeros(3,1);transpose(phatf(:,k+1))]);
                            ptilf(:,k+1)=conj(giv(2,k))*ptilf(:,k) + conj(giv(1,k))* ptilf(:,k+1);
                        end
                    end
                    uf1=sqrt(abs(an(n-1))^2 + abs(bn(n-1))^2);
                    giv(1,n-1)=-conj(an(n-1))/uf1;
                    giv(2,n-1)=conj(bn(n-1))/uf1;
                    zz=transpose(ptilf(:,n-1))*zhatf(:,n);
                    an(n)=conj(giv(2,n-1))*zz+conj(giv(1,n-1))*an(n);
                    an(n)=an(n)/abs(an(n));
                    %%% VERIFY THE REDUCTION %%
                    %                     qm=eye(n);
                    %                     for k=1:n-1
                    %                         givm=[-giv(1,k), giv(2,k); conj(giv(2,k)), conj(giv(1,k))];
                    %                         qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                    %                                 zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                    %                         aprova=qma*aprova;
                    %                         %aprova
                    %                         qm=qma*qm;
                    %                     end
                    %                     qma=eye(n);
                    %                     qma(n,n)=conj(an(n));
                    %                     qm=qma*qm;
                    %                     qm3=qm;
                    %                     aprova=qma*aprova;
                    %                     aprova
                    %                     disp('pippo');
                    %                     qm3*qm2*qm1*hh
                    %verify the Hessenberg reconstruction%%%%
                    %                     pp
                    %                     qm3*pp
                    %                     pp=zeros(n);
                    %                     for k=1:n-1
                    %                         pp(k+1,k)=bn(k);
                    %                     end
                    %                     for k=1:n
                    %                         saux=1;
                    %                         for j=k:n
                    %                             if(k==1)
                    %                                 pp(k,j)=saux*an(j);
                    %                                 if(j<n)
                    %                                     saux=saux*conj(bn(j));
                    %                                 end
                    %                             else
                    %                                 pp(k,j)=-conj(an(k-1))*saux*an(j);
                    %                                 if(j<n)
                    %                                     saux=saux*conj(bn(j));
                    %                                 end
                    %                             end
                    %                         end
                    %                     end
                    %                     pp
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%% RE-COMPUTING THE QUASISEPARABLE STRUCTURE %%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %cputime
                    if (n>=3)
                        for k=1:n-2
                            b(1, 2*k-1)=giv(2,k);
                            b(1, 2*k)=0;
                            b(2, 2*k-1)=conj(giv2(1,k))*giv(1,k);
                            b(2, 2*k)=giv2(2,k);
                        end
                    end
                    if(n>2)
                        for k=1:n-2
                            recss(k)=conj(giv(2,k))*giv2(2,k);
                        end
                    end
                    recs(1)=-conj(giv(2,1))*conj(giv2(1,1));
                    if (n>=4)
                        for k=1:n-3
                            uf1=+conj(giv2(1,k))*giv(1,k)*conj(giv(1,k+1))-giv2(2,k)*conj(giv2(1,k+1))*conj(giv(2,k+1));
                            t(:,k)=[giv(2,k)*conj(giv(1,k+1)); uf1];
                            recs(k+1)=-giv2(2,k)*giv(1,k)*conj(giv(1,k+1))-giv2(1,k)*conj(giv2(1,k+1))*conj(giv(2,k+1));
                        end
                    end
                    if(n>=3)
                        uf1=+conj(giv2(1,n-2))*giv(1,n-2)*conj(giv(1,n-1))+giv2(2,n-2)*conj(giv(2,n-1));
                        t(:,n-2)=[giv(2,n-2)*conj(giv(1,n-1)); uf1];
                        recs(n-1)=-giv2(2,n-2)*giv(1,n-2)*conj(giv(1,n-1))+giv2(1,n-2)*conj(giv(2,n-1));
                    end
                    t(:,n-1)=[-giv(2,n-1)*an(n); giv(1,n-1)* an(n)];
                    q(:,1)=[-1; 0];
                    da(1)=-conj(giv(1,1));
                    da(2)=transpose([0;1])* t(:,1);
                    if(n>=3)
                        q(:,2)=transpose([0,1]*b(:, [1,2]));
                        qhat=[-giv2(2,1)*giv(1,1); giv2(1,1)];
                        if(abs(faipr(1))>0)
                            w(1)=w(1)*conj(recss(1))/conj(faipr(1));
                        end
                        if(n>=4)
                            for k=3:n-1
                                zz=transpose(qhat)*b(:, [2*k-3, 2*k-2]);
                                q(:,k)=-conj(giv3(1,k-2))*transpose(zz) + giv3(2,k-2)*[-giv2(2,k-1)*giv(1,k-1); giv2(1,k-1)];
                                da(k)=-conj(giv3(1,k-2))*transpose(qhat)* t(:,k-1)+ giv3(2,k-2)*recs(k);
                                recs(k)=giv3(2,k-2)*transpose(qhat)* t(:,k-1) + giv3(1,k-2)*recs(k);
                                aaa=recss(k-1);
                                recss(k-1)=giv3(2,k-2)*recs(k-1) + giv3(1,k-2)*recss(k-1);
                                if(abs(faipr(k-1))>0)
                                    w(k-1)=w(k-1)*conj(recss(k-1))/conj(faipr(k-1));

                                end
                                qhat=giv3(2,k-2)*transpose(zz) + giv3(1,k-2)*[-giv2(2,k-1)*giv(1,k-1); giv2(1,k-1)];
                                recs(k-1)=-conj(giv3(1,k-2))*recs(k-1) + giv3(2,k-2)* aaa;
                            end
                        end
                        da(n)=transpose(qhat)* t(:,n-1);
                    end
                    %                     for k=1:n-1
                    %                         subda(k)=recs(k)+u(k+1)*conj(w(k));
                    %                     end
                    %                      subda(1:n-1)
                    %                      recs(1:n-1)
                    %                      u(2:n)
                    %                      w(1:n-1)
                    subda(1:n-1)=recs(1:n-1) +u(2:n).*conj(transpose(w(1:n-1)));
                    %                     for k=1:n
                    %                         da(k)=da(k)+u(k)*conj(w(k));
                    %                     end
                    da(1:n)=da(1:n)+u(1:n).*conj(transpose(w(1:n)));
                else
                    %cputime
                    if (n>=3)
                        if(n==m-1)
                            % uf=ucum;
                            faipr(n-1)=-ucum*wtil(n-1);
                            %[r,ae,be]=giv1(u(n),ucum);
                            [Gf,Yf]=planerot([u(n); ucum]);
                            dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                            Gf=dddf*Gf;
                            ae=-Gf(1,1);
                            be=real(Gf(1,2));
                            Yf=dddf*Yf;
                            r=Yf(1);
                            giv3(1,n-2)=ae; giv3(2,n-2)=be;
                            uf=r;
                            gf=-ae* ggtil(n-1)- be* ucum*wtil(n-1);
                            hf=-ae*hhtil(n)-be* (ucum*wtil(n) -subda(n));
                            subgf=be*ggtil(n-1) - conj(ae)* ucum*wtil(n-1);
                            subhf=be*hhtil(n) - conj(ae)* (ucum*wtil(n)-subda(n));
                        else
                            %uf=ucum;
                            faipr(n-1)=-ucum*wtil(n-1);
                            %[r,ae,be]=giv1(u(n),ucum);
                            [Gf,Yf]=planerot([u(n); ucum]);
                            dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                            Gf=dddf*Gf;
                            ae=-Gf(1,1);
                            be=real(Gf(1,2));
                            Yf=dddf*Yf;
                            r=Yf(1);
                            giv3(1,n-2)=ae; giv3(2,n-2)=be;
                            uf=r;
                            gf=-ae* ggtil(n-1)- be* ucum*wtil(n-1);
                            hf=-ae*hhtil(n)-be* (ucum*wtil(n) -hhfaux);
                            subgf=be*ggtil(n-1) - conj(ae)* ucum*wtil(n-1);
                            subhf=be*hhtil(n) - conj(ae)* (ucum*wtil(n)-hhfaux);
                        end
                        if(n>3)
                            for k=n-1:-1:3
                                %[r,ae,be]=giv1(u(k),uf);
                                %[Gf,Yf]=planerot([u(k); uf]);
                                xf=[u(k); uf];
                                if xf(2) ~= 0
                                    Yf1 = norm(xf);
%                                     Gf = [xf'; -xf(2) xf(1)]/Yf1;
%                                     %Yf=[r; 0];
%                                     dddf=diag([xf(2), -conj(xf(2))])/abs(xf(2));
%                                     Gf=dddf*Gf;
%                                     %Gf
%                                     Yf1=dddf(1,1)*Yf1;
%                                     %x = [r; 0];
                                     ccvv=(xf(2)/abs(xf(2)));
                             ae=-(conj(xf(1))/Yf1)* ccvv; 
                             be=real((conj(xf(2))/Yf1)*ccvv);  
                             Yf1=ccvv*Yf1; 
                                else
                                    Gf = eye(2,class(xf));
                                    Yf1=-xf(1);
                                    ae=-1.; 
                                    be=0.; 
                                end
                                %                                dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                                %   Gf=dddf*Gf;
                                %ae=-Gf(1,1);
                                %be=real(Gf(1,2));
                                %  Yf=dddf*Yf;
                                r=Yf1;
                                giv3(1,k-2)=ae;
                                giv3(2,k-2)=be;
                                qf(k)=be*hhtil(k) +conj(ae)*gf;
                                uf1=uf*wtil(k-1);
                                faipr(k-1)=-uf1;
                                rf(k-1)=be*ggtil(k-1)-conj(ae)*uf1;
                                zhatf(:, k+1)=[ztil(:,k); hf];
                                hf=-ae*hhtil(k) +be *gf;
                                gf=-ae*ggtil(k-1) -be*uf1;
                                uf=r;
                                ptilf(:, k+1)=[be*qhatn(:, k); conj(ae)];
                                phatf(:,k)=[-ae*qhatn(:, k); be];
                            end
                        end
                        qf(2)=gf;
                        uf1=uf*wtil(1);
                        rf(1)=-uf1;
                        faipr(1)=-uf1;
                        zhatf(:, 3)=[ztil(:,2); hf];
                        ptilf(:,3)=[zeros(3,1); 1];
                    else
                        subgf=- ucum*wtil(n-1);
                        faipr(1)=-uf*wtil(1);
                        subhf=- ucum*wtil(n) +hhfaux;
                    end
                    hf=-hhtil(2);
                    gf=-ggtil(1);
                    phatf(:,2)=[-qhatn(:, 2); 0];
                    qf(1)=-gf;
                    zhatf(:, 2)=[ztil(:,1); hf];
                    ptilf(:,2)=[zeros(3,1); -1];
                    phatf(:,1)=[qhatn(:, 1); 0];
                    zhatf(:, 1)=[zeros(3,1); hhtil(1)];
                    ptilf(:,1)=[zeros(3,1); 1];
                    %%%%VERIFY THE REDUCTION INTO A BANDED FORM  %%%%%%%%%%%%
                    %                     pp=zeros(n);
                    %                     if(n>=3)
                    %                         for k=1:n-2
                    %                             pp(k+2,k)=rf(k);
                    %                         end
                    %                     end
                    %                     for k=1:n-1
                    %                         pp(k+1, k)=qf(k);
                    %                     end
                    %                     for k=1:n
                    %                         E=eye(4);
                    %                         for j=k:n
                    %                             pp(k,j)=transpose(ptilf(:,k))*E*zhatf(:,j);
                    %                             if(j<m)
                    %                                 if(j==1)
                    %                                     E=E*[eye(3), zeros(3,1);transpose(phatf(:,1))];
                    %                                 else
                    %                                     E=E*[f(:, [3*j-5, 3*j-4, 3*j-3]), zeros(3,1);transpose(phatf(:,j))];
                    %                                 end
                    %                             end
                    %                         end
                    %                     end
                    %                     pp
                    %                     qm=eye(n);
                    %                     aprova=hh;
                    %                     if(n>3)
                    %                         for k=n-1:-1:3
                    %                             givm=[-giv3(1,k-2), giv3(2,k-2); giv3(2,k-2), conj(giv3(1,k-2))];
                    %                             qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                    %                                     zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                    %                             aprova=qma*aprova;
                    %                             qm=qma*qm;
                    %                         end
                    %                     end
                    %                     qm1=qm;
                    %                     aprova
                    %                     aprova -pp
                    %%% REDUCTION INTO A HESSSENBERG FORM %%%%
                    an(1)=transpose(ptilf(:,1))*zhatf(:,1);
                    zz=transpose( ptilf(:,1))* [eye(3), zeros(3,1);transpose(phatf(:,1))];
                    ptilf(:,1)=transpose(zz);
                    pp1=zeros(4); 
                    qf(:)=qf(:); 
                    rf(:)=rf(:); 
%                    Gf=zeros(2); 
%                    dddf=zeros(2); 
%                     f(:,:)=f(:,:); 
%                     phatf(:,:)=phatf(:,:); 
%                     an(:)=an(:); 
%                     bn(:)=bn(:); 
                    ffm=zeros(4, 4*(n-2)); 
                    f(:,:)=f(:,:); 
                    phatf(:,:)=phatf(:,:); 
                    for k=1:n-2
                        ffm(:, [4*k-3, 4*k-2, 4*k-1, 4*k])= [f(:, [3*k-2, 3*k-1, 3*k]), zeros(3,1);transpose(phatf(:,k+1))];
                    end

                    if(n>=3)
                        for k=2:n-1
                            %[r,ae,be]=giv1(qf(k-1),rf(k-1));
                            %[Gf,Yf]=planerot([qf(k-1); rf(k-1)]);
                            xf=[qf(k-1); rf(k-1)];
                            if xf(2) ~= 0
                                Yf1 = norm(xf);
%                                 Gf = [xf'; -xf(2) xf(1)]/Yf1;
%                                 %Yf1=[r; 0];
%                                 dddf=diag([xf(2), -conj(xf(2))])/abs(xf(2));
%                                 Gf=dddf*Gf;
%                                 %Gf
%                                 Yf1=dddf(1,1)*Yf1;
%                                 %x = [r; 0];
                                ccvv=(xf(2)/abs(xf(2)));
                             ae=-(conj(xf(1))/Yf1)* ccvv; 
                             be=real((conj(xf(2))/Yf1)*ccvv);  
                             Yf1=ccvv*Yf1; 
                            else
                                Gf = eye(2,class(xf));
                                Yf1=xf(1);
                                ae=-1.; 
                                be=0.; 
                            end
                            %   dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                            %   Gf=dddf*Gf;
                            %ae=-Gf(1,1);
                            %be=real(Gf(1,2));
                            %  Yf=dddf*Yf;
                            r=Yf1;
                            giv2(1, k-1)=ae;
                            giv2(2,k-1)=be;
                            bn(k-1)=r;
                            uf1=transpose(ptilf(:,k))*zhatf(:,k);
                            an(k)=-ae*uf1+be*qf(k);
                            qf(k)=be*uf1+conj(ae)*qf(k); 
                            pp1=ffm(:, [4*k-7, 4*k-6, 4*k-5, 4*k-4]); 
                            %pp1=[f(:, [3*k-5, 3*k-4, 3*k-3]), zeros(3,1);transpose(phatf(:,k))]; 
                            zz=transpose( ptilf(:,k))*pp1;
                            ppp1=[-ae, be; be, conj(ae)]; 
                            ppp1=[transpose(zz)  ptilf(:,k+1)]*ppp1; 
                            ptilf(:,k)=ppp1( :, 1); 
                            ptilf(:,k+1)=ppp1(:, 2); 
%                             ptilf(:,k)=-ae*transpose(zz)+ be*ptilf(:,k+1);
%                             ptilf(:,k+1)=be*transpose(zz)+ conj(ae)*ptilf(:,k+1);
                        end
                    end
                    an(n)=transpose(ptilf(:,n))*zhatf(:,n); bn(n-1)=qf(n-1);
                    %[r,ae,be]=giv1(qf(n-1),subgf);
                    [Gf,Yf]=planerot([qf(n-1); subgf]);
                    dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                    Gf=dddf*Gf;
                    ae=-Gf(1,1);
                    be=real(Gf(1,2));
                    Yf=dddf*Yf;
                    r=Yf(1);
                    giv2(1, n-1)=ae; zz=an(n);
                    giv2(2,n-1)=be; an(n)=-ae*an(n)+be*subhf; bn(n-1)=r;
                    subhf=be*zz+ conj(ae)*subhf;
                    if (n<=m-2)
                        %[r,ae,be]=giv1(subhf,hhf);
                        [Gf,Yf]=planerot([subhf; hhf]);
                        dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                        Gf=dddf*Gf;
                        ae=-Gf(1,1);
                        be=real(Gf(1,2));
                        %Yf=dddf*Yf;
                        %r=Yf(1);
                        giv2m=ae;
                        subhf=-ae*subhf+be*hhf;
                    else
                        giv2m=-1;
                    end
                    %% VERIFY THE REDUCTION %%%%
                    %                     pp=zeros(n);
                    %                     for k=1:n-1
                    %                         pp(k+1,k)=bn(k);
                    %                     end
                    %                     for k=1:n
                    %                         pp(k, k)=an(k);
                    %                     end
                    %                     for k=1:n-1
                    %                         E=eye(4);
                    %                         for j=k+1:n
                    %                             pp(k,j)=transpose(ptilf(:,k))*E*zhatf(:,j);
                    %                             if(j<m)
                    %                                 E=E*[f(:, [3*j-5, 3*j-4, 3*j-3]), zeros(3,1);transpose(phatf(:,j))];
                    %                             end
                    %                         end
                    %                     end
                    %                     pp
                    %                     %aprova
                    %                     qm=eye(n);
                    %                     if(n>=3)
                    %                         for k=2:n-1
                    %                             givm=[-giv2(1,k-1), giv2(2,k-1); giv2(2,k-1), conj(giv2(1,k-1))];
                    %                             qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                    %                                     zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                    %                             aprova=qma*aprova;
                    %                             %aprova
                    %                             qm=qma*qm;
                    %                         end
                    %                     end
                    %                     qm2=qm;
                    %                     aprova
                    %%% REDUCTION INTO A DIAGONAL FORM %%%%%
                    if(n>=3)
                        pp1=zeros(4); 
                        for k=1:n-2
                            uf1=sqrt(abs(an(k))^2 + abs(bn(k))^2);
                            giv(1,k)=-conj(an(k))/uf1;
                            giv(2,k)=conj(bn(k))/uf1;
                            zz=transpose(ptilf(:,k))*zhatf(:,k+1);
                            an(k+1)=conj(giv(2,k))*zz+conj(giv(1,k))*an(k+1);
                            pp1=ffm(:, [4*k-3, 4*k-2, 4*k-1, 4*k]); 
                            %pp1=[f(:, [3*k-2, 3*k-1, 3*k]), zeros(3,1);transpose(phatf(:,k+1))]; 
                            ptilf(:,k)=transpose(transpose(ptilf(:,k))* pp1);
                            ptilf(:,k+1)=conj(giv(2,k))*ptilf(:,k) + conj(giv(1,k))* ptilf(:,k+1);
                        end
                    end
                    uf1=sqrt(abs(an(n-1))^2 + abs(bn(n-1))^2);
                    giv(1,n-1)=-conj(an(n-1))/uf1;
                    giv(2,n-1)=conj(bn(n-1))/uf1;
                    zz=transpose(ptilf(:,n-1))*zhatf(:,n);
                    an(n)=conj(giv(2,n-1))*zz+conj(giv(1,n-1))*an(n);
                    uf1=sqrt(abs(an(n))^2 + abs(subhf)^2);
                    giv(1,n)=-conj(an(n))/uf1;
                    giv(2,n)=conj(subhf)/uf1;
                    %%% VERIFY THE REDUCTION %%%%%%%%%
                    %                     qm=eye(n);
                    %                     for k=1:n-1
                    %                         givm=[-giv(1,k), giv(2,k); conj(giv(2,k)), conj(giv(1,k))];
                    %                         qma=[eye(k-1), zeros(k-1, 2), zeros(k-1, n-k-1); zeros(2,k-1), givm, zeros(2, n-k-1); ...................
                    %                                 zeros(n-k-1, k-1), zeros(n-k-1,2), eye(n-k-1, n-k-1)];
                    %                         aprova=qma*aprova;
                    %                         %aprova
                    %                         qm=qma*qm;
                    %                     end
                    %                     qma=eye(n);
                    %                     qma(n,n)=conj(an(n));
                    %                     qm=qma*qm;
                    %                     qm3=qm;
                    %                     aprova=qma*aprova;
                    %                     aprova
                    %                     disp('pippo');
                    %                     qm3*qm2*qm1*hh
                    %verify the Hessenberg reconstruction%%%%
                    %                     pp
                    %                     qm3*pp
                    %                     pp=zeros(n);
                    %                     for k=1:n-1
                    %                         pp(k+1,k)=bn(k);
                    %                     end
                    %                     for k=1:n
                    %                         saux=1;
                    %                         for j=k:n
                    %                             if(k==1)
                    %                                 pp(k,j)=saux*an(j);
                    %                                 if(j<n)
                    %                                     saux=saux*conj(bn(j));
                    %                                 end
                    %                             else
                    %                                 pp(k,j)=-conj(an(k-1))*saux*an(j);
                    %                                 if(j<n)
                    %                                     saux=saux*conj(bn(j));
                    %                                 end
                    %                             end
                    %                         end
                    %                     end
                    %                     pp
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%% RE-COMPUTING THE QUASISEPARABLE STRUCTURE %%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %cputime
                    giv(:,:)=giv(:,:); 
                    giv3(:,:)=giv3(:,:); 
                    giv2(:,:)=giv2(:,:); 
                    b(:,:)=b(:,:); 
                    t(:,:)=t(:,:); 
                    da(:)=da(:); 
                    q(:,:)=q(:,:); 
                    recs(:)=recs(:); 
                    qhat(:)=qhat(:);
                    xf=zeros(2,1); 
                    recss(:)=recss(:);
                    if (n>=3)
                        for k=1:n-2
                            b(1, 2*k-1)=giv(2,k);
                            b(1, 2*k)=0;
                            b(2, 2*k-1)=conj(giv2(1,k))*giv(1,k);
                            b(2, 2*k)=giv2(2,k);
                        end
                    end
                    if(n>=2)
                        for k=1:n-1
                            recss(k)=conj(giv(2,k))*giv2(2,k);
                        end
                    end
                    recs(1)=-conj(giv(2,1))*conj(giv2(1,1));
                    if (n>=3)
                        for k=1:n-2
                            uf1=+conj(giv2(1,k))*giv(1,k)*conj(giv(1,k+1))-giv2(2,k)*conj(giv2(1,k+1))*conj(giv(2,k+1));
                            t(:,k)=[giv(2,k)*conj(giv(1,k+1)); uf1];
                            recs(k+1)=-giv2(2,k)*giv(1,k)*conj(giv(1,k+1))-giv2(1,k)*conj(giv2(1,k+1))*conj(giv(2,k+1));
                        end
                    end
                    if(n>=2)
                        uf1=+conj(giv2(1,n-1))*giv(1,n-1)*conj(giv(1,n))-giv2(2,n-1)*conj(giv(2,n))*conj(giv2m);
                        recs(n)=-giv2(2,n-1)*giv(1,n-1)*conj(giv(1,n))-giv2(1,n-1)*conj(giv(2,n))*conj(giv2m);
                        t(:,n-1)=[giv(2,n-1)*conj(giv(1,n)); uf1];
                    end
                    q(:,1)=[-1; 0];
                    da(1)=-conj(giv(1,1));
                    da(2)=transpose([0;1])* t(:,1);
                    if(n>=3)
                        if(abs(faipr(1))>0)
                            w(1)=w(1)*conj(recss(1))/conj(faipr(1));
                        end
                        q(:,2)=transpose([0,1]*b(:, [1,2]));
                        qhat=[-giv2(2,1)*giv(1,1); giv2(1,1)];
                        for k=3:n
                            zz=transpose(qhat)*b(:, [2*k-3, 2*k-2]);
                            xf=[-giv2(2,k-1)*giv(1,k-1); giv2(1,k-1)]; 
                            qqq=transpose(qhat)* t(:,k-1); 
                            sspp=[zz, qqq; transpose(xf), recs(k)]; 
                            sspp=[-conj(giv3(1,k-2)), giv3(2,k-2)]*sspp; 
                            q(:,k)=transpose(sspp(:, 1:2)); 
                            da(k)=sspp(:,3); 
%                             q(:,k)=-conj(giv3(1,k-2))*transpose(zz) + giv3(2,k-2)*xf;
%                             da(k)=-conj(giv3(1,k-2))*transpose(qhat)* t(:,k-1)+ giv3(2,k-2)*recs(k);
                            aaa=recss(k-1);
                            sspp=[recs(k), aaa; qqq, recs(k-1)]; 
                            sspp=transpose(giv3(:, k-2))*sspp; 
                            recs(k)=sspp(:,1); 
                            recss(k-1)=sspp(:,2); 
%                             recs(k)=giv3(2,k-2)*transpose(qhat)* t(:,k-1) + giv3(1,k-2)*recs(k);
%                             aaa=recss(k-1);
%                             recss(k-1)=giv3(2,k-2)*recs(k-1) + giv3(1,k-2)*recss(k-1);
                            if(abs(faipr(k-1))>0)
                                w(k-1)=w(k-1)*conj(recss(k-1))/conj(faipr(k-1));
                            end
                            xf=[-giv2(2,k-1)*giv(1,k-1); giv2(1,k-1)]; 
                            qhat=[xf transpose(zz)]*giv3(:,k-2);
                            %qhat=transpose(qhat); 
                            %qhat=giv3(2,k-2)*transpose(zz) + giv3(1,k-2)*xf;
                            recs(k-1)=-conj(giv3(1,k-2))*recs(k-1) + giv3(2,k-2)* aaa;
                        end
                    end
                    %                     for k=1:n-1
                    %                         subda(k)=recs(k)+u(k+1)*conj(w(k));
                    %                     end
                    subda(1:n-1)=recs(1:n-1) +u(2:n).*conj(transpose(w(1:n-1)));
                    %                     for k=1:n
                    %                         da(k)=da(k)+u(k)*conj(w(k));
                    %                     end
                    da(1:n)=da(1:n)+u(1:n).*conj(transpose(w(1:n)));
                end
                iter=iter+1;
                alastold=alastnew;
                alastnew=da(n);
            end
            dapp(n)=alastnew;
            %disp('eigenvalue is found');
            %dapp(n)
            %cputime
            if (n==m)
                ucum=u(n);
            else
                %[r,ae,be]=giv1(u(n), ucum);
                [Gf,Yf]=planerot([u(n); ucum]);
                dddf=diag([exp(-sqrt(-1)*angle(Gf(1,2))), exp(sqrt(-1)*angle(Gf(1,2)))]);
                Gf=dddf*Gf;
                ae=-Gf(1,1);
                be=real(Gf(1,2));
                Yf=dddf*Yf;
                r=Yf(1);
                ucum=r;
                hhf=be*subda(n);
                hhfaux=-ae*subda(n);
            end
            err=min(abs(eigtrue-dapp(n)*ones(m,1)));
            %disp('error')
            %err
            %disp('number')
            %n
            n=n-1;
            numit(nmdim)=numit(nmdim)+iter;
        end
        %cputime
        dapp(1)=da(1);
        %disp('maximum error');
        for j=1:m
            err(j)=min(abs(eigtrue-dapp(j)*ones(m,1)));
            err1(j)=min(abs(transpose(dapp)-eigtrue(j)*ones(m,1)));
        end
        errs=max(max(err), max(err1));
        %err
        timev(nmdim)=cputime-tttime;
        timev(nmdim)
        errt(nmdim)=errs+errt(nmdim);
        numit(nmdim)=numit(nmdim)/msize;
    end
end
errt=errt/maxit;
condt=condt/maxit;
timev=timev/maxit;
for k=1:kstart
    condt(k)=max(condt(k),1);
end
condt=1./condt;
disp('errt')
for k=1:kstart
    errt(k)
end
disp('condt')
for k=1:kstart
    condt(k)
end
disp('time')
for k=1:kstart
    timev(k)
end
numit