Non-contiguous sparse matrix views

[Santiago Peñate Vera]

Hi,

I am programming a circuit simulation program, and I need to have submatrices that do not match the matrix block logic.

submatrix = sparse_matrix[row_indices, col_indices]

For Armadillo, I made the following routine:

/**
    * @brief sp_submatrix Function to extract columns and rows from a sparse matrix
    * @param A Sparse matrix pointer
    * @param rows vector of the rown indices to keep (must be sorted)
    * @param cols vector of the clumn indices to keep (must be sorted)
    * @return Sparse matrix of the indicated indices
    */
    template<typename T> arma::SpMat<T> sp_submatrix(arma::SpMat<T> *A, arma::uvec *rows, arma::uvec *cols) {

        // declare reduced matrix
        std::size_t n_rows = rows->size();
        std::size_t n_cols = cols->size();

        bool found = false;
        std::size_t n = 0;
        std::size_t p = 0;
        std::size_t found_idx = 0;

        arma::Col<T> new_val(A->n_nonzero);
        arma::uvec new_row_ind(A->n_nonzero);
        arma::uvec new_col_ptr(n_cols + 1);

        new_col_ptr(p) = 0;

        for (auto const& j: *cols) { // for every column in the cols vector

            for (std::size_t k = A->col_ptrs[j]; k < A->col_ptrs[j + 1]; k++) {  // k is the index of the "values" and "row_indices" that corresponds to the column j

                // search row_ind[k] in rows
                found = false;
                found_idx = 0;
                while (!found && found_idx < n_rows) {
                    if (A->row_indices[k] == rows->at(found_idx))
                        found = true;
                    found_idx++;
                }

                // store the values if the row was found in rows
                if (found) { // if the row index is in the designated rows...
                    new_val(n) = A->values[k]; // store the value
                    new_row_ind(n) = found_idx - 1;  // store the index where the original index was found inside "rows"
                    n++;
                }
            }

            p++;
            new_col_ptr(p) = n;
        }
        new_col_ptr(p) = n ;

        // reshape the vectors to the actual number of elements
        new_val.reshape(n, 1);
        new_row_ind.reshape(n, 1);

        // return new object
        return arma::SpMat<T>(new_row_ind, new_col_ptr, new_val, n_rows, n_cols);
    }

This routine needs the sparse matrix to have a CSC structure. For Blaze, I don't see how I can access the internal elements so I can mimic the function above.

This is both a feature suggestion and a question about how to access the internal elements of a CSC (column major) sparse matrix in Blaze.

[Klaus Iglberger]

Hi Santiago!

In Blaze you can solve this via row selections and column selections. The following example gives you an impression how to select the first and third row and the first and second column of a 4x5 column-major sparse matrix:

blaze::CompressedMatrix<int,columnMajor> A{ { 1, 0, 2, 0, 3 }
                                          , { 0, 4, 0, 5, 0 }
                                          , { 6, 7, 0, 0, 0 }
                                          , { 0, 0, 8, 9, 0 } };

const std::vector<size_t> rowsToKeep{ 1UL, 3UL };
const std::vector<size_t> colsToKeep{ 1UL, 2UL };

blaze::CompressedMatrix<int,columnMajor> B( columns( rows( A, rowsToKeep ), colsToKeep ) );

std::cout << "\n B =\n" << B << "\n\n";

If you know the row and column indices at compile time, you can even omit the vectors and pass the indices as template parameters:

blaze::CompressedMatrix<int,columnMajor> A{ { 1, 0, 2, 0, 3 }
                                          , { 0, 4, 0, 5, 0 }
                                          , { 6, 7, 0, 0, 0 }
                                          , { 0, 0, 8, 9, 0 } };

blaze::CompressedMatrix<int,columnMajor> B( columns<1UL,2UL>( rows<1UL,3UL>( A ) ) );

std::cout << "\n B =\n" << B << "\n\n";

I hope this is exactly the feature you are looking for. If not, please give an example of what you would like to achieve.

Best regards,

Klaus!

[Santiago Peñate Vera]

Hi Klaus,

Thank you for the prompt reply!

Indeed this is what I need.

[Klaus Iglberger]

Duplicate of #48.