function kernels = vkrnsvd(F,G,b,c,N,T,m,P,msk)
%VKRNSVD computes triangular Volterra kernels from the Carleman
%Bilinearization of a system
%
% KERNELS = VKRNSVD(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. The kernels are returned in vector
% form.
%
% KERNELS = VKRN(F,G,B,C,N,T,M,P,msk) uses a cell of P pre-calculated masks
% (see util/mask_P*_N*.m)

% See also CARLBILIN and CARLBL
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

warning('cbvk:vkrnsvd','Consider the use of vkrnsvdx, which is usually faster');


% Vector initialization
h1 = zeros(1,N);
if P > 1
    h2 = zeros(1,nchoosek(N+1,2));
    if P > 2
        h3 = zeros(1,nchoosek(N+2,3));
        if P > 3
            h4 = zeros(1,nchoosek(N+3,4));
            if P > 4
                h5 = zeros(1,nchoosek(N+4,5));
            end
        end
    end
end    


% 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=cell(1,N);
eV=cell(1,N);
eUV=cell(1,N);

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

% Kernel indices
i2=1;i3=1;i4=1;i5=1;
% Kernel masks
m2=1;m3=1;m4=1;m5=1;
if nargin <9
    % Main loop with mask calculation
    if P>1
        for k1=1:N
            % Repetition counter for k indices
            rk = 1;
            for k2=k1:N
                if(k2==k1)
                    rk = rk + 1;
                else
                    rk = 1;
                    % Mask for h2
                    m2 = 1/rk;
                end
                h2(i2) = m2 * eU{k1}*eV{k2-k1+1};
                i2=i2+1;
                if P>2
                    v2=eU{k1}*eUV{k2-k1+1};
                    for k3=k2:N
                        if(k3==k2)
                            rk = rk + 1;
                        else
                            rk = 1;
                            % Mask for h3
                            m3 = m2/rk;
                        end
                        h3(i3) = m3 * v2*eV{k3-k2+1};
                        i3=i3+1;
                        if P>3
                            v3=v2*eUV{k3-k2+1};
                            for k4=k3:N
                                if(k4==k3)
                                    rk = rk + 1;
                                else
                                    rk = 1;
                                    m4 = m3/rk;
                                end
                                h4(i4) = m4 * v3*eV{k4-k3+1};
                                i4=i4+1;
                                if P > 4
                                    v4=v3*eUV{k4-k3+1};
                                    for k5=k4:N
                                        if(k5==k4)
                                            rk = rk + 1;
                                        else
                                            rk = 1;
                                            m5 = m4/rk;
                                        end
                                        h5(i5) = m5 * v4*eV{k5-k4+1};
                                        i5=i5+1;
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    if P > 1
        kernels{2} = h2;
        if P > 2
            kernels{3} = h3;
            if P > 3
                kernels{4} = h4;
                if P > 4
                    kernels{5} = h5;
                end
            end
        end
    end

else
    % Main loop without mask calculation
    if P>1
        for k1=1:N
            for k2=k1:N
                h2(i2) = eU{k1}*eV{k2-k1+1};
                i2=i2+1;
                if P>2
                    v2=eU{k1}*eUV{k2-k1+1};
                    for k3=k2:N
                        h3(i3) = v2*eV{k3-k2+1};
                        i3=i3+1;
                        if P>3
                            v3=v2*eUV{k3-k2+1};
                            for k4=k3:N
                                h4(i4) = v3*eV{k4-k3+1};
                                i4=i4+1;
                                if P > 4
                                    v4=v3*eUV{k4-k3+1};
                                    for k5=k4:N
                                        h5(i5) = v4*eV{k5-k4+1};
                                        i5=i5+1;
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    if P > 1
        kernels{2} = h2.*msk{2};
        if P > 2
            kernels{3} = h3.*msk{3};
            if P > 3
                kernels{4} = h4.*msk{4};
                if P > 4
                    kernels{5} = h5.*msk{5};
                end
            end
        end
    end

end


