# (FEM) Efficient CRS vectors evaluation using elements connectivities

What is an efficient way of evaluating the column (col_ind) and the row (row_ptr) vectors for the CRS (Compressed Row Storage) storage format using the Connectivity Array? The Connectivity Array is a matrix with "E" rows (E = number of elements at the mesh) and "N" columns (N = number of nodes per element) where each row contains the node's labels for each element. I've already done that but the time complexity is of O(n²) and I'm trying to find out a faster way. The notation "col_ind" and "row_ptr" are presented at the following link: http://www.netlib.org/utk/people/JackDongarra/etemplates/node373.html

• What is $n$? The number of nodes per element? – Wolfgang Bangerth Nov 25 '19 at 15:39
• Professor Bangerth, it's the total number of unconstrained degrees of freedom at the mesh (the order of the original square matrix after imposition of boundary conditions). – Gustavo Costa Nov 25 '19 at 15:48
• So $n=N$? in your question? – Wolfgang Bangerth Nov 25 '19 at 21:58
• No. N = number of nodes per element (N = 4 for a linear tetrahedra, for example). On the other hand, n = total number of unconstrained degrees of freedom at the mesh (for example, in a thermal problem it would be the nodes with Neumann BC) – Gustavo Costa Nov 25 '19 at 22:12
• John Burkhardt's online repository has good examples doing this (look for fem2D... codes). You collect in pairs interacting dofs, you are creating sparse triplet or ij or coordinate storage format sparsity pattern. Afterwards sort it with efficient sorting algorithm (heapsort..) then you can 'condensate' and easily create CRS with columns (length nnz) and offset (length nrows+1) arrays. – Johntra Volta Nov 26 '19 at 9:55

After struggling a little bit, I developed my own code for the CSR vectors evaluation on Fortran. I believe to be very efficient.

! COMPRESSED STORAGE ROW SUBROUTINE (CSR)

! PROVIDED DATA:
! conec: elements connectivity array. conec(elm,1) is the Element label and
! conec(elm,2) to conec(elm,ndf*nnoelm) are the nodes labels
! nelem: number of elements at the mesh
! nnode: number of nodes at the mesh
! nnoelm: number of nodes per element
! ndf: number of degrees of freedom per node
! ndof: total number of unconstrained (free) degrees of freedom
! vdf: Explained below
! For the following vector, degrees of freedom 1 and 4 are constrained (equal to 0)
! and degrees of freedom 2 and 3 are unconstrained. The unconstrained must be sorted,
! thus vdf(2) = 2 and vdf(3) = 1 would be wrong.
! vdf(1) = 0
! vdf(2) = 1
! vdf(3) = 2
! vdf(4) = 0

! OBTAINED DATA:
! values: non-zero values vector
! ja: column vector
! ia: row vector
! nnz: number of non-zero entries
! nlin: ndof + 1

subroutine csr(conec,nelem,nnode,nnoelm,ndof,ndf,vdf,values,ja,ia,nnz,nlin)
implicit none
integer*4 nnode,nelem,ndf,nnoelm,ndof
integer*4 a,i,j,k,m,n,s,elm,nnz,nlin,conec(nelem,1+nnoelm),vdf(nnode*ndf)
integer*4 temp(ndf*nnoelm)
integer*4,dimension(:),allocatable::vec,list1,ja,ia
double precision,dimension(:),allocatable::values

if (allocated(values)) then
deallocate(values)
end if
if (allocated(ja)) then
deallocate(ja)
end if
if (allocated(ia)) then
deallocate(ia)
end if
if (allocated(vec)) then
deallocate(vec)
end if

allocate(vec(ndof+1))

do i = 1, ndof+1
vec(i) = 0
end do

do elm = 1, nelem
s = 0
do i = 1, nnoelm
do j = 1, ndf
s = s + 1
temp(s) = ndf*conec(elm,i+1)-(ndf-j)
end do
end do
do i = 1, ndf*nnoelm
if (vdf(temp(i)) .ne. 0) then
s = 0
do j = 1, ndf*nnoelm
s = s + 3*vdf(temp(j))/(3*vdf(temp(j))-1)
end do
vec(vdf(temp(i))) = vec(vdf(temp(i))) + s
end if
end do
end do

s = 0
do i = 1, ndof
s = s + vec(i)
end do
vec(ndof+1) = s+1

do i = ndof, 1, -1
vec(i) = vec(i+1) - vec(i)
end do

allocate(list1(vec(ndof+1)-1))

do elm = 1, nelem
s = 0
do i = 1, nnoelm
do j = 1, ndf
s = s + 1
temp(s) = ndf*conec(elm,i+1)-(ndf-j)
end do
end do
do i = 1, ndf*nnoelm
s = 0
do j = 1, ndf*nnoelm
if (vdf(temp(i))*vdf(temp(j)) .ne. 0) then
s = s + 3*vdf(temp(j))/(3*vdf(temp(j))-1)
list1(((3*vdf(temp(i))/(3*vdf(temp(i))-1))*(3*vdf(temp(j))/(3*vdf(temp(j))-1)))*(vec(vdf(temp(i)))+s-1)) = temp(j)
end if
end do
vec(vdf(temp(i))) = vec(vdf(temp(i))) + s
end do
end do

do i = ndof, 2, -1
vec(i) = vec(i-1)
end do
vec(1) = 1

do i = 1, ndof
do j = vec(i)+1, vec(i+1)-1
a = list1(j)
do m = j-1, vec(i), -1
if (list1(m) .le. a) goto 10
list1(m+1) = list1(m)
end do
m = vec(i)-1
10 continue
list1(m+1) = a
end do
end do

nnz = 0
do i = 1, ndof
nnz = nnz+1
do j = vec(i)+1, vec(i+1)-1
nnz = nnz + 3*(list1(j)-list1(j-1))/(3*(list1(j)-list1(j-1))-1)
end do
end do
nlin = ndof+1

allocate(values(nnz),ja(nnz),ia(nlin))
do i = 1, nnz
values(i) = 0.d0
end do

s = 0
k = 0
do i = 1, ndof
s = s+1
k = k + 1
ja(s) = vdf(list1(k))
ia(i) = s
do j = vec(i)+1, vec(i+1)-1
k = k + 1
s = s + 3*(list1(j)-list1(j-1))/(3*(list1(j)-list1(j-1))-1)
ja(s) = vdf(list1(k))
end do
end do
ia(ndof+1) = s+1

end subroutine csr