Оптимизация произведения двух разреженных матриц в CRS формате

// Разреженное матричное умножение двух разреженных матриц.
// Алгоритм Ф. Густавсона IBM 1978.
template <typename doublerealT>
void my_parallel8_sparse_matrix_by_matrix_multiplication_RA(Ak2& Amat,
    Ak1*& P, integer& istartAnew, integer*& istartAnew_m,
    integer*& row_ind_SR,
    integer*& row_ind_ER,
    integer*& row_startA,
    integer numberofcoarcenodes,
    const int iKnumber_thread,
    bool**& hash_table_m,
    integer**& index_visit_m,
    doublerealT**& vector_sum_m,
    integer*& index_size_m,
    integer& nsizeA,
    integer n, integer nnz, integer ilevel,
    bool bprint_mesage_diagnostic,
    integer*& n_a, doublereal AccumulqtorA_m_SIZE8,
    Ak1**& AccumulqtorA_m
    //, integer istartAnew_mem
) {
 
#ifdef _OPENMP
    int i_my_num_core_parallelesation = omp_get_max_threads();
    
 
    integer inum_core = number_cores();
    omp_set_num_threads(inum_core);
    
 
    const integer iSIZE =  n + 2;
 
    // Данные используемые для частичного формирователя суммы.
 
    for (int i_9 = 0; i_9 < iKnumber_thread; ++i_9) {
 
#pragma omp parallel for schedule (static)
        for (integer i_91 = 0; i_91 < iSIZE; ++i_91) {
            hash_table_m[i_9][i_91] = false;// inicialization
        }
        index_size_m[i_9] = 0;
        istartAnew_m[i_9] = 0;
    }
 
    const integer LIMIT_ANEW_m = (integer)((1.0/(1.0* iKnumber_thread)) * AccumulqtorA_m_SIZE8 * nnz + 1);
 
    //omp_set_num_threads(iKnumber_thread); // только 8 потоков.
 
    // Сканируем первый операнд построчно.
    // глобальные переменные не перечисляются.
#pragma omp parallel for schedule (static) 
    for (integer istr = 1; istr <= numberofcoarcenodes; ++istr) {
 
#ifdef _OPENMP 
        const int tid = omp_get_thread_num();
#else
        const int tid = 0;
#endif
 
        // на основе хеш-таблицы. 
 
        // Сканируем текущую i-ую строку поэлементно
        for (integer ii = row_ind_SR[istr]; ii <= row_ind_ER[istr]; ++ii) {
            const integer col_ind = P[ii].j;
            // Сканируем col_ind строку второго операнда
 
            // Общую переменную объявим на уровень выше.
            const doublerealT left_operand = P[ii].aij;
 
            const integer i_1start = row_startA[col_ind];
            const integer i_1end = row_startA[col_ind + 1] - 1;
            for (integer i_1 = i_1start; i_1 <= i_1end; ++i_1) {
 
                const doublerealT right_operand = Amat.aij[i_1];
                const integer iaddind = Amat.j[i_1];
 
                if (hash_table_m[tid][iaddind]) {
 
                    vector_sum_m[tid][iaddind] += left_operand * right_operand;
                }
                else {
                    // Первое добавление.
 
                    index_size_m[tid]++;
 
                    index_visit_m[tid][index_size_m[tid]] = iaddind;
 
                    hash_table_m[tid][iaddind] = true;
 
                    vector_sum_m[tid][iaddind] = left_operand * right_operand;
 
                }
            }
        }
        
        
 
        const integer ism = index_size_m[tid];
        if (nsizeA > istartAnew + ism) {
            if (istartAnew_m[tid] + ism < LIMIT_ANEW_m) {
                for (integer i_6 = 1; i_6 <= ism; ++i_6) {
                    const integer jstr = index_visit_m[tid][i_6];
 
                    hash_table_m[tid][jstr] = false; // инициализируем хеш-таблицу для следующих проходов.
 
                    const doublerealT vs1 = vector_sum_m[tid][jstr];
 
                    // 7 ноября 2016 игнорируем чистые нули:
                    if (fabs(vs1) > 1.0e-37) {
                        // Мы не записываем в матрицу чистый ноль.
                        Ak1 Atemp;
                        Atemp.aij = vs1;
                        Atemp.i = istr;
                        Atemp.j = jstr;
 
 
                        AccumulqtorA_m[tid][istartAnew_m[tid]++] = Atemp;
                    }
                }
            }
            else {
                printf("AccumulqtorA_m index overflow\n");
                system("pause");
            }
        }
        else {
            // Слишком мало памяти выделено под матрицу А и в неё не умещаются все данные.
            // Нужно увеличить объём выделяемой памяти для А и перезапустить приложение.
            printf("Amat lack of memory\n");
            printf("yuo mast increase the size of the matrix Amat and restart solver\n");
            printf("please, press any key to exit.\n");
            system("pause");
            exit(1);
        }
 
 
 
        index_size_m[tid] = 0; // Сброс индикатора, строка обработана.           
 
    }
 
 
    omp_set_num_threads(i_my_num_core_parallelesation);
 
#endif
 
}

