Я работаю с алгоритмами для матрицы строк сжатого хранилища. Некоторое время я использую Eigen .. Однако я понял, что Eigen не использует истинные основные матрицы строк (по крайней мере, не после makeCompressed или setFromTriplets), поскольку память в строках не выровнена (даже с опцией RowMajor); или я что-то упустил? коэффициенты в каждой строке. Может ли кто-нибудь помочь мне сгенерировать более быстрый код или указать мне статью, в которой объясняется эффективная реализация (я уже пробовал использовать двоичный поиск вместо поиска с переходом и не получил улучшений)?
Моя реализация:
#include <vector>
#include <iostream>
#include <algorithm>
//--------------------------------------------------------------------------------
//method for sorting:
//--------------------------------------------------------------------------------
//mirrored quicksort algorithm
template <typename T>
int mirror_partition(std::vector<unsigned>& arr, const int& low, const int& high, std::vector<T>& mirror)
{
int pivot = arr[high]; // pivot
int i = (low - 1);
for (int j = low; j <= high- 1; j++)
{
//if current element is smaller than pivot, increment the low element
//swap elements at i and j
if (arr[j] <= pivot)
{
i++; // increment index of smaller element
std::swap(arr[i], arr[j]);
std::swap(mirror[i], mirror[j]);
}
}
std::swap(arr[i + 1], arr[high]);
std::swap(mirror[i + 1], mirror[high]);
return (i + 1);
}
template <typename T>
void mirror_quicksort(std::vector<unsigned>& arr, const int& low, const int& high, std::vector<T>& mirror)
{
if (low < high)
{
//partition the array
int pivot = mirror_partition(arr, low, high, mirror);
//sort the sub arrays independently
mirror_quicksort(arr, low, pivot - 1, mirror);
mirror_quicksort(arr, pivot + 1, high, mirror);
}
}
//--------------------------------------------------------------------------------
//method for iterating:
//--------------------------------------------------------------------------------
// jump search method
inline int min_int(const int& i, const int& j){
if (i > j) return j;
else return i;
}
inline int max_int(const int& i, const int& j){
if (i > j) return i;
else return j;
}
template<typename T>
inline int jump_search(std::vector<T>& vec, const int& start, const int& end,
const int& step, const T& val){
int a = start;
while(vec[min_int(a,end-1)] < val){
if (a > end) return -1;
a += step;
}
for(int k = max_int(a-step,start); k < min_int(end,a+1); ++k){
if(vec[k] == val) return k;
}
return -1;
}
//--------------------------------------------------------------------------------
//CSR matrix implementation:
//--------------------------------------------------------------------------------
template<typename T>
class CSR_matrix{
public:
std::vector<T> data; //data vec
std::vector<unsigned> p_rows, cols_id, row_nnz; //row pointers, column IDs, number of nnz per row
unsigned nnz, rows, cols, step;
bool use_jump;
public:
//constructor
explicit CSR_matrix(const unsigned& rows, const unsigned& cols) : rows(rows), cols(cols), nnz(0), use_jump(false){
p_rows.resize(rows);
row_nnz.resize(rows);
}
//reserve number of cols in each row
void reserve(const unsigned& cols_per_row){
data.resize(rows*cols_per_row);
cols_id.resize(rows*cols_per_row);
for(unsigned i = 0; i < rows; ++i){
p_rows[i] = i*cols_per_row;
}
}
//get a ref to a coeff
T& coeff_ref(const unsigned& row, const unsigned& col){
if(use_jump){
int pos = jump_search(cols_id, p_rows[row], row_nnz[row] + p_rows[row], step, col);
if (pos != -1) {
return data[pos];
}
}
else{
for (unsigned i = p_rows[row]; i < row_nnz[row] + p_rows[row]; ++i) {
if (cols_id[i] == col) {
return data[i];
}
}
}
return insert(row, col);
}
//get coeff
T coeff(const unsigned& row, const unsigned& col){
for(unsigned i = p_rows[row]; i < row_nnz[row]+p_rows[row]; ++i){
if(cols_id[i] == col){
return data[i];
}
}
return 0.0;
}
//get a reference for insertion at row,col
T& insert(const unsigned& row, const unsigned& col){
cols_id[p_rows[row] + row_nnz[row]] = col;
row_nnz[row]++;
++nnz;
return data[p_rows[row] + row_nnz[row] - 1];
}
//insert the value in at row,col
template<typename S>
void insert(const unsigned& row, const unsigned& col, const S& in){
cols_id[p_rows[row] + row_nnz[row]] = col;
data[p_rows[row] + row_nnz[row]] = in;
row_nnz[row]++;
++nnz;
}
//print out the matrix
void print(){
for(unsigned i = 0; i < rows; ++i){
for(unsigned j = 0; j < cols; ++j){
std::cout << coeff(i,j) << " ";
}
std::cout << std::endl;
}
}
//rewrites all values to range from the vector's start
void compress(){
unsigned c = p_rows[0]+row_nnz[0];
for(unsigned i = 0; i < rows-1; ++i){
auto ind = p_rows[i+1];
for(unsigned j = 0; j < row_nnz[i+1]; ++j, ++ind){
data[c] = data[ind];
cols_id[c] = cols_id[ind];
++c;
}
p_rows[i+1] = p_rows[i] + row_nnz[i];
}
data.resize(nnz);
cols_id.resize(nnz);
}
void truncate() {
data.shrink_to_fit();
cols_id.shrink_to_fit();
}
//nullify all none zero entries
void nullify(const bool& par = true){
if(par){
std::for_each(std::execution::par_unseq, data.begin(), data.end(), [&](auto& i){
i = 0.;
});
} else{
std::for_each(std::execution::seq, data.begin(), data.end(), [&](auto& i){
i = 0.;
});
}
}
void sort(){
for(unsigned i = 0; i < p_rows.size(); ++i){
mirror_quicksort(cols_id,p_rows[i],row_nnz[i]+p_rows[i]-1,data);
}
if (row_nnz[1] >= 10){
use_jump = true; //switch to jump search
step = sqrt(row_nnz[1]);
}
}
//debug method for viewing data
void print_data(){
unsigned c = 0;
for(auto& i : data){
std::cout << i << " ";
++c;
if(c >= 11){
std::cout << std::endl;
c = 0;
}
}
}
void print_indices(){
for(unsigned i = 0; i < p_rows.size(); ++i){
std::cout << i << ": ";
for(int j = p_rows[i]; j < p_rows[i]+row_nnz[i]; ++j){
std::cout << cols_id[j] << " ";
}
std::cout << std::endl;
}
}
};