MPI域的分解和块Tridiagonal矩阵

Question

我需要一些建议，即如何构建用MPI分解2D域的块三基因矩阵。 让我解释一下，我需要在2D中求解像方程一样的热扩散，然后使用虚拟的笛卡尔拓扑结构在用MPI的块中分解了域。 报道了分解的域 在这里，其中单元格中的数字代表全局索引中的单元格坐标。

例如，当我使用有限差异分散方程式时，我获得了一个块tridiagonal矩阵，例如矩阵中的数字再次是与以前相同的全局索引中的单元格坐标。

腮腺

• 作为
• 我使用
• 求解器

， 1D域（导致Tridiagonal矩阵）也很容易。 问题在于，我无法在每个过程中找到与DO循环分配的清晰的索引序列，因为离散化后计算域（第一个图像）的位置导致矩阵中的非连续元素（第二张图像）。

我知道图形重新排序是由诸如Parmetis之类的库执行的，但是据我所知，它并不能保留原始域分解模式，因为它在矩阵级别上进行了重新定位。 我认为这是一个非常普遍的问题，这应该是微不足道的，但是我看不出一种清晰的方法。

Answer 1

您的分区域基于顺序排序。将其转换为本地订购很棘手。如果您首先在Proc 0上编号所有点会更容易，则在Proc 1上的所有点等。这为您提供了一个具有略有不同结构的矩阵，但是相同的解决方案和很大相同的属性。

Answer 2

在维克多的建议下，我在“本地”专栏的主要订购后重新排序了矩阵的全局索引。
现在，物理域被映射为：
单击
其中4个左上元素属于该过程，等级为0，左下4个元素，该过程具有等级1等的过程。
在将其移植到Fortran之前，我创建了一个草稿代码来组装矩阵。如果有人遇到同样的问题，我会留在这里。该代码可能会经常简化。

clc 
clear
% Global number of rows and columns 
gRow = 4;
gCol = 4;
% Number of processor along rows and columns
rProc = 2;
cProc=2;
% END OF USER INPUTS
rLoc = gRow/rProc;
cLoc = gCol/cProc;

%Global Indexing following column major
M = zeros(gRow,gCol);
gidx = 1;
for l = 1:cProc
    for i=1:rProc
        for j=1:cLoc
            for k=1:rLoc
                M((i-1)*rLoc+k,(l-1)*cLoc+j) = gidx;
                gidx = gidx+1;
            end
        end
    end
end
rowidx = []; % rows of non zero elements 
colidx=[]; % Columns of non zero elements
gidx = 1;
for l = 1:cProc
    for i=1:rProc
        for j=1:cLoc
            for k=1:rLoc
                if (i-1)*rLoc+k ~=1 && (i-1)*rLoc+k ~=gRow && (l-1)*cLoc+j~=1 && (l-1)*cLoc+j~=gCol
                    A =  [M((i-1)*rLoc+k-1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j+1)...
                        M((i-1)*rLoc+k+1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j-1)];
                    A = sort(A);

                elseif (i-1)*rLoc+k ==1 && (i-1)*rLoc+k ~=gRow && (l-1)*cLoc+j~=1 && (l-1)*cLoc+j~=gCol
                    A =  [ M((i-1)*rLoc+k,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j+1)...
                        M((i-1)*rLoc+k+1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j-1)];
                    A = sort(A);
                elseif (i-1)*rLoc+k ~=1 && (i-1)*rLoc+k ==gRow && (l-1)*cLoc+j~=1 && (l-1)*cLoc+j~=gCol
                    A =  [M((i-1)*rLoc+k-1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j+1)...
                        M((i-1)*rLoc+k,(l-1)*cLoc+j-1)];
                    A = sort(A);
                elseif (i-1)*rLoc+k ~=1 && (i-1)*rLoc+k ~=gRow && (l-1)*cLoc+j==1 && (l-1)*cLoc+j~=gCol
                    A =  [M((i-1)*rLoc+k-1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j+1)...
                        M((i-1)*rLoc+k+1,(l-1)*cLoc+j) ];
                    A = sort(A);
                elseif (i-1)*rLoc+k ~=1 && (i-1)*rLoc+k ~=gRow && (l-1)*cLoc+j~=1 && (l-1)*cLoc+j==gCol
                    A =  [M((i-1)*rLoc+k-1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j)...
                        M((i-1)*rLoc+k+1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j-1)];
                    A = sort(A);
                elseif (i-1)*rLoc+k ==1 && (l-1)*cLoc+j==1
                    A =  [ M((i-1)*rLoc+k,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j+1)...
                        M((i-1)*rLoc+k+1,(l-1)*cLoc+j) ];
                    A = sort(A);
                elseif (i-1)*rLoc+k ==1 && (l-1)*cLoc+j==gCol
                    A =  [ M((i-1)*rLoc+k,(l-1)*cLoc+j)...
                        M((i-1)*rLoc+k+1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j-1)];
                    A = sort(A);
                elseif (i-1)*rLoc+k ==gRow && (l-1)*cLoc+j==1
                    A =  [M((i-1)*rLoc+k-1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j+1)];
                    A = sort(A);
                elseif (i-1)*rLoc+k ==gRow && (l-1)*cLoc+j==gCol
                    A =  [M((i-1)*rLoc+k-1,(l-1)*cLoc+j) M((i-1)*rLoc+k,(l-1)*cLoc+j)...
                        M((i-1)*rLoc+k,(l-1)*cLoc+j-1)];
                    A = sort(A);
                end
                colidx = cat(2,colidx, ones(1,length(A))*gidx);
                rowidx = cat(2,rowidx,A);
                gidx = gidx+1;
        
            end
        end
        %UNCOMMENT IF YOU WANT TO PLOT DATA OWNED BY EACH PROCESS
%         S = sparse(rowidx,colidx,ones(1,length(colidx)),gRow*gCol,gRow*gCol);
%         figure
%         spy(S)
%         rowidx = [];
%         colidx = [];
    end
end
        S = sparse(rowidx,colidx,ones(1,length(colidx)),gRow*gCol,gRow*gCol);
        figure
        spy(S)

唯一要传递的输入是全局行和列的数量以及行沿行和列的过程数。该代码已使用不同数量的全局行和列以及沿两个方向进行不同数量的过程进行了测试。它应该是灵活的。
但是，请务必检查是否正常工作。
通过删除末端附近的部分，您可以绘制矩阵中每个过程所拥有的元素。

如果绘制了整个矩阵，它将看起来像 this （无线或文本框）