How to get a lower-right part from sparse matrix?

1
votes

What's the best way to get a lower-right quarter of sparse matrix in "coordinates format"? I need it for Schur complement.

Coordinates are indexed from 1 (I'm using MUMPS).

EXAMPLE:

int order_of_matrix = 4

int row_coords[] = {1,1,2,3,3,4};
int col_coords[] = {1,3,3,3,4,4};
int vals[]       = {1,2,3,4,5,6};

So my matrix should look like this:

|1 0 2 0|
|0 0 3 0|
|0 0 4 5|
|0 0 0 6|

And now, I need to get values from the lower-right corner like this:

int schur_complement_size = 1
|6|

int schur_complement_size = 2
|4 5|
|0 6|

int schur_complement_size = 3
|0 3 0|
|0 4 5|
|0 0 6|

int schur_complement_size = 4
|1 0 2 0|
|0 0 3 0|
|0 0 4 5|
|0 0 0 6|

Is there any elegant way to achieve this without converting matrix to the dense format? The dense format is very problematic for me, because matrices in real usage of my program are going to be very large.

3
c++cmatrix
First thing I'd do is to create a struct so you don't have to keep three vectors of equal size for the different parts of the matrix. Then, consider which parts of the existing matrix for the new matrix and how the elements differ. Lastly, decide for one language, as the answers will differ. – Ulrich Eckhardt
I must use several arrays, it's the way it was used it this project even before me. I've inserted 'C++' into the title, so now it's more specific. And I don't understand the part with "consideration" - I simply want to extract the "lower-right" corner of the matrix... – Eenoku
Are the row and column vectors indexing the nonzero values guaranteed to be sorted? If so - how? That would definitely change the "most elegant" solution over a brute force approach because you could do a binary search to find the starting point. For huge sparse matrices the speed up could be significant. – Floris
Well, if you want to keep the trouble of using three different, dynamically allocated vectors and the ensuing memory management burden, so be it. Anyhow, logically you still have triples of row, columns and value, regardless of how you store them. What I wanted you to think about is how to determine whether a triplet is part of the matrix to extract from. Then, think about which new triplet will be used in the new matrix. – Ulrich Eckhardt

3 Answers

2
votes

Assuming that you cannot change how to model your sparse matrix, you can try this:

  • Make sure that row_coords is ordered (as in the example you provided). Note that you can choose col_coords as well, what you need is that one of the two arrays is ordered.
  • calculate q = order_of_matrix - order_of_submatrix + 1
  • run a binary search in the sorted array (row_coords) to find out the index of the first element row_coords[i]=q

  • now you can create the submatrix using 3 new arrays (you know their maximum size) and adding an element only if col_coords[i]>=q.

    sub_col_coords[j] = col_coords[i] - q + 1; sub_row_coords[j] = col_coords[i] - q + 1; sub_vals[j] = vals[i];

0
votes

For your task, I have written a working example and have used std::vector for outputs if you don't mind because memory management with dynamic allocation may get complicated:

// Example program
#include <iostream>
#include <vector>

void GetSubSparseMatrixLowerRight(int order_of_matrix , int pointCount, int * row_coords, int *col_coords, int * vals, int schur_complement_size, std::vector<int> &row_coordsOut, std::vector<int> &col_coordsOut, std::vector<int> &valsOut)
{
    int upperLeftRow = order_of_matrix - schur_complement_size;
    int upperLeftCol = order_of_matrix - schur_complement_size;

    for(int i=0; i<pointCount; i++)
    {
        if(row_coords[i] > upperLeftRow && col_coords[i] > upperLeftCol)
        {
            row_coordsOut.push_back(row_coords[i] - upperLeftRow);
            col_coordsOut.push_back(col_coords[i] - upperLeftCol);
            valsOut.push_back(vals[i]);
        }
    }
}
int main()
{
  int order_of_matrix = 4;

  int row_coords[] = {1,1,2,3,3,4};
  int col_coords[] = {1,3,3,3,4,4};
  int vals[]       = {1,2,3,4,5,6};

  int schur_complement_size = 4;

  std::vector<int> row_coordsOut;
  std::vector<int> col_coordsOut;
  std::vector<int> valsOut;

  GetSubSparseMatrixLowerRight(order_of_matrix , sizeof(row_coords)/sizeof(row_coords[0]), (int *)row_coords, (int *)col_coords, (int *) vals, schur_complement_size, row_coordsOut,col_coordsOut, valsOut);

  std::cout<<valsOut[0]<<" "<<row_coordsOut[0]<<" "<<col_coordsOut[0]<<std::endl; //to control the output
}

To simplify the function arguments, you an define a struct or a class to represent a sparse matrix with member row coordinates, column coordinates and values, as mentioned in the comments. Nevertheless, this code gets the job done.

0
votes

If there is no special structure to how you're storing the nonzero entries, then essentially the only algorithm is to simply iterate over all of the entries, copying the ones you want to a new matrix.

Proof sketch: You have to inspect all of the entries, otherwise you might miss a relevant entry; the lack of structure means that you can't rule any position out.

(loophole: if your format contains no repeated entries, and you've read enough entries to completely fill a dense matrix, then you can abort)

Since you have to read every entry anyways, it's hard to imagine any algorithm can be better than just reading through them and taking the ones you want.

Examples of "special structure" that could potentially allow faster algorithms are: the entries are stored in sorted order (lexicographically by coordinates), or maybe in some other clever order (e.g. along space filling curve) that makes it easier to extract blocks, or you keep around an auxiliary data structure that aids this operation.

There could be an interesting question along the lines of "What data structure can I use to make this operation easier?" But that needs to have clearly delineated requirements. (and should be asked in a different posting)