%CODISTRIBUTOR2DBC   2D block-cyclic codistributor object
%   The 2D block-cyclic codistributor can only be used for two dimensional
%   matrices. It distributes matrices along two subscripts over a rectangular
%   computational grid of workers in a blocked, cyclic manner. The 2D block-cyclic
%   codistributor is employed by the parallel matrix computation software
%   library called ScaLAPACK.
%
%   Example: Look at how a small matrix is distributed using the 2D block-cyclic
%   distribution scheme.  For the purposes of the example, we choose the matrix
%   size such that each worker stores exactly one block.
%       spmd
%           sz = codistributor2dbc.defaultBlockSize*codistributor2dbc.defaultWorkerGrid;
%           codistr = codistributor2dbc(codistributor2dbc.defaultWorkerGrid, ...
%                                       codistributor2dbc.defaultBlockSize, ...
%                                       codistributor2dbc.defaultOrientation, ...
%                                       sz);
%           [e1, f1] = codistr.globalIndices(1);
%           [e2, f2] = codistr.globalIndices(2);
%           fprintf('Worker %d stores elements [%d:%d, %d:%d].\n', ...
%                   spmdIndex, e1, f1, e2, f2);
%       end
%   When run on 4 workers, the matrix size is [128, 128], and this shows:
%       Worker 1 stores elements [1:64, 1:64].
%       Worker 2 stores elements [1:64, 65:128].
%       Worker 3 stores elements [65:128, 1:64].
%       Worker 4 stores elements [65:128, 65:128].
%   That is, the 128-by-128 matrix is split into 4 blocks, and they are divided
%   amongst the workers as:
%      |   Block stored on worker 1   |    Block stored on worker 2   |
%      |   Block stored on worker 3   |    Block stored on worker 4   |
%
%   Example: We can see the cyclic nature of the distribution scheme when we
%   work with a larger matrix than in the previous example.
%       spmd
%           sz = 3*codistributor2dbc.defaultBlockSize*codistributor2dbc.defaultWorkerGrid;
%           codistr = codistributor2dbc(codistributor2dbc.defaultWorkerGrid, ...
%                                       codistributor2dbc.defaultBlockSize, ...
%                                       codistributor2dbc.defaultOrientation, ...
%                                       sz);
%           [e1, f1] = codistr.globalIndices(1);
%           rows = strtrim(sprintf('%d:%d ', [e1; f1]));
%           [e2, f2] = codistr.globalIndices(2);
%           cols = strtrim(sprintf('%d:%d ', [e2; f2]));
%           fprintf('Worker %d stores elements [%s, %s].\n', ...
%                    spmdIndex, rows, cols);
%       end
%   When run on 4 workers, the matrix is of size 384-by-384, and this shows:
%       Worker 1 stores elements [1:64 129:192 257:320, 1:64 129:192 257:320].
%       Worker 2 stores elements [1:64 129:192 257:320, 65:128 193:256 321:384].
%       Worker 3 stores elements [65:128 193:256 321:384, 1:64 129:192 257:320].
%       Worker 4 stores elements [65:128 193:256 321:384, 65:128 193:256 321:384].
%
%   codistributor2dbc methods:
%     codistributor2dbc/codistributor2dbc  - Create 2D block-cyclic codistributor object
%     globalIndices    - Global indices for the local part corresponding to codistributor
%     isComplete       - Return true if codistributor has all of its information set
%     CELL    - Create codistributed cell array using codistributor
%     COLON   - Build codistributed arrays of the form j:d:k using codistributor
%     EYE     - Identity codistributed matrix using codistributor
%     FALSE   - False codistributed array using codistributor
%     INF     - Infinity codistributed array using codistributor
%     NAN     - Build codistributed array containing Not-a-Number using codistributor
%     ONES    - Ones codistributed array using codistributor
%     RAND    - codistributed array of uniformly distributed pseudorandom numbers using codistributor
%     RANDN   - codistributed array of normally distributed pseudorandom numbers using codistributor
%     SPALLOC - Allocate space for sparse codistributed matrix using codistributor
%     SPARSE  - Create sparse codistributed matrix using codistributor
%     SPEYE   - Overloaded to create a codistributed array
%     SPRAND  - Sparse uniformly distributed random codistributed matrix using codistributor
%     SPRANDN - Sparse normally distributed random codistributed matrix using codistributor
%     TRUE    - True codistributed array using codistributor
%     ZEROS   - Zeros codistributed array using codistributor
%
%   codistributor2dbc static methods and properties:
%     defaultBlockSize   - MATLAB's choice for ScaLAPACK block size
%     defaultWorkerGrid  - MATLAB's choice for computational grid
%     defaultOrientation - MATLAB's choice for ScaLAPACK orientation
%
%   codistributor2dbc object properties:
%     BlockSize   - Block size of a 2D block-cyclic codistributor
%     WorkerGrid  - Worker grid of a 2D block-cyclic codistributor
%     Orientation - The orientation of worker grid of a 2D block-cyclic codistributor
%
%   See also CODISTRIBUTOR, CODISTRIBUTED, CODISTRIBUTOR1D.

