0

总之，我的问题是关于如何在Scalapack（BLACS）中的两个不同进程网格上的两个块循环分布矩阵之间实现矩阵复制。我正在尝试使用 pdgemr2d_ 来实现这一点，我经常在其他情况下使用它，我在同一进程网格上的两个矩阵之间进行复制。

下面是对我遇到的问题状态的相当技术性的讨论。我已经把它归结为一个基本问题，但在我看来没有解决方案......虽然必须有，因为 Scalapack 特别指出我正在尝试的操作类型是可能的。我在任何地方都找不到足够的例子。

在 C 中使用 MPI 运行的 Scalapack 中的 1x1 计算网格的初始化通常如下所示：

[...]
int NUM_TASKS, RANK;

int MPI_STARTUP = MPI_Init (&argc, &argv);
if (MPI_STARTUP != MPI_SUCCESS)
  MPI_Abort (MPI_COMM_WORLD, MPI_STARTUP);

MPI_Comm_size (MPI_COMM_WORLD, &NUM_TASKS);
MPI_Comm_rank (MPI_COMM_WORLD, &RANK);

int CONTEXT;

Cblacs_pinfo (&RANK, &NUM_TASKS);
Cblacs_get (-1, 0, &CONTEXT);
Cblacs_gridinit (&CONTEXT, "Row", 1, 1);
Cblacs_gridinfo (CONTEXT, &NPROW, &NPCOL, &MYPROW, &MYPCOL);
[...]

此代码将生成一个 1x1 网格，无论 MPI 知道多少处理器（{1, 1}，网格的大小，被传递给 Cblacs_gridinit）。在这里，CONTEXT 向 Scalapack 函数指示我们正在处理哪个网格（可以同时使用多个网格，并且由 Cblacs_get 生成）。Cblacs_gridinfo 将 NPROW 和 NPCOL 设置为处理器行数和列数（在本例中为 {1, 1}）。MYPROW 和 MYPCOL 向每个处理器指示哪个网格块属于它。在这种情况下，在 1x1 网格上，只有一个处理器参与，其网格 ID 为 {0, 0}。

为简单的块循环分布式 100x100 矩阵初始化矩阵描述符通常也很简单：

int info;
int desc[9];
int block_size = 32;
int zero = 0; int one = 1;
int num_rows = 100; int num_cols = 100;

int num_cols_local = numroc_ (&num_cols, &block_size, &mypcol, &zero, &npcol);
int num_cols_local_protect = MAX (1, num_cols_local);
int num_rows_local = numroc_ (&num_rows, &block_size, &myprow, &zero, &nprow);
int num_rows_local_protect = MAX (1, num_rows_local);
descinit_ (desc, &num_rows, &num_cols, &block_size, &block_size, &zero, &zero, \
  &CONTEXT, &num_rows_local_protect, &info);

/* now allocate array with per-processor size num_cols_local_protect * num_rows_local_protect */

（稍后我们将看到为什么需要“保护”变量，因为在某些处理器上 num_cols_local 或 num_rows_local 将非常正确地返回为负整数。）

上面的大部分内容都是不言自明的，除了传递给 descinit_ 的 &zeros，它表示矩阵的第一行所在的处理器行，以及第一列所在的处理器列。这些值在 descinit_ 函数中使用时具有非常明确的界限。从 Fortran 函数本身来看，

[...]
ELSE IF( IRSRC.LT.0 .OR. IRSRC.GE.NPROW ) THEN
  INFO = -6
ELSE IF( ICSRC.LT.0 .OR. ICSRC.GE.NPCOL ) THEN
  INFO = -7
[...]

我们在这里将 IRSRC 和 ICSRC 作为零传递，因为 {0,0} 是我们单个网格块的正确索引。即使网格大得多，我们可能仍会传递 {0,0}，因为第一个处理器块可能会存储第一行和第一列的值。

在一个处理器上运行时，效果很好。NPROW、NPCOL、MYPROW 和 MYPCOL 的唯一处理器 RANK 0 上的值分别为 1、1、0 和 0。在这种情况下，CONTEXT 为 0，它的非负性表明它所指的网格在此 RANK 上是活动的。这些值表示存在 1x1 进程网格，并且第一个处理器具有正确的表示属于 RANK 0 的正确进程网格块。在这种情况下，它是唯一的块。

然而，当在两个处理器上运行时，事情会发生故障，而且它们不应该正式出现。在第一个和第二个 RANK 中，我们有 CONTEXT、NPROW、NPCOL、MYPROW 和 MYCOL：

RANK 0: 0, 1, 1, 0, 0
RANK 1: -1, -1, -1, -1, -1

