function F = interpSpline(X, localX, V, localV, dimInterp)
%INTERPSPLINE calculates the local parts of the piecewise polynomial
%representation of a cubic spline fit on a distributed grid.
%   Inputs: X         distributed grid.
%           localX    local part of X padded by 2 elements on either side.
%           V         distributed value array.
%           localV    local part of V padded by 2 elements on either side.
%           dimInterp the dimension of X and V on which interpolation is
%                     done (also the distribution dimension).
%   Output: F piecewise polynomial structure for the local part of the fit.

%   Copyright 2021-2022 The MathWorks, Inc.

% Do processing that would have been done in chckxy.
% If V is ND, reshape to a matrix by combining all dimensins but the
% distribution dimension, which must be last.
sizev = size(V);
dims = 1:numel(sizev);
dims(dimInterp) = [];
dims = [dimInterp dims];  % Put interpolation dimension first
localV = permute(localV, dims);
lvn = size(localV, 1);
sizev(dimInterp) = [];
vd = prod(sizev);

if length(sizev) > 1  % V is not a vector
    localV = reshape(localV, lvn, vd);
end

localX = localX(:);
dLocalX = diff(localX);
divdif = diff(localV, 1, 1)./dLocalX;

% First find out what grid points each worker has.
codX = getCodistributor(X);
lnX = length(X);
padSz = 2;
partitionX = codX.Partition;
ixLast = cumsum(partitionX);  % Last element of X on each worker.
diagsOffset = [0 ixLast(1:end-1)];  % displacement of distributed tridiagonals
leftPad = [0 min(ixLast, padSz)];
leftPad(end) = [];
rightPad = min(padSz, lnX - ixLast);
dLocalXstart = leftPad + 1;
dLocalXend = size(localV, 1) - min(padSz, rightPad);
firstUsedWorker = find(ixLast, 1);                % First worker that has data.
lastUsedWorker = find(ixLast == ixLast(end), 1);  % Last worker with data.

if lnX == 2
    F = mkpp(localX, [divdif; localV(1, :)].', sizev);
elseif lnX == 3
    if length(localX) == 2
        F = mkpp(localX, [divdif; localV(1, :)].', sizev);
    else
        localV(3, :) = diff(divdif)/(localX(3) - localX(1));
        localV(2, :) = divdif(1, :) - localV(3, :)*dLocalX(1);
        F = mkpp(localX([1, 3]), localV([3 2 1], :).', sizev);
    end
else  % lnX > 3
    % Set up the sparse, tridiagonal, linear system s = C\B for the slopes.
    % Build B, the matrix of rhs's.
    localB = zeros(size(localV, 1), vd);
    if partitionX(spmdIndex) > 0 % Computation unnecessary on workers without X
        localB(2:end-1, :) = 3*(dLocalX(2:end).*divdif(1:end-1, :)...
            + dLocalX(1:end-1).*divdif(2:end, :));
        if spmdIndex == firstUsedWorker  % worker calculates first element of B
            x31 = localX(3) - localX(1);
            localB(1, :) = ((dLocalX(1)+2*x31)*dLocalX(2)*divdif(1, :)...
                + dLocalX(1)^2*divdif(2, :))/x31;
        end
        if spmdIndex == lastUsedWorker % worker calculates last B
            xn = localX(end) - localX(end-2);
            localB(end, :) = (dLocalX(end)^2*divdif(end-1, :)...
                + (2*xn + dLocalX(end))*dLocalX(end-1)*divdif(end, :))/xn;
        end
    end
    Bcodist = codistributor1d(1, partitionX, [lnX, vd]);
    B = codistributed.pDoBuildFromLocalPart(localB(dLocalXstart(spmdIndex):...
        dLocalXend(spmdIndex), :), Bcodist);
    
    % Build C, the tridiagonal matrix
    if partitionX(spmdIndex) > 0 % Computation unnecessary on workers without X
        tridiags = [-1 0 1] + diagsOffset(spmdIndex);
        localRPad = rightPad(spmdIndex);
        localLPad = leftPad(spmdIndex);
        if spmdIndex == firstUsedWorker  % First block of C.
            localDiag = [0 dLocalX(2) (localX(3) - localX(1))];
            if length(localX) > 3
                localDiag = [localDiag; [dLocalX(2:end+1-localRPad) ...
                    2*(dLocalX(2:end+1-localRPad) + dLocalX(1:end-localRPad))...
                    dLocalX(1:end-localRPad)]]; 
            end
        elseif spmdIndex == lastUsedWorker  % Last block of C.
            localDiag = [localX(end) - localX(end-2) dLocalX(end-1) 0];
            if length(localX) > 3
                localDiag = [[dLocalX(localLPad+1:end) ...
                    2*(dLocalX(localLPad+1:end) + dLocalX(localLPad:end-1))...
                    dLocalX(localLPad:end-1)]; localDiag];
            end
        else   % Rest of C
            localDiag = [dLocalX(localLPad+1:end+1-localRPad)...
                2*(dLocalX(localLPad+1:end+1-localRPad) + dLocalX(localLPad:end-localRPad))...
                dLocalX(localLPad:end-localRPad)];
        end
        localC = spdiags(localDiag, tridiags, partitionX(spmdIndex), lnX);
    else
        localC = zeros(0, lnX);
    end
    
    Ccodist = codistributor1d(1, partitionX, [lnX lnX]);
    C = codistributed.pDoBuildFromLocalPart(localC, Ccodist);
    
    % Solve for slopes
    S = scalaTridiagSolve(distributedutil.TransposeFlags.NO_TRANSPOSE, C, B);
    
    % Construct local piecewise cubic Hermite interpolant.
    % Get padded localS to match localX and localV.
    [~, localS] = padLocalGrid({X}, S, padSz);
    % pwch expects matrices to be oriented as rows.
    F = pwch(localX.', localV.', localS.', dLocalX.', divdif.');
end
F.dim = sizev;
F.GridVectors = {localX};  % Needed at evaluation time.
F.Values = localV;
F.dims = dims;
end
