MPI. Умножение матрицы на вектор.

Приветствую. Давно не писал, но тут встала задача по MPI, а именно умножение матрицы на матрицу, а потом еще и на вектор. Разумеется матрицы могут быть огромными, по этому не так все просто.

И так, как это сделал я, смотрим ниже.

Как же? Что же?
Начиная решать поставленную задачу, пришел к мысли, что нужно описать предметную область, т.к. работать будем с матрицами, надо описать помощника, который бы и создавал, производил махинации и освобождал матрицу. Кто впервые слышит о MPI, то писать будем на небезызвестном языке C.

Матрица
Если подумать, что бывает иногда сложно, матрица с одной строкой или с одним столбцом есть вектор. Поэтому нас устроит данная постановка матрицы. Что есть матрица? Разумеется это массив однотипных элементов, ширина (столбцы), высота (строки) и иногда еще нужно знать размер, по сути это ширина * высота. Пусть структура матрицы будет MATRIX_STRUCT, а сама матрица MATRIX, указатель на соответствующую структуру.

Что нам понадобиться?
Для работы с нашей матрицей, нам понадобиться: создание матрицы, инициализация случайными числами, заполнение одним числом, удаление, форматированный вывод в консоль, умножение матрицы на матрицу и вывод результата, как новая матрица. В принципе этого достаточно, но если что-то потребуется еще, ничего тяжелого нет в том, чтобы добавить свои функции.

И так, какой у нас вышел Matrix.h:

#ifndef MATRIX_H
#define MATRIX_H

#define MATRIX_TYPE double

struct MATRIX_STRUCT {
    MATRIX_TYPE *items;
    int columns;
    int rows;
    int size;
};

typedef struct MATRIX_STRUCT *MATRIX;

#define MATRIX_SUCCESS 1
#define MATRIX_ERROR_INVALID_ARGUMENTS -1
#define MATRIX_ERROR_NOT_ENOUGH_MEMORY -2
#define MATRIX_ERROR_NOT_SUPPORTED_DEMENSION -3

#define MATRIX_RANDOM_INITIALIZE() srand(time(NULL))
#define MATRIX_RANDOM(min, max) ((MATRIX_TYPE)(min + (max - min) * (float)rand() / RAND_MAX))

#define MATRIX_FIELD(matrix, column, row) (matrix)->items[(matrix)->columns * row + column]

int createMatrix(MATRIX *matrix, int columns, int rows);
int initMatrix(MATRIX matrix, int minValue, int maxValue);
int fillMatrix(MATRIX matrix, MATRIX_TYPE value);
int freeMatrix(MATRIX *matrix);

#define MATRIX_PRINT_FORMAT "%4.2f  "

int printMatrix(MATRIX matrix, char *title);

int multiplyMatrix(MATRIX matrix1, MATRIX matrix2, MATRIX *result);

#endif /* MATRIX_H */

Реализация функций с матрицами
В принципе, все что нам нужно описано выше, так что не будем тянуть ризину, сразу перейдем к делу, а именно к Matrix.c:

#include "matrix.h"
#include 
#include 

int createMatrix(MATRIX *matrix, int columns, int rows) {
    if (columns <= 0 || rows <= 0 || !matrix) {
        return MATRIX_ERROR_INVALID_ARGUMENTS;
    }
    
    (*matrix) = (MATRIX)malloc(sizeof(struct MATRIX_STRUCT)); 
    if (!(*matrix)) {
        return MATRIX_ERROR_NOT_ENOUGH_MEMORY;
    }
    
    (*matrix)->items = (MATRIX_TYPE*)malloc(columns * rows * sizeof(MATRIX_TYPE));
    if (!(*matrix)->items) {
        free((*matrix));
        return MATRIX_ERROR_NOT_ENOUGH_MEMORY;
    }
    
    (*matrix)->columns = columns;
    (*matrix)->rows = rows;
    (*matrix)->size = columns * rows;
    
    return MATRIX_SUCCESS;
}

int initMatrix(MATRIX matrix, int minValue, int maxValue) {
    int i, j;
    
    if (!matrix || minValue < 0 || maxValue < minValue) {
        return MATRIX_ERROR_INVALID_ARGUMENTS;
    }
    
    for (i = 0; i < matrix->rows; i++) {
        for (j = 0; j < matrix->columns; j++) {
            MATRIX_FIELD(matrix, j, i) = MATRIX_RANDOM(minValue, maxValue);
        }
    }
    
    return MATRIX_SUCCESS;
}

int freeMatrix(MATRIX *matrix) {
    if (!matrix || !(*matrix)) {
        return MATRIX_ERROR_INVALID_ARGUMENTS;
    }
    
    free((*matrix)->items);
    free((*matrix));
    (*matrix) = NULL;
    
    return MATRIX_SUCCESS;
}

