Parallellizable Algorithms to traverse a 2D matrix being aware of both col/row-wise neighborhood

Question

I have a fairly large N*N integer matrix Matrix2D (assume sufficient memory),

1, within each row/column, I need to record the col/row index of an element if it's value is different than it's right/lower neighbour.

2, I want to find an optimal algorithms that is parallellizable, ideally by OMP.

So, finally I will have some data structures like,

std::vector<std::vector<int>>   RowWiseDiscontinuity(N);// N= #of rows
std::vector<std::vector<int>>   ColWiseDiscontinuity(N);// N= #of cols

where inner std::vector<int> records the row/col indices.

I put my serial version here but find it difficult to be OMP parallelized... Someone could provide some idea how to implement the traversing over this 2D matrix with omp?

code snippet,

std::vector<std::vector<int>>   RowWiseDiscontinuity(N);// N= #of rows
std::vector<std::vector<int>>   ColWiseDiscontinuity(N);// N= #of cols
std::vector<int> TempX1;
std::vector<int> TempX2;
for (int y=0; y<N; ++y)
{
    TempX1.clear();
    for (int x =0; x<N; ++x)
    {
        int value = Matrix2D(x,y);
        TempX1.push_back(value);
    }

    auto iter1 = TempX1.begin();
    auto iter2 = TempX2.begin();

    if (y>0)
    for (int x =0; x<N; ++x)
    {
         if (*iter1 !=*(iter1+1))
         {
             RowWiseDiscontinuity[y].push_back(x); //Critical for OMP
         }
         ++iter1;
         ++iter2;
         if (*iter1 != *iter2)
         {
             ColWiseDiscontinuity[x].push_back(y); //Critical for OMP
         }
     }

     TempX2.swap(TempX1); // proceed to next row, remember previous

}

Answer 1

I would create another array that holds the nearest neighbor both column and row-wise. This would have to be done as a first pass, obviously. I recommend creating a 2d array of pairs (pair) that holds the indices you want. Instead of two vectors, I would do a vector of pairs. Pairs are parallellizable and easily sorted.

vector<vector<pair<int, int>>> elements(N);

Answer 2

Perform two passes (which can be performed on different threads) on the matrix, one for the row wise discontinuities and the other for the column wise discontinuities.

The row pass looks as follows:

for (int y = 0; y < N; ++y)    // Can be parallelized
{
    for (int x = 0; x < N - 1; ++x)
    {
        if(Matrix(x, y) != Matrix(x + 1, y))
            RowWiseDiscontinuity[y].push_back(x);
    }
}

The column pass is similar:

for (int x = 0; x < N; ++x)    // Can be parallelized
{
    for (int y = 0; y < N - 1; ++y)
    {
        if(Matrix(x, y) != Matrix(x, y + 1))
            ColWiseDiscontinuity[x].push_back(y);
    }
}

The outer loop in both cases can be parallelized. A different element of Row / ColWiseDiscontinuity is mutated in each iteration of the outer loop, thus preventing a data race. The passes themselves can be performed on different threads.

As a side note, you can optimize this algorithm further by reducing cache misses (at the expense of memory) by storing the matrix in both row-major and column-major order, and using each order when appropriate. In a row-major order, element (x + 1, y) is always next to (x, y). The same holds true for element (x, y + 1) in a column-major order.

Answer 3

Here is an algorithm that does a base test of finding the adjacent diagonal neighbor and records the results using a 4x4 identity matrix. This does not include any use of OMP or parallel computing. However this is a generic class template of a MxN matrix that is simple enough to use. Instead of storing the contents in a vector of vectors; I've flattened out the data to a single 1D vector and the amount of memory is already preserved upon instantiation of the template. I'm using a function template to compare elements from within a matrix passing back the indexes (M,N) or (x,y) as well as if the result was true or false. I use a struct here to contain the relationship of the x-y indexes & the bool result. The heuristics of checking the neighbors avoids looking in the last column & last row of the matrix since there will not be any elements father to the right nor farther down: this can be seen within the main function. This may be of help to you to where you can try to apply the class, struct & function to the OMP library.

template<unsigned Col, unsigned Row>
class Matrix2D {
public:
    const unsigned col_size = Col;
    const unsigned row_size = Row;
    const unsigned stride_ = col_size;
    const unsigned matrix_size = col_size * row_size;

private:
    std::vector<int> data_;

public:
    Matrix2D() {
        data_.resize( matrix_size );
    }

    void addElement( unsigned x, unsigned y, int val ) {
        data_[(x * col_size + y)] = val;
    }

    /*int getElement( unsigned x, unsigned y ) {
        int value = data_[(x * col_size + y)];
        return value;
    }*/

    int getElement( unsigned idx ) {
        return data_[idx];
    }
};

struct Neighbor {
    unsigned indexCol;
    unsigned indexRow;
    bool     notSame;
};


template<unsigned Col, unsigned Row>
void compareMatrixDiagonals( Matrix2D<Col, Row>& mat, Neighbor& n, unsigned colIdx, unsigned rowIdx );

int main() {

    Matrix2D<4, 4> mat4x4;
    mat4x4.addElement( 0, 0, 1 );
    mat4x4.addElement( 0, 1, 0 );
    mat4x4.addElement( 0, 2, 0 );
    mat4x4.addElement( 0, 3, 0 );

    mat4x4.addElement( 1, 0, 0 );
    mat4x4.addElement( 1, 1, 1 );
    mat4x4.addElement( 1, 2, 0 );
    mat4x4.addElement( 1, 3, 0 );

    mat4x4.addElement( 2, 0, 0 );
    mat4x4.addElement( 2, 1, 0 );
    mat4x4.addElement( 2, 2, 1 );
    mat4x4.addElement( 2, 3, 0 );

    mat4x4.addElement( 3, 0, 0 );
    mat4x4.addElement( 3, 1, 0 );
    mat4x4.addElement( 3, 2, 0 );
    mat4x4.addElement( 3, 3, 1 );

    unsigned idx = 0;
    for ( unsigned i = 0; i < mat4x4.matrix_size; i++ ) {
        std::cout << mat4x4.getElement( i ) << " ";
        idx++;

        if ( idx == 4 ) {
            std::cout << "\n";
            idx = 0;
        }
    }
    std::cout << "\n";    

    unsigned colIdx = 0;
    unsigned rowIdx = 0;
    std::vector<Neighbor> neighbors;
    Neighbor n;

    // If we are in the last col or row we can ignore
    // (0,3),(1,3),(2,3),(3,3),(3,0),(3,1),(3,2), {*(3,3)* already excluded}
    // This is with a 4x4 matrix: we can substitute and use LastCol - LastRow 
    // for any size MxN Matrix.
    const unsigned LastCol = mat4x4.col_size - 1;
    const unsigned LastRow = mat4x4.row_size - 1;

    for ( unsigned i = 0; i < LastCol; i++ ) {
        for ( unsigned j = 0; j < LastRow; j++ ) {
            compareMatrixDiagonals( mat4x4, n, i, j );
            neighbors.push_back( n );
        }
    }

    for ( unsigned i = 0; i < neighbors.size(); i++ ) {
        std::cout << "(" << neighbors[i].indexCol
            << "," << neighbors[i].indexRow
            << ") " << neighbors[i].notSame
            << "\n";
    }

    std::cout << "\nPress any key & enter to quit." << std::endl;
    char c;
    std::cin >> c;

    return 0;
}

template<unsigned Col, unsigned Row>
void compareMatrixDiagonals( Matrix2D<Col, Row>& mat, Neighbor& N, unsigned colIdx, unsigned rowIdx ) {
    unsigned firstIdx = (colIdx * mat.col_size + rowIdx);
    unsigned nextIdx  = ((colIdx + 1) * mat.col_size +  (rowIdx + 1));
    if ( mat.getElement( firstIdx ) != mat.getElement( nextIdx ) ) {
        N.indexCol = colIdx;
        N.indexRow = rowIdx;
        N.notSame  = true;          
    } else {
        N.indexCol = colIdx;
        N.indexRow = rowIdx;
        N.notSame  = false;     
    }
}