Discussion:
Code for sparse stiffness matrix assembly
(too old to reply)
rusi_pathan
2010-09-04 03:46:08 UTC
Permalink
I _need_ to explicitly assemble the sparse stiffness matrix resulting
from finite element discretization.

Does anyone know of any Fortran code/subroutine which can do this
using linked lists etc?

Any sparse format like CSR, COO would do.

Thanks in advance.
rusi_pathan
2010-09-04 03:55:41 UTC
Permalink
Post by rusi_pathan
I _need_ to explicitly assemble the sparse stiffness matrix resulting
from finite element discretization.
Does anyone know of any Fortran code/subroutine which can do this
using linked lists etc?
Any sparse format like CSR, COO would do.
Thanks in advance.
Just want to clarify that I already have the element stiffness matrix
and simply want to assemble the global sparse matrix in the
appropriate format.
dumashu
2010-09-04 07:57:08 UTC
Permalink
Post by rusi_pathan
Post by rusi_pathan
I _need_ to explicitly assemble the sparse stiffness matrix resulting
from finite element discretization.
Does anyone know of any Fortran code/subroutine which can do this
using linked lists etc?
Any sparse format like CSR, COO would do.
Thanks in advance.
Just want to clarify that I already have the element stiffness matrix
and simply want to assemble the global sparse matrix in the
appropriate format.
first, you should get the location that the element of the matrix in
global sparse matrix

if you use CSR format,there are 3 array , val(*),ia(*),ja(*)
now i show that:
ek:element stiffness matrix
MK:global sparse matrix
we assemble ek(i,j) to MK
id1 = l_g(i) !l_g store the degree of freedom
id2 = l_g(j)
id2 is between ja(ia(id1)) and ja(ia(id1+1)-1)
get the length id2 to ia(id1),
so the location of ek(i,j) is in val(ia(id1) + length)

there is the code.
sorry for my poor english!

subroutine assemble_ek(ek,l_g)
use usertype,only : wp
use geometry_mod,only : MK=>stiffmatrix
implicit none
INTEGER :: i,j,l,m,lm,kk,tt,len1
real(wp),INTENT(IN) :: ek(:,:)
INTEGER,INTENT(IN) :: l_g(:)
! integer,external :: BINARY_SEARCH
integer :: isize
isize = size(ek,dim=1)
do i=1,isize
m = l_g(i)
if(m<=0) cycle
tt = MK%ia(m); len1 = MK%ia(m+1) - tt
do j=1,isize
l = l_g(j)
if(l<=0) cycle
if(l < m) cycle !only the upper matrix
!"BINARY_SEARCH(MK%ja(tt), len1, l) - 1" is the length
kk = BINARY_SEARCH(MK%ja(tt), len1, l) - 1 + tt
MK%val(kk) = MK%val(kk) + ek(j,i)
end do
end do
end subroutine
rusi_pathan
2010-09-07 16:17:00 UTC
Permalink
Post by dumashu
Post by rusi_pathan
Post by rusi_pathan
I _need_ to explicitly assemble the sparse stiffness matrix resulting
from finite element discretization.
Does anyone know of any Fortran code/subroutine which can do this
using linked lists etc?
Any sparse format like CSR, COO would do.
Thanks in advance.
Just want to clarify that I already have the element stiffness matrix
and simply want to assemble the global sparse matrix in the
appropriate format.
first, you should get the location that the element of the matrix in
global sparse matrix
if you use CSR format,there are 3 array , val(*),ia(*),ja(*)
ek:element stiffness matrix
MK:global sparse matrix
we assemble ek(i,j) to MK
id1 = l_g(i) !l_g store the degree of freedom
id2 = l_g(j)
id2 is between ja(ia(id1)) and ja(ia(id1+1)-1)
get the length id2 to ia(id1),
so the location of ek(i,j) is in val(ia(id1) + length)
there is the code.
sorry for my poor english!
subroutine assemble_ek(ek,l_g)
     use usertype,only : wp
     use geometry_mod,only : MK=>stiffmatrix
     implicit none
     INTEGER :: i,j,l,m,lm,kk,tt,len1
     real(wp),INTENT(IN) :: ek(:,:)
     INTEGER,INTENT(IN)  :: l_g(:)
     !    integer,external    :: BINARY_SEARCH
     integer             :: isize
     isize = size(ek,dim=1)
     do i=1,isize
         m = l_g(i)
         if(m<=0) cycle
         tt = MK%ia(m); len1 = MK%ia(m+1) - tt
         do j=1,isize
             l = l_g(j)
             if(l<=0) cycle
             if(l < m) cycle !only the upper matrix
!"BINARY_SEARCH(MK%ja(tt), len1, l) - 1" is the length
             kk = BINARY_SEARCH(MK%ja(tt), len1, l) - 1 + tt
             MK%val(kk) = MK%val(kk) + ek(j,i)
         end do
     end do
end subroutine
Dumashu:
Thanks for your note but how are you calculating the size of the "val"
array

Jared:
Yes I am aware of PETSc and Sparskit. There is even an example in
Sparskit but I need to know the number of nonzeros in the matrix first.
Jared Ahern
2010-09-07 19:09:24 UTC
Permalink
... but I need to know the number of nonzeros in the matrix first.
Two points:

1) I assume that you are assembling the stiffness matrix on an element-
by-element basis; inserting each element matrix into the global matrix
in a loop. If so, you would like to directly insert into a global
sparse matrix, correct? (Since the alternative would be to make a
full dense matrix and then copy it into a sparse duplicate...) This
was not entirely clear in your original post.