int fillMatrix(MATRIX matrix, MATRIX_TYPE value) {
    if (!matrix) {
        return MATRIX_ERROR_INVALID_ARGUMENTS;
    }
    
    MATRIX_TYPE *p = matrix->items;
    int i;
    
    for (i = 0; i < matrix->size; i++, p++) {
        (*p) = value;
    }
    
    return MATRIX_SUCCESS;
}

int printMatrix(MATRIX matrix, char *title) {
    if (!matrix) {
        return MATRIX_ERROR_INVALID_ARGUMENTS;
    }
    
    int i, j;
    
    if (title) {
        printf("%s:\n", title);
    }
    
    for (i = 0; i < matrix->rows; i++) {
        printf("\t");
        for (j = 0; j < matrix->columns; j++) {
            printf(MATRIX_PRINT_FORMAT, MATRIX_FIELD(matrix, j, i));
        }
        printf("\n");
    }
    printf("\n");
    
    return MATRIX_SUCCESS;
}

int multiplyMatrix(MATRIX matrix1, MATRIX matrix2, MATRIX *result) {
    if (!matrix1 || !matrix2 || !result || matrix1->columns != matrix2->rows) {
        return MATRIX_ERROR_INVALID_ARGUMENTS;
    }
    
    int ret, i, j, k;
    MATRIX_TYPE n;
    
    ret = createMatrix(result, matrix2->columns, matrix1->rows);
    if (ret != MATRIX_SUCCESS) {
        return ret;
    }
    
    for (i = 0; i < matrix1->rows; i++) {
        for (j = 0; j < matrix2->columns; j++) {
            n = 0;
            for (k = 0; k < matrix1->columns; k++) {
                n += MATRIX_FIELD(matrix1, k, i) * MATRIX_FIELD(matrix2, j, k);
            }
            MATRIX_FIELD(*result, j, i) = n;
        }
    }
    
    return MATRIX_SUCCESS;
}

К делу
Если вы все еще здесь, я рад, спасибо. И так, не стану описывать все сразу, пойдем по ходу дела, и описание будет в комментариях. И так, main.c:

#include 
#include 
#include 
#include "matrix.h"

#define RANDOM_MIN 0
#define RANDOM_MAX 9.99

#define MATRIX_MPI_TYPE MPI_DOUBLE

int main(int argc, char *argv[]) {
    ...
}