auto& Atemp=AccumulqtorA_m[tid][istartAnew_m[tid]++];
Atemp.aij = vs1;
Atemp.i = istr;
Atemp.j = jstr;

/// C++17
///----------------------------------------------------------------------------|
/// Стенд
/// для анализа возможностей оптимизации умножения двух разряжённых матриц.
/// version(2021.08.29)
///----------------------------------------------------------------------------:
#include <iostream>
#include <vector>
#include <string>
#include <iomanip>
 
typedef double T;
typedef std::vector<T>     vec_t;
typedef std::vector<vec_t> mat_t;
 
///----------------------------------------|
/// Конфиг.                                |
///----------------------------------------:
///--------------------|
/// Размерность.       |
///--------------------:
const size_t Ca = 14   ;
const size_t Ra = 9    ;
 
const size_t Cb = 12   ;
const size_t Rb = Ca   ; /// НЕ ТРОГАТЬ! assert(Ca == Rb)
 
///--------------------|
/// Как инициализируем.|
///--------------------:
const int    SID       = 2021; /// Набор рандома.
const size_t DISCHARGE = 96  ; /// Разряжение(в процентах): 96%
const size_t PREC      =  2  ; /// Точность вывода.
 
///------------------------|
/// Критерий разряженности.|
///------------------------:
inline bool IS_DISCHARGE(T a){return -0.0000000001 < a && a < 0.0000000001;}
 
#define l(v)      std::cout << std::setw(12) << #v << " = " << (v) << "\n";
#define ASSERT(v) if(!(v)) {std::cout << "ERROR: line = " << __LINE__ << '\n';\
                            std::cin.get();}
 
///--------------------------------------|
/// Пока тут умножаем по классике...     |
///--------------------------------------:
static mat_t mult_classic(const mat_t& a, const mat_t& b)
{   mat_t c(a.size(), vec_t(b[0].size(), 0));
 
    if (a[0].size() != b.size())
    {   std::cout << "ERROR: bad data: Cb != Ra\n";
        return mat_t(0, vec_t(0));
    }
 
    for         (size_t i = 0; i < a   .size(); i++)
    {   for     (size_t j = 0; j < b[0].size(); j++)
        {   for (size_t k = 0; k < b   .size(); k++)
            {
                c[i][j] += a[i][k] * b[k][j];
            }
        }
    }
 
    return c;
}
 
///--------------------------------------|
/// Умножаем по СУПЕР-ПУПЕР методе.     <---: Объект исследования ...
///--------------------------------------:
static mat_t mult_super_rocket(const mat_t& a, const mat_t& b)
{   mat_t c(a.size(), vec_t(b[0].size(), 0));
 
    if (a[0].size() == b.size())
    {   std::cout << "ERROR: bad data: Ca == Rb\n";
        return mat_t(0, vec_t(0));
    }
 
    ///-----------------|
    /// ...             |
    ///-----------------:
    /// TODO.
 
    return c;
}
 
