Block sparse matrix with sparse matrix

[Li__Yapeng]

I have a sparse matrix M, and a sparse matrix K,  both of the sparse matrices obtain the same non-zeros location, how to form a sparse matrix with following format:

A= [ M,  0.2K;  0.1M, 0.1K];

This was used to solve linear equation with Ax = b.

Thanks a lot.

[Pamela_H_Intel]

Yapeng,

Please explain "obtain the same non-zeros location". Obtain them from where? Or do you mean the matrices have non-zero values in the same row,column positions?

Also, please identify or explain the format you are using.

A= [ M,  0.2K;  0.1M, 0.1K];

Pamela

[Spencer_P_Intel]

So, this is a sort of a tricky case to do with existing functionality,  you could probably use the mkl_sparse_?_add functionality, but that would require you to call it 3 times and adjust some arrays as you go with their offsets.  I would probably instead write if from scratch.  Maybe something like the following:

Two steps to do full computation.  I will assume you are using CSR format to represent the matrices. For this to make sense, M and K must have the same number of rows, but columns can be different, so suppose M is nxm  and K is nxk., then your new sparse matrix has global dimensions (nrows=2n) x (ncols=m+k).  Let M be represented by Mrowptr, Mcolumns and Mvalues arrays, and K be represented by Krowptr, Kcolumns, Kvalues arrays.  Finally let rowptr, columns and values arrays be the CSR representation for the new A matrix.

1.  Construct new rowptr array and initialize memory for new matrix.  

  Let us consider the number of sparse elements of first row (same as the m+1 row).  We will use 0 based (c-style) indexing in our arrays.    The the first row of new matrix has Mrowptr[1]-Mrowptr[0] + Krowptr[1]-Krowptr[0] nonzero elements.   then rowptr[0] = 0  and rowptr[1] = Mrowptr[1]-Mrowptr[0] + Krowptr[1]-Krowptr[0]

You can go through and calculate new row lengths independently of each other.  So if using parallelism, this is a simple parallel for loop (OMP or TBB) otherwise it is just a for loop.  I would recommend filling in each value with the number of nonzero elements on that row.  Then you do a second pass to calculate the prefix sum (partial sums increasing up to the end) of this rowptr array.  Once this is done, you know how big the columns arrays need to be and where each row can start. Now allocate memory for the new columns and values arrays.

2.  Now you can fill in proper values of the column and values arrays. 

Again this can be done independently on each row so a simple (parallel) for loop will handle the computations. For columns of M being copied, there are no offsets, it can be copied as is, but for the columns of K, you need to offset (read this as add) by m (the number of columns in M) to shift it to the right blocks.  The values for each row can be copied and scaled as desired to match your 4 blocks,  one row at a time...

I hope that was detailed enough for you to see how to do it.  Best of luck!

Spencer