Mesh Mess

FEM, CFD, CAE, and all that.

Tag: Fortran

Size of deallocated pointer in Fortran

Something I just found out. A deallocated pointer retains its size in terms of the return  of the SIZE function. Specifically, if the pointer is first allocated with size n and later deallocated, SIZE function still returns n. I’ve checked this on GNU and PGI Fortran and got same results, but not sure whether this is compiler dependent. Check standards later.


Store a polynomial in non-symbolic language

Sometimes one wants to store a polynomial in a non-symbolic language, such as Fotran or C. A easy way to achieve this is to take advantage of the isomorphism between polynomial space P(n) and space R^{n+1}. Namely, we can build an array with entries as the coefficients of the polynomial to be stored. Furthermore, considering these coefficients are just the coordinates of the polynomial in the space spanned by the canonical basis 1,x,x^2,\dots, we can express the polynomial using any basis expansion.

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(:)


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 ?
write(*,*) b1                 ! b1 deallocated along with a ?
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


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