@ --------------------------- CoRRAM_4.src --------------------------

  Alfred Maußner
  15 September 2006
  Last modification: 28 December 2010

  This file collects several routines that compute the matrices of the linear part
  of the policy rules that approximate the solution of DSGE models.
  
  List of procedures:

    SolveLA
    SolveLA2
    SolveLA3
    SolveLA4
    MySchur
    MyGSchur
    SchurS
    L2Norm
    InterChange
    Givens
    
-------------------------------------------------------------- @

/* SolveLA:
**
** usage: {Lxx,Lxz,Llx,Llz,Lux,Luz}=SolveLA(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,rho,nx)
**
** Input:   the matrices from (2.36), where Pi=rho, and nx=n(x)
** Output:  the policy matrices explained in Section 2.2.2
*/

proc (6)=SolveLA(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,rho,nx);

   local Lxx,Lxz,Llx,Llz,Lux,Luz,Cui,temp,w,r,v,vi,t,
         toli,lambda,q1,q2,phi,i,j,tll,vlxi,
         vlli,vllii,itest,nl,nz,nu,phi2, phi3;

   clear Lxx,Lxz,Llx,Llz,Lux,Luz;

    /* Tolerence in case of numerical imprecision producing
    ** complex policy matrices   */
    toli=1e-10;

    /* get dimension */
 
   nl=cols(Dxl)-nx;
   if ismiss(Dz) or ismiss(Fz) or ismiss(rho);
     nz=0; 
   else;
     nz=cols(Cz);
   endif;

/* Compute the dynamic equation
**
** E_t(x[t+1]|l[t+1]) = w*(x[t]|l[t]) + r*Et(z[t+1]) + q*z[t]
**
*/

if ((rows(Cu)==1) and (cols(Cu)==1));
   if (Cu==0);
       Cui=0;
       nu=0;
   else;
       Cui=1/Cu;
       nu=cols(Cu);
   endif;
else;
   nu=cols(Cu);
   Cui=inv(Cu);
endif;

temp=inv(Dxl-Du*Cui*Cxl);

w=-temp*(Fxl-Fu*Cui*Cxl);
r=temp*( (Dz+Du*Cui*Cz)*Rho + (Fz+Fu*Cui*Cz) );

/*  Schur decomposition of w, test of solvability, and separation of variables */

{t,v}=SchurS(w); @ was SchurO, but now with different sorting mechanism @
output on;
?"Eigenvalues of W:";
lambda=diag(t);
?lambda[sortind(abs(lambda))];
output off;
lambda=abs(lambda);
if nx>0;
    if sumc(abs(lambda[1:nx]) .> 1)>0;
        ?"Not all of the " ftos(nx,"%*.*lf",2,0)
        " eigenvalues are within the unit circle!";
        ?lambda;?"";
        ?"The program terminates! Press any key!";wait;stop;
        retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
    endif;
endif;

if sumc( abs(lambda[nx+1:rows(lambda)]) .< 1) > 0;
  ?"Not all of the " ftos(nl,"%*.*lf",2,0)
   " eigenvalues are out of the unit circle!"; 
  ?lambda;?"";
  ?"The program terminates. Press any key!";wait;stop;
  retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;

vi=v';
if nx>0;
    vlxi=vi[nx+1:rows(vi),1:nx];
endif;
vlli=vi[nx+1:rows(vi),nx+1:cols(vi)];
vllii=inv(vlli);

temp=vi*r;
if nx>0;
    q1=temp[1:nx,.];
endif;
q2=temp[nx+1:rows(temp),.];

/* Solution for lambda */
if nz>0;
    phi=zeros(nl,nz);
    tll=t[nx+1:rows(t),nx+1:cols(t)];
    i=nl;
    do while i>0;
        j=i+1;
        temp=zeros(1,nz);
            do until j>nl;
              temp=temp+tll[i,j]*phi[j,.];
              j=j+1;
           endo;
        phi[i,.]=(q2[i,.]+temp)*inv((Rho-tll[i,i]*eye(nz)));
        i=i-1;
    endo;
    /* alternative code to compute phi
    **    phi2=rho';
    **   phi2=phi2.*.eye(nl);
    **   phi2=phi2-(eye(nz).*tll);
    **   phi2=inv(phi2);
    **   phi2=phi2*vec(q2);
    **   phi2=reshape(phi2',nl,nz);
    */
    Llz=vllii*phi;
else;
    Llz=zeros(nl,1);
endif;
if nx>0;
    Llx=-vllii*vlxi;
endif;

if nx>0;
    Lxx=v[1:nx,1:nx]*t[1:nx,1:nx]*inv(v[1:nx,1:nx]);
    if nz>0;  Lxz=w[1:nx,nx+1:cols(w)]*Llz+r[1:nx,.]; else; Lxz=zeros(nx,1); endif;
else;
    Lxx=0;
    Lxz=0;
endif;
if nx>0;
   if nu>0; Lux=Cui*Cxl*(eye(nx)|Llx); else; Lux=0; endif;
else;
    Lux=0;
endif;
if nx>0;
    if nz>0;
         Luz=Cui*Cxl*(zeros(nx,nz)|Llz)+Cui*Cz;
    else;
        if nu>0; Luz=zeros(nu,1); else; Luz=0; endif;
    endif;
else;
    Luz=Cui*Cxl*Llz+Cui*Cz;
endif;

/* Test for complex coefficients */
itest=maxc(maxc(imag(Llx)));

if itest > toli;
   ?"Warning:There are complex coefficients in Llx";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
else;
   Llx=real(Llx);
endif;

itest=maxc(maxc(imag(Llz)));

if itest > toli;
   ?"Warning:There are complex coefficients in Llz";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
else;
   Llz=real(Llz);
endif;

itest=maxc(maxc(imag(Lxx)));

if itest > toli;
   ?"Warning:There are complex coefficients in Lxx";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
else;
   Lxx=real(Lxx);
endif;

itest=maxc(maxc(imag(Lxz)));

if itest > toli;
   ?"Warning:There are complex coefficients in Lxz";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
else;
   Lxz=real(Lxz);
endif;

itest=maxc(maxc(imag(Lux)));

if itest > toli;
   ?"Warning:There are complex coefficients in Lux";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
else;
   Lux=real(Lux);
endif;

itest=maxc(maxc(imag(Luz)));

if itest > toli;
   ?"Warning:There are complex coefficients in Luz";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
else;
   Luz=real(Luz);
endif;


retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endp;


/* SolveLA2:
**
** usage: {Lxx,Lxz,Llx,Llz,Lux,Luz}=SolveLA2(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,rho,nx)
**
** Input:   the matrices from (2.36), where Pi=rho, and nx=n(x)
** Output:  the policy matrices explained in Section 2.2.2
*/

proc (6)=SolveLA2(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,rho,nx);

   local Lxx,Lxz,Llx,Llz,Lux,Luz,Cui,temp,w,r,v,vi,t,
         toli,lambda,q1,q2,phi,i,j,tll,vlxi,
         vlli,vllii,itest,nl,nz;

   clear Lxx,Lxz,Llx,Llz,Lux,Luz;

    /* Tolerence in case of numerical imprecision producing
    ** complex policy matrices   */
    toli=1e-10;

    /* get dimension */
 
   nl=cols(Dxl)-nx;
   nz=cols(Cz);

/* Compute the dynamic equation
**
** E_t(x[t+1]|l[t+1]) = w*(x[t]|l[t]) + r*Et(z[t+1]) + q*z[t]
**
*/

if ((rows(Cu)==1) and (cols(Cu)==1));
   if (Cu==0);
       Cui=0;
   else;
       Cui=1/Cu;
   endif;
else;
   Cui=inv(Cu);
endif;

temp=inv(Dxl-Du*Cui*Cxl);

w=-temp*(Fxl-Fu*Cui*Cxl);
r=temp*( (Dz+Du*Cui*Cz)*Rho + (Fz+Fu*Cui*Cz) );

/*  Schur decomposition of w, test of solvability, and separation of variables */
{t,v}=MySchur(w);
output on;
?"Eigenvalues of W:";
lambda=diag(t);
?lambda[sortind(abs(lambda))];
output off;
lambda=abs(lambda);
if sumc(abs(lambda[1:nx]) .> 1)>0;
  ?"Not all of the " ftos(nx,"%*.*lf",2,0)
   " eigenvalues are within the unit circle!";
  ?lambda;?"";
  ?"The program terminates! Press any key!";wait;stop;
  retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;

if sumc( abs(lambda[nx+1:rows(lambda)]) .< 1) > 0;
  ?"Not all of the " ftos(nl,"%*.*lf",2,0)
   " eigenvalues are out of the unit circle!"; 
  ?lambda;?"";
  ?"The program terminates. Press any key!";wait;stop;
  retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;

vi=v';
vlxi=vi[nx+1:rows(vi),1:nx];
vlli=vi[nx+1:rows(vi),nx+1:cols(vi)];
vllii=inv(vlli);

temp=vi*r;
q1=temp[1:nx,.];
q2=temp[nx+1:rows(temp),.];

/* Solution for lambda */
phi=zeros(nl,nz);
tll=t[nx+1:rows(t),nx+1:cols(t)];
i=nl;
do while i>0;
   j=i+1;
   temp=zeros(1,nz);
   do until j>nl;
      temp=temp+tll[i,j]*phi[j,.];
         j=j+1;
   endo;
   phi[i,.]=(q2[i,.]+temp)*inv((Rho-tll[i,i]*eye(nz)));
          i=i-1;
endo;

Llx=-vllii*vlxi;
Llz=vllii*phi;
Lxx=v[1:nx,1:nx]*t[1:nx,1:nx]*inv(v[1:nx,1:nx]);
Lxz=w[1:nx,nx+1:cols(w)]*Llz+r[1:nx,.];

/* Test for complex coefficients */
itest=maxc(maxc(imag(Llx)));

if itest > toli;
   ?"Warning:There are complex coefficients in Llx";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Llx=real(Llx);

itest=maxc(maxc(imag(Llz)));

if itest > toli;
   ?"Warning:There are complex coefficients in Llz";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Llz=real(Llz);

itest=maxc(maxc(imag(Lxx)));

if itest > toli;
   ?"Warning:There are complex coefficients in Lxx";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Lxx=real(Lxx);

itest=maxc(maxc(imag(Lxz)));

if itest > toli;
   ?"Warning:There are complex coefficients in Lxz";
   ?"itest=" ftos(itest,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Lxz=real(Lxz);

/* Compute Lux and Luz */
Lux=Cui*Cxl*(eye(nx)|Llx);
Luz=Cui*Cxl*(zeros(nx,nz)|Llz)+Cui*Cz;


retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endp;

/* SolveLA3:
**
** usage: {Lxx,Lxz,Llx,Llz,Lux,Luz}=SolveLA3(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,rho,nx)
**
** Input:   the matrices from (2.36), where Pi=rho, and nx=n(x)
** Output:  the policy matrices explained in Section 2.2.2
**
** Remarks: Different from SolveLA2 it is not necessary
**          that B=(Dxl-Du*Cui*Cxl) is invertible.
**          Here we use the generalized Schur factorization
**          to get new variables in which the model is solved.
**
*/

proc (6)=SolveLA3(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,rho,nx);

   local Lxx,Lxz,Llx,Llz,Lux,Luz,Cui,
         nl, nz, nxl, vt,
         a, b, c, d, u, v, s, t, lambda, test, toli,        
         phi, temp1, temp2, dx, dl, sll, tll, vxxi, sxxi;

   clear Lxx,Lxz,Llx,Llz,Lux,Luz;

    /* Tolerence in case of numerical imprecision producing
    ** complex policy matrices   */
    toli=1e-10;

    /* get dimension */
 
   nxl=cols(Dxl);
   nl=nxl-nx;
   nz=cols(Cz);
   
/* Compute the dynamic equation
**
** B E_t(x[t+1]|l[t+1]) = A*(x[t]|l[t]) + C*z[t]
**
*/

if ((rows(Cu)==1) and (cols(Cu)==1));
   if (Cu==0);
       Cui=0;
   else;
       Cui=1/Cu;
   endif;
else;
   Cui=inv(Cu);
endif;

b=Dxl-Du*Cui*Cxl;
a=-(Fxl-Fu*Cui*Cxl);
c=(Dz+Du*Cui*Cz)*Rho + (Fz+Fu*Cui*Cz);

/*  generalized Schur decomposition of (a,b) */
{u,v,s,t}=MyGSchur(a,b);
output on;
?"Eigenvalues of W:";
test=0;
for i (1,nxl,1);
  if abs(S[i,i]) /= 0;
        lambda=t[i,i]/s[i,i];
        ?abs(lambda);
        if (abs(lambda)<1.0);
            test=test+1;
        endif;
  else;
        ?"Infinite Eigenvalue"; @ this will happen, if B is not invertible @
  endif;
endfor;
output off;
if test<nx;
    ?"Not all of the " ftos(nx,"%*.*lf",2,0)
     " eigenvalues are within the unit circle!";
    ?"The program terminates! Press any key!";wait;stop;
    retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;
if test>nx;
    ?"Not all of the " ftos(nl,"%*.*lf",2,0)
     " eigenvalues are out of the unit circle!";   
    ?"The program terminates. Press any key!";wait;stop;
    retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;

/* Compute phi */
Temp1=(u')*c;
dl=Temp1[nx+1:nx+nl,1:nz];
dx=Temp1[1:nx,1:nz];
sll=s[nx+1:nx+nl,nx+1:nx+nl];
tll=t[nx+1:nx+nl,nx+1:nx+nl];
phi=zeros(nl,nz);

	
for i (nl,1,-1);
    Temp2=zeros(1,nz);
    for j (i+1,nl,1);
    	Temp2[1,1:nz]=Temp2[1,1:nz]+(Tll[i,j]*Phi[j,1:nz]-Sll[i,j]*(Phi[j,1:nz]*Rho));
	endfor;
	Temp1=Sll[i,i]*Rho-Tll[i,i]*Eye(nz);
	Temp1=inv(Temp1);
	Phi[i,1:nz]=(Dl[i,1:nz]+Temp2[1,1:nz])*Temp1;
        
endfor;
Vxxi=V[1:nx,1:nx];
Vxxi=inv(Vxxi);
Temp1=V[nx+1:nxl,nx+1:nxl];
Temp1=Temp1-(V[nx+1:nxl,1:nx]*(Vxxi*V[1:nx,nx+1:nxl]));
Llz=Temp1*Phi; @ the shock part @
Llx=V[nx+1:nxl,1:nx]*Vxxi; @ the x-part @
Sxxi=inv(S[1:nx,1:nx]);

/* The policy function for xt+1 */
Lxx=(V[1:nx,1:nx]*Sxxi)*(T[1:nx,1:nx]*Vxxi);
Lxz=V[1:nx,nx+1:nxl]*(Phi*Rho); @ Vxl*Phi*Rho @
Temp1=Dx;
Temp1=Temp1+(T[1:nx,nx+1:nxl]*Phi); 
Temp1=Temp1-(S[1:nx,nx+1:nxl]*(Phi*Rho));
Temp1=Temp1-((T[1:nx,1:nx]*Vxxi)*(V[1:nx,nx+1:nxl]*Phi));
Lxz=Lxz+((V[1:nx,1:nx]*Sxxi)*Temp1);

/* Test for complex coefficients */
test=maxc(maxc(imag(Llx)));

if test > toli;
   ?"Warning:There are complex coefficients in Llx";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Llx=real(Llx);

test=maxc(maxc(imag(Llz)));

if test > toli;
   ?"Warning:There are complex coefficients in Llz";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Llz=real(Llz);


test=maxc(maxc(imag(Lxx)));

if test > toli;
   ?"Warning:There are complex coefficients in Lxx";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Lxx=real(Lxx);

test=maxc(maxc(imag(Lxz)));

if test > toli;
   ?"Warning:There are complex coefficients in Lxz";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Lxz=real(Lxz);

/* Compute Lux and Luz */

Lux=Cui*Cxl*(eye(nx)|Llx);
Luz=Cui*Cxl*(zeros(nx,nz)|Llz)+Cui*Cz;


retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endp;

/* SolveLA4:
**
** usage: {Lxx,Lxz,Llx,Llz}=SolveLA4(A,B,C,rho,nx)
**
** Input:   the matrices from (2) and (1c), respectively,
**          where Pi=rho, and nx=n(x)
** Output:  the policy matrices for the states (Lxx and Lxz) and
**          for the non predetermined variables in the vector lambda
**          (Llx and Llz).
**
**
*/

proc (4)=SolveLA4(A,B,C,Rho,nx);

   local Lxx,Lxz,Llx,Llz,Lux,Luz,Cui,
         nl, nz, nxl, vt,
         d, u, v, s, t, lambda, test, toli,        
         phi, temp1, temp2, dx, dl, sll, tll, vxxi, sxxi;

   clear Lxx,Lxz,Llx,Llz;

    /* Tolerence in case of numerical imprecision producing
    ** complex policy matrices   */
    toli=1e-10;

    /* get dimension */
 
   nxl=cols(B);
   nl=nxl-nx;
   nz=cols(rho);
  

/*  generalized Schur decomposition of (a,b) */
{u,v,s,t}=MyGSchur(a,b);

@ test factorization @
test=(u')*b*v;
test=test-s;
test=maxc(maxc(abs(test)));
output on; ?"Test of factorization: " ftos(test,"*.*lf",12,9); output off;
if test>1.e-6;
    ?"Incorrect factorization of b. The solution may be (very) missleading! Press any key to continue."; wait;
endif;
test=(u')*a*v-t;
if (maxc(maxc(abs(test))))>1.e-6;
    ?"Incorrect factorization of a. The solution may be (very) missleading! Press any key to continue."; wait;
endif;
output on;
?"Eigenvalues of W:";
test=0;
for i (1,nxl,1);
  if abs(S[i,i]) /= 0;
        lambda=t[i,i]/s[i,i];
        ?abs(lambda);
        if (abs(lambda)<1.0);
            test=test+1;
        endif;
  else;
        ?"Infinite Eigenvalue"; @ this will happen, if B is not invertible @
  endif;
endfor;
output off;

if test<nx;
    ?"Not all of the " ftos(nx,"%*.*lf",2,0)
     " eigenvalues are within the unit circle!";
    ?"The program terminates! Press any key!";wait;stop;
    retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;
if test>nx;
    ?"Not all of the " ftos(nl,"%*.*lf",2,0)
     " eigenvalues are out of the unit circle!";   
    ?"The program terminates. Press any key!";wait;stop;
    retp(Lxx,Lxz,Llx,Llz,Lux,Luz);
endif;

/* Compute phi */
Temp1=(u')*c;
dl=Temp1[nx+1:nx+nl,1:nz];
dx=Temp1[1:nx,1:nz];
sll=s[nx+1:nx+nl,nx+1:nx+nl];
tll=t[nx+1:nx+nl,nx+1:nx+nl];
phi=zeros(nl,nz);

	
for i (nl,1,-1);
    Temp2=zeros(1,nz);
    for j (i+1,nl,1);
    	Temp2[1,1:nz]=Temp2[1,1:nz]+(Tll[i,j]*Phi[j,1:nz]-Sll[i,j]*(Phi[j,1:nz]*Rho));
	endfor;
	Temp1=Sll[i,i]*Rho-Tll[i,i]*Eye(nz);
	Temp1=inv(Temp1);
	Phi[i,1:nz]=(Dl[i,1:nz]+Temp2[1,1:nz])*Temp1;
        
endfor;
Vxxi=V[1:nx,1:nx];
Vxxi=inv(Vxxi);
Temp1=V[nx+1:nxl,nx+1:nxl];
Temp1=Temp1-(V[nx+1:nxl,1:nx]*(Vxxi*V[1:nx,nx+1:nxl]));
Llz=Temp1*Phi; @ the shock part @
Llx=V[nx+1:nxl,1:nx]*Vxxi; @ the x-part @
Sxxi=inv(S[1:nx,1:nx]);

/* The policy function for xt+1 */
Lxx=(V[1:nx,1:nx]*Sxxi)*(T[1:nx,1:nx]*Vxxi);
Lxz=V[1:nx,nx+1:nxl]*(Phi*Rho); @ Vxl*Phi*Rho @
Temp1=Dx;
Temp1=Temp1+(T[1:nx,nx+1:nxl]*Phi); 
Temp1=Temp1-(S[1:nx,nx+1:nxl]*(Phi*Rho));
Temp1=Temp1-((T[1:nx,1:nx]*Vxxi)*(V[1:nx,nx+1:nxl]*Phi));
Lxz=Lxz+((V[1:nx,1:nx]*Sxxi)*Temp1);

/* Test for complex coefficients */
test=maxc(maxc(imag(Llx)));

if test > toli;
   ?"Warning:There are complex coefficients in Llx";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Llx=real(Llx);

test=maxc(maxc(imag(Llz)));

if test > toli;
   ?"Warning:There are complex coefficients in Llz";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Llz=real(Llz);


test=maxc(maxc(imag(Lxx)));

if test > toli;
   ?"Warning:There are complex coefficients in Lxx";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Lxx=real(Lxx);

test=maxc(maxc(imag(Lxz)));

if test > toli;
   ?"Warning:There are complex coefficients in Lxz";
   ?"itest=" ftos(test,"*.*lf",15,10);
   ?"Press any key...";wait;
endif;
Lxz=real(Lxz);

retp(Lxx,Lxz,Llx,Llz);
endp;


/* MySchur: obtain the Schur factorization of w via a call
** to the dll-function schur_f in the dynamic link library tools.dll
*/

proc(2)=MySchur(w);

    local n, sr, si, tr, ti, t, v;
    n=rows(w);
    sr=zeros(n,n);
    si=zeros(n,n);
    tr=zeros(n,n);
    ti=zeros(n,n);
    dlibrary -a Tools;
    dllcall -r schur_f(n,w,sr,si,tr,ti);
    t=complex(sr,si);
    v=complex(tr,ti);
    dlibrary -d;
retp(t,v);

endp;

/* MyGSchur: obtain the generalized Schur factorization of (a,b) via a call
** to the dll-function gschur_f in the dynamic link library tools.dll
*/

proc(4)=MyGSchur(a,b);

    local n, ur,ui,vr,vi,sr, si, tr, ti, u,v,s,t;
    n=rows(a);
    ur=zeros(n,n);
    ui=zeros(n,n);
    vr=zeros(n,n);
    vi=zeros(n,n);
    sr=zeros(n,n);
    si=zeros(n,n);
    tr=zeros(n,n);
    ti=zeros(n,n);
    dlibrary -a Tools; 
    dllcall  gschur_f(n,a,b,ur,ui,vr,vi,sr,si,tr,ti);
    u=complex(ur,ui);
    v=complex(vr,vi);
    s=complex(sr,si);
    t=complex(tr,ti);
    dlibrary -d;
retp(u,v,s,t);

endp;

/* SchurS */
proc(2)=SchurS(W);

   local T, S, i1, i2, i3, n, S1;

  if rows(W) == 1;
    S=W;
    T=1;
  else;

/* Step 1: Get the real Schur factorization */
{S,T}=schur(W);

/* Step 2: Get the possibly complex Schur factorization */

{S,T}=schtoc(S,T);

/* Step 3:  Sort diagonal elements of T in ascending order  */
n=rows(S); @?"Start Sorting";wait;@
i1=1;
@output off;@
do while i1 < n;
locate 5,5; @?"i1= " ftos(i1,"*.*lf",6,0);@
i2=minindc(abs(diag(S[i1:n,i1:n])));
   i3=i2+i1-1;
locate 5,25;@?"imin= " ftos(i3,"*.*lf",6,0);   @
   do until i3<(i1+1);
      {S,T}=InterChange(S,T,i3-1);
      i3 = i3-1;
   endo;
 i1=i1+1;
endo;
@output on;@
endif;

/* Normalize T */
T=L2Norm(T);
retp(S,T);
endp;

/* InterChange:
**
** Usage: {S,T}=Interchange(S,T,k)
**
** Purpose: Interchange two neighboring elements on the main diagonal of
**          the Schur matrix S and modify the transformation matrix T accordingly.
** 
** Input:   S Schur matrix
**          T tranformation matrix
**          k index of the element being interchanged with element k+1
**
** Remarks: The computations are based on algorithm 7.6.1 in the book of Golub and van Loan (1996), p. 365f.
**          See also Appendix 1 of Heer and Maussner.
*/

proc(2)=InterChange(S,T,k);

   local n, c, s1, gmat;

   n=cols(S);

   { c,s1 } = givens(S[k,k+1],S[k+1,k+1]-S[k,k]);

   gmat=(c~s1)|(-s1'~c'); 

   S[k:k+1,k:n]  =gmat'S[k:k+1,k:n];
   S[1:k+1,k:k+1]=S[1:k+1,k:k+1]*gmat;   
   T[1:n,k:k+1]=T[1:n,k:k+1]*gmat;

retp(S,T);
endp;

/* Givens:
**
** Usage {c,s}=Givens(a,b);
**
** Purpose: Get the coefficients of a Givens rotation that interchanges
**          the element a[i,i] with the element a[i+1,i+1] of a (complex) matrix A.
**
**
** Input: a=a[i,i+1]
**        b=a[i+1,i+1]-a[i,i]
**
** Output: c and s as defined by algorithm 5.1.3 in the book by Golub and van Loan (1996), p. 216.
*/


proc(2) = givens(a,b); 

   local tau, s, c;

if b == 0;

   c=1; s=0;

elseif abs(b) > abs(a);
   tau=-a/b; s=1/(sqrt(1+tau*tau')); c=s'*tau;

else;
   tau=-b/a; c=1/(sqrt(1+tau*tau')); s=c*tau;s=s';

endif;

retp(c,s);
endp;


/* L2Norm:
**
** usage: y=L2Norm(x)
**
** purpose: normalizes the vector x so that it has unit length
**
** input:  x, column vector
** output: y, column vector with unit length
*/

proc(1)=L2Norm(x);

  local i, n, norm;

  n=cols(x);
  i=1;
  do until i>n;
     norm=sqrt(sumc(abs(x[.,i])^2));
     x[.,i]=x[.,i]/norm;
     i=i+1;
  endo;

retp(x);

endp;