///--------------------------------------|
/// Индикация либо на экран, либо в файл.|
/// Если control_big = true,             |
///     то вывод с ограничениями ...     |
///--------------------------------------:
template<class OUT>
void debug(OUT& out, const mat_t& m, const char* s ="", bool control_big = true)
{   out << "\nmat_t " << s
        << "[ " << m.size() << ", " << m[0].size() << " ] --------------------:"
        << "\n";
    int rr = 0;
    for(    const auto& r : m)
    {       if( control_big && ++rr%12 == 0){ out << "  ..."; break;}
        int cc = 0;
        for(const auto& c : r)
        {   if( control_big && ++cc%12 == 0){ out << "  ..."; break;}
            out << std::setprecision(PREC) << std::setw(7) << c << ' ' ;
        }   out                                                 << '\n';
    }       out                                                 << '\n';
}
 
///--------------------------------------|
/// Сравниваем два контейнера.           |
///--------------------------------------:
template<class C>
bool is_equal(const C& a, const C& b)
{   if(a.size() != b.size() || a[0].size() != b[0].size()) return false;
    for(auto ia = a.begin(),
             ib = b.begin(); ia != a.end() && ib != b.end(); ++ia, ++ib)
    {   if(*ia != *ib) return false;
    }
    return true;
}
 
///--------------------------------------|
/// ГПСЧ.                                |
///--------------------------------------:
#include <ctime>
#include <cstdlib>
struct Rand
{   Rand(            ){   srand(static_cast<unsigned>(time(nullptr)));}
    Rand(unsigned sid){   srand(sid)                                 ;}
    int operator(    )(int range_min, int range_max) const
    {   return rand() % (range_max - range_min) + range_min;
    }
    static int get(int range_min, int range_max)
    {   static Rand r(SID); return r(range_min, range_max);
    }
};
 
///--------------------------------------|
/// Рандомно заполняем матрицу.          |
/// с заданным разряжением DISCHARGE.    |
///--------------------------------------:
static void fillrand(mat_t& m)
{   auto R = m   .size();
    auto C = m[0].size();
 
    size_t AMVAL = R * C * (100 - DISCHARGE) / DISCHARGE;
    for(size_t i = AMVAL; i;)
    {   size_t c = static_cast<size_t>(Rand::get(0, static_cast<int>(C)));
        size_t r = static_cast<size_t>(Rand::get(0, static_cast<int>(R)));
 
        if(0.0 == m[r][c]) {m[r][c] = static_cast<T>(r * c * 1234) / 100; --i;}
    }
}
 
///----------------------------------------------------------------------------|
/// Старт тест.
///----------------------------------------------------------------------------:
int main()
{   ///--------|
    /// Конфиг.|
    ///--------:
    l(Ra*Ca)
    l(Rb*Cb)
    l(DISCHARGE)
 
    ///----------|
    /// На входе.|
    ///----------:
    mat_t a(Ra, vec_t(Ca, 0.)); fillrand(a);
    mat_t b(Rb, vec_t(Cb, 0.)); fillrand(b);
 
    ///-------------------------|
    /// Классика умножения.     |
    ///-------------------------:
    mat_t c = mult_classic(a, b);
    mat_t d = mult_classic(a, b);
 
    ASSERT(is_equal(c, d)) // ок
 
    mat_t e = mult_super_rocket(a, b);
  //ASSERT(is_equal(c, e)) /// тут пока ожидается срабатывание ... (TODO)
 
    debug(std::cout, a, "a");
    debug(std::cout, b, "b");
    debug(std::cout, c, "c");
 
    std::cout << "\nTEST FINISH!\n";
    std::cin.get();
}