%   Copyright 2007-2022 The MathWorks, Inc.

classdef codistributor2dbc < TensorProductCodistributor
    properties(SetAccess = private)
        %BlockSize  Block size of a 2D block-cyclic codistributor
        %   codistr.BlockSize is the ScaLAPACK block size associated
        %   with the 2D block-cyclic codistributor codistr.
        %
        %   See also codistributor2dbc, codistributor2dbc.defaultBlockSize.
        BlockSize

        %WorkerGrid Worker grid of a 2D block-cyclic codistributor
        %   codistr.WorkerGrid is the worker grid of a 2D block-cyclic codistributor
        %   codistr. The worker grid is the row vector of length 2, [nprow,npcol], 
        %   used by the ScaLAPACK library. The product of codistr.WorkerGrid must 
        %   equal spmdSize. he default is given by 
        %   codistributor2dbc.defaultWorkerGrid.
        %
        %   See also codistributor2dbc, codistributor2dbc.defaultWorkerGrid.
        WorkerGrid

        %Orientation The orientation of worker grid of a 2D block-cyclic codistributor
        %   codistr.Orientation is the orientation of the worker grid of a 2D
        %   block-cyclic codistributor codistr.  The supported values are
        %   'row' or 'col', which means that the workers are laid out in the
        %   worker grid in either a row-wise or column-wise manner, respectively.
        %
        %   Example:
        %      With Orientation equal to 'row' and WorkerGrid equal to [2, 3], the
        %      workers are laid out in the worker grid as follows:
        %      |   worker 1  |  worker 2   |  worker 3 |
        %      |   worker 4  |  worker 5   |  worker 6 |
        %
        %   See also codistributor2dbc, codistributor2dbc.defaultOrientation.
        Orientation
    end

    properties (Hidden, Dependent)
        LabGrid % Discouraged - use WorkerGrid
    end

    methods(Access = private, Static = true)
        function ex = notImplemented(fcnstr)
            ex = MException(message(...
                'parallel:codistributors:c2dbcMethodNotImplemented', fcnstr));
        end
    end

    methods (Access = private)
        function col = pSpmdIndexToProcessorCol(codistr, workerIdx)
        % Map workerIdx of workers into the column in the processor grid.
            switch codistr.Orientation
              case 'row'
                col = mod(workerIdx - 1, codistr.WorkerGrid(2)) + 1;
              case 'col'
                col = ceil(workerIdx/codistr.WorkerGrid(1));
              otherwise
                error(message('parallel:codistributors:c2dbcSpmdIndexToProcessorColOrientUnsupported', codistr.Orientation));
            end
        end

        function row = pSpmdIndexToProcessorRow(codistr, workerIdx)
        % Map workerIdx of the workers into the row in the processor grid.
            switch codistr.Orientation
              case 'row'
                row = ceil(workerIdx/codistr.WorkerGrid(2));
              case 'col'
                row = mod(workerIdx - 1, codistr.WorkerGrid(1)) + 1;
              otherwise
                error(message('parallel:codistributors:c2dbcSpmdIndexToProcessorRowOrientUnsupported', codistr.Orientation));
            end
        end
        [data, isOnEdge] = pCatToEdgeOfWorkerGrid(codistr, data, dim)
        isOnEdge = pIsOnWorkerGridEdge(codistr, dim)
    end

    methods(Access = protected)
        function codistr = pTransposeCodistributor(codistr)
            sz = codistr.Cached.GlobalSize;
            wkrGrid = codistr.WorkerGrid;

            % (c)transpose without communication relies on transposing
            % the Worker Grid and using the opposite orientation to set up
            % the new codistributor.
            orient = codistr.Orientation;
            if strcmp(orient, 'row')
                altOrient = 'col';
            else
                altOrient = 'row';
            end

            codistr = codistributor2dbc([wkrGrid(2) wkrGrid(1)], ...
                                        codistr.BlockSize, ...
                                        altOrient, [sz(2), sz(1)]);
        end % End pTransposeCodistributor
        
        function attrs = pGetCodistAttributes(codistr)
            attrs = {class(codistr), ...
                codistr.WorkerGrid, codistr.BlockSize, codistr.Orientation, ...
                codistr.Cached.GlobalSize};
        end

    end

    % Implementation of hidden, abstract methods.
    methods(Hidden = true)
        [LP, codistr] = hBuildFromReplicatedImpl(codistr, srcWorker, X)
        tf = hCatCheck(codistr, catdim, cellOfCodistr, cellOfLPs)
        [LP, codistr] = hCatImpl(codistr, catdim, codistrs, ...
                                 LPs, templ, targetSize)

        [LP, codistr] = hCell2StructImpl(codistr, LP, fields, dim)

        [LP, codistr] = hColonImpl(codistr, a, d, b)
        function tf = hDoGatherBsxfunInput(codistr, szInput)
            tf = any(szInput == 1) && ~isequal(szInput, codistr.Cached.GlobalSize);
        end

        function tf = hFftnCheck(codistr, fftSize)
            tf = (max(fftSize)<=codistr.BlockSize && numel(fftSize)<3);
        end
        [LP, codistr] = hFftnImpl(codistr, LP, fftSize, sym, isFFT)

        localLinInd = hFindDiagElementsInLocalPart(codistr)
        X = hGatherImpl(codistr, LP, destWorker)
        function LP = hGetBsxfunLocalPartFromReplicated(codistr, replicatedInput)
            szIn = size(replicatedInput);
            szOut = codistr.Cached.GlobalSize;
            ind = {1,1};
            for dim = 1:2
                % The test below on the size of the output codistributor
                % in the distribution dimension is necessary because 
                % if both inputs are singleton in the distribution 
                % dimension they should not be gathered.   Since arrays are 
                % distributed in both dimensions by codistributor2d, 
                % the test must be applied to both dimensions.
                if szIn(dim) ~= 1 || szOut(dim) == 1
                    ind{dim} = codistr.globalIndices(dim);
                end
            end
            LP = replicatedInput(ind{:});
        end
        
        % hGetColumn Return the column of a 2D matrix with global index 
        % globalColIndex
        column = hGetColumn(codistr, LP, globalColIndex)
        
        function codistr = hGetCompleteForSize(codistr, wantedGlobalSize)
        % Return a codistributor that is based on the input codistributor and matches
        % the desired global size.
            if ~codistr.isComplete() || ~isequal(codistr.Cached.GlobalSize, wantedGlobalSize)
                codistr = codistributor2dbc(codistr.WorkerGrid, codistr.BlockSize, ...
                                            codistr.Orientation, wantedGlobalSize);
            end
        end % End of hGetCompleteForSize.
        function codistr = hGetNewForSize(codistr, wantedGlobalSize)
        % Return a codistributor that is based on the input codistributor
        % and matches the desired global size.
            codistr = hGetCompleteForSize(codistr, wantedGlobalSize);
        end
        function dims = hGetDimensions(~)
            dims = [1, 2];
        end
        function [LP, codistr] = hGetRepmatCatDimensionOrdering(codistr, LP, sz, replFcn)
            [~, index] = sort(sz);
            for d = index
                if sz(d) ~= 1
                    LP = replFcn(LP, d, sz(d));
                end
            end
        end

        % hGetRow Return the row of a 2D matrix with global index 
        % globalRowIndex
        row = hGetRow(codistr, LP, globalRowIndex)
        
        varargout = hGlobalIndicesImpl(codistr, dim, workerIdx)
        localIdx = hGlobalToLocalIndices(codistr, globalIdx, dim, workerIdx)

        tf = hIsGlobalIndexOnWorker(codistr, dim, gIndexInDim, workerIdx)
        function col = hSpmdIndexToProcessorCol(codistr, workerIdx)
            col = codistr.pSpmdIndexToProcessorCol(workerIdx);
        end
        function row = hSpmdIndexToProcessorRow(codistr, workerIdx)
            row = codistr.pSpmdIndexToProcessorRow(workerIdx);
        end
        szs = hLocalSize(codistr, workerIdx)
        globalIdx = hLocalToGlobalIndices(codistr, localIdx, dim, workerIdx)
        globalLinIdx = hLocalToGlobalLinearIndices(codistr, localLinIdx, workerIdx)

        [LPY, LPI, codistr] = hMinMaxImpl(codistr, fcnMinMax, LP, dim, wantI, linearI)
        function tf = hMtimesCheck(codistr, areBothMatricesDense) %#ok<INUSL>
            tf = areBothMatricesDense;
        end

        [LP, codistr] = hMtimesImpl(codistr, codistrA, LPA, transA, ...
            codistrB, LPB, transB)

        function tf = hMtimesReplicatedCheck(codistr, areBothMatricesDense) %#ok<INUSL>
            tf = areBothMatricesDense;
        end

        [LP, codistr] = hMtimesReplicatedImpl(codistr, codistrA, LPA, transA, ...
            codistrB, LPB, transB)

        function numWorkers = hNumWorkers(codistr)
            numWorkers = prod(codistr.WorkerGrid);
        end

        function tf = hPermuteCheck(codistr, lenOrder) %#ok<INUSL>
            tf = lenOrder <= 2;
        end

        function [LP, codistr] = hPermuteImpl(codistr, LP, order)
            if isequal(order, [2 1])
                [LP, codistr] = codistr.hTransposeTemplateImpl(LP, @transpose);
            end
        end

        [LP, codistr] = hReductionOpAlongDimImpl(codistr, fcn, LP, dim)
        
        % 2D block-cyclic codistributors do not support sparse arrays, so
        % we always throw an error for h(Spalloc/Sparsify)Impl.
        function [LP, codistr] = hSpallocImpl(codistr, m, n, nzmx) %#ok<INUSD,STOUT>
            ex = MException(message('parallel:codistributors:c2dbcUnsupportedSparse'));
            throwAsCaller(ex);
        end
        function [LP, codistr] = hSparsifyImpl(codistr, fcn, LP) %#ok<INUSD>
            ex = MException(message('parallel:codistributors:c2dbcUnsupportedSparse'));
            throwAsCaller(ex);
        end

        function tf = hSupportsDimensionality(codistr, ndim) %#ok<INUSL>
            tf = ndim <= 2;
        end
        
        B = hSubsrefParensImpl(codistr, A, aidx, varargin)
        
        [LP, codistr] = hTypecastImpl(codistr, LP, datatype)

        function hVerifySupportsSparse(~)
        % 2D block-cyclic codistributors do not support sparse arrays, so 
        % we always throw an error.
            ex = MException(message('parallel:codistributors:c2dbcUnsupportedSparse'));
            throwAsCaller(ex);
        end
        
        function hVerifySupportsTable(~)
            % Tables cannot be distributed using a 2DBC scheme
            ex = MException(message('parallel:codistributors:c2dbcUnsupportedTable'));
            throwAsCaller(ex);
        end

    end % End of hidden, abstract methods.

    properties (Constant = true)
        %codistributor2dbc.defaultBlockSize    MATLAB's choice for ScaLAPACK block size
        %   BLKSIZE = codistributor2dbc.defaultBlockSize
        %
        %   Example:
        %      codistr = codistributor2dbc()
        %   returns a distribution scheme with codistr.BlockSize set to
        %   codistributor2dbc.defaultBlockSize.
        %
        %   See also codistributor2dbc, codistributor2dbc/defaultWorkerGrid,
        %   codistributor2dbc/defaultOrientation
        defaultBlockSize = 64;

        %codistributor2dbc.defaultOrientation    MATLAB's choice for ScaLAPACK orientation
        %    orient = codistributor2dbc.defaultOrientation
        %
        %   Example:
        %      codistr = codistributor2dbc()
        %   returns a distribution scheme with codistr.Orientation set to
        %   codistributor2dbc.defaultOrientation.
        %
        %   See also codistributor2dbc, codistributor2dbc/defaultWorkerGrid,
        %   codistributor2dbc/defaultBlockSize
        defaultOrientation = 'row';
    end % End of constant properties.

    methods(Static = true)
        grd = defaultWorkerGrid();
    end

    methods(Static=true, Hidden=true)
        % Hidden as discouraged
        function dgrid = defaultLabGrid()
            dgrid = codistributor2dbc.defaultWorkerGrid();
        end
    end


    methods
        function A = codistributor2dbc(wkgrid, blksize, orientation, siz)
        %CODISTRIBUTOR2DBC   Create 2D block-cyclic codistributor object
        %   The 2D block-cyclic codistributor can be used only for two dimensional
        %   arrays. It distributes arrays along two subscripts over a rectangular
        %   computational grid of workers in a blocked, cyclic manner. The 2D block-cyclic
        %   codistributor is employed by the ScaLAPACK parallel matrix computation
        %   software library.
        %
        %   Each of the following forms a 2D block-cyclic codistributor object with
        %   its global size unspecified.  When provided, the worker grid, block size,
        %   and orientation are set to WKGRID, BLKSIZE, and ORIENT, respectively,
        %   otherwise they are set to their default values.  The resulting
        %   codistributor is incomplete as its global size is not specified.
        %   DIST = CODISTRIBUTOR2DBC()
        %   DIST = CODISTRIBUTOR2DBC(WKGRID)
        %   DIST = CODISTRIBUTOR2DBC(WKGRID, BLKSIZE)
        %   DIST = CODISTRIBUTOR2DBC(WKGRID, BLKSIZE, ORIENT)
        %   A codistributor constructed in this manner can then be used as an argument
        %   to other functions as a template codistributor when creating codistributed
        %   arrays.
        %   The default values for the worker grid, block size, and orientation are given
        %   by codistributor2dbc.defaultWorkerGrid, codistributor2dbc.defaultBlockSize,
        %   and codistributor2dbc.defaultOrientation, respectively.
        %
        %   DIST = CODISTRIBUTOR2DBC(WKGRID, BLKSIZE, ORIENT, GSIZE) forms a
        %   codistributor object that distributes arrays over worker grid WKGRID with
        %   block size BLKSIZE and the global size of the codistributed arrays being
        %   GSIZE.  The resulting codistributor object is complete and can therefore
        %   be used to build a codistributed array from its local parts with
        %   codistributed.build.
        %
        %   Example: Use a codistributor2dbc object to create an N-by-N matrix of ones.
        %   N = 1000;
        %   spmd
        %       codistr = codistributor2dbc();
        %       D = codistributed.ones(N, codistr);
        %   end
        %
        %   Example: Use a fully specified codistributor2dbc object to create a trivial
        %   N-by-N codistributed matrix from its local parts.  Then visualize which
        %   elements are stored on worker 2.
        %   N = 1000;
        %   spmd
        %       codistr = codistributor2dbc(codistributor2dbc.defaultWorkerGrid, ...
        %                                   codistributor2dbc.defaultBlockSize, ...
        %                                   codistributor2dbc.defaultOrientation, ...
        %                                   [N, N]);
        %       localSize = [length(codistr.globalIndices(1)), length(codistr.globalIndices(2))];
        %       localPart = spmdIndex*ones(localSize);
        %       D = codistributed.build(localPart, codistr);
        %   end
        %   spy(D == 2);
        %
        %   See also CODISTRIBUTOR, CODISTRIBUTED, CODISTRIBUTOR1D.

        % Error check wkgrid or provide default value.
            if nargin >= 1
                wkgrid = distributedutil.CodistParser.gatherIfCodistributed(wkgrid);
                if isValidWorkerGrid(wkgrid)
                    wkgrid = double(wkgrid);
                else
                    error(message('parallel:codistributors:c2dbcWorkerGridInput'));
                end
            else
                wkgrid = codistributor2dbc.defaultWorkerGrid();
            end

            % Error check blksize or provide default value.
            if nargin >= 2
                blksize = distributedutil.CodistParser.gatherIfCodistributed(blksize);
                if isValidBlockSize(blksize)
                    blksize = double(blksize);
                else
                    error(message('parallel:codistributors:c2dbcBlockSizeInput'));
                end
            else
                blksize = codistributor2dbc.defaultBlockSize;
            end

            % Error check orientation or provide default value.
            if nargin >= 3
                orientation = distributedutil.CodistParser.gatherIfCodistributed(orientation);
                if ischar(orientation) || (isstring(orientation) && isscalar(orientation))
                    if any( strcmpi( orientation, {'row', 'col'} ) )
                        orientation = char(lower(orientation));
                    else
                        error(message('parallel:codistributors:c2dbcOrientationInput', orientation));
                    end
                else
                    error(message('parallel:codistributors:c2dbcNonCharOrientation', class( orientation )))
                end
            else
                orientation = codistributor2dbc.defaultOrientation;
            end

            % Error check siz or provide default value.
            if nargin >= 4 && ~isempty(siz)
                siz = distributedutil.CodistParser.gatherIfCodistributed(siz);
                % Set the property so that we can benefit from all the error checking
                % performed in the parent class.
                A.Cached.GlobalSize = siz;
                siz = A.Cached.GlobalSize;
                if ~isValidSize(siz)
                    error(message('parallel:codistributors:c2dbcSizeInput'));
                end
            else
                siz = [];
            end

            A.WorkerGrid = wkgrid;
            A.BlockSize = blksize;
            A.Orientation = orientation;
            A.Cached.GlobalSize = siz;
        end % End of codistributor2dbc.

        function g = get.LabGrid(codistr)
            g = codistr.WorkerGrid;
        end
        function codistr = set.LabGrid(codistr, g)
            codistr.WorkerGrid = g;
        end
    end % End of public methods block.
end