int multiplyMatrixVector(MATRIX matrix, MATRIX vector, MATRIX resultMatrix, int column) {
    // понадобиться не мало переменных
    int i, j, matrixSize, result, blockSize, blockWidth, blockHeight, reduceRank, rank, numberOfProcesses;
    MATRIX blockA, blockB, blockC, resultVector;
    MPI_Datatype matrixBlockType, vectorBlockType;
    MPI_Request requestBlocks[2];
    MPI_Status statusRequestBlocks[2];
    MPI_Comm commReduce, commGather;
    
    // текущий номер процесса в MPI
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    // количество работающих над задачей процессов
    MPI_Comm_size(MPI_COMM_WORLD, &numberOfProcesses);
    
    // ведущим я зделал нулевой (он же первый) процесс, он знает все :)
    if (rank == 0) {
        matrixSize = matrix->size;
    }
    
    // сообщим всем процессам о размере исходной матрицы
    MPI_Bcast(&matrixSize, 1, MPI_INT, 0, MPI_COMM_WORLD);
    
    // разобъем входящую матрицу на блоки, размер, ширина и высота
    blockSize = matrixSize / numberOfProcesses;
    blockWidth = sqrt(blockSize);
    blockHeight = blockSize / blockWidth;
    
    // временный локальный блок для получения кусочка от входящей матрицы
    // само собой, только ведущий процесс (нулевой) знает всю входящую матрицу,
    // остальные процессы получают на вход данной функции matrix == NULL
    result = createMatrix(&blockA, blockWidth, blockHeight);
    if (result != MATRIX_SUCCESS) {
        return result;
    }
    
    // временный локальный блок для получения кусочка от входящего вектора
    // на счет видимости, аналогично matrix
    result = createMatrix(&blockB, 1, blockWidth);
    if (result != MATRIX_SUCCESS) {
        freeMatrix(&blockA);
        return result;
    }
    
    // если ведущий процесс, создадим результирующий вектор, т.к. 
    // результатом умножением матрицы на вектор, есть вектор
    if (rank == 0) {
        result = createMatrix(&resultVector, 1, vector->rows);
        if (result != MATRIX_SUCCESS) {
            freeMatrix(&blockA);
            freeMatrix(&blockB);
            return result;
        }
    }
    
    // запрашиваем получение блоков, асинхронно
    MPI_Irecv(&blockA->items[0], blockA->size, MATRIX_MPI_TYPE, 0, 0, MPI_COMM_WORLD, &requestBlocks[0]);
    MPI_Irecv(&blockB->items[0], blockB->size, MATRIX_MPI_TYPE, 0, 0, MPI_COMM_WORLD, &requestBlocks[1]);
    
    // если ведущий процесс начнем раздачу блоков
    if (rank == 0) {
        // типы данных MPI для формирования и раздачи строк блоков для матрицы
        MPI_Type_vector(blockHeight, blockWidth, matrix->columns, MATRIX_MPI_TYPE, &matrixBlockType);
        MPI_Type_commit(&matrixBlockType);
        
        // ... для вектора
        MPI_Type_vector(blockWidth, 1, vector->columns, MATRIX_MPI_TYPE, &vectorBlockType);
        MPI_Type_commit(&vectorBlockType);
        
        // бегаем по матрице и вектору, и раздаем блоки
        // даже не спрашивайте как подобраны индексы, долго считал и подбирал
        for (i = 0; i < matrix->columns / blockWidth; i++) {
            for (j = 0; j < matrix->rows / blockHeight; j++) {
                MPI_Send(&matrix->items[i * blockWidth + j * blockHeight * matrix->columns], 1, matrixBlockType, i + j * blockHeight, 0, MPI_COMM_WORLD);
                MPI_Send(&vector->items[column + i * blockWidth * vector->columns], 1, vectorBlockType, i + j * blockHeight, 0, MPI_COMM_WORLD);
            }
        }
    }
    
    // все процессы ждут своих блоков
    MPI_Waitall(2, requestBlocks, statusRequestBlocks);
    
    // желаете посмотреть как блоки делятся? прошу
    // printf("%d: ", rank); printMatrix(blockA, "Ma");
    // printf("%d: ", rank); printMatrix(blockB, "Mb");
    
    // наша функция из Matrix.h/.c для перемножения матриц, результат есть вектор blockC
    result = multiplyMatrix(blockA, blockB, &blockC);
    
    // более не нуждаемся в блоках
    freeMatrix(&blockA);
    freeMatrix(&blockB);
    
    // думаю уже заметили, везде проверяем результат, в int main()... если результат плох, то
    // MPI_Abort(MPI_COMM_WORLD, result);
    if (result != MATRIX_SUCCESS) {
        return result;
    }
    
    // посмотреть локальный результат
    // printf("%d: ", rank); printMatrix(blockC, "Mc");
    
    // надо разбить на строки блоков всю рабочую область MPI
    MPI_Comm_split(MPI_COMM_WORLD, rank / blockWidth, 0, &commReduce);
    MPI_Comm_rank(commReduce, &reduceRank); // и получить номер относительно разбитой области
    
    // если ведущий в разбитой области, создадим временный блок
    // для сложения всех блоков в соответствующей строке разбитой области
    // так получим строки (кол-во равное строкам блоков) результата
    if (reduceRank == 0) {
        result = createMatrix(&blockA, 1, blockWidth);
        if (result != MATRIX_SUCCESS) {
            freeMatrix(&blockC);
            return result;
        }
    }
    
    // MPI_SUM решает
    MPI_Reduce(&blockC->items[0], reduceRank == 0 ? &blockA->items[0] : NULL, blockWidth, MATRIX_MPI_TYPE, MPI_SUM, 0, commReduce);
    freeMatrix(&blockC);
    
    // все как обычно, хотим смотрим, хотим не смотрим
    // if (reduceRank == 0) {
    //     printf("%d: ", rank); printMatrix(blockA, "Ma");
    // }
    
    // разбиваем рабочую область на столбцы, уже известных ведущих по строкам
    // сложно сказал ведь :)
    MPI_Comm_split(MPI_COMM_WORLD, reduceRank, 0, &commGather);
    
    // если ведущий по строкам, то раздадим ведущему
    if (reduceRank == 0) {
        MPI_Gather(&blockA->items[0], blockWidth, MATRIX_MPI_TYPE, rank == 0 ? &resultVector->items[0] : NULL, blockWidth, MATRIX_MPI_TYPE, 0, commGather);
        freeMatrix(&blockA);
    }
    
    // заполним нужный столбец в результирующей матрице используя наш локальный вектор-результат
    if (rank == 0) {
        // printMatrix(resultVector, "Me");
        
        for (i = 0; i < resultMatrix->rows; i++) {
            MATRIX_FIELD(resultMatrix, column, i) = MATRIX_FIELD(resultVector, 0, i);
        }
        
        freeMatrix(&resultVector);
    }
    
    return MATRIX_SUCCESS;
}

Ну что же, удачи в ваших начинаниях, скоро напишу о том, как разбитые блоки собрать со всех процессов обратно в единую матрицу, в соответствии с номерами процессов и блоков.

Спасибо и до встречи.