omp/matrix.cpp - master - NenNil/vmk2025-parallel-computing

omp/matrix.cpp
274 строки10 KB
Нил Ненахов
axpy fix
08 ноя 2025, 13:09
08 ноя 2025, 13:0912253d4
#include "matrix.h"

#include <vector>
#include <unordered_map>
#include <algorithm> // Для std::sort
#include <iostream>
#include <cmath>
#include <unordered_map>
#include <iomanip> // Для setw()
#include <chrono>

using namespace std;

// Создание новой матрицы
CSRMatrix* new_matrix(int x, int y) {
    auto matrix = new CSRMatrix();
    matrix->rows = x*y;
    matrix->cols = x*y;

    int amount = 5*y*x - 2*y - 2*x;

    matrix->nns = amount;
    matrix->A.resize(amount, 0.0);
    matrix->JA.resize(amount, 0);
    matrix->IA.resize(matrix->rows + 1, 0);

    return matrix;
}

// Быстрое добавление значения в матрицу
void fast_add(CSRMatrix* matrix, double value, int colIndex, int* num)
{
    matrix->A[*num] = value;
    matrix->JA[*num] = colIndex;
    (*num)++;
}

// Заполнение вектора значением
void vfill(int n, std::vector<double>& vector, double value) {

    std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
    double duration;   
    start = std::chrono::high_resolution_clock::now();

    #pragma omp parallel for
    for (int i = 0; i < n; ++i) {
        vector[i] = value;
    }
    end = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
    duration = duration/1000000;
    std::cout << "5:" << duration << ";";
    
}

// Функция печати матрицы в формате CSR
void print_csr_matrix(const CSRMatrix* matrix) {
    if (!matrix || matrix->rows <= 0 || matrix->cols <= 0) return;
    
    cout << "Rows: " << matrix->rows << ", Columns: " << matrix->cols << endl;
    cout << "Non-zero elements count: " << matrix->nns << endl;

    // Печать массивов A, JA и IA
    cout << "Values (A):\n";
    for (double val : matrix->A) {
        cout << fixed << setprecision(2) << val << ' ';
    }
    cout << "\nColumn indices (JA):\n";
    for (int idx : matrix->JA) {
        cout << idx << ' ';
    }
    cout << "\nRow pointers (IA):\n";
    for (int ptr : matrix->IA) {
        cout << ptr << ' ';
    }
    cout << endl;

    // Полноценный вывод самой матрицы
    cout << "\nFull Matrix Representation:\n";
    for (int row = 0; row < matrix->rows; ++row) {
        cout << "Row " << row << ": ";
        
        // Получаем диапазон индексов для текущей строки
        int start = matrix->IA[row];       // Начало диапазона для текущей строки
        int end = (row + 1 < matrix->rows ? matrix->IA[row + 1] : matrix->nns); // Конец диапазона (учитываем последнюю строку)
        
        // Проходим по каждому столбцу строки
        for (int col = 0; col < matrix->cols; ++col) {
            bool found = false;
            
            // Проверяем, является ли значение в текущей строке ненулевым
            for (int k = start; k < end && !found; ++k) {
                if (matrix->JA[k] == col) { // Если нашли ненулевое значение в данном столбце
                    cout << fixed << setprecision(2) << matrix->A[k] << " "; // Вывести значение
                    found = true;
                }
            }
            
            if (!found) {
                cout << "0.00 "; // Иначе вывести ноль
            }
        }
        cout << endl;
    }
}

// Умножение разреженной матрицы на вектор
void SpMV(const CSRMatrix* A, const std::vector<double>& v, std::vector<double>& result) {

    std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
    double duration;   
    start = std::chrono::high_resolution_clock::now();

    int n = A->IA.size() - 1;

    #pragma omp parallel for
    // #pragma omp parallel for schedule(static)
    for (int i = 0; i < n; ++i) { // строки матрицы параллельно обрабатываются
        double sum = 0.0;
        int line_start = A->IA[i];     // начало строки
        int line_end = A->IA[i + 1];   // конец строки
        
        for (int j = line_start; j < line_end; ++j) {
            sum += A->A[j] * v[A->JA[j]]; // накопляем сумму
        }
    
        result[i] = sum; // сохраняем результат вычислений для каждой строки
    }

    end = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
    duration = duration/1000000;
    std::cout << "6:" << duration << ";";
    
}

