I am writing a finite element solver in C++. The main bottle neck is assembling the global stiffness matrix in sparse compressed row storage (so far I am only solving steady problems). Because I don't know how many nonzero entries exist in each row, I am currently assuming a constant upper bound on the number of non-zeros per row. My codes currently looks like:
int nrow[fem.Np+1]; //row vector in CRS format initialized to 0
int ncol[fem.Np*MAX_NNZ]; //column vector in CRS format initialized to -1
double A[fem.Np*MAX_NNZ]; //values vector in CRS format initialized to 0.0
//update nrow, ncol, and A
//MAX_NNZ is the constant upper bound on the number of nonzeros per row
//Np=='number of points (nodes)
//Npe=='number of points per element e.g. 3 for linear triangles
//kmat[Npe][Npe] is the previously computed local element stiffness matrix
for(int i=0;i<e.Npe;i++)
{
for(int j=0;j<e.Npe;j++)
{
r = e.nodes[i];
c = e.nodes[j];
for(int p=MAX_NNZ*r-MAX_NNZ;p<MAX_NNZ*r;p++)
{
if(ncol[p]==-1)
{
for(int l=r;l<fem.Np+1;l++)
nrow[l]++;
ncol[p] = c-1;
A[p] = e.kmat[i][j];
break;
}
else if(ncol[p]==c-1)
{
A[p] = A[p] + e.kmat[i][j];
break;
}
}
}
Is this the most efficient way to assemble the global matrix in CRS format? Are there improvements that I can make? Thanks.