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(j3) 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=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(j3) 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(j3) % 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=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=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(j3) % 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=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=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