[C++ 성능 최적화 및 고급 테크닉] Day 26: 프로젝트: 고성능 매트릭스 라이브러리 개발 (3)

프로젝트 목표

이 단계에서는 매트릭스 라이브러리의 성능을 더욱 향상시키기 위해 추가적인 최적화 기법을 적용하고, 병렬 프로그래밍을 도입합니다. 특히, 다음과 같은 부분을 다룹니다:

1. 병렬 프로그래밍 도입: OpenMP를 사용하여 병렬로 행렬 연산을 수행합니다.
2. 매트릭스 라이브러리 완성: 라이브러리의 다양한 기능을 통합하고 최적화합니다.

 

Step 1: 병렬 프로그래밍 도입

OpenMP를 사용하여 행렬 덧셈과 행렬 곱셈을 병렬로 수행하도록 수정합니다.

 

Matrix.cpp: 병렬 행렬 덧셈

#include "Matrix.h"
#include <omp.h>

Matrix operator+(const Matrix& lhs, const Matrix& rhs) {
    if (lhs.getRows() != rhs.getRows() || lhs.getCols() != rhs.getCols()) {
        throw std::invalid_argument("Matrices must have the same dimensions for addition");
    }

    Matrix result(lhs.getRows(), lhs.getCols());
    size_t rows = lhs.getRows();
    size_t cols = lhs.getCols();

    #pragma omp parallel for
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < cols; ++j) {
            result(i, j) = lhs(i, j) + rhs(i, j);
        }
    }

    return result;
}

 

Matrix.cpp: 병렬 행렬 곱셈

#include "Matrix.h"
#include <omp.h>

Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
    if (lhs.getCols() != rhs.getRows()) {
        throw std::invalid_argument("Number of columns of the first matrix must be equal to the number of rows of the second matrix");
    }

    Matrix result(lhs.getRows(), rhs.getCols());
    size_t lhsCols = lhs.getCols();
    size_t rhsCols = rhs.getCols();
    size_t rows = lhs.getRows();

    #pragma omp parallel for
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < rhsCols; ++j) {
            double sum = 0.0;
            for (size_t k = 0; k < lhsCols; k += 4) {
                sum += lhs(i, k) * rhs(k, j)
                     + lhs(i, k + 1) * rhs(k + 1, j)
                     + lhs(i, k + 2) * rhs(k + 2, j)
                     + lhs(i, k + 3) * rhs(k + 3, j);
            }
            result(i, j) = sum;
        }
    }

    return result;
}

 

Step 2: 라이브러리 완성

최종적으로 매트릭스 라이브러리의 다양한 기능을 통합하고, 성능을 최적화합니다.

 

Matrix.h

#ifndef MATRIX_H
#define MATRIX_H

#include <vector>
#include <iostream>
#include <stdexcept>

class Matrix {
private:
    std::vector<std::vector<double>> data;
    size_t rows, cols;

public:
    Matrix(size_t rows, size_t cols);
    double& operator()(size_t row, size_t col);
    const double& operator()(size_t row, size_t col) const;
    size_t getRows() const;
    size_t getCols() const;
    void print() const;

    friend Matrix operator+(const Matrix& lhs, const Matrix& rhs);
    friend Matrix operator*(const Matrix& lhs, const Matrix& rhs);
    friend Matrix transpose(const Matrix& matrix);
    friend Matrix inverse(const Matrix& matrix);
    friend double determinant(const Matrix& matrix);
};

#endif // MATRIX_H

 

Matrix.cpp

#include "Matrix.h"
#include <omp.h>

// 생성자
Matrix::Matrix(size_t rows, size_t cols) : rows(rows), cols(cols), data(rows, std::vector<double>(cols, 0.0)) {}

// 연산자 오버로딩
double& Matrix::operator()(size_t row, size_t col) {
    if (row >= rows || col >= cols) {
        throw std::out_of_range("Index out of range");
    }
    return data[row][col];
}

const double& Matrix::operator()(size_t row, size_t col) const {
    if (row >= rows || col >= cols) {
        throw std::out_of_range("Index out of range");
    }
    return data[row][col];
}

size_t Matrix::getRows() const {
    return rows;
}

size_t Matrix::getCols() const {
    return cols;
}

void Matrix::print() const {
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < cols; ++j) {
            std::cout << data[i][j] << " ";
        }
        std::cout << std::endl;
    }
}

Matrix operator+(const Matrix& lhs, const Matrix& rhs) {
    if (lhs.getRows() != rhs.getRows() || lhs.getCols() != rhs.getCols()) {
        throw std::invalid_argument("Matrices must have the same dimensions for addition");
    }

    Matrix result(lhs.getRows(), lhs.getCols());
    size_t rows = lhs.getRows();
    size_t cols = lhs.getCols();

    #pragma omp parallel for
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < cols; ++j) {
            result(i, j) = lhs(i, j) + rhs(i, j);
        }
    }

    return result;
}

Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
    if (lhs.getCols() != rhs.getRows()) {
        throw std::invalid_argument("Number of columns of the first matrix must be equal to the number of rows of the second matrix");
    }

    Matrix result(lhs.getRows(), rhs.getCols());
    size_t lhsCols = lhs.getCols();
    size_t rhsCols = rhs.getCols();
    size_t rows = lhs.getRows();

    #pragma omp parallel for
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < rhsCols; ++j) {
            double sum = 0.0;
            for (size_t k = 0; k < lhsCols; k += 4) {
                sum += lhs(i, k) * rhs(k, j)
                     + lhs(i, k + 1) * rhs(k + 1, j)
                     + lhs(i, k + 2) * rhs(k + 2, j)
                     + lhs(i, k + 3) * rhs(k + 3, j);
            }
            result(i, j) = sum;
        }
    }

    return result;
}

