function[F]=hesseigstru(m, itype)
% This routine computes the structural representation  of the matrix M(H),
% where  H is a unitary Hessenberg matrix H of size n=50*m.   
% The variable itype specifies the type (I,II or  III) of H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
rand('state',sum(100*clock)); 
kk=m; 
% Computing the matrix H explicitly
if (itype==1)  
    n=50*kk +1;
    sizeh(kk)=n-1;
    a=zeros(1, n); 
    a(1)=1.; 
    for k=2:n-1
        a(k)=rand *exp(sqrt(-1)*pi*2* rand);
    end
    a(n)=exp(sqrt(-1)*pi*2*rand);
    aold=a;
    n=n-1;
    h=zeros(n,n);
    for i=1:n-1
        h(i+1,i)=sqrt(1-abs(a(i+1))^2);
    end
    for i=1:n-1
        b=1;
        for j=i:n
            h(i,j)=-conj(a(i))*b*a(j+1);
            if (j<n)
                b=b*h(j+1,j);
            end
        end
    end
    h(n,n)=-conj(a(n))*a(n+1);
else 
    if (itype==2)   
        n=50*kk;
        h=zeros(n-1,n-1);
        sizeh(kk)=n-1;
        a=zeros(1,n); 
        a(1)=1.; 
        for k=2:n-1
            th=rand;
            a(k)=cos(pi*th);
            h(k, k-1)=sin(pi*th);
            be(k-1)=h(k, k-1);
        end
        a(n)=1;
        aold=a;
        n=n-1;
        for i=1:n-1
            b=1;
            for j=i:n
                h(i,j)=-conj(a(i))*b*a(j+1);
                if (j<n)
                    b=b*h(j+1,j);
                end
            end
        end
        h(n,n)=-conj(a(n))*a(n+1);
    else 
        if (itype==3) 
            n=50*kk+1;
            np=5;
            h=zeros(n-1,n-1);
            sizeh(kk)=n-1;
            a=zeros(1, n); 
            a(1)=1.; 
            for k=2:np
                a(k)=rand *exp(sqrt(-1)*pi*2/rand);
                h(k, k-1)=sqrt(1-abs(a(k))^2);
                be(k-1)=h(k, k-1);
            end
            be(np)=10^(-3);
            a(np+1)=sqrt(1-(be(np))^2);
            h(np+1, np)=be(np);
            for jj=1: 10*kk-2
                for kkk=1:np
                    a(np*jj+ kkk+1)=a(kkk+1);
                    be(np*jj+ kkk)=be(kkk);
                    h(np*jj+ kkk+1, np*jj+ kkk)=be(kkk);
                end
            end
            jj=10*kk-1;
            for kkk=1:np-1
                a(np*jj+ kkk+1)=a(kkk+1);
                be(np*jj+ kkk)=be(kkk);
                h(np*jj+ kkk+1, np*jj+ kkk)=be(kkk);
            end
            a(n)=exp(sqrt(-1)*pi*2/rand);
            aold=a;
            n=n-1;
            for i=1:n-1
                b=1;
                for j=i:n
                    h(i,j)=-conj(a(i))*b*a(j+1);
                    if (j<n)
                        b=b*h(j+1,j);
                    end
                end
            end
            h(n,n)=-conj(a(n))*a(n+1);
        end
    end
end
disp( 'size of the matrix'); 
size(h)
eigm=eig(h);
% Computing the parameters which define the Moebius transformation.
% The selection is repeated four times by choosing the best selection
for i=1:n-1
    be(i)=sqrt(1-abs(a(i+1))^2);
end
for i=1:n
    apr(i)=conj(a(i))*a(i+1) ; 