2) Your global stiffness matrices likely contains rows that correspond
to nodal equations (possibly displacement or velocity) or elemental
equations (possibly pressure). For nodal equations, the number of
nonzeros often corresponds to the number of variables employed,
multiplied by the number of nodes or element points connected to the
node in question. It is similar for non-nodal equations.

For example, if you had a 2D mesh consisting of 3-node triangular
elements, to get the number of nonzeros for an equation at node "A",
you could multiply the number of nodal variables per node (say 2: u_x
and u_y) by the total number of unique nodes contained by elements
which themselves contain "A". It is possible to calculate this
exactly for node "A", but it may be more efficient just to estimate
that you will have a max of, say, 10 elements intersecting at any
node, so the max number of nonzeros per row would be approximately
2*(10*(3-2)+1)=22. Note that the (3-2) term accounts for nodes
duplicated by adjacent elements. This assumption could be false, as
it is mesh-dependent, but you could just estimate slightly more, too.

Remember, you can always set the number of nonzeros in each row to to
be up to the length of the entire row itself; you just end up with a
"sparse" matrix that is much slower (and larger) than a simple dense
rank-2 array. So eliminate as many as you can, but you can generally
have some extra zero entries floating around without too much cause
for concern.
rusi_pathan
2010-09-07 19:47:36 UTC
Permalink
Thanks for the reply
Post by Jared Ahern
... but I need to know the number of nonzeros in the matrix first.
1) I assume that you are assembling the stiffness matrix on an element-
by-element basis; inserting each element matrix into the global matrix
in a loop.  If so, you would like to directly insert into a global
sparse matrix, correct?  (Since the alternative would be to make a
full dense matrix and then copy it into a sparse duplicate...)  This
was not entirely clear in your original post.
Yes that is correct.
Post by Jared Ahern
It is possible to calculate this
exactly for node "A", but it may be more efficient just to estimate
that you will have a max of, say, 10 elements intersecting at any
node, so the max number of nonzeros per row would be approximately
2*(10*(3-2)+1)=22. Note that the (3-2) term accounts for nodes
duplicated by adjacent elements.  This assumption could be false, as
it is mesh-dependent, but you could just estimate slightly more, too.
I dont think this is a good approach because
- Memory is wasted and for large 3D meshes it can be substantial
- There really is no upper limit on the number of elements a node can
be in (specially in unstructured meshes). Think of a circle with 1000
sides and all nodes directly connected to the center node giving you
1000 triangular elements (though aspect ratios will be bad).

I think finding the exact number of nnzs is also recommended in the
PETSc manual.

Since this is such a common operation I was hoping that there are many
such routines out there but I cant seem to find any.
Jared Ahern
2010-09-07 21:45:25 UTC
Permalink
Post by rusi_pathan
I dont think this is a good approach because
- Memory is wasted and for large 3D meshes it can be substantial
- There really is no upper limit on the number of elements a node can
be in (specially in unstructured meshes). Think of a circle with 1000
sides and all nodes directly connected to the center node giving you
1000 triangular elements (though aspect ratios will be bad).
Well, as I mentioned, you can calculate the number of nonzeros exactly
for each node if you desire. Just loop through the nodes; for each
node, loop over each element and append the node numbers to a running
list. This approach uses little memory, but a relatively large amount
of computation. You could alternatively create a logical (or bit)
matrix (or sparse matrix), and loop over the elements marking off the
entries used. Then just count up the true values for each row at the
end. This takes more memory, but little processing. There are other
approaches, but these are the simple methods that come to mind. If
you generate the mesh (not read in from an external program), you may
be able to couple the generation to the nonzero calculation.

And yes, there would be some wasted space for the simple approach, but
if you end up with 100+ elements touching a corner node you'll
probably have major mesh problems anyways.
rusi_pathan
2010-09-08 08:33:41 UTC
Permalink
Post by Jared Ahern
Post by rusi_pathan
I dont think this is a good approach because
- Memory is wasted and for large 3D meshes it can be substantial
- There really is no upper limit on the number of elements a node can
be in (specially in unstructured meshes). Think of a circle with 1000
sides and all nodes directly connected to the center node giving you
1000 triangular elements (though aspect ratios will be bad).
Well, as I mentioned, you can calculate the number of nonzeros exactly
for each node if you desire. Just loop through the nodes; for each
node, loop over each element and append the node numbers to a running
list.  This approach uses little memory, but a relatively large amount
of computation. You could alternatively create a logical (or bit)
matrix (or sparse matrix), and loop over the elements marking off the
entries used.  Then just count up the true values for each row at the
end.  This takes more memory, but little processing.  There are other
approaches, but these are the simple methods that come to mind.  If
you generate the mesh (not read in from an external program), you may
be able to couple the generation to the nonzero calculation.
Thanks for your note. I think I will have to write my own routine
using a running list as I cant find much on this online.
Post by Jared Ahern
And yes, there would be some wasted space for the simple approach, but
if you end up with 100+ elements touching a corner node you'll
probably have major mesh problems anyways.
Actually this code will also be used by other people so I cannot make
any assumptions about the element quality.

Jared Ahern
2010-09-04 13:27:55 UTC
Permalink
Post by rusi_pathan
Does anyone know of any Fortran code/subroutine which can do this
using linked lists etc?
I've been using PETSc to do this same task for a while now. It is not
a native Fortran library (C I believe), but there are Fortran
bindings. It can also be used with MPI.

http://www.mcs.anl.gov/petsc/petsc-as/

SPARSKIT would be another option.

http://www-users.cs.umn.edu/~saad/software/SPARSKIT/sparskit.html
Loading...