function kernels = vkrnsvdx(F,G,b,c,N,T,m,P)
%VKRNSVDX computes triangular Volterra kernels from the Carleman
%Bilinearization of a system
%
% KERNELS = VKRNSVDX(F,G,B,C,N,T,M,P) computes T-spaced samples up to time
% N*T of the first P<=5 Volterra triangular kernels of the bilinear system
% [F,G,b,c] obtained by Pth order Carleman Bilinearization of a
% linear-analytic sistem whose state vector has dimension M. SVD is used to
% exploit the rank defficiency of G. A mex function is used for greater
% speed. The kernels are returned in vector form
%
% See also CARLBILIN and CARLBL
%
%   CBVK Toolbox
%   JHMG, PMSB 2016



% Rank of G
mp=calcnumterms(P,m);
rG=sum(mp(2:end));

% SVD and order reduction
[u,s,v]=svd(G);
S=sqrt(diag(s));
S=S(1:rG);
u=u(:,1:rG)*diag(S);
v=v(:,1:rG)*diag(S);

eT = expm(F*T);


%---------------------------------------------------
% This is for greater precision as P grows
eF1=expm(F(1:m,1:m)*T);
if P>2
    M2=m+nchoosek(m+2-1,2);
    eF2=expm(F(1:M2,1:M2)*T);
    if P>3
        M3=M2+nchoosek(m+3-1,3);
        eF3=expm(F(1:M3,1:M3)*T);
        if P>4
            M4=M3+nchoosek(m+4-1,4);
            eF4=expm(F(1:M4,1:M4)*T);
        end
    end
end
if P==2
    eT(1:m,1:m)=eF1;
elseif P==3
    eF2(1:m,1:m)=eF1;
    eT(1:M2,1:M2)=eF2;
elseif P==4
    eF2(1:m,1:m)=eF1;
    eF3(1:M2,1:M2)=eF2;
    eT(1:M3,1:M3)=eF3;
elseif P==5
    eF2(1:m,1:m)=eF1;
    eF3(1:M2,1:M2)=eF2;
    eF4(1:M3,1:M3)=eF3;
    eT(1:M4,1:M4)=eF4;
end
%------------------------------------

eU=zeros(N*rG,1);
eV=zeros(N*rG,1);
eUV=zeros(N*(rG^2),1);

h1=zeros(1,N);

% Pre-computed matrices
ce = c;									 
be = b;									 
eF = v';								 
  for k=1:N								 
    h1(k)=ce*b;								 
    if P>1								 
      eU((k-1)*rG+1:k*rG) = (ce*u)';
      eV((k-1)*rG+1:k*rG)= v'*be;
      et=eF*u;
      eUV((k-1)*rG*rG+1:k*rG*rG)=et(:);								 
    end									 
    ce = ce*eT;								 
    be = eT*be;								 
    eF = eF*eT;								 
  end 
kernels=cell(1,P);
kernels{1}=h1;

% Size of vectorized kernels
lp=zeros(1,4);
for p=2:5
  lp(p-1)=nchoosek(N+p-1,p);
end

% Call C routine
if (P==2)
  kernels{2}=voltkrn(eU,eV,eUV,rG,N,P,lp);
elseif(P==3)
  [kernels{2},kernels{3}]=voltkrn(eU,eV,eUV,rG,N,P,lp);
elseif(P==4)
  [kernels{2},kernels{3},kernels{4}]=voltkrn(eU,eV,eUV,rG,N,P,lp);
elseif(P==5)
  [kernels{2},kernels{3},kernels{4},kernels{5}]=voltkrn(eU,eV,eUV,rG,N,P,lp);
end