// Получение элемента матрицы
double get_elem(const CSRMatrix* matrix, int row, int col) {

    int row_start = matrix->IA[row];
    int row_end = matrix->IA[row + 1];

    for (int i = row_start; i < row_end; ++i) {
        if (matrix->JA[i] == col) {
            return matrix->A[i];
        }
    }
    return 0.0;
}

double dot(int n, const std::vector<double>& a, const std::vector<double>& b) {

    std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
    double duration;   
    start = std::chrono::high_resolution_clock::now();

    // Объявляем локальные переменные вне области видимости OMP-директив
    double result = 0.0;
    
    // Используем pragma omp parallel for для распараллеливания цикла,
    // а также используем специальную секцию reduction для суммирования
    #pragma omp parallel for reduction(+ : result)
    for (int i = 0; i < n; ++i) {
        result += a[i] * b[i];   // Параллельно суммируем элементы произведений
    }

    end = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
    duration = duration/1000000;
    std::cout << "7:" << duration << ";";
    
    return result;
}

// Линейная комбинация векторов: result = x + alpha * y
void axpy(int n, const std::vector<double>& x, const std::vector<double>& y, double alpha, std::vector<double>& result) {

    std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
    double duration;   
    start = std::chrono::high_resolution_clock::now();


    #pragma omp parallel for 
    for (int i = 0; i < n; ++i) {
        result[i] = x[i] + alpha * y[i];
    }

    end = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
    duration = duration/1000000;
    std::cout << "8:" << duration << ";";
    
}


// Копирование вектора
void vcopy(int n, const std::vector<double>& src, std::vector<double>& dst) {

    std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
    double duration;   
    start = std::chrono::high_resolution_clock::now();

    #pragma omp parallel for
    for (int i = 0; i < n; ++i) {
        dst[i] = src[i];
    }

    end = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
    duration = duration/1000000;
    std::cout << "9:" << duration << ";";
    
}

  void copy_diag(CSRMatrix* M1, const CSRMatrix* M2)
{
    int n = M2->rows;
    M1->rows = n;
    M1->cols = n;
    M1->nns  = n;
    M1->A.assign(n, 0.0);
    M1->JA.resize(n);
    M1->IA.resize(n+1);
    #pragma omp for
    for (int i = 0; i < n; ++i) {
        M1->JA[i] = i;
        M1->IA[i] = i;
        double di = get_elem(M2, i, i);
        if (di == 0.0) throw std::runtime_error("Нулевой диагональный элемент для предобуславливателя");
        M1->A[i] = 1.0 / di;
    }
    M1->IA[n] = n;
}

// Функция, создающая подматрицу CSR из исходной матрицы
CSRMatrix submatrix(const CSRMatrix& mat, int m1, int m2, int n1, int n2) {

    // Формируем новую матрицу
    CSRMatrix result = {n2-n1+1, m2-m1+1, 0, {}, {}, {}};

    if (n1 > n2 || m1 > m2 || n1 >= mat.cols || n2 >= mat.cols || m1 >= mat.rows || m2 >= mat.rows)
        return result;

    for(int i = n1; i <= n2; ++i) {
        // Получаем начало и конец диапазона текущего ряда
        int start = mat.IA[i];
        int end = mat.IA[i + 1];  // Следующая строка начинается там же, где заканчивается текущая
        bool first_item = true;
        
        // Проходим по элементам строки
        for(int j = start; j < end; ++j) {
            // Проверяем попадание столбца в нужный диапазон
            if(mat.JA[j] >= m1 && mat.JA[j] <= m2) {
                // Добавляем элемент в новую матрицу
                result.A.push_back(mat.A[j]);
                result.JA.push_back(mat.JA[j]-m1);  // Относительный номер столбца
                result.nns += 1;
                
                // Заполняем вектор смещений
                if(first_item) {
                    result.IA.push_back(result.A.size()-1);
                    first_item = false;
                }
            }
        }
        if (first_item) result.IA.push_back(result.A.size());
    }

    // Закрываем последнюю строку
    result.IA.push_back(result.A.size());
    result.nns = result.A.size();

    return result;
}