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.

On Sep 3, 10:46=A0pm, rusi_pathan <tabrez...@gmail.com> wrote:
> 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.

于 2010/9/4 11:55, rusi_pathan 写道:
> On Sep 3, 10:46 pm, rusi_pathan<tabrez...@gmail.com>  wrote:
>> 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

> 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

On Sep 4, 2:57=C2=A0am, dumashu <duma...@gmail.com> wrote:
> =E4=BA=8E 2010/9/4 11:55, rusi_pathan =E5=86=99=E9=81=93:> On Sep 3, 10:4=
6 pm, rusi_pathan<tabrez...@gmail.com> =C2=A0wrote:
> >> 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=EF=BC=8Cthere 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 =3D l_g(i) !l_g store the degree of freedom
> id2 =3D 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)
> =C2=A0 =C2=A0 =C2=A0use usertype,only : wp
> =C2=A0 =C2=A0 =C2=A0use geometry_mod,only : MK=3D>stiffmatrix
> =C2=A0 =C2=A0 =C2=A0implicit none
> =C2=A0 =C2=A0 =C2=A0INTEGER :: i,j,l,m,lm,kk,tt,len1
> =C2=A0 =C2=A0 =C2=A0real(wp),INTENT(IN) :: ek(:,:)
> =C2=A0 =C2=A0 =C2=A0INTEGER,INTENT(IN) =C2=A0:: l_g(:)
> =C2=A0 =C2=A0 =C2=A0! =C2=A0 =C2=A0integer,external =C2=A0 =C2=A0:: BINAR=
Y_SEARCH
> =C2=A0 =C2=A0 =C2=A0integer =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 :: =
isize
> =C2=A0 =C2=A0 =C2=A0isize =3D size(ek,dim=3D1)
> =C2=A0 =C2=A0 =C2=A0do i=3D1,isize
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0m =3D l_g(i)
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if(m<=3D0) cycle
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0tt =3D MK%ia(m); len1 =3D MK%ia(m+1) - =
tt
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0do j=3D1,isize
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0l =3D l_g(j)
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if(l<=3D0) cycle
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if(l < m) cycle !only the=
 upper matrix
> !"BINARY_SEARCH(MK%ja(tt), len1, l) - 1" is the length
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0kk =3D BINARY_SEARCH(MK%j=
a(tt), len1, l) - 1 + tt
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0MK%val(kk) =3D MK%val(kk)=
 + ek(j,i)
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0end do
> =C2=A0 =C2=A0 =C2=A0end 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.

> ... 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.

Thanks for the reply

On Sep 7, 2:09=A0pm, Jared Ahern <jared.ah...@gmail.com> wrote:
> > ... 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. =A0If so, you would like to directly insert into a global
> sparse matrix, correct? =A0(Since the alternative would be to make a
> full dense matrix and then copy it into a sparse duplicate...) =A0This
> was not entirely clear in your original post.
Yes that is correct.

> 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)=3D22. Note that the (3-2) term accounts for nodes
> duplicated by adjacent elements. =A0This 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.

> 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.

On Sep 7, 4:45=A0pm, Jared Ahern <jared.ah...@gmail.com> wrote:
> > 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. =A0This 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. =A0Then just count up the true values for each row at the
> end. =A0This takes more memory, but little processing. =A0There are other
> approaches, but these are the simple methods that come to mind. =A0If
> 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.
>
> 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.