mtx_avx2_arr_f64.cpp

Go to the documentation of this file.

/*************************************************************************
 *
 *  Project
 *                         _____ _____  __  __ _____
 *                        / ____|  __ \|  \/  |  __ \
 *  ___  _ __   ___ _ __ | |  __| |__) | \  / | |__) |
 * / _ \| '_ \ / _ \ '_ \| | |_ |  ___/| |\/| |  ___/
 *| (_) | |_) |  __/ | | | |__| | |    | |  | | |
 * \___/| .__/ \___|_| |_|\_____|_|    |_|  |_|_|
 *      | |
 *      |_|
 *
 * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
 *
 * This software is licensed as described in the file LICENSE, which
 * you should have received as part of this distribution. The terms
 * among other details are referenced in the official documentation
 * seen here : https://akielaries.github.io/openGPMP/ along with
 * important files seen in this project.
 *
 * You may opt to use, copy, modify, merge, publish, distribute
 * and/or sell copies of the Software, and permit persons to whom
 * the Software is furnished to do so, under the terms of the
 * LICENSE file. As this is an Open Source effort, all implementations
 * must be of the same methodology.
 *
 *
 *
 * This software is distributed on an AS IS basis, WITHOUT
 * WARRANTY OF ANY KIND, either express or implied.
 *
 ************************************************************************/
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <iostream>
#include <openGPMP/linalg/mtx.hpp>
#include <vector>
 
#if defined(__x86_64__) || defined(__amd64__) || defined(__amd64)
 
/************************************************************************
 *
 * Matrix Operations for AVX ISA
 *
 ************************************************************************/
#if defined(__AVX2__)
 
// AVX family intrinsics
#include <immintrin.h>
 
/************************************************************************
 *
 * Matrix Operations on Arrays
 *
 ************************************************************************/
// matrix addition using Intel intrinsics, accepts double arrays as matrices
void gpmp::linalg::Mtx::mtx_add(const double *A,
                                const double *B,
                                double *C,
                                int rows,
                                int cols) {
    if (rows > 8) {
        for (int i = 0; i < rows; ++i) {
            int j = 0;
            // requires at least size 4x4 size matrices
            for (; j < cols - 3; j += 4) {
                // load 4 elements from A, B, and C matrices using SIMD
                __m256d a = _mm256_loadu_pd(&A[i * cols + j]);
                __m256d b = _mm256_loadu_pd(&B[i * cols + j]);
                __m256d c = _mm256_loadu_pd(&C[i * cols + j]);
                // perform vectorized addition and accumulate the result
                c = _mm256_add_pd(a, b);
 
                // store the result back to the C matrix
                _mm256_storeu_pd(&C[i * cols + j], c);
            }
 
            // handle the remaining elements that are not multiples of 8
            for (; j < cols; ++j) {
                C[i * cols + j] = A[i * cols + j] + B[i * cols + j];
            }
        }
    } else {
        // use standard matrix addition
        std_mtx_add(A, B, C, rows, cols);
    }
}
 
void gpmp::linalg::Mtx::mtx_mult(const double *A,
                                 const double *B,
                                 double *C,
                                 int rows_a,
                                 int cols_a,
                                 int cols_b) {
    if (cols_a != rows_a) {
        // Matrix dimensions don't match for multiplication
        std::cerr << "Matching error";
        return;
    }
 
    if (rows_a > 8) {
 
        for (int i = 0; i < rows_a; ++i) {
            for (int j = 0; j < cols_b - 3; j += 4) {
                // creat result vector of zeros
                __m256d sum_vec = _mm256_setzero_pd();
 
                for (int k = 0; k < cols_a; ++k) {
                    __m256d a_vec = _mm256_set1_pd(A[i * cols_a + k]);
 
                    __m256d b_vec = _mm256_loadu_pd(&B[k * cols_b + j]);
 
                    __m256d prod = _mm256_mul_pd(a_vec, b_vec);
 
                    sum_vec = _mm256_add_pd(sum_vec, prod);
                }
                _mm256_storeu_pd(&C[i * cols_b + j], sum_vec);
            }
 
            // handle remaining elements not multiples of 4
            for (int j = cols_b - cols_b % 4; j < cols_b; ++j) {
                double sum = 0.0;
 
                for (int k = 0; k < cols_a; ++k) {
                    sum += A[i * cols_a + k] * B[k * cols_b + j];
                }
                C[i * cols_b + j] = sum;
            }
        }
 
    }
 
    else {
        std_mtx_mult(A, B, C, rows_a, cols_a, cols_b);
    }
}
 
#endif
 
// x86
#endif