1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ What is $n$? The number of nodes per element? $\endgroup$ – Wolfgang Bangerth Nov 25 '19 at 15:39
  • $\begingroup$ 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). $\endgroup$ – Gustavo Costa Nov 25 '19 at 15:48
  • $\begingroup$ So $n=N$? in your question? $\endgroup$ – Wolfgang Bangerth Nov 25 '19 at 21:58
  • $\begingroup$ 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) $\endgroup$ – Gustavo Costa Nov 25 '19 at 22:12
  • 1
    $\begingroup$ 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. $\endgroup$ – Johntra Volta Nov 26 '19 at 9:55
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.