FEM Sparse Matrix Assembly in CSR format

[elmeliegy__abdelrahm]

Hi;

I have C++ FEM code wrote long time ago with its own functions. Now, i am trying to upgrade this code to be faster. i am trying to convert the existing linear system function to the mkl linear system functions. the bottelneck here is the assembly of the global stiffness matrix in csr format. i want to use PARDISO solver, so i have to assemble the global matrix in csr format taking in consideration that we do not know the No. of nonzero elements (NNZ) so we have to maximize NNZ to the max possible number. 

Is there any exsiting function/s that can do this?

Any Ideas will be appreciated.

[mecej4]

You may have to make two passes over the mesh. In the first pass, only the connectivity is of interest, and the calculations to find the element contributions to the global matrix can be skipped. During this process, there can be multiple contributions to A_ij coming from adjacent elements, and book-keeping work is needed to keep track of and consolidate these contributions. At the end of this pass, the row index array IA will have been formed. The final value in this array will be nnz, and that value is used to allocate the column index and value arrays JA and A.

During the second pass, the element contributions are calculated and inserted into their proper places in A.

Alternatively, you can assemble the matrix in coordinate oriented format (COO), and call mkl_csrcoo() to obtain the CSR data.

[elmeliegy__abdelrahm]

Thank you for your response

Is that very tricky to be implemented directly?, i mean , i spent some time searching but i did not find a way to implement it. even, i did not find a way to assemble in coo format.

[mecej4 · Accepted solution]

There is nothing tricky about assembling FEM global matrices. For each element, you have to determine the mapping between local node numbers (e.g., 1-2-3-4 for a four-node element) and global node numbers. If you write A(i,j) = expr. when assembling a dense A matrix, instead you write:

   nnz = nnz+1

   row(nnz) = i

   col(nnz) = j

   val(nnz) = expr

When you have done this for all elements in the mesh, you have the COO representation of the global matrix. There will be some multiple entries, that is, more than one entry in VAL may be present for the same values of ROW and COL. to consolidate them, you have to sort by row and column, and sum up entries that share row and column, thus obtaining the CSR data. MKL_?CSRCOO does this conversion, if you wish to use it.

[elmeliegy__abdelrahm]

I am starting learning C++ and how to link it with mkl. I tried this simple/poor code to assemble 2-bars stiffness matrix in coo and then convert it using mkl_dcsrcoo. The assembly to coo is ok. Now, i do not understand what is the problem when i am trying to use mkl_csrcoo?.

this may dumb questions, sorry for that and thanks in advance.

Here is the code:

//coo format variables

std::vector<double> val; //non-zero values

std::vector<MKL_INT> ia; // row indices

std::vector<MKL_INT> ja; // column indices

//csr format variables

double Acsr;

MKL_INT AI; 

MKL_INT AJ; 

//local stiffness matrix for element 1

double lhs1[2][2] = { {1,-1},{-1,1} };

//local stiffness matrix for element 2

double lhs2[2][2] = { { 2,-2 },{ -2,2 } };

//connectivity array for element 1

int connect1[] = { 0,1 };

//connectivity array for element 2

int connect2[] = { 1,2 };

//mumber of nodes/element

int n = 2;

int nnz = 0;

//for element 1

for (int i = 0; i < n; i++) {

for (int j = 0; j < n; j++) {

ia.push_back(connect1);

ja.push_back(connect1);

val.push_back(lhs1);

nnz++;

}

}

//for element 2

for (int i = 0; i < n; i++) {

for (int j = 0; j < n; j++) {

ia.push_back(connect2);

ja.push_back(connect2);

val.push_back(lhs2);

nnz++;

}

}

//At the end of this loop , val, ia and ja filled in with correct numbers

job[0] = 1;

job[5] = 2;

info = INFO;

//Here, i want to convert the val,ia,ja (COO format) to CSR format (AI,AJ,Acsr)

mkl_dcsrcoo(job, &n, Acsr, AJ, AI, &nnz, val, ia, ja, &info);

[mecej4]

You have shown some bits of code. Depending on what is in the other lines of code in the program, and on what was left out but should have been present, the program may or may not work. We would not write a program to handle a problem with exactly two elements. With a small number of elements, a dense matrix solver, such as one of the Lapack routines, would be easier to use.

Your job[] values are definitely wrong. Consult the documentation, look at the examples provided with MKL, and correct your mistakes.

[elmeliegy__abdelrahm]

This is a separate trial code i have made for myself to try how things work. I know that i should make a general code for any element/number of element. I am sorry for no attaching the error message. But i received this messages while trying to compile the code. 

Error C2664 'void mkl_dcsrcoo(const int *,const int *,double *,int *,int *,int *,double *,int *,int *,int *)': cannot convert argument 3 from 'double' to 'double *' COOtoCSR c:\users\aaelmeli\source\repos\cootocsr\cootocsr.cpp 90

Error (active) E0167 argument of type "double" is incompatible with parameter of type "double *" COOtoCSR c:\Users\aaelmeli\source\repos\COOtoCSR\COOtoCSR.cpp 90

Error (active) E0167 argument of type "int" is incompatible with parameter of type "int *" COOtoCSR c:\Users\aaelmeli\source\repos\COOtoCSR\COOtoCSR.cpp 90

Error (active) E0167 argument of type "int" is incompatible with parameter of type "int *" COOtoCSR c:\Users\aaelmeli\source\repos\COOtoCSR\COOtoCSR.cpp 90

Error (active) E0413 no suitable conversion function from "std::vector<double, std::allocator<double>>" to "double *" exists COOtoCSR c:\Users\aaelmeli\source\repos\COOtoCSR\COOtoCSR.cpp 90

Error (active) E0413 no suitable conversion function from "std::vector<int, std::allocator<int>>" to "int *" exists COOtoCSR c:\Users\aaelmeli\source\repos\COOtoCSR\COOtoCSR.cpp 90

Error (active) E0413 no suitable conversion function from "std::vector<int, std::allocator<int>>" to "int *" exists COOtoCSR c:\Users\aaelmeli\source\repos\COOtoCSR\COOtoCSR.cpp 90

[mecej4]

The error messages are perfectly clear. For instance, take the third argument to mkl_dcsrcoo(). What should its type be, according to the documentation? What is the type that you are passing? What should its value be? What is the value that you are passing?

These are matters that you must learn to handle by yourself.