Efficiently searching arrays in FORTRAN

Question

I am trying to store the stiffness matrix in FORTRAN in sparse format to save memory, i.e. I am using three vectors of non-zero elements (irows, icols, A). After finding out the size of these arrays the next step is to insert the values in them. So I am using gauss points, i.e. for each gauss point I am going to find out the local stiffness matrix and then insert this local stiffness matrix in the Global (irows, icols, A) one.

The main problem with this insertion is that every time we have to check that either the new value is exists in the global array or not, so if the value exists add the new to the old but if not append to the end. i.e. we have to search the whole array to find that either the value exists or not. If the size of these arrays (irows, icols, A) are large so this search is computationally very expensive.

Can any one suggest a better way of insertion of the local stiffness matrix for each gauss point the global stiffness matrix.

Answer 1

I am fairly sure that this is a well known problem in FEM analysis - I found reference to it in this scipy documentation, but of course the principals are language independent. Basically what you should do is create your matrix in the format you have, but instead of searching the matrix to see whether an entry already exists, just assume that it doesn't. This means that you will end up with duplicate entries which need to be added together to get the correct value.

Once you have constructed your matrix, you will usually convert it to some more efficient form for solving it (e.g. CSR etc.) - the exact format may be determined by the sparse solver you are using. During this conversion process duplicate entries should get added together - and some sparse matrix libraries will do this for you. I know that scipy does this, and many of its internal routines are written in fortran, so you might be able to use one of them (they are all open source). Or you could check if anything suitable is on netlib.

Answer 2

If you use a data structure that is pre-sorted it would be very efficient to search it. Either as your primary data structure or as an auxiliary data structure. You want one that you can insert another entry into the middle. For example, a binary search tree (http://en.wikipedia.org/wiki/Binary_search_tree).