Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <vector>
- #include <algorithm>
- #include "omp.h"
- double RandomDouble(double min = -1, double max = 1)
- {
- double f = (double)std::rand() / RAND_MAX;
- return min + f * (max - min);
- }
- std::vector<double> RandomDoubleVector(size_t size, double min = -1, double max = 1)
- {
- std::vector<double> vector = std::vector<double>(size);
- for (int i = 0; i < size; ++i) {
- vector[i] = RandomDouble(min, max);
- }
- return vector;
- }
- #define MATRIX_SIZE 100
- class Matrix {
- public:
- Matrix(size_t size = MATRIX_SIZE) : _size(size), _data(size * size) {}
- Matrix(size_t size, std::vector<double> data) : _size(size), _data(size * size) {
- _data = data;
- }
- Matrix(const Matrix& other)
- {
- _size = other._data.size();
- _data.resize(_size * _size);
- for (size_t i = 0; i < _size; ++i)
- for (size_t j = 0; j < _size; ++j)
- this->operator()(i, j) = other(i, j);
- }
- Matrix& operator=(const Matrix& other) {
- for (size_t i = 0; i < _size; ++i)
- for (size_t j = 0; j < _size; ++j)
- this->operator()(i, j) = other(i, j);
- return *this;
- }
- double& operator()(size_t i, size_t j) {
- return _data[i * _size + j];
- }
- double operator()(size_t i, size_t j) const {
- return _data[i * _size + j];
- }
- static Matrix I(size_t size) {
- Matrix matrix = Matrix(size);
- for (size_t i = 0; i < size; ++i)
- matrix(i, i) = 1;
- return matrix;
- }
- static Matrix Zero(size_t size) {
- Matrix matrix = Matrix(size);
- return matrix;
- }
- static Matrix RandomMatrix(size_t size, double min = -1, double max = 1) {
- Matrix matrix = Matrix(size);
- for (size_t i = 0; i < size; ++i)
- for (size_t j = 0; j < size; ++j)
- matrix(i, j) = RandomDouble(min, max);
- return matrix;
- }
- Matrix operator+(const Matrix& b) {
- Matrix sum = Matrix(_size, _data);
- for (size_t i = 0; i < _size; ++i)
- for (size_t j = 0; j < _size; ++j)
- sum(i, j) += b(i, j);
- return sum;
- }
- Matrix operator+(const double & d) {
- Matrix sum = Matrix(_size, _data);
- for (size_t i = 0; i < _size; ++i)
- sum(i, i) += d;
- return sum;
- }
- Matrix operator+=(const Matrix& b) {
- for (size_t i = 0; i < _size; ++i)
- for (size_t j = 0; j < _size; ++j)
- operator()(i, j) += b(i, j);
- return *this;
- }
- Matrix operator*(const Matrix& b) {
- Matrix product = Matrix(_size);
- for (size_t i = 0; i < _size; ++i)
- for (size_t j = 0; j < _size; ++j)
- for (size_t k = 0; k < _size; ++k) {
- product(i, j) += operator()(i, k) * b(k, j);
- }
- return product;
- }
- Matrix operator*(const double& d) {
- Matrix scalarMult = Matrix(_size, _data);
- for (size_t i = 0; i < _size; ++i)
- for (size_t j = 0; j < _size; ++j)
- scalarMult(i, j) *= d;
- return scalarMult;
- }
- Matrix Pow(size_t n) {
- Matrix result = Matrix(_size);
- for (size_t i = 0; i < _size; i++)
- {
- result(i, i) = 1;
- }
- for (size_t i = 0; i < n; i++) {
- result = result * *this;
- }
- return result;
- }
- void Show() {
- for (size_t i = 0; i < _size; ++i) {
- for (size_t j = 0; j < _size; ++j)
- std::cout << operator()(i, j) << " ";
- std::cout << std::endl;
- }
- std::cout << std::endl;
- }
- const size_t GetSize() {
- return _size;
- }
- double* GetData() {
- return _data.data();
- }
- private:
- size_t _size;
- std::vector<double> _data;
- };
- Matrix ForwardSum(std::vector<double> polinomialCoefs, Matrix& A)
- {
- auto result = Matrix(A.GetSize());
- for (size_t n = 0; n < polinomialCoefs.size(); n++) {
- result += A.Pow(n) * polinomialCoefs[n];
- }
- return result;
- }
- #pragma omp declare reduction( + : Matrix: \
- omp_out = omp_out + omp_in)
- int main(int argc, char* argv[]) {
- size_t matrixSize = MATRIX_SIZE;
- Matrix A = Matrix::RandomMatrix(matrixSize);
- std::vector<double> polinomialCoefs = RandomDoubleVector(20, -5, 5);
- std::cout << omp_get_max_threads();
- int threadsCount;
- clock_t begin = clock();
- Matrix result = Matrix::Zero(matrixSize);
- #pragma omp parallel for reduction(+ : result)
- for (int currentThread = 0; currentThread < omp_get_num_threads(); currentThread++) {
- threadsCount = omp_get_num_threads();
- std::vector<size_t> indexes(0);
- for (int i = currentThread; i < polinomialCoefs.size(); i += threadsCount) {
- indexes.push_back(i);
- }
- Matrix APowered = A.Pow(threadsCount);
- Matrix I = Matrix::I(matrixSize);
- Matrix part = Matrix::I(matrixSize) * polinomialCoefs[indexes[indexes.size() - 1]];
- for (int n = indexes.size() - 2; n >= 0; n--) {
- part = part * APowered + I * polinomialCoefs[indexes[n]];
- }
- part = A.Pow(currentThread) * part;
- result += part;
- }
- #pragma omp barrier
- clock_t end = clock();
- if (true) {
- std::cout << "\nDone in " << double(end - begin) / CLOCKS_PER_SEC << " sec." << std::endl;
- //result.Show();
- //ForwardSum(polinomialCoefs, A).Show();
- }
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement