How to efficiently "remove" non-zero elements of a matrix of C ++ matrices

I worked on an algorithm for removing nonzero elements from a large sparse (ish) matrix as a way to reduce the computation time of a numerical solution for a large system of connected ODEs of the form of a logical equation .

The basic idea is that each T-timesteps you evaluate which elements in the state matrix, B, are greater than a certain threshold. If those deemed insignificant (i.e. below the threshold), you remove (set to zero) the corresponding row / column in the sparse interaction matrix A so that fewer flops are required to compute AxB. (Note that the diagonal remains non-zero so that the associated elements of the state matrix do not explode)

The algorithm looks like this:

• Scan through vector B for items that are less than create and write space

• Iterating through the nonzero elements of matrix A corresponding to the locations written in 1.and writing the locations

• Set the A elements corresponding to the locations written from 2. to zero.

I need to separate steps 2. and 3. because if you set the non-zero elements of the sparse matrix (arma :: sp_mat) to zero, you will destroy the iterators and therefore cannot achieve the result in one pass.

Can anyone suggest any suggestions on how to improve this algorithm?

(Note, I am not a computer scientist, and so there may be an obvious improvement that I am not aware of. MMMM Many thanks for trying !!)

#include <iostream>
#include <armadillo>


using namespace std;
using namespace arma;
int main() {

mat A = {{1, 0.2, 0, 0.2, 0, 0, 0.2, 0, 0.2, 0},
         {0, 1, 0.2, 0, 0.2, 0, 0, 0.2, 0 , 0},
         {0, 0, 1, 0.2, 0, 0, 0.2, 0, 0, 0.2},
         {0.2, 0, 0, 1, 0.2, 0, 0.2, 0.2, 0, 0.2},
         {0, 0.2, 0, 0.2, 1, 0.2, 0, 0, 0, 0.2},
         {0, 0, 0.2, 0, 0.2, 1, 0, 0.2, 0, 0},
         {0.2, 0, 0, 0.2, 0.2, 0, 1, 0, 0, 0.2},
         {0, 0.2, 0, 0.2, 0, 0, 0, 1, 0.2, 0},
         {0, 0, 0.2, 0, 0.2, 0.2, 0, 0, 1, 0},
         {0, 0, 0, 0.2, 0, 0.2, 0, 0, 0.2, 1}};
sp_mat Asp(A); // Interaction matrix in sparse representation

vec B = {1, 1, 0.1, 1, 0.1, 1, 1, 0.1, 1, 1};

double thresh = 0.2;

vec index(Asp.n_rows); // obj for storing locations of insignificant elements b

int count = 0;
for (int i=0; i<B.n_rows; i++) {
    if (B(i) < thresh) {
        index(count) = i;
        count ++;
    }
}

index.resize(count);

sp_mat tmp = Asp; // temp copy of interaction matrix

for (int i=0; i<index.n_rows; i++) {
    // iterate through non-zero elements in col corresponding to each element of index
    sp_mat::iterator bc = tmp.begin_col(index(i));
    sp_mat::iterator ec = tmp.end_col(index(i));

    // iterate through non-zero elements in row corresponding to each element of index
    sp_mat::row_iterator br = tmp.begin_row(index(i));
    sp_mat::row_iterator er = tmp.end_row(index(i));

    // mat for storing locations of non-zero elements to be removed
    mat loc(2, 2 * tmp.n_cols);

    int count2 = 0;
    for (auto j = bc; j != ec; j++) {
        loc(0, count2) = j.row();
        loc(1, count2) = j.col();
        count2++;
    }

    for (auto j = br; j != er; j++) {
        loc(0, count2) = j.row();
        loc(1, count2) = j.col();
        count2++;
    }

    loc.resize(2, count2);

    for (int j = 0; j < loc.n_cols; j++) {
        tmp(loc(0, j), loc(1, j)) = 0.0;
    }
}

tmp.diag().ones();

cout << "Number of non-zero terms interaction matrix = " << Asp.n_nonzero << endl;
cout << "Number of non-zero terms reduced matrix = " << tmp.n_nonzero << endl;
}

      

        
        
        
      

    

Output:

Number of non-zero terms interaction matrix = 45
Number of non-zero terms reduced matrix = 24