Common block in Domain Decomposition implementation

Fortran 77 essentially does not have global variables for subroutines, instead, common blocks are used to passing data across subroutines when arguments are not preferred. This is conceptually unsafe, though widely used in some legacy F77 codes. The code I’ve been working on adopts this “feature” extensively. The problem is, in Dirichlet-Neumann-like domain decomposition implementation, the solver will be called twice for each iteration, each time with different stiffness matrix setup(due to different subdomains) and boundary conditions. On the other hand, the common blocks make those two calls undistinguishable.

To reuse the legacy solver, my first try was to construct two copies of the solver by putting it into two modules, hoping the explicit interface of the modules would identify the common blocks in different copies. This turns out to be wrong: the subroutines from the two modules still share the same common blocks.

What I finally decided to do is to build a copy of the original solver by renaming all the common blocks. Sed & AWK make this  not to difficult. But still I don’t think that’s the most efficient way to tackle this. Consider this is a lesson that not all sequential solvers are easy to parallelize.

A subtle point in Fortran 90 pointer

Though I found out it today, what’s discussed below is already pointed out at here.

First, a few toy lines of f90 ‘s pointer feature.

program main

integer, ALLOCATABLE, TARGET :: a(:), c(:)
integer, POINTER :: b1(:)

allocate(a(3))
allocate(c(3))
allocate(b1(3))

a = [3,4,893]
c = [4,5,90]
b1 = a
write(*,*) b1                 ! b1 pointed to a
a = c
write(*,*) b1                 ! b1 pointed to a or c ?
deallocate(a)
write(*,*) b1                 ! b1 deallocated along with a ?
deallocate(c)
deallocate(b1)                ! b1 deallocated with a
write(*,*) b1

end program main

The output is:

% a.out
           3           4         893
           3           4         893
           3           4         893

Explanation:

The first and last output are trivial: just the content of b1 and null (blank line).  What’s interesting is the 2nd and 3rd output. At the 2nd output, “a” is reassigned with “c”, but the pointer “b1″ still gives “a” ‘s initial value, the same happens at 3rd output, when “a” is already deallocated. This is what pointed out by the link mentioned above: assignment operator gives an temporary copy, which persists even after the original version is destroyed (DEALLOCATED). This gives a subtle way of assigning pointers: by updating the same array and using assignment operator, different values of the array could be assigned to different pointers.

From F77 to F95/2003

No, this is not a metaphor you found on tutorial texts of F95 or F2003. This is literally the real thing I am working on now. In a middle of a project trying to coupling two solvers, I am working on one of them, whose origin gives it considerable amount of unsecured variables due to F77 tradition.  What I am trying to do is to modulize it and pack & tag those variables, functions, subroutines. And, during this process, indeed some naming pitfalls showed up, apparently not caught during past coding or running. Even though somehow they do not contribute any malfunction of the solver, it indeed shows me why people make those rules of new versions of Fortran. However, believe me, it’s not the best job in the world….

FYI, what I am doing is absolutely more than this:

via Terece Tao’s GBuzz

Pitfalls in fortran

I am glad to run into this page because some of those mistakes were honorably made. In particular, the variable passing with array is always a place to make things go south in high performance computing, and I found that module is usually the easy way to go, assuming assumed-shape arrays.