MPI 域分解和分块三对角矩阵

Question

我需要一些关于如何使用 MPI 分解二维域来构建分块三对角矩阵的建议。让我解释一下，我需要求解类似热扩散的二维方程，为此我使用虚拟笛卡尔拓扑通过 MPI 将域分解为块。报告分解域 here，其中单元格中的数字表示全局索引中的单元格坐标。

例如，当我用有限差分离散化我的方程时，我得到了一个像 this 这样的块三对角矩阵，其中矩阵中的数字再次是与以前相同的全局索引中的单元格坐标。

我使用 MUMPS 作为求解器，它需要为每个过程输入具有三个向量的 spars 矩阵，包含：

全局行索引
全局列索引
对应元素条目（题项系数）

在序列中，如何做到这一点很简单，在具有一维域（导致三对角矩阵）的 MPI 中也很容易。问题是我无法在每个进程上找到要分配给 do 循环的清晰索引序列，因为离散化后计算域（第一张图像）中的位置导致矩阵中的非连续元素（第二张图像).

我知道图形重新排序是由 Parmetis 等库执行的，但据我了解，它不会保留原始域分解模式，因为它在矩阵级别重新排序。我认为这应该是微不足道的，因为这是一个非常普遍的问题，但我看不出一个明确的方法来做到这一点。

Answer 1

您的分区域是基于顺序排列的。将其转换为本地订单有点棘手。如果您首先对 proc 0 上的所有点进行编号，然后对 proc 1 上的所有点进行编号，等等，这样会更容易。这为您提供了一个结构略有不同但解相同且属性基本相同的矩阵。

Answer 2

根据 Victor 的建议，我按照“本地”列主要排序对矩阵的全局索引进行了重新排序。现在物理域映射为： click 其中左上角的 4 个元素属于等级为 0 的过程，左下角的 4 个元素属于等级为 1 的过程，依此类推。在将它移植到 Fortran 之前，我在 Matlab 上为 assemble 矩阵创建了一个代码草案。如果有人有同样的问题，我会把它留在这里。代码大概可以简化很多。

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（没有线条或文本框）

MPI 域分解和分块三对角矩阵

MPI Domain Decomposition and Block Tridiagonal Matrix

matrix

fortran

mpi