function y = sumMatrixIncludeNaN(b,identity,x,dim,vlen,vstart)
%MATLAB Code Generation Library Function

%   Sum of matrix (x is 2D) columns in column-major mode or rows in
%   row-major mode with includenan semantics. The x array can be ND when
%   dim = 1. It will work, but note that higher dimensions will be
%   flattened during processing, which is inconsistent with the preserve
%   array dimensions option. The function will still work when the majority
%   doesn't match the direction of summation, but the code will be slower
%   in that case because of the striding.
%
%   NOTE: x must not be empty.

%   Copyright 2002-2020 The MathWorks, Inc.
%#codegen

coder.inline('always');
coder.internal.prefer_const(b,identity,dim,vlen,vstart);
coder.internal.allowEnumInputs;
coder.internal.allowHalfInputs;
if dim == 1
    y = sumMatrixColumns(b,identity,x,vlen,vstart);
else
    y = sumMatrixRows(b,identity,x,vlen,vstart);
end

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

function y = sumMatrixColumns(b,identity,x,vlen,vstart)
% x can be ND. Optimized for column-major mode.
coder.internal.prefer_const(b,identity,vlen,vstart);

sz = size(x);
sz(1) = 1;
y = coder.nullcopy(eml_expand(identity,sz));
ncols = coder.internal.prodsize(x,'above',1);
b4 = coder.const(b*4);
AVOID_SUMCOLUMN4 = tooSmallToBlock(b4,x,1,vlen,vstart);
if AVOID_SUMCOLUMN4 || vlen < b4
    % This loop could be parallelized if there are multiple columns.
    for col = 1:ncols
        y(1,col) = sumColumnB(b,identity,x,col,vlen,vstart);
    end
else
    nfb = coder.internal.indexDivide(vlen,b4); % number of full blocks
    inb = nfb*b4; % number of elements in full blocks
    lidx = vstart + inb; % first index of leftovers
    nleft = vlen - inb; % number of leftovers
    % This loop could be parallelized if there are multiple columns.
    for col = 1:ncols
        s = sumColumnB4(b,identity,x,col,vstart);
        % This loop could be parallelized if there is only one column.
        for ib = 2:nfb
            s = s + sumColumnB4(b,identity,x,col,vstart + (ib - 1)*b4);
        end
        if nleft > 0
            s = s + sumColumnB(b,identity,x,col,nleft,lidx);
        end
        y(1,col) = s;
    end
end

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

function y = sumMatrixRows(b,identity,x,vlen,vstart)
% x must be 2D. Optimized for row-major mode.
coder.inline('always');
coder.internal.prefer_const(b,identity,vlen,vstart);

sz = size(x);
sz(2) = 1;
y = coder.nullcopy(eml_expand(identity,sz));
nrows = coder.internal.indexInt(size(x,1));
b4 = coder.const(b*4);
AVOID_SUMROW4 = tooSmallToBlock(b4,x,2,vlen,vstart);
if AVOID_SUMROW4 || vlen < b4
    % This loop could be parallelized if there are multiple rows.
    for row = 1:nrows
        y(row,1) = sumRowB(b,identity,x,row,vlen,vstart);
    end
else
    nfb = coder.internal.indexDivide(vlen,b4); % number of full blocks
    inb = nfb*b4; % number of elements in full blocks
    lidx = vstart + inb; % first index of leftovers
    nleft = vlen - inb; % number of leftovers
    % This loop could be parallelized if there are multiple rows.
    for row = 1:nrows
        s = sumRowB4(b,identity,x,row,vstart);
        % This loop could be parallelized if there is only one row.
        for ib = 2:nfb
            s = s + sumRowB4(b,identity,x,row,vstart + (ib - 1)*b4);
        end
        if nleft > 0
            s = s + sumRowB(b,identity,x,row,nleft,lidx);
        end
        y(row,1) = s;
    end
end

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

function y = sumColumnB4(b,identity,x,col,vstart)
% Column-major blocked and partially-unrolled column summation.
% y = sum(v,1) for vlen = 4*b, b >= 1.
coder.inline('never'); % to avoid loop fusion
coder.internal.prefer_const(b,identity);

b2 = 2*b;
b3 = 3*b;
if coder.internal.avoidArrayFlattening
    i1 = vstart;
    i2 = vstart + b;
    i3 = vstart + b2;
    i4 = vstart + b3;
    psum1 = cast(x(i1,col),'like',identity);
    psum2 = cast(x(i2,col),'like',identity);
    psum3 = cast(x(i3,col),'like',identity);
    psum4 = cast(x(i4,col),'like',identity);
    for k = 1:b-1
        psum1 = psum1 + cast(x(i1 + k,col),'like',identity);
        psum2 = psum2 + cast(x(i2 + k,col),'like',identity);
        psum3 = psum3 + cast(x(i3 + k,col),'like',identity);
        psum4 = psum4 + cast(x(i4 + k,col),'like',identity);
    end
else
    i1 = vstart + (col - 1)*coder.internal.indexInt(size(x,1));
    i2 = i1 + b;
    i3 = i1 + b2;
    i4 = i1 + b3;
    psum1 = cast(x(i1),'like',identity);
    psum2 = cast(x(i2),'like',identity);
    psum3 = cast(x(i3),'like',identity);
    psum4 = cast(x(i4),'like',identity);
    for k = 1:b-1
        psum1 = psum1 + cast(x(i1 + k),'like',identity);
        psum2 = psum2 + cast(x(i2 + k),'like',identity);
        psum3 = psum3 + cast(x(i3 + k),'like',identity);
        psum4 = psum4 + cast(x(i4 + k),'like',identity);
    end