Matrix transpose(const Matrix& matrix) {
    Matrix result(matrix.getCols(), matrix.getRows());
    size_t blockSize = 16; // 캐시 친화적 접근을 위한 블록 크기 설정

    for (size_t i = 0; i < matrix.getRows(); i += blockSize) {
        for (size_t j = 0; j < matrix.getCols(); j += blockSize) {
            for (size_t k = i; k < i + blockSize && k < matrix.getRows(); ++k) {
                for (size_t l = j; l < j + blockSize && l < matrix.getCols(); ++l) {
                    result(l, k) = matrix(k, l);
                }
            }
        }
    }

    return result;
}

Matrix inverse(const Matrix& matrix) {
    if (matrix.getRows() != matrix.getCols()) {
        throw std::invalid_argument("Matrix must be square to find its inverse");
    }

    size_t n = matrix.getRows();
    Matrix augmented(matrix.getRows(), matrix.getCols() * 2);

    // Augment the matrix with the identity matrix
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            augmented(i, j) = matrix(i, j);
            augmented(i, j + n) = (i == j) ? 1.0 : 0.0;
        }
    }

    // Perform Gaussian elimination
    for (size_t i = 0; i < n; ++i) {
        if (augmented(i, i) == 0) {
            throw std::runtime_error("Matrix is singular and cannot be inverted");
        }

        for (size_t j = 0; j < 2 * n; ++j) {
            augmented(i, j) /= augmented(i, i);
        }

        for (size_t k = 0; k < n; ++k) {
            if (k != i) {
                double factor = augmented(k, i);
                for (size_t j = 0; j < 2 * n; ++

j) {
                    augmented(k, j) -= factor * augmented(i, j);
                }
            }
        }
    }

    // Extract the inverse matrix
    Matrix result(n, n);
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            result(i, j) = augmented(i, j + n);
        }
    }

    return result;
}

double determinant(const Matrix& matrix) {
    if (matrix.getRows() != matrix.getCols()) {
        throw std::invalid_argument("Matrix must be square to calculate its determinant");
    }

    size_t n = matrix.getRows();

    if (n == 1) {
        return matrix(0, 0);
    }

    if (n == 2) {
        return matrix(0, 0) * matrix(1, 1) - matrix(0, 1) * matrix(1, 0);
    }

    double det = 0.0;
    for (size_t i = 0; i < n; ++i) {
        Matrix subMatrix(n - 1, n - 1);
        for (size_t j = 1; j < n; ++j) {
            size_t colIndex = 0;
            for (size_t k = 0; k < n; ++k) {
                if (k == i) continue;
                subMatrix(j - 1, colIndex++) = matrix(j, k);
            }
        }
        det += ((i % 2 == 0) ? 1 : -1) * matrix(0, i) * determinant(subMatrix);
    }

    return det;
}

 

테스트 코드

main.cpp

#include <iostream>
#include "Matrix.h"

int main() {
    Matrix mat1(3, 3);

    mat1(0, 0) = 1; mat1(0, 1) = 2; mat1(0, 2) = 3;
    mat1(1, 0) = 0; mat1(1, 1) = 1; mat1(1, 2) = 4;
    mat1(2, 0) = 5; mat1(2, 1) = 6; mat1(2, 2) = 0;

    std::cout << "Matrix 1:" << std::endl;
    mat1.print();

    Matrix inv = inverse(mat1);
    std::cout << "Inverse of Matrix 1:" << std::endl;
    inv.print();

    double det = determinant(mat1);
    std::cout << "Determinant of Matrix 1: " << det << std::endl;

    Matrix transposed = transpose(mat1);
    std::cout << "Transposed Matrix 1:" << std::endl;
    transposed.print();

    Matrix mat2(3, 3);
    mat2(0, 0) = 7; mat2(0, 1) = 8; mat2(0, 2) = 9;
    mat2(1, 0) = 4; mat2(1, 1) = 5; mat2(1, 2) = 6;
    mat2(2, 0) = 1; mat2(2, 1) = 2; mat2(2, 2) = 3;

    std::cout << "Matrix 2:" << std::endl;
    mat2.print();

    Matrix sum = mat1 + mat2;
    std::cout << "Sum of Matrix 1 and Matrix 2:" << std::endl;
    sum.print();

    Matrix product = mat1 * mat2;
    std::cout << "Product of Matrix 1 and Matrix 2:" << std::endl;
    product.print();

    return 0;
}

 

이제 스물여섯 번째 날의 학습을 마쳤습니다. 고성능 매트릭스 라이브러리 개발 프로젝트의 세 번째 단계를 통해 병렬 프로그래밍을 도입하고, 라이브러리의 다양한 기능을 통합하여 최적화했습니다.

질문이나 피드백이 있으면 언제든지 댓글로 남겨 주세요. 내일은 "프로젝트: 실시간 데이터 처리 시스템 개발 (1)"에 대해 학습하겠습니다.