function Z = pagemtimes(varargin)
%MATLAB Code Generation Library Function

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

narginchk(2, 4);
coder.internal.errorIf(nargin == 3, 'MATLAB:pagemtimes:threeInputNotSupported');

transpX = 'N';
transpY = 'N';
if nargin == 2
    X = varargin{1};
    Y = varargin{2};
else  % 4 args
    coder.internal.assert(coder.internal.isConst(varargin{2}), 'Coder:builtins:ConstantRequired');
    coder.internal.assert(coder.internal.isConst(varargin{4}), 'Coder:builtins:ConstantRequired');
    X = varargin{1};
    transpX = processTransposeValue(varargin{2});
    Y = varargin{3};
    transpY = processTransposeValue(varargin{4});
end
checkArrayType(X);
checkArrayType(Y);
outputSize = processArraySizes(X, Y, transpX, transpY);

pageCount = coder.internal.indexInt(prod(outputSize(3:end)));
scalarPagesX = (size(X,1) == 1 && size(X,2) == 1);
scalarPagesY = (size(Y,1) == 1 && size(Y,2) == 1);

if coder.gpu.internal.isGpuEnabled && ~scalarPagesX && ~scalarPagesY
    % When using gpucoder.stridedMatrixMultiply, the two input collections
    % of pages must have the same dimensions, unless one input is a matrix.
    % The two lines below explicitly expand the input arrays (if necessary)
    % so that they can be passed to gpucoder.stridedMatrixMultiply.
    X_expanded = expandArray(X, outputSize);
    
    if isa(X, 'single') && isa(Y, 'double')
        Y_expanded = single(expandArray(Y, outputSize));
    else
        Y_expanded = expandArray(Y, outputSize);
    end
    
    % The class of the output of gpucoder.stridedMatrixMultiply is the
    % class of the first input. In the case X is double but Y is single,
    % X is converted to single to ensure the output is single.
    if isa(X,'double') && isa(Y,'single')
        X_expanded2 = single(X_expanded);
    else        
        X_expanded2 = X_expanded;
    end
    
    transpValue = [transpX transpY];
    Z = gpucoder.stridedMatrixMultiply(X_expanded2, Y_expanded, 'transpose', transpValue);
else
    Z = coder.nullcopy(zeros(outputSize, 'like', coder.internal.scalarEg(X,Y)));

    sizeX = size(X);
    sizeY = size(Y);
    sameDimsAfterFirstTwo = isequal(sizeX(3:end), sizeY(3:end));
    for ii = 1 : pageCount
        if sameDimsAfterFirstTwo
            X_idx = ii;
            Y_idx = ii;
        else
            X_idx = computeMultiplicationIndex(X, outputSize, ii);
            Y_idx = computeMultiplicationIndex(Y, outputSize, ii);
        end
        
        pageX = X(:,:,X_idx);
        pageY = Y(:,:,Y_idx);
        if scalarPagesX || scalarPagesY
            transformedPageX = transformPage(pageX, transpX);
            transformedPageY = transformPage(pageY, transpY);
            if scalarPagesX
                alpha = transformedPageX(1);
                Z(:,:,ii) = alpha * transformedPageY;
            else
                alpha = transformedPageY(1);
                Z(:,:,ii) = transformedPageX * alpha;
            end
        else
            Z(:,:,ii) = coder.internal.blas.mtimes(pageX, pageY, transpX, transpY);
        end
    end
end



function pageOut = transformPage(pageIn, transp)
% Takes a page sliced from an input array and applies a transposition
% transformation on it. "transp" is assumed to take one of the following
% values:
%   'N' = no transposition (returns pageIn)
%   'T' = nonconjugate transposition
%   'C' = complex conjugate transposition

coder.inline("always");
if transp == 'N'
    pageOut = pageIn;
elseif transp == 'T'
    pageOut = pageIn.';
else
    pageOut = pageIn';
end



