[petsc-users] Writing a domain decomposition code with PETSc

[Åsmund Ervik]

Jed,

On ma. 07. okt. 2013 kl. 17.14 +0200, Jed Brown wrote:
> Åsmund Ervik <asmund.ervik at ntnu.no> writes:
>>   subroutine transform_petsc_us(array,f,dof)
>>     !Note: assumed shape-array does the "coordinate transformation"
>>     real, intent(in), dimension(dof,1-stw:,1-stw:,1-stw:) :: array
>>     real,intent(inout),dimension(1-stw:,1-stw:1-stw:,dof) :: f
>
> This is a really bad ordering for cache-based architectures because as
> you traverse the mesh, you need to pull values from segments of the
> array that are far apart.  We have yet to encounter a practical project
> in which indexing the first way (dofs adjacent in memory) is not faster,
> sometimes by large factors.  (For example, Clawpack was recently
> convinced to change and saw a 50% speedup.  The PETSc-FUN3D papers had
> many comparisons.)  I recommend transforming your code to use the former
> ordering, or if you must, use CPP macros so that you can change it at
> compile time.

Thanks for this comment; I have been suspecting this is a bad idea, 
nice to get confirmation. I will follow your advice.

Other than that, does the mapping from local to global using Fortran 
assumed-shape arrays make sense, or is there a better way of doing it?

Åsmund