end
a=[transpose(a), transpose(a), transpose(a), transpose(a)];
for itt=1:4
    gpt=2*rand;
    gth=rand;
    gamma1(itt)=gpt*(cos(2*pi*gth) + sqrt(-1)*sin(2*pi*gth));
    dpt=2*rand;
    dth=rand;
    delta1(itt)=dpt*(cos(2*pi*dth) + sqrt(-1)*sin(2*pi*dth));
    th=rand;
    alpha1(itt)=gpt*(cos(2*pi*th) + sqrt(-1)*sin(2*pi*th));
    th=gth+th-dth;
    beta1(itt)=dpt*(cos(2*pi*th) + sqrt(-1)*sin(2*pi*th));
    gamma1(itt)=-gamma1(itt);
    th=delta1(itt);
    delta1(itt)=-alpha1(itt);
    alpha1(itt)=th;
    for i=1:n
        dd(i, itt)=-gamma1(itt)*apr(i) -delta1(itt);
    end
    for i=1:n-1
        subdd(i)=gamma1(itt)*be(i);
    end
    aaus(n+1, itt)=a(n+1, itt); 
    for k=1:n-1;
        [b,c,s]=giv(subdd(n-k), dd(n-k+1, itt));
        aa=conj(s);
        bb=c;
        av(n-k, itt)=aa;
        bv(n-k, itt)=bb;
        dd(n-k+1, itt)=b;
        pqp=be(n-k)*a(n-k+2, itt);
        pqp1=pqp*gamma1(itt)*conj(a(n-k, itt)); 
        ddnew=-aa*dd(n-k, itt)-pqp1*bb;
        spnew=bb*dd(n-k, itt)-pqp1*conj(aa);
        anew=-aa*a(n-k+1, itt)+bb*pqp;
        a(n-k+2, itt)=a(n-k+1, itt)*bb+conj(aa)*pqp;
        a(n-k+1, itt)=anew;
        aaus(n-k+1, itt)=anew; 
        dd(n-k, itt)=ddnew;
        superdd(n-k, itt)=spnew;
    end
    rmat=zeros(n,n); 
    for i=1:n 
        rmat(i, i)=dd(i,itt); 
    end
    for i=1:n-1
        rmat(i, i+1)=superdd(i,itt); 
    end
    for i=1:n-2
        b=be(i);
        for j=i+2:n
            rmat(i,j)=-gamma1(itt)*conj(aold(i))*b*a(j+1, itt);
            if (j<n)
                b=b*be(j-1);
            end
        end
    end
    [nm(itt)]=norm1est(n, aold, be, a(:, itt), dd(:, itt), superdd(:, itt), gamma1(itt)); 
end
[y,itt]=min(nm); 
for i=1:n
    ds(i)=-alpha1(itt)*apr(i) -beta1(itt);
end
for i=1:n-1
    subds(i)=alpha1(itt)*be(i);
end
for k=1:n-1;
    aa=av(n-k, itt); 
    bb=bv(n-k, itt); 
    sbd=subds(n-k);
    subds(n-k)=-aa*sbd+bb*ds(n-k+1);
    ds(n-k+1)=bb*sbd +conj(aa) *ds(n-k+1);
    pqp=alpha1(itt)*conj(aold(n-k))*be(n-k)*aaus(n-k+2, itt); 
    dsnew=-aa*ds(n-k)-pqp*bb;
    ssnew=bb*ds(n-k)-pqp*conj(aa);
    ds(n-k)=dsnew;
    superds(n-k)=ssnew;
end
dsep(1)=real(ds(1)/dd(1, itt));
v(1)=1/dd(1,  itt);
dsep(2)=real((ds(2)-subds(1)*v(1)*superdd(1, itt))/dd(2, itt));
v(2)=(conj(av(1, itt)) - v(1)*superdd(1, itt)*bv(1, itt))/dd(2, itt);
s1=-conj(aold(1))* v(1)*be(1); 
s2=s1*bv(1, itt);
dsep(3)=real((ds(3)-subds(2)*v(2)*superdd(2, itt)-gamma1(itt)*a(4, itt)*subds(2)*s2)/dd(3, itt));
s4=s2*bv(2, itt);
v(3)=(conj(av(2, itt))-v(2)*bv(2, itt)*superdd(2, itt)-gamma1(itt)*a(4, itt)*s4)/dd(3, itt);
for k=4:n
    s2=-conj(aold(k-2))* v(k-2)*be(k-2)*bv(k-2, itt)+ bv(k-2, itt)*be(k-2)*s2;  
    pqp=v(k-1)*superdd(k-1, itt); 
    pqp1=gamma1(itt)*a(k+1, itt); 
    dsep(k)=real((ds(k)-subds(k-1)*pqp-pqp1*subds(k-1)*s2)/dd(k, itt));
    if (k<n) 
        s4=s2*bv(k-1, itt);
        v(k)=(conj(av(k-1, itt))-bv(k-1, itt)*pqp-pqp1*s4)/dd(k, itt);
    end
end
%%Definition of the quasiseparable structure
u(1)=0;
t(1)=0;
for k=2:n
    u(k)=subds(k-1);
    v(k-1)=conj(v(k-1));
    tt(k)=bv(k-1, itt);
end
for j=1:n
    z(j)=0;
    w(j)=0;
end
for j=1:n-1
    t(j)=tt(j);
end
[F]=matrixex(n, dsep, u,v,t,z,w);