Sparse addition

[Amar_K_1]

Hello!

My question is regarding the addition of a sparse matrix on to itself. For e.g in the dense sense, performing an operation like A = A + B, where A and B are dense matrices.

Speaking in the sparse CSR sense, A may be represented by ia,ja and aval and B may be represented by ib,jb,bval and C may be represented by ic,jc and cval.

1. operation: A = B + C

Since the number of non-zero entries in A is unknown, mkl_dcsradd is called first with request = 1. After this call is executed, the number of non-zero entries will be known and then mkl_dcsradd is called with request = 2. At this point, the result of the addition is known and stored in ia,ja and aval.

Problem 1: I don't clearly understand the use of pointers here, but what I have seen is that though ia, ja and aval as per the CSR format are arrays, mkl_?csradd accepts and outputs ia, ja and aval as scalars. Printing them gives the notion of them being arrays to the end user, but they actually are scalars.

(code attached with output )

2. Considering operation: A = A + B

Problem 2: Suppose after the operation 1 is performed, I need to perform the operation A = A + B, it wouldn't be possible to do this straight away because the input arguments to the call of mkl_?csradd require the CSR representation of A (3 1D arrays). On the other hand, the data I would have is just 3 scalars (not arrays!). Is there a solution to this problem? Suppose we can circumvent this problem, will it be a problem to then perform A = A + B, because there is an inconsistency in the size of A on LHS and RHS, as mkl would see? From my point of view, I'm trying to just make the matrix A bigger and am wanting to just add all the entries in B to the same number of entries in A, hence no inconsistency, physically speaking.

I hope I have presented the problem clearly. The code and output attached should hopefully make things clear!

Many Thanks,

Amar

[mecej4]

I found your post rather confusing, since you mix dense and sparse matrices in your discussion and the program. Your code also manipulates a global matrix G, which added to the confusion.

Let us focus on what I think is the main task: adding two sparse matrices A and B. You seem to have adapted the code in the MKL example file dcsr_addition.f, and the addition is done correctly in your code. The matrix C = A + B is output correctly by the code. Your code initializes A and B as dense matrices represented as 2-D arrays, and goes through the conversion to CSR representation, but this is a secondary task.

Now, the POINTER statement in the code, which seems to be a stumbling block, is a declaration of what are often called CRAY POINTERs. These pointers are quite different from (and predate) Fortran 9X pointers. You can read about them in the Intel Fortran documentation. Cray pointers are used in the MKL example code that I cited and in your code, but are best avoided in new code. You can easily avoid them and allocate the ic, jc, vc arrays for matrix C if you do not mind a bit of over-allocation for jc and vc, noting that nnz_C cannot exceed nnz_A + nnz_B . Therefore, you can allocate jc and vc of size nnz_A + nnz_B , call mkl_?csradd to populate ic, jc and vc, and trim the resulting arrays using REALLOC if you wish.

If, on the other hand, you wish to use Cray pointers, you will need to read about their properties and use them accordingly.

[SergeyKostrov]

>>...My question is regarding the addition of a sparse matrix on to itself... It looks like a multiplication of the sparse matrix by a scalar value, that is, by 2 in your case. Am I wrong?

[Amar_K_1]

mecej4 wrote:

I found your post rather confusing, since you mix dense and sparse matrices in your discussion and the program. Your code also manipulates a global matrix G, which added to the confusion.

Thank you for your reply!

The code might have seemed confusing because the initial bit performs the addition I intend to do in the dense sense. The global matrix is what I need. Since the problem in hand will involve millions of degrees of freedom, I would like to construct the global matrix in the sparse sense, to save memory. The code is not complete, but in the end the aim is to construct the global matrix sparsely and densely and check the results for validation purposes.

I'm still wondering if A = A + B will be a legal operation?

[Amar_K_1]

Sergey Kostrov wrote:

It looks like a multiplication of the sparse matrix by a scalar value, that is, by 2 in your case. Am I wrong?

I intend to perform A = A + B.

mkl_?dcsradd on the other hand is designed to perform A = A + (scalar_value * B). I just set the scalar_value to 1.

[mecej4]

Amar,

From your description, it appears that there is neither need for nor benefit from using a sparse matrix representation for the individual (FEM) element matrices, which are probably small, square, of fixed size, and dense. Instead, you should concentrate on:

• scanning all the small dense element matrices to obtain the number of nonzero entries in the global matrix -- this involves only the mapping of node/element numbers between the local and global matrices, involving only integer arithmetic
• process all the elements again, this time computing the element integrals and accumulating their values into the global matrix using the CSR representation built in the previous step.

Note that in this suggested scheme the element matrices are dense, and the global matrix is sparse. No attempt is to be made to store a sparse version of the element matrices (not needed, waste of time) and no attempt is made to store a dense version of the global matrix (waste of memory to hold zeroes, or not enough memory to store full matrix including zero entries).

I intend to perform A = A + B.

The compact mathematical notation obscures quite a bit of detail at the programming level. In the familiar dense matrix notation, the two A-s and the B are all rectangular arrays of the same dimensions. In a sparse matrix representation, each matrix is represented by two (or three) integer arrays and a real array. The number of nonzero entries of A on the left side can be more than the number of nonzero entries of the A on the right, so in many ways the A on the left is a new and different matrix, and the arrays may have different dimensions. If there is no overlap between the nonzero entries of A and B on the right, the number nnz_A of the nonzeroes of the A on the left will equal nnz_A + nnz_B from the right.

[Todd_R_Intel]

In his first post above, mecej4, mentioned use of allocate instead of CRAY POINTERS, but because of the copyright notices on the source files he decided not to post the files that he'd created. Because these files might be useful to the community, I'd like to share them on this thread. 

Thanks mecej4, for your many contributions on the Intel MKL user forum.

Todd