end
y = (psum1 + psum2) + (psum3 + psum4);

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

function y = sumRowB4(b,identity,x,row,vstart)
% Row-major blocked and partially-unrolled row summation.
% y = sum(v,2) for vlen = 4*b, b >= 1.
coder.inline('never'); % to avoid loop fusion
coder.internal.prefer_const(b,identity);

b2 = 2*b;
b3 = 3*b;
i1 = vstart;
i2 = i1 + b;
i3 = i1 + b2;
i4 = i1 + b3;
psum1 = cast(x(row,i1),'like',identity);
psum2 = cast(x(row,i2),'like',identity);
psum3 = cast(x(row,i3),'like',identity);
psum4 = cast(x(row,i4),'like',identity);
for k = 1:b-1
    psum1 = psum1 + cast(x(row,i1 + k),'like',identity);
    psum2 = psum2 + cast(x(row,i2 + k),'like',identity);
    psum3 = psum3 + cast(x(row,i3 + k),'like',identity);
    psum4 = psum4 + cast(x(row,i4 + k),'like',identity);
end
y = (psum1 + psum2) + (psum3 + psum4);

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

function y = sumColumnB(b,identity,x,col,vlen,vstart)
% Column-major blocked column summation sum(x(vstart:vstart+vlen-1,col)).
coder.inline('never');
coder.internal.prefer_const(b,identity);

AVOID_BLOCKING = tooSmallToBlock(b,x,1,vlen,vstart);
if AVOID_BLOCKING || vlen <= b
    y = sumColumn(identity,x,col,vlen,vstart);
else
    nfb = coder.internal.indexDivide(vlen,b); % number of full blocks
    inb = nfb*b; % number of elements in full blocks
    y = sumColumn(identity,x,col,b,vstart);
    for k = 2:nfb
        y = y + sumColumn(identity,x,col,b,vstart + (k - 1)*b);
    end
    if vlen > inb
        nleft = vlen - inb; % number of leftovers
        lidx = vstart + inb; % first index of leftovers
        y = y + sumColumn(identity,x,col,nleft,lidx);
    end
end

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

function y = sumColumn(identity,x,col,vlen,vstart)
% Column-major unblocked column summation sum(x(vstart:vstart+vlen-1,col)).
coder.inline('always');
coder.internal.prefer_const(identity);

if coder.internal.avoidArrayFlattening
    y = cast(x(vstart,col),'like',identity);
    i1 = vstart + 1;
    i2 = vstart + vlen - 1;
    for k = i1:i2
        y = y + cast(x(k,col),'like',identity);
    end
else
    i0 = vstart + (col - 1)*coder.internal.indexInt(size(x,1));
    y = cast(x(i0),'like',identity);
    for k = 1:vlen-1
        y = y + cast(x(i0 + k),'like',identity);
    end
end

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

function y = sumRowB(b,identity,x,col,vlen,vstart)
% Row-major blocked row summation sum(x(row,vstart:vstart+vlen-1)).
coder.inline('never');
coder.internal.prefer_const(b,identity);

AVOID_BLOCKING = tooSmallToBlock(b,x,2,vlen,vstart);
if AVOID_BLOCKING || vlen <= b
    y = sumRow(identity,x,col,vlen,vstart);
else
    nfb = coder.internal.indexDivide(vlen,b); % number of full blocks
    inb = nfb*b; % number of elements in full blocks
    y = sumRow(identity,x,col,b,vstart);
    for k = 2:nfb
        y = y + sumRow(identity,x,col,b,vstart + (k - 1)*b);
    end
    if vlen > inb
        nleft = vlen - inb; % number of leftovers
        lidx = vstart + inb; % first index of leftovers
        y = y + sumRow(identity,x,col,nleft,lidx);
    end
end

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

function y = sumRow(identity,x,row,vlen,vstart)
% Row-major unblocked row summation sum(x(row,vstart:vstart+vlen-1)).
coder.inline('always');
coder.internal.prefer_const(identity);

y = cast(x(row,vstart),'like',identity);
i1 = vstart + 1;
i2 = vstart + vlen - 1;
for k = i1:i2
    y = y + cast(x(row,k),'like',identity);
end

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

function p = tooSmallToBlock(b,x,dim,vlen,vstart)
% We assume vstart+vlen-1 <= size(x,dim), but there are cases where the
% compiler may be able to generate dead code that indexes out of range.
% This constant predicate can be used to avoid generating blocking code,
% with its fixed-length loops, when the code will never be called. What it
% does, in other words, is connect all the available information to the
% question of whether the sums will be blocked, whereas the code generator
% does not always put all the pieces together.
coder.inline('always');
coder.internal.prefer_const(b,dim,vlen,vstart);

SMALL_VLEN = coder.internal.isConstTrue(vlen <= b);
SMALL_X = coder.internal.isConstTrue(size(x,dim) <= b);
lastIdx = vstart + (b - 1);
SHORT_VECTOR = coder.internal.isConstTrue(lastIdx > size(x,dim));
p = coder.const(SMALL_VLEN || SMALL_X || SHORT_VECTOR);

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