所有值都是负数。最重要的是，RANK 1 上的 CONTEXT 为负数，表明此 RANK 不参与我们的 1x1 处理器网格。现在调用 descinit_ 立即成为所有处理器的问题。如果我们从 descinit_ 引用 Fortran 代码，我们有（为了清楚起见从上面重复）：

[...]
ELSE IF( IRSRC.LT.0 .OR. IRSRC.GE.NPROW ) THEN
  INFO = -6
ELSE IF( ICSRC.LT.0 .OR. ICSRC.GE.NPCOL ) THEN
  INFO = -7
[...]

只要每个处理器都参与网格，这些限制就有意义。索引不能为负数或大于或等于进程网格中的总行数或列数，因为这样的网格块不存在！

然后在 RANK 1 上，IRSRC 作为零传递，但 NPROW 和 NPCOL 从网格初始化返回为 -1，因此 descinit_ 将始终失败。

通过简单地将矩阵描述符的初始化和所有后续操作限制为参与当前网格的处理器，可以轻松地克服上述所有问题。就像是：

if (CONTEXT > -1) {
[...]

但是，我不只有一个处理器网格，我有两个，我需要它们通过 pdgemr2d_ 函数进行通信。此函数的目的是将一个网格上的分布式矩阵 A 的子集复制到另一个网格上的分布式矩阵 B。网格不需要以任何方式相互关联，并且可以部分或完全不相交。这应该是一个微不足道的操作。例如，我想将具有上下文 CONTEXT_A 的处理器网格中的完整矩阵复制到具有上下文 CONTEXT_B 的处理器网格中。每个上下文中矩阵的描述符以 desc_A 和 desc_B 的形式给出。

pdgemr2d_ (&num_rows, &num_cols, matrix_A, &one, &one, \
  desc_A, matrix_B, &one, &one, desc_B, &CONTEXT_B);

这也是相当不言自明的。它必须在任一上下文具有任何网格成员的所有处理器上运行。在我的例子中，CONTEXT_A 有一个跨越 MPI 知道的所有处理器的网格，而 CONTEXT_B 是一个 1x1 单处理器网格。

pdgemr2d_ 必须提供一个上下文标识符，至少包含包含在 CONTEXT_A 和 CONTEXT_B 中的所有处理器，并且对于那些不属于 CONTEXT_A 或 CONTEXT_B 的处理器，元素 desc_A[CTXT] 或 desc_B[CTXT] 必须分别设置为 - 1 在那个处理器上。

理论上，descinit_ 可以优雅地做到这一点，因为 Cblacs_gridinit 返回的 CONTEXT 值在任何不参与该上下文网格的处理器上都是 -1。但是，由于上面详述的 NPROW 和 NPCOL 负值的限制，descinit_ 不会在任何不参与网格的处理器上生成正确的矩阵描述符。

为了进行适当的不相交网格通信，必须在参与任一上下文的所有处理器上定义这样的矩阵描述符。

显然，pdgemr2d_ 不能用这个作为一个不可克服的缺陷来编写，因为代码中的函数描述特别指出：

PDGEMR2D 将 A 的子矩阵复制到 B 的子矩阵上。A 和 B 可以有不同的分布：它们可以在不同的处理器网格上，它们可以有不同的块大小，要复制的区域的开头可以在 A 上的不同位置和 B。

非常感谢您的帮助，我知道这是一个相当专业的问题。

4

1 回答 1

2

我在尝试弄清楚如何使用 PDGEMR2D 时遇到了类似的困难，我想在这里分享我的结论。

简而言之，如果您尝试设置多个上下文，请不要使用提供的 DESCINIT 子例程。它的错误检查假定所有进程都参与了正在初始化的数组描述符的上下文，如果您尝试使用 PDGEMR2D，情况并非如此。

您可以轻松地初始化自己的描述符，而无需使用 DESCINIT，因为它只是一个长度为 9 的整数数组。前 8 个字段（dtype、ctxt、m、n、mb、nb、csrc 和 rsrc）是全局的，并且应该具有所有进程的值相同。只有第 9 个字段 (lld) 是本地的，并且在不是定义数组的上下文成员的进程上忽略它的值。

ScaLAPACK 源代码（此处为在线版本）中的示例程序 pdgemrdrv.c 在我试图弄清楚这些东西时对我很有用。它包含很多不必要的复杂性，但你可以推断出以下关键点：

  1. 数组描述符应该全局初始化（在发送和接收上下文中）。
  2. 数组的内存只需要在数组所在的上下文中分配。
  3. PDGEMR2D 应该被全局调用（在发送和接收上下文中）。只有作为发送或接收上下文成员的进程才会真正参与。

希望这可以帮助。干杯。

于 2012-02-20T20:20:47.757 回答