function in_idx = computeMultiplicationIndex(X, outputSize, out_idx)
% When multiplying two pages, we need to find their indices in their arrays
% of pages. Some calculation is needed here since an input array of pages
% need not have the same size as the output array of pages due to implicit
% expansion.
%
% - out_idx: The index of the page in the output array to be calculated
% - in_idx: The index of the page in the input array to retrieve

coder.inline("always");
in_idx = coder.internal.indexInt(1);
v = out_idx - 1;
p = coder.internal.indexInt(1);

coder.unroll();
for ii = 3 : numel(outputSize)
    if size(X,ii) ~= 1
        in_idx = in_idx + p * mod(v, outputSize(ii));
        p = p * outputSize(ii);
    end
    v = idivide(v, outputSize(ii));
end



function X_expanded = expandArray(X, outputSize)

coder.internal.prefer_const(outputSize);
coder.inline("always")
sizeX = size(X);
if ismatrix(X) || isequal(sizeX(3:end), outputSize(3:end))
    X_expanded = X;
else
    % Expand X so that the dimensions of X beyond the first two match
    % outputSize
    reps = ones(1, numel(outputSize), 'coder.internal.indexInt');
    coder.unroll();
    for ii = 3 : numel(outputSize)
        if size(X,ii) == 1
            reps(ii) = outputSize(ii);
        end
    end
    X_expanded = repmat(X, reps);
end



function outputSize = processArraySizes(X, Y, transpX, transpY)

coder.internal.prefer_const(transpX);
coder.internal.prefer_const(transpY);
coder.inline("always");

[rowsX, colsX] = getRowsAndCols(X, transpX);
[rowsY, colsY] = getRowsAndCols(Y, transpY);

if rowsX == 1 && colsX == 1
    % Scalar expansion case
    rowsZ = rowsY;
    colsZ = colsY;
elseif rowsY == 1 && colsY == 1
    % Scalar expansion case
    rowsZ = rowsX;
    colsZ = colsX;
else
    coder.internal.errorIf(colsX ~= rowsY, 'MATLAB:pagemtimes:mismatchPageDimensions');
    rowsZ = rowsX;
    colsZ = colsY;
end

numOutputDims = max(coder.internal.ndims(X), coder.internal.ndims(Y));
outputSize = coder.nullcopy(zeros(1, numOutputDims, 'coder.internal.indexInt'));
outputSize(1:2) = [rowsZ colsZ];
coder.unroll();
for ii = 3 : numOutputDims
    sizeX = coder.internal.indexInt(size(X, ii));
    sizeY = coder.internal.indexInt(size(Y, ii));
    if sizeX ~= 1
        if sizeY ~= 1
            coder.internal.errorIf(sizeX ~= sizeY, 'MATLAB:pagemtimes:mismatchImplicitExpansion');
        end
        outputSize(ii) = sizeX;
    else
        outputSize(ii) = sizeY;
    end
end



function [rows, cols] = getRowsAndCols(X, transp)

coder.internal.prefer_const(transp);
coder.inline("always");
if transp == 'N'
    rows = coder.internal.indexInt(size(X, 1));
    cols = coder.internal.indexInt(size(X, 2));
else
    rows = coder.internal.indexInt(size(X, 2));
    cols = coder.internal.indexInt(size(X, 1));
end



function checkArrayType(X)

coder.inline("always");
% Input must be full and of type single or double. Subclasses of single or
% double are not allowed, so can't just call "isfloat"
coder.internal.assert(strcmp(class(X),'double') || strcmp(class(X),'single'), ...
    'MATLAB:pagemtimes:invalidArray'); %#ok<STISA>
coder.internal.errorIf(issparse(X), 'MATLAB:pagemtimes:sparseNotSupported');



function valOut = processTransposeValue(valIn)

coder.inline("always");

valTmp = validatestring(valIn, {'none', 'transpose', 'ctranspose'}, "pagemtimes");
if strcmp(valTmp, 'none')
    valOut = 'N';
elseif strcmp(valTmp, 'transpose')
    valOut = 'T';
else
    valOut = 'C';
end

