function y = ltitr(a,b,u,x0)
%MATLAB Code Generation Library Function

%   Copyright 1984-2017 The MathWorks, Inc.
%#codegen

coder.columnMajor;
eml_invariant(nargin >= 3, ...
    'MATLAB:minrhs');
check_class(a);
check_class(b);
check_class(u);
eml_invariant(ismatrix(a) && ismatrix(b) && ismatrix(u), ...
    'Coder:toolbox:ltitr_inputsMustBe2D');
am = cast(size(a,1),coder.internal.indexIntClass);
n = cast(size(a,2),coder.internal.indexIntClass);
eml_invariant(am == n, ...
    'Coder:MATLAB:square');
bm = cast(size(b,1),coder.internal.indexIntClass);
bn = cast(size(b,2),coder.internal.indexIntClass);
um = cast(size(u,1),coder.internal.indexIntClass);
un = cast(size(u,2),coder.internal.indexIntClass);
eml_invariant(n == bm && bn == un, ...
    'Coder:MATLAB:dimagree');
abu0 = coder.internal.scalarEg(a,b,u);
ONE = ones(coder.internal.indexIntClass);
if nargin < 4
    y = eml_expand(abu0,[um,n]);
    if isempty(u)
        return
    end
else
    check_class(x0);
    eml_invariant(isvector(x0) && numel(x0) == n, ...
        'Coder:MATLAB:ltitr_invalidInitialCondition');
    y = eml_expand(coder.internal.scalarEg(abu0,x0),[um,n]);
    for k = ONE:n
        y(1,k) = x0(k);
    end
end
CZERO = coder.internal.scalarEg(y);
CONE = CZERO + 1;
for i = ONE:coder.internal.indexMinus(um,ONE)
    ip1 = coder.internal.indexPlus(i,ONE);
    y = coder.internal.blas.xgemv('N',am, n,1,a,ONE,am,y,i,um,CZERO,y,ip1,um);
    y = coder.internal.blas.xgemv('N',bm,bn,1,b,ONE,bm,u,i,um,CONE,y,ip1,um);
end

%--------------------------------------------------------------------------

function check_class(x)
% Note that the algorithm above would work with single as well, but MATLAB 
% requires double, so we require double.
eml_invariant(isa(x,'double'), ...
    'Coder:toolbox:unsupportedClass','ltitr',class(x));

%--------------------------------------------------------------------------
