As an illustration, we will consider a program which calculates properties of a cube. Some of the properties depend on the edge length of the cube, some others depend on the Cartesian coordinates of the vertices.
The Implicit Reference to Parameters method
In scientific programming, a program can almost always be seen as a pure function of its data:
output = program(input)
In this functional view, a program can be represented as an acyclic graph, where:
- the vertices of the graph are the entities of interest
- the edges of the graph represent the relation {needs/needed by}
We call this graph the production tree.
Consider the example of the electronic energy of a molecular system:
$$ E = E_{\text{nucl}} + E_{\text{elec}} $$
decomposed as the nuclear repulsion energy and the electronic energy. The electronic energy can be expressed as the sum of the potential and kinetic electronic energies:
$$ E_{\text{elec}} = E_{\text{pot}} + E_{\text{kin}} $$
The production tree of E is represented as:
Usual programming
The calculation of E could be done in Fortran using subroutine calls to imperatively ask the program to realize the needed calculations, such as
program compute_Energy
double precision :: E_nucl, E_pot, E_kin, E_elec, E
call compute_E_nucl(E_nucl)
call compute_E_pot(E_pot)
call compute_E_kin(E_kin)
call compute_E_elec(E_elec,E_pot,E_kin)
call compute_E(E,E_nucl,E_elec)
print *, 'Energy = ', E
end
In this example, it is not clear which are input and output arguments in the subroutine calls. A clearer way to express the same code would be using functions:
program compute_Energy
double precision :: E_nucl, E_pot, E_kin, E_elec, E
double precision :: compute_E_nucl, compute_E_pot, compute_E_kin, &
compute_E_elec, compute_E
E_nucl = compute_E_nucl()
E_pot = compute_E_pot()
E_kin = compute_E_kin()
E_elec = compute_E_elec(E_pot,E_kin)
E = compute_E(E_nucl,E_elec)
print *, 'Energy = ', E
end
In both example, the programmer needs to know the production tree of E, in order to be position the statements in the correct order. For instance, the line
E_nucl = compute_E_nucl()
has to be positioned before
E = compute_E(E_nucl,E_elec)
otherwise the program is wrong.
If the code is written in this way:
program compute_Energy
double precision :: compute_E_nucl, compute_E_pot, compute_E_kin, &
compute_E_elec, compute_E
print *, 'Energy = ', compute_E (compute_E_nucl(), &
compute_E_elec (compute_E_pot(), compute_E_kin()) )
end
the production tree is expressed in the list of arguments, and the programmer
does not need to worry about the sequence of the function calls.
The call sequence is now handled by the compiler. Some compilers will call
compute_E_nucl
first, some others will call compute_E_elec
first.
However, for a real code, this practice is impossible because
the code would be unreadable by a programmer.
IRPF90 programming
If the parameters of function calls could be automatically inserted by a tool, the program would be as simple as
program compute_E
print *, 'Energy = ', E
end
The role of the IRPF90 tool is to allow programmers to write code with this style. IRPF90 is a program which generates Fortran90 code, from a language which is an extension of Fortran. In IRPF90, the concept of ‘‘entity of interest’’ or ‘‘IRP entity’’ is introduced. An IRP entity is a node of the production tree. It is defined by a provider function, whose role is to build the value of the entity. In this example, the provider function of E would be
BEGIN_PROVIDER [ double precision, E ]
E = E_nucl + E_elec
END_PROVIDER
where E_nucl
and E_elec
are also IRP entities. When an entity is used, if
the entity has already been computed, the value of the last evaluation is
returned. Otherwise, the code in the provider function is executed, and the
result is returned. However, an entity can be invalidated (by modifying the
value of a needed entity, for example), it will be recomputed. This mechanism
guarantees that the same quantity is never computed more than once, and that
when it is used it is always valid.
The main benefit of using IRPF90 is that the programmer never worries about the calling sequence of the code. As soon as he uses an entity, it is guaranteed that this entity is computed and valid. At first sight, with simple examples it is difficult to realize to what extent IRP programming makes things simpler. For a real code with thousands (millions) of lines, the code written with IRPF90 is as easy to read and modify as a code with a few hundreds of lines. Most of the complexity of a large code is now handled by the computer, and not the programmer.
A simple example
Create a new directory named cube for example. Go into this directory and run the command
$ irpf90 --init
Two temporary directories are created:
$ ls
IRPF90_man IRPF90_temp Makefile
and a standard Makefile
is build, using the gfortran compiler:
IRPF90 = irpf90 #-a -d
FC = gfortran
FCFLAGS= -ffree-line-length-none -O2
SRC=
OBJ=
LIB=
include irpf90.make
irpf90.make: $(wildcard *.irp.f)
$(IRPF90)
We can now start to write the code.
The main program
First, we write a very simple program which prints the surface and the
volume of a cube. In a file named cube_example.irp.f
, write:
program cube_example
implicit none
print *, 'Surface Area :', surface
print *, 'Volume :', volume
end
Remark that there is no explicit directive to run a computation. Hence, this code is clear as it expresses the intention of the programmer with a very epurated style : the goal of the main program is to print the surface and the volume. The way these quantities are computed should not appear here.
The properties of the cube
The IRPF90 environment guarantees that when an entity is used, it is valid.
In the main program, the action to print the variable surface
will
automatically request the value of the surface before the print
statement : the value of the entity surface
will be provided.
We can now write in the file named properties.irp.f
how the properties are
computed. The calculation of a property appears in the provider function
of the property.
BEGIN_PROVIDER [ real, surface ]
BEGIN_DOC
! Surface of the cube
END_DOC
surface = 6. * edge**2
END_PROVIDER
BEGIN_PROVIDER [ real, volume ]
BEGIN_DOC
! Volume of the cube
END_DOC
volume = edge**3
END_PROVIDER
These two properties depend on the value of the edge
variable which can be
given in the standard input. The edge is the way we define a cube, so we
write in a file named cube.irp.f
:
BEGIN_PROVIDER [ real, edge ]
BEGIN_DOC
! Value of the edge of the cube
END_DOC
print *, "Value of the edge of the cube"
read(*,*) edge
ASSERT (edge > 0.)
END_PROVIDER
Here the ASSERT
keyword guarantees that when edge
is provided, its value
is positive. Otherwise, if the -a
option of the irpf90
command is used,
the program will fail will an error message:
Stack trace: 0
-------------------------
cube_example
provide_center
provide_vertex
provide_edge
bld_edge
-------------------------
bld_edge: Assert failed:
file: cube.irp.f, line: 24
(edge > 0.)
edge = -1.000000
The assertion ensures that the value of edge
will always be positive
everywhere in the code.
The BEGIN_DOC ... END_DOC
groups contain the documentation of the IRP entities,
which helps to write dynamically the technical documentation of the program.
The documentation of the IRP entities can be viewed through man pages
using the irpman tool (shown in next section).
Code compilation
To understand what happens at the execution, turn on the -a
and -d
options
of irpf90
in the Makefile
by removing the #
in the first line. Then, run make
.
$ make
Makefile:9: irpf90.make: No such file or directory
irpf90 -a -d
IRPF90_temp/cube.irp.F90
IRPF90_temp/properties.irp.F90
IRPF90_temp/cube_example.irp.F90
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/cube.irp.F90 -o IRPF90_temp/cube.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/properties.irp.F90 -o IRPF90_temp/properties.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/cube_example.irp.F90 -o IRPF90_temp/cube_example.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/irp_stack.irp.F90 -o IRPF90_temp/irp_stack.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/irp_touches.irp.F90 -o IRPF90_temp/irp_touches.irp.o
gfortran -o cube_example IRPF90_temp/cube_example.irp.o IRPF90_temp/irp_stack.irp.o IRPF90_temp/cube.irp.o
IRPF90_temp/properties.irp.o IRPF90_temp/irp_touches.irp.o
Many files have been created:
$ ls
cube_example cube.irp.f irpf90.make IRPF90_temp properties.irp.f
cube_example.irp.f irpf90_entities IRPF90_man Makefile
For each *.irp.f
file, a corresponding Fortran90 module has been built. The
file irpf90.make
was generated, and contains the dependencies between the *.irp.f
files. The file irpf90_entities
contains the list of the IRP entities,
their type and the file in which they are defined:
$ cat irpf90_entities
cube.irp.f : real :: edge
properties.irp.f : real :: surface
properties.irp.f : real :: volume
Now, if you run:
$ irpman surface
you obtain a man page describing the entity surface and its dependencies:
IRPF90 entities(l) surface IRPF90 entities(l)
Declaration
real :: surface
Description
Surface of the cube
File
properties.irp.f
Needs
edge
IRPF90 entities surface IRPF90 entities(l)
Code execution
Recall that the -d
and -a
options were activated in the Makefile
.
Run the program, and choose 2 for the value of the edge of the cube:
$ ./cube_example
0 : -> provide_volume
0 : -> provide_edge
0 : -> edge
Value of the edge of the cube
2.
0 : <- edge 0.0000000000000000
0 : <- provide_edge 0.0000000000000000
0 : -> volume
0 : <- volume 0.0000000000000000
0 : <- provide_volume 0.0000000000000000
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.0000000000000000
0 : <- provide_surface 0.0000000000000000
0 : -> cube_example
Surface Area : 24.000000
Volume : 8.0000000
0 : <- cube_example 0.0000000000000000
The debug option -d
prints a lot of output. It corresponds to the
exploration of the production trees of surface
and volume
, which are
needed in the cube_example
program. A right arrow tells us that we enter
inside a subroutine, and the left arrow tells us that we leave the subroutine.
When we leave a subroutine, the CPU time spent in the subroutine is printed
(here, it is always 0 seconds).
In this example, we can see that the print statements of the main program appear at the bottom of the output. Many things happen before the code we wrote.
We first need to provide surface
-> provide_surface
As surface
is not valid, we have to build it, but as surface
needs
edge
, we first need to provide edge
before computing the value of
surface
.
-> provide_edge
As edge
is not valid, it has to be built.
-> bld_edge
Building edge
asks the user to enter the value with the standard
input.
Value of the edge of the cube
Now edge
is valid
<- bld_edge 0.00000000000000
<- provide_edge 0.00000000000000
and we can build the value of surface
:
-> bld_surface
<- bld_surface 0.00000000000000
We are now back into the main program with a valid value for the surface.
The value of volume
is also requested in the main program, so the
exploration of the production tree of volume
starts:
-> provide_volume
volume
needs edge
. As edge
is valid, it is not re-built at its
value is used to calculate volume
-> bld_volume
<- bld_volume 0.00000000000000
<- provide_volume 0.00000000000000
We can now execute the main program with valid values of surface
and
volume
.
-> cube_example
These values can now be used to be printed:
Surface Area : 24.00000
Volume : 8.000000
The last line tells us we leave the program
<- cube_example 0.00000000000000
From this simple example, we can notice that a lot of energy has been saved for the programmer. First, the makefiles have been automatically generated. Then, the sequence of execution of the code is absolutely not controlled by the programmer. The only thing he had to worry about was the correctness of the definition of the entities of interest. Using IRPF90, there are questions that programmers will never ask anymore, for example:
At this particular place of the code I need x
. Is it already calculated?
The programmer will just use x
without worrying if it has
been calculated or not, and he will be sure that its value is valid if the
ASSERT
statements have been properly inserted.
Improvement of the example
Every day, developers are improving already existing codes, which may or may not be written by them. Let’s try to add a new functionality to the code of the previous section: we now want to print the Cartesian coordinates of the center of the faces of a cube. To add this new feature, we have to introduce the coordinates of the vertices of the cube, find the groups of 4 vertices which constitute the faces, and calculate the centers.
The main program
In the main code, just add the information that you want to print the coordinates of the centers of the faces of the cube:
program cube_example
implicit none
print *, 'Surface Area :', surface
print *, 'Volume :', volume
print *, ''
print *, 'Centers of the faces:'
integer :: i
do i=1,face_num
integer :: j
print *, (center(j,i), j=1,3)
end do
end
In this code, the declarations of the integers ‘‘i’’ and ‘‘j’’ have been introduced just before their first use. In IRPF90, the declarations can appear anywhere.
The addition of a new property
Now, we introduce the ‘‘center’’ entity, which contains the values of the
centers of the faces of the cube, in the properties.irp.f
file:
BEGIN_PROVIDER [ real, center, (3,face_num) ]
implicit none
BEGIN_DOC
! Coordinates of the center of the faces cube
END_DOC
integer :: i
do i=1,face_num
integer :: k
do k=1,3
center(k,i) = 0.
integer :: j
do j=1,4
integer :: l
l = face(j,i)
center(k,i) = center(k,i) + vertex(k,l)
enddo
center(k,i) = center(k,i) / 4.
enddo
enddo
END_PROVIDER
The center
entity is an array of dimension (3,face_num)
, where
face_num
is the number of faces in a cube. For each center, the
coordinates are computed as the average of the coordinates of the 4 vertices
constituting the face. The coordinates of the vertices of the cube are present
in the array vertex
. The array face
contains, for each face, the
indices of the array vertex
corresponding to the vertices of the face:
these last arrays are defined in the cube.irp.f
:
BEGIN_PROVIDER [ integer, vertex_num ]
implicit none
BEGIN_DOC
! Number of vertices
END_DOC
vertex_num = 8
END_PROVIDER
BEGIN_PROVIDER [ real, vertex, (3,vertex_num) ]
implicit none
BEGIN_DOC
! Coordinates of the vertices of the cube
END_DOC
integer :: k, i
! Initialize the array
do i=1,vertex_num
do k=1,3
vertex(k,i) = 0.
enddo
enddo
! The 1st point is the origin
! Build the 3 points on the axes
do k=1,3
vertex(k,k+1) = edge
enddo
! Build the 3 points in the xy, yz and zx planes
integer :: knew(3) = (/2, 3, 1 /)
do k=1,3
vertex(k,k+4) = edge
vertex(knew(k),k+4) = edge
enddo
! The last point
do k=1,3
vertex(k,8) = edge
enddo
END_PROVIDER
For simplicity, the cube is chosen to have one point at the origin, and faces in the xy, yz and xz planes. The faces are computed using the squared distance matrix of the vertices:
BEGIN_PROVIDER [ integer, face_num ]
implicit none
BEGIN_DOC
! Number of face of a cube
END_DOC
face_num = 6
END_PROVIDER
BEGIN_PROVIDER [ integer, face, (4,face_num) ]
implicit none
BEGIN_DOC
! Indices of the vertices for each face of the cube
END_DOC
integer :: i, ii
i=1
do ii = 0,1
! Pick the 1st point and find the 3 points at 'edge' distance of it
integer :: j, ifound(3), inext
inext = 1
do j=1,vertex_num
if (distance2(j,i) == edge2) then
ifound ( inext ) = j
inext = inext + 1
endif
enddo
ASSERT ( inext == 4 )
integer :: istart
istart = 3*ii
integer :: k, knew(3)
knew(:) = (/ 2, 3, 1 /)
do k=1,3
face(1,k+istart) = i
face(2,k+istart) = ifound(k)
face(3,knew(k)+istart) = ifound(k)
enddo
! Find the 4th point of those 3 faces
do j=1,vertex_num
if (distance2(j,i) == 2.*edge2) then
do k=1, 3
if ( (distance2(j,face(2,k+istart)) == edge2) .and. &
(distance2(j,face(3,k+istart)) == edge2) ) then
face(4,k+istart) = j
endif
enddo
endif
enddo
! Find the point opposite to the point i for 2nd iteration
integer :: inew
do j=1,vertex_num
if (distance2(j,i) == 3.*edge2) then
inew = j
endif
enddo
i = inew
end do
END_PROVIDER
where the squared distance matrix is defined with:
BEGIN_PROVIDER [ real, distance2, (vertex_num,vertex_num) ]
implicit none
BEGIN_DOC
! Squared distance matrix of the vertices
END_DOC
integer :: i, j, k
do i=1,vertex_num
do j=1,vertex_num
distance2(j,i) = 0.
do k=1,3
distance2(j,i) = distance2(j,i) + &
(vertex(k,i) - vertex(k,j))**2
enddo
enddo
enddo
END_PROVIDER
and the squared value of the edge is computed only once using:
BEGIN_PROVIDER [ real, edge2 ]
BEGIN_DOC
! edge2 : Square of the value of the edge
END_DOC
edge2 = edge**2
END_PROVIDER
Code execution
Compile the program
$ make
irpf90 -a -d
IRPF90_temp/cube.irp.F90
IRPF90_temp/properties.irp.F90
IRPF90_temp/cube_example.irp.F90
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/cube.irp.F90 -o IRPF90_temp/cube.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/properties.irp.F90 -o IRPF90_temp/properties.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/cube_example.irp.F90 -o IRPF90_temp/cube_example.irp.o
gfortran -ffree-line-length-none -O2 -c IRPF90_temp/irp_touches.irp.F90 -o IRPF90_temp/irp_touches.irp.o
gfortran -o cube_example IRPF90_temp/cube_example.irp.o IRPF90_temp/irp_stack.irp.o IRPF90_temp/cube.irp.o
IRPF90_temp/properties.irp.o IRPF90_temp/irp_touches.irp.o
and run it:
$ ./cube_example
0 : -> provide_volume
0 : -> provide_edge
0 : -> edge
Value of the edge of the cube
2
0 : <- edge 1.00000000000000002E-003
0 : <- provide_edge 1.00000000000000002E-003
0 : -> volume
0 : <- volume 0.0000000000000000
0 : <- provide_volume 1.00000000000000002E-003
0 : -> provide_center
0 : -> provide_vertex
0 : -> provide_vertex_num
0 : -> vertex_num
0 : <- vertex_num 0.0000000000000000
0 : <- provide_vertex_num 0.0000000000000000
0 : -> vertex
0 : <- vertex 0.0000000000000000
0 : <- provide_vertex 1.00000000000000002E-003
0 : -> provide_face_num
0 : -> face_num
0 : <- face_num 0.0000000000000000
0 : <- provide_face_num 0.0000000000000000
0 : -> provide_face
0 : -> provide_distance2
0 : -> distance2
0 : <- distance2 0.0000000000000000
0 : <- provide_distance2 0.0000000000000000
0 : -> provide_edge2
0 : -> edge2
0 : <- edge2 0.0000000000000000
0 : <- provide_edge2 0.0000000000000000
0 : -> face
0 : <- face 0.0000000000000000
0 : <- provide_face 0.0000000000000000
0 : -> center
0 : <- center 0.0000000000000000
0 : <- provide_center 1.00000000000000002E-003
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.0000000000000000
0 : <- provide_surface 0.0000000000000000
0 : -> cube_example
Surface Area : 24.000000
Volume : 8.0000000
Centers of the faces:
1.0000000 0.0000000 1.0000000
1.0000000 1.0000000 0.0000000
0.0000000 1.0000000 1.0000000
2.0000000 1.0000000 1.0000000
1.0000000 2.0000000 1.0000000
1.0000000 1.0000000 2.0000000
0 : <- cube_example 0.0000000000000000
Modification of the core of the program
Modification of the input data
In this section, we propose to modify the design of the cube. We want the user
to give in input the Cartesian coordinates of the vertices, and the entity
edge
will be computed. We will show that this modification has minor
impact on the rest of the code.
In the cube.irp.f
file, we first modify the provider of the edge
entity.
Considering that the vertices constitute a cube, the edge is the minimum value
of the distance matrix
BEGIN_PROVIDER [ real, edge ]
BEGIN_DOC
! Value of the edge of the cube
END_DOC
edge = sqrt(edge2)
END_PROVIDER
BEGIN_PROVIDER [ real, edge2 ]
BEGIN_DOC
! edge2 : Square of the value of the edge
END_DOC
integer :: i
edge2 = huge(1.)
do i=1,vertex_num
do j=i+1,vertex_num
edge2 = min(edge2,distance2(j,i))
enddo
enddo
END_PROVIDER
Then, we modify the provider of the vertex
entity to read 8 points from
the input file.
BEGIN_PROVIDER [ real, vertex, (3,vertex_num) ]
implicit none
BEGIN_DOC
! Coordinates of the vertices of the cube
END_DOC
integer :: k, i
! Initialize the array
print *, 'Vertices of the cube:'
do i=1,vertex_num
read(*,*) (vertex(k,i), k=1,3)
enddo
logical :: is_a_cube
ASSERT ( is_a_cube(vertex) )
END_PROVIDER
To be sure that the data read by the code is valid, we write the function
is_a_cube
which returns .True.
when the points given are vertices of a
cube. This function is a standard Fortran function:
logical function is_a_cube(v)
implicit none
real, intent(in) :: v(3,vertex_num)
is_a_cube = .True.
integer :: j
do j=1,vertex_num-1
! Choose a vector, then compute the dot product with all other vectors
real :: dot(vertex_num)
integer :: i
do i=1,vertex_num
dot(i) = 0.
integer :: k
do k=1,3
dot(i) = dot(i) + (v(k,i) - v(k,j))*(v(k,j+1) - v(k,j))
enddo
enddo
! Sort the dot array
integer :: pos(1)
real :: temp
do i=1,vertex_num
pos = minloc(dot(i:))
pos(1) = pos(1) + i - 1
temp = dot(i)
dot(i) = dot(pos(1))
dot(pos(1)) = temp
enddo
! Normalize to unity
real :: norm
norm = dot(vertex_num)
do i=1,vertex_num
dot(i) = dot(i) / norm
enddo
! Check the values of the dot products
real :: ref(vertex_num)
ref = (/ 0., 0., 0., 0., 1., 1., 1., 1. /)
do i=1,vertex_num
is_a_cube = is_a_cube .and. (dot(i) == ref(i))
enddo
if (.not.is_a_cube) then
is_a_cube = .True.
ref = (/ 0., 0., 1., 1., 1., 1., 2., 2. /)
ref = ref/2.
do i=1,vertex_num
is_a_cube = is_a_cube .and. (dot(i) == ref(i))
enddo
endif
if (.not.is_a_cube) then
is_a_cube = .True.
ref = (/ 0., 1., 1., 1., 2., 2., 2., 3. /)
ref = ref/3.
do i=1,vertex_num
is_a_cube = is_a_cube .and. (dot(i) == ref(i))
enddo
endif
if (.not.is_a_cube) then
return
endif
enddo
end function
Code execution
Build the new executable, and run the program :
$ make
irpf90 -a -d
IRPF90_temp/cube.irp.F90
gfortran -ffree-line-length-none -O0 -g -c IRPF90_temp/cube.irp.F90 -o IRPF90_temp/cube.irp.o
gfortran -ffree-line-length-none -O0 -g -c IRPF90_temp/properties.irp.F90 -o IRPF90_temp/properties.irp.o
gfortran -ffree-line-length-none -O0 -g -c IRPF90_temp/debug.irp.F90 -o IRPF90_temp/debug.irp.o
gfortran -ffree-line-length-none -O0 -g -c IRPF90_temp/cube_example.irp.F90 -o IRPF90_temp/cube_example.irp.o
gfortran -o debug IRPF90_temp/debug.irp.o IRPF90_temp/irp_stack.irp.o IRPF90_temp/properties.irp.o IRPF90_temp/cube.irp.o
gfortran -o cube_example IRPF90_temp/cube_example.irp.o IRPF90_temp/irp_stack.irp.o IRPF90_temp/properties.irp.o IRPF90_temp/cube.irp.o
$ ./cube_example
0 : -> provide_volume
0 : -> provide_edge
0 : -> provide_edge2
0 : -> provide_distance2
0 : -> provide_vertex_num
0 : -> vertex_num
0 : <- vertex_num 0.0000000000000000
0 : <- provide_vertex_num 0.0000000000000000
0 : -> provide_vertex
0 : -> vertex
Vertices of the cube:
0 0 0
2 2 2
2 0 0
2 2 0
0 2 0
0 0 2
0 2 2
2 0 2
0 : -> is_a_cube
0 : <- is_a_cube 0.0000000000000000
0 : <- vertex 0.0000000000000000
0 : <- provide_vertex 0.0000000000000000
0 : -> distance2
0 : <- distance2 0.0000000000000000
0 : <- provide_distance2 0.0000000000000000
0 : -> edge2
0 : <- edge2 0.0000000000000000
0 : <- provide_edge2 0.0000000000000000
0 : -> edge
0 : <- edge 0.0000000000000000
0 : <- provide_edge 0.0000000000000000
0 : -> volume
0 : <- volume 0.0000000000000000
0 : <- provide_volume 0.0000000000000000
0 : -> provide_center
0 : -> provide_face_num
0 : -> face_num
0 : <- face_num 0.0000000000000000
0 : <- provide_face_num 0.0000000000000000
0 : -> provide_face
0 : -> face
0 : <- face 0.0000000000000000
0 : <- provide_face 0.0000000000000000
0 : -> center
0 : <- center 0.0000000000000000
0 : <- provide_center 0.0000000000000000
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.0000000000000000
0 : <- provide_surface 0.0000000000000000
0 : -> cube_example
Surface Area : 24.000000
Volume : 8.0000000
Centers of the faces:
1.0000000 0.0000000 1.0000000
1.0000000 1.0000000 0.0000000
0.0000000 1.0000000 1.0000000
2.0000000 1.0000000 1.0000000
1.0000000 2.0000000 1.0000000
1.0000000 1.0000000 2.0000000
0 : <- cube_example 0.0000000000000000
The result is the same as what was obtained in the previous section, but now we can see that the sequence of the code is different. The entities of interest are computed in a different order.
Changing value of entities
In real applications, iterative processes are often used, and the values of entities change. In this section, we show how modification of entities is realized.
The iterative program
We write a new program which prints the value of the surface of the cube, as
long as the value of the surface is below a threshold. At the end of one
iteration, the length of the edge is incremented by the value increment
.
In file iterative_test.irp.f
, write:
program iterative_test
implicit none
do while (surface < threshold)
print *, surface
edge = edge + increment
TOUCH edge
enddo
end program
In this iterative process, the IRP entity edge
is modified. The TOUCH keyword
is mandatory informs the IRPF90 that the value of edge
is valid, but all
the values of the entities which depend on edge
are not valid anymore, and
have to be provided again.
The threshold
and increment
values are given in input, in file
control.irp.f
:
BEGIN_PROVIDER [ real, threshold ]
BEGIN_DOC
! Threshold for the value of the surface
END_DOC
print *, 'Threshold for the surface:'
read (*,*) threshold
ASSERT (threshold >= 0.)
END_PROVIDER
BEGIN_PROVIDER [ real, increment ]
BEGIN_DOC
! The increment of the value of the edge at each iteration
END_DOC
print *, 'Increment of the edge'
read(*,*) increment
ASSERT (increment > 0.)
END_PROVIDER
Now, if you build the program
$ make
irpf90 -a -d
gfortran -ffree-line-length-none -O0 -g -c IRPF90_temp/control.irp.F90 -o IRPF90_temp/control.irp.o
gfortran -o cube_example IRPF90_temp/cube_example.irp.o IRPF90_temp/irp_stack.irp.o
IRPF90_temp/properties.irp.o IRPF90_temp/control.irp.o IRPF90_temp/cube.irp.o
gfortran -o debug IRPF90_temp/debug.irp.o IRPF90_temp/irp_stack.irp.o IRPF90_temp/properties.irp.o
IRPF90_temp/control.irp.o IRPF90_temp/cube.irp.o
gfortran -ffree-line-length-none -O0 -g -c IRPF90_temp/iterative_test.irp.F90 -o IRPF90_temp/iterative_test.irp.o
gfortran -o iterative_test IRPF90_temp/iterative_test.irp.o IRPF90_temp/irp_stack.irp.o
IRPF90_temp/properties.irp.o IRPF90_temp/control.irp.o IRPF90_temp/cube.irp.o
you can remark that the executable cube_example
of the previous section is still built, and a new executable iterative_test
is created.
Running the program gives:
$ ./iterative_test
0 : -> provide_threshold
0 : -> threshold
Threshold for the surface:
100.
0 : <- threshold 0.00000000000000
0 : <- provide_threshold 0.00000000000000
0 : -> provide_edge
0 : -> provide_edge2
0 : -> provide_distance2
0 : -> provide_vertex_num
0 : -> vertex_num
0 : <- vertex_num 0.00000000000000
0 : <- provide_vertex_num 0.00000000000000
0 : -> provide_vertex
0 : -> vertex
Vertices of the cube:
0 0 0
2 2 2
2 0 0
2 2 0
0 2 0
0 0 2
0 2 2
2 0 2
0 : -> is_a_cube
0 : <- is_a_cube 0.00000000000000
0 : <- vertex 0.00000000000000
0 : <- provide_vertex 0.00000000000000
0 : -> distance2
0 : <- distance2 0.00000000000000
0 : <- provide_distance2 0.00000000000000
0 : -> edge2
0 : <- edge2 0.00000000000000
0 : <- provide_edge2 0.00000000000000
0 : -> edge
0 : <- edge 0.00000000000000
0 : <- provide_edge 0.00000000000000
0 : -> provide_increment
0 : -> increment
Increment of the edge
1.
0 : <- increment 0.00000000000000
0 : <- provide_increment 0.00000000000000
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.00000000000000
0 : <- provide_surface 0.00000000000000
0 : -> iterative_test
24.00000
0 : -> touch_edge
0 : <- touch_edge 0.00000000000000
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.00000000000000
0 : <- provide_surface 0.00000000000000
54.00000
0 : -> touch_edge
0 : <- touch_edge 0.00000000000000
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.00000000000000
0 : <- provide_surface 0.00000000000000
96.00000
0 : -> touch_edge
0 : <- touch_edge 0.00000000000000
0 : -> provide_surface
0 : -> surface
0 : <- surface 0.00000000000000
0 : <- provide_surface 0.00000000000000
0 : <- iterative_test 0.00000000000000
In this example, it appears clearly that when the edge value is modified, only the surface is computed again. The volume is invalid, but as it is not requested, it is not computed.
Embedding shell scripts
It is common practice to use shell scripts to gather data at compilation time.
In the IRPF90 environment, shell scripts can be directly introduced in the
code. Let us write a simple header for the code, which returns the name of the
user who compiled the code, and the date of compilation. We now remove the
-d
option of irpf90
for shorter outputs.
program iterative_test
implicit none
print *, 'Program ', irp_here
BEGIN_SHELL [ /bin/sh ]
echo "print *, \'Compiled by : $USER \'"
echo "print *, \'Compilation date: `date`\'"
END_SHELL
do while (surface < threshold)
print *, surface
edge = edge + increment
TOUCH edge
enddo
end program
Running this program returns the following output:
Threshold for the surface:
100
Vertices of the cube:
0 0 0
2 2 2
2 0 0
2 2 0
0 2 0
0 0 2
0 2 2
2 0 2
Increment of the edge
1
Program iterative_test
Compiled by : scemama
Compilation date: Fri Sep 25 16:02:22 CEST 2009
24.00000
54.00000
96.00000
First, note the use of the irp_here
variable. This is a character*(*)
variable which takes as value the name of the subroutine or function in which
it is used. Then, you can remark that the print statements were executed after
gathering information from the input. This is due to the fact that in IRPF90,
the entities are provided as soon as possible. If you really want to print the
header at the beginning of the program, you can use the following trick:
program iterative_test
implicit none
print *, 'Program ', irp_here
BEGIN_SHELL [ /bin/sh ]
echo "print *, \'Compiled by : $USER \'"
echo "print *, \'Compilation date: `date`\'"
END_SHELL
call run_iterative_process
end program
subroutine run_iterative_process
implicit none
do while (surface < threshold)
print *, surface
edge = edge + increment
TOUCH edge
enddo
end program
which gives the output:
Program iterative_test
Compiled by : scemama
Compilation date: Fri Sep 25 16:03:58 CEST 2009
Threshold for the surface:
100
Vertices of the cube:
0 0 0
2 2 2
2 0 0
2 2 0
0 2 0
0 0 2
0 2 2
2 0 2
Increment of the edge
1
24.00000
54.00000
96.00000
Shell scripts can also be used to write templates. For example, if you need to sort arrays of real, double precision or integers, you can use the following Python script:
BEGIN_SHELL [ /usr/bin/python ]
for i in [('' ,'real'), \
('d','double precision'), \
('i','integer'), \
]:
print "subroutine "+i[0]+"sort (x,iorder,isize)"
print " implicit none"
print " "+i[1]+" :: x(*), xtmp"
print " integer :: iorder(*)"
print " integer :: isize"
print " integer :: i, i0, j, jmax"
print ""
print " do i=1,isize"
print " xtmp = x(i)"
print " i0 = iorder(i)"
print " do j=i-1,1,-1"
print " if ( x(j) > xtmp ) then "
print " x(j+1) = x(j)"
print " iorder(j+1) = iorder(j)"
print " else"
print " exit"
print " endif"
print " enddo"
print " x(j+1) = xtmp"
print " iorder(j+1) = i0"
print " enddo"
print ""
print "end"
print ""
which builds in one shot three subroutines:
isort
for integersdsort
for double precisionsort
for reals
Entities of interest can also be generated by scripts. The following example
(documentation.irp.f
) builds a character*(*)
IRP entity for each entity
which contains its documentation:
BEGIN_SHELL [ /usr/bin/python ]
import os
doc = {}
for filename in os.listdir('.'):
if filename.endswith('.irp.f'):
file = open(filename,'r')
inside_doc = False
for line in file:
if line.strip().lower().startswith('begin_provider'):
name = line.split(',')[1].split(']')[0].strip()
doc[name] = ""
elif line.strip().lower().startswith('begin_doc'):
inside_doc = True
elif line.strip().lower().startswith('end_doc'):
inside_doc = False
elif inside_doc:
doc[name] += line[1:].strip()+" "
file.close()
lenmax = 0
for e in doc.keys():
lenmax = max(len(e),lenmax)
print "BEGIN_PROVIDER [ character*(%d), entities, (%d) ]"%(lenmax,len(doc))
print " BEGIN_DOC"
print "! List of IRP entities"
print " END_DOC"
for i,e in enumerate(doc.keys()):
print "entities(%d) = '%s'"%(i+1, e)
print "END_PROVIDER"
for e in doc.keys():
print "BEGIN_PROVIDER [ character*(%d), %s_doc ]"%(len(doc[e]),e)
print " BEGIN_DOC"
print "! Documentation of variable %s"%(e,)
print " END_DOC"
print " %s_doc = '%s'"%(e,doc[e])
print "END_PROVIDER"
END_SHELL
and a new main program is created (get_doc.irp.f
) to print the documentation of
a variable if it is present in the command line:
program get_doc
integer :: iargc
character*(32) :: arg
integer :: i, j
if (iargc() == 0) then
print *, 'List of IRP entities'
do j=1,size(entities)
print *, entities(j)
enddo
return
endif
do i=1,iargc()
call getarg(i,arg)
BEGIN_SHELL [ /usr/bin/python ]
import os
entities = []
for filename in os.listdir('.'):
if filename.endswith('.irp.f'):
file = open(filename,'r')
for line in file:
if line.strip().lower().startswith('begin_provider'):
name = line.split(',')[1].split(']')[0].strip()
entities.append(name)
file.close()
for e in entities:
print " if (arg == '%s') then"%(e,)
print " print *, %s_doc"%(e,)
print " endif"
END_SHELL
enddo
end
Execution of this code gives:
$ ./get_doc
List of IRP entities
distance2
center
vertex
surface
face
volume
vertex_num
edge2
edge
increment
threshold
face_num
$ ./get_doc volume
Volume of the cube
Other features
Freeing memory
Memory of an IRP entity can be freed using the FREE
keyword
FREE x
where x
is an array entity. The memory occupied by x
will be freed,
and its status will be tagged as invalid. If x
is needed later, the memory
will be re-allocated, and x
will be re-built.
If it is not an array, the FREE
keyword will only mark it as invalid.
Conditional compilation
Conditional compilation is possible using the IRP_IF ... IRP_ELSE .. IRP_ENDIF
directives.
IRP_IF MPI
include 'mpif.h'
print *, 'Multiprocessor code'
IRP_ELSE
print *, 'Monoprocessor code'
IRP_ENDIF
Compiling the previous code with irpf90 -DMPI
will compile code suitable
for the use of the MPI library, otherwise, the mono-processor code will be
compiled.
Conclusion
Many observations can be made from this simple example.
In the first section, we started to write a simple code. At the time we wrote code, we did not build the design for future improvement: only the edge length of the cube was needed.
Then, we wanted to compute another property depending on the coordinates of the vertices. This quantity was easily introduced, without any interference with the code which was written before.
Then, we chose to change the internal representation of the cube. The user gave in input the coordinates of the vertices and the edge value was computed from them. Making this modification did not interfere at all with the rest of the program.
From this example, we can conclude that, using the IRPF90 environment, scientific programming is simpler. This simplicity of writing code lets the scientific programmer focus on science instead of focusing on memory allocation or makefiles.
The code is also clearer, since the information is very localized. There is only one way to compute a quantity, and it is located in the provider of this quantity.
The use of embedded scripts allows the programmer to reduce considerably the number of lines of code, and also to automatically update certain parts of the code upon modification. For instance, in the presented example, if the programmer adds a new IRP entity, the documentation program will be automatically updated.
The resulting code is usually faster than code written in Fortran. Indeed, as the programmer is forced to partition his code in small functions (providers), very few memory locations are used en each function and the compiler can memory accesses are well optimized. Moreover, programmers will write (unintentionally) code which will be easier for the compilers to optimize.