Как посчитать ненулевые значения на строку в диагональном блоке матрицы CSR - PullRequest
Новая
       6

Как посчитать ненулевые значения на строку в диагональном блоке матрицы CSR

0 голосов
/ 02 сентября 2018

У меня есть матрица в формате CSR, и мне нужен вектор c ++, который содержит количество ненулевых записей (количество) на строку, ограниченное квадратом, диагональных блоков разных размеров. 

 // The matrix in CSR format
    std::vector<int> row_idx = {0,2,4,6,10,13}; // size n+1 where 0-n are the idx of the row starts in values and column_idx and n+1 the TOTAL number of values
    std::vector<int> values = {1,6,2,7,3,8,10,11,4,9,12,13,5}; // nonzero matrix values
    std::vector<int> column_idx = {0,3,1,3,2,4,0,1,3,4,2,3,4};  // column indices of the values

В приведенном ниже примере есть два блока разных размеров A и B (интересующие блоки всегда квадратные и по диагонали).

Желаемым результатом для этого примера будет nnz_in_ranges [n] = {1,1,2,2,3}, но так как он должен быть встроен в другую подпрограмму, я в первую очередь ищу подпрограмму для вычисления отдельных блоков с использованием C ++. Как то так: 

// block A
int rangeStart = 0;
int rangeEnd = 2;

// block B
//int rangeStart = 2;
//int rangeEnd = n;

for (int i = rangeStart; i<rangeEND; ++i)
{
    nnz_in_ranges[n] = ...
}

// desired result for block A: nnz_in_ranges[n] = {1,1,0,0,0}
// desired result for block B: nnz_in_ranges[n] = {0,0,2,2,3}

Я пытался решить ее с помощью функции std :: count ..., но не смог расширить приведенный ниже код, считая ненулевые значения на строку, так как я не смог ввести диапазон столбцов.

У кого-нибудь есть идеи, как подойти к этой проблеме? 

#include <stdio.h>
#include <iostream>
#include <vector>
#include <algorithm>

int main()
{

// NxN matrix example

/*
index     0   1    2   3   4
        ______________________
     0  | 1   0  | 0   6   0 |
        |   A    |           |
     1  | 0   2  | 0   7   0 |
        |--------------------|
     2  | 0   0  | 3   0   8 |
        |        |   B       |
     3  | 10  11 | 0   4   9 |  expected result: nnz_in_ranges[n] = {1,1,2,2,3} 
        |        |           |  here ranges are A and B
     4  | 0   0  |12  13   5 |
        ----------------------

*/


// matrix in CSR format

    int n = 5; // matrix size
    int nnz = 13; // number of nonzero values

    // The matrix in CSR format
    std::vector<int> row_idx = {0,2,4,6,10,13}; // size n+1 where 0-n are the idx of the row starts in values and column_idx and n+1 the TOTAL number of values
    std::vector<int> values = {1,6,2,7,3,8,10,11,4,9,12,13,5}; // nonzero matrix values
    std::vector<int> column_idx = {0,3,1,3,2,4,0,1,3,4,2,3,4};  // column indices of the values

    std::vector<int> tmp = {0,0,0,0,0,0,0,0,0,0,0,0,0};

    std::vector<int> sum(n);

    // count nonzeros per row  sum[] = {2,2,2,4,3}
    for(size_t i = 0; i < row_idx.size()-1; ++i) {
    sum[i] = std::count(tmp.begin() + row_idx[i], tmp.begin() + row_idx[i + 1], 0);
    }


    std::cout << "nnz_in_range = " << std::endl;
    for (int i=0; i<n; i++)
    {
    std::cout << ' ' << sum[i];
    }

return 0;

}

1 Ответ

0 голосов
/ 02 сентября 2018

Разреженные матрицы могут явно содержать ноль записей. Из вашего вопроса не ясно, хотите ли вы посчитать позиции или фактические значения. Я буду считать первое, потому что ваш счетный код не использует values.

Тогда мы просто следуем определению формата CSR и получаем что-то вроде этого: 

std::vector<int> count_positions_in_block(int block_begin, int block_size, 
    const std::vector<int>& row_idx, const std::vector<int>& column_idx)
{
    std::vector<int> cnt(block_size, 0);

    const auto block_end = block_begin + block_size;
    assert(block_end < row_idx.size());

    for (auto row = block_begin; row < block_end; ++row)
    {
        auto first = row_idx[row];
        auto last = row_idx[row + 1];
        assert(first <= last);

        for (auto i = first; i < last; ++i)
            if (column_idx[i] >= block_begin && column_idx[i] < block_end)
                ++cnt[row - block_begin];
    }

    return cnt;
}

auto nnz1 = count_positions_in_block(0, 2, row_idx, column_idx);
// nnz1 = [1, 1]
auto nnz2 = count_positions_in_block(2, 3, row_idx, column_idx);
// nnz2 = [2, 2, 3]

Его можно переписать с помощью std::count_if: 

std::vector<int> count_positions_in_block(int block_begin, int block_size, 
    const std::vector<int>& row_idx, const std::vector<int>& column_idx)
{
    std::vector<int> cnt(block_size, 0);

    const auto block_end = block_begin + block_size;
    assert(block_end < row_idx.size());

    const auto is_in_block = [&](auto col)
        { return (col >= block_begin && col < block_end); };

    for (auto row = block_begin; row < block_end; ++row)
    {
        auto first = row_idx[row];
        auto last = row_idx[row + 1];
        assert(first <= last);

        const auto cib = column_idx.begin();
        cnt[row - block_begin] = std::count_if(
            cib + first, cib + last, is_in_block);
    }

    return cnt;
}
...