function [Fi,Fc,Gi,Gc] = carlbidirsp(F1i,F1c,G1i,G1c,G11,P)
%CARLBIDIRSP computes the direct Carleman Bilinearization of a
%linear-analytic system, using sparsity
%
% [Fi,Fc,Gi,Gc] = CARLBIDIRSP(F1i,F1c,G1i,G1c,G11,P) computes the Pth order
% direct Carleman Bilinearization v' = Fv + G [1; v]u of the
% linear-analytic system x'= f(x) + g(x)u from its Taylor expansion
% coeficients F1 and G1. The procedure takes advantage of the
% sparsity of the matrices involved. Inputs fi, fc, gi and gc and outputs
% are a sparsity oriented representation.

%   CBVK Toolbox
%   JHMG, PMSB 2016

warning('cbvk:carlbidirsp','Consider the use of carlbilinsp, which is possibly faster');

m = length(F1i);

% Conversion from cell to matrix format
F1im=zeros(m,1);
F1cm=F1im;
G1im=F1im;
G1cm=F1im;

for i=1:m
    F1im(i,1:length(F1i{i}))=F1i{i};
    F1cm(i,1:length(F1i{i}))=F1c{i};
    
    G1im(i,1:length(G1i{i}))=G1i{i}+1;
    G1cm(i,1:length(G1i{i}))=G1c{i};
    
    if G11(i) ~=0
        G1im(i,1:length(G1im(i))+1)=[1 G1im(i)];
        G1cm(i,1:length(G1cm(i))+1)=[G11(i) G1cm(i)];
    end
     
end

% Bilinearization

allterms = m.^(1:P);
MP = sum(allterms);
numtermsF = cumsum([1 allterms]);
numtermsG = cumsum([1 1 allterms]);


[Fi,Fc,Gi,Gc]=cbdsp(F1im',F1cm',G1im',G1cm',m,P,MP,allterms,numtermsF,numtermsG);


end


