Таким образом, мой вопрос о том, как реализовать матричную копию между двумя блочно-циклически распределенными матрицами в двух разных технологических сетях в Scalapack (BLACS). Я пытаюсь реализовать это с помощью pdgemr2d_, который я часто использую в других случаях, когда копирую между двумя матрицами в одной и той же технологической сетке.
Ниже приведено довольно техническое обсуждение состояния проблемы, с которой я сталкиваюсь. Я прибил это к основной проблеме, но мне кажется, что нет решения ... хотя оно должно быть, поскольку Скалапак определенно заявляет, что тип операции, для которой я пытаюсь, возможен. Я не нахожу адекватных примеров этого нигде.
Инициализация вычислительной сетки 1x1 в Scalapack, работающей с MPI в C, обычно происходит примерно так:
[...]
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 участвует только один процессор, и его идентификаторы сетки: {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 */
(Мы увидим позже, почему необходимы переменные "protect", так как на некоторых процессорах num_cols_local или num_rows_local будут возвращаться, совершенно правильно, как отрицательные целые числа.)
Большинство вышесказанного говорит само за себя, за исключением нулей, переданных в descinit_, которые указывают строку процессора, в которой распределена первая строка матрицы, и столбец процессора, в котором распределен первый столбец. Эти значения имеют очень четкие границы при использовании в функции descinit_. Из самой функции Фортрана,
[...]
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}, поскольку первый блок процессора, скорее всего, будет хранить значения первой строки и столбца.
При запуске на одном процессоре это работает очень хорошо. Значения на единственном процессоре RANK 0 для NPROW, NPCOL, MYPROW и MYPCOL равны 1, 1, 0 и 0 соответственно. CONTEXT в этом случае равен 0, его неотрицательность указывает на то, что сетка, на которую он ссылается, активна на этом RANK. Эти значения указывают на наличие сетки процесса 1x1, а первый процессор имеет правильное значение, указывает на правильный блок сетки процесса, принадлежащий RANK 0. В этом случае это единственный блок.
Однако при работе на двух процессорах все выходит из строя и формально не должно. На первом и втором рангах у нас есть для КОНТЕКСТ, NPROW, NPCOL, MYPROW и MYCOL:
RANK 0: 0, 1, 1, 0, 0
RANK 1: -1, -1, -1, -1, -1
Все значения отрицательные. Самое главное, что CONTEXT на RANK 1 отрицательный, что указывает на то, что этот RANK не участвует в нашей сетке процессоров 1x1. Вызов descinit_ сейчас сразу же становится проблемой для всех процессоров. Если мы ссылаемся на код Фортрана из descinit_, мы имеем (повторяется сверху для ясности):
[...]
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_, теоретически, делает это элегантно, поскольку значения CONTEXT, возвращаемые Cblacs_gridinit, равны -1 на любом процессоре, не участвующем в сетке этого контекста.Однако descinit_ не будет генерировать правильный дескриптор матрицы на любом процессоре, который не участвует в сетке из-за ограничения, детализированного выше для отрицательных значений NPROW и NPCOL.
Для обеспечения правильной связи по непересекающейся сетке такой матричный дескриптор должен быть определен на всех процессорах, которые участвуют в любом контексте.
Очевидно, что pdgemr2d_ не может быть записан с этим как непреодолимый недостаток, так как описание функции в коде специально гласит:
PDGEMR2D копирует подматрицу A на подматрицу B. A и B могут иметь разные распределения: они могут быть на разных процессорных сетках, они могут иметь разные размеры блоков, начало области для копирования может находиться в разных местах на A и B.
Большое спасибо за любую помощь, я знаю, что это довольно специализированный вопрос.