% file = irks.m
%
% This function computes explicit and implicit
% irks methods.
%
% lambda=0 explicit, with error constant cp1
% lambda>0 Lstable, cp1 is ignored
% 
% function glm=irks(lambda,cp1,c,beta,T,dbeta,dT)
%
% Date: 17 March 2006
%
% -----------------------------------------------

function glm=irks(lambda,cp1,c,beta,T,dbeta,dT)

s=length(c); p=s-1;

if lambda>1e-5 % compute the error constant for L-stabiliy
    pp=factorial(p)./(factorial(0:p).^2.*factorial(p:-1:0));
    cp1=polyval([pp,1/factorial(s)],-lambda);
end

colx=[0;0;-beta'];
K=diag(ones(s-1,1),1);
%E=polyvalm(1./gamma(s:-1:1),K);
E=toeplitz([1;zeros(s-1,1)],1./gamma(1:s));

Psi=zeros(s,s); Psi(s,s)=1;

for i=s-1:-1:1
    Psi(i,s)=lambda*Psi(i+1,s)-colx(i+1);
end

for i=s-1:-1:1
    jj=s-1:-1:1;
    Psi(i,jj)=lambda*Psi(i+1,jj)+Psi(i+1,jj+1);
end
Psii=inv(Psi);

X=lambda*eye(s)+Psi*K'*Psii;
C=repmat(c',1,s).^repmat(0:(s-1),s,1)/diag(gamma(1:s));

es=[zeros(s-1,1);1]; Xes=X(:,s); Ees=E(:,s);



% Generate the linear systems for Btilde
rhs=-[ Psii(2:end,:)*(E*Xes-X*Ees) 
          zeroV(E,T)
         1/factorial(s)-cp1
          zeros(s-1,1)];
Bt=zeros(s,s); G=[];

for i=1:s
    for j=1:i
        Bt(i,j)=1; 
        eq=[
          Bt(2:end,:)*C*(es-K*Xes)+Psii(2:end,:)*(E*Xes-X*Ees) % IRKS
          zeroV(E-Psi*Bt*C*K,T)                      % ZeroStab
          1/factorial(s)-cp1-Bt(1,1)*(C(1,:)-lambda*C(1,:)*K)*Psi(:,s) % Error Const
          (Psi(2,:)*Bt(:,1:s-1))'                    % FSAL condition
        ];              
        G=[G,eq+rhs]; Bt(i,j)=0;
    end
end


bt=G\rhs;
for i=1:s
    l=1+i*(i-1)/2;
    jj=l:l+i-1;
    Bt(i,1:i)=bt(jj)'; 
end
glm.btilde=Bt;

glm.B=Psi*Bt;
glm.A=(Bt\(Psi\X))*glm.B;
glm.V=E-Psi*Bt*C*K;
glm.U=C-glm.A*C*K;


if nargin<=5, return, end % ====================================

dcolx=[0;0;-dbeta'];
dPsi=zeros(s,s);

for i=s-1:-1:1
    dPsi(i,s)=lambda*dPsi(i+1,s)-dcolx(i+1);
end

for i=s-1:-1:1
    jj=s-1:-1:1;
    dPsi(i,jj)=lambda*dPsi(i+1,jj)+dPsi(i+1,jj+1);
end

dPsii=-inv(Psi)*dPsi*inv(Psi);

dX=dPsi*K'*inv(Psi)+Psi*K'*dPsii;

dXes=dX(:,s);

row2end=@(r)reshape(r(2:end),s-1,1);


dEQ1=@(Bt)row2end(...
   Bt*C*(-K*dXes)+dPsii*(E*Xes-X*Ees)+Psi\(E*dXes-dX*Ees));  % IRKS condition
dEQ2=@(Bt)zerodV(E-Psi*Bt*C*K,-dPsi*Bt*C*K,T,dT);     % Zero stability
dEQ3=@(Bt)Bt(1,1)*(C(1,:)-lambda*C(1,:)*K)*Psi(:,s); % Error Const
dEQ4=@(Bt)(dPsi(2,:)*Bt(:,1:s-1))';                
dEQ=@(Bt)[dEQ1(Bt);dEQ2(Bt);dEQ3(Bt);dEQ4(Bt)];       

dBt=zeros(s,s); drhs=-dEQ(dBt); dG=[];
for i=1:s
    for j=1:i
        dBt(i,j)=1; dG=[dG,dEQ(dBt)+drhs]; dBt(i,j)=0;
    end
end

dGi=-inv(G)*dG*inv(G);

bt=G\drhs+dGi*rhs;
for i=1:s
    l=1+i*(i-1)/2;
    jj=l:l+i-1;
    dBt(i,1:i)=bt(jj); 
end

dbtilde=dBt;
dB=dPsi*Bt+Psi*dBt;
dBi=-inv(glm.B)*dB*inv(glm.B);

dA=dBi*X*glm.B+inv(glm.B)*dX*glm.B+inv(glm.B)*X*dB;
dV=-dPsi*Bt*C*K-Psi*dBt*C*K;
dU=-dA*C*K;

glm.btilde=dbtilde;
glm.A=dA; glm.B=dB; glm.U=dU; glm.V=dV;

%-------- some helpers .......................

function elem=zeroV(V,T)
elem=[trilpart(T\V(3:end,3:end)*T)];

function elem=zerodV(V,dV,T,dT)
dTi=-inv(T)*dT*inv(T);
elem=[trilpart(dTi*V(3:end,3:end)*T + ...
               T\dV(3:end,3:end)*T   + ...
               T\V(3:end,3:end)*dT)];

function elem=trilpart(T)
s=length(T);
elem=zeros(s*(s+1)/2,1);
for i=1:s
      jj=1:i;
      kk=i*(i-1)/2+jj;
      elem(kk)=T(i,jj);
end