[petsc-users] Preallocation woes again

[Tabrez Ali]

PETSc Gurus

First I want to thank you for patiently answering questions that I have 
asked in past regarding preallocation. Unfortunately I am still having 
problems.

I have a small unstructured FE (elasticity) code that uses PETSc.

Unfortunately I am not yet able to find an _efficient_ way of 
calculating the non-zero structure so I simply overestimate the o_nz and 
d_nz values (= dimension x max_number_of nodes _a_node_contributes_to) 
which means that my stiffness matrix consumes at least 2X more memory. 
For example, for a perfectly structured mesh of linear hexes this number 
is 3*27=81=o_nz=d_nz. In complicated 3D linear tet/hex meshes that I 
have generated I rarely need to set a value greater than 150. In general 
I have found that as long as I keep DOF/core between 100-200K (assuming 
1GB/per core) there is enough local RAM left even when memory for 
stiffness matrix is overestimated by 2X-6X.

In any case now I do want to preallocate exactly for better memory 
performance and the ambiguity involved in choosing a reasonable 
o_nz/d_nz. The way I am trying to do it involves loops like ...

do i=1, num_local_elements
    do j=1, num_total_nodes
       ...
    end do
end do

or using lists that involve searches.

This unfortunately takes much time (much more than assembly and solve) 
due to the second loop. I am aware of many posts/slides by PETSc authors 
that mention that the non-zero structure can be found simply by looping 
once through elements but what to do next is not quite clear.

I am also aware of this post by Barry [ 
http://lists.mcs.anl.gov/pipermail/petsc-users/2008-May/003020.html ] 
but I cannot get it to work for a simple 4 element problem. Can some one 
please expand on the strategy being suggested by Barry specially where 
the loops are discussed (I understand the arrays and the vecscatter 
part). I am also not sure how does using a t of 0.5 and 1 prevents 
double counting. Shouldn't it be 0 and 1?

I am also aware of some capabilities in DMMESH but right now I dont 
understand it well enough to utilize it.

Alternatively does PETSc support MATMPIAIJ of size=integer(1)/logical 
instead of real(8) which I can use to put zeros/ones (when I loop over 
elements the first time) to get the non-zero structure?

Thanks in advance.

Tabrez