irpf90/website/content/post/2010-01-18-tutorial.md

---
title: Simple tutorial
subtitle: First steps with IRPF90
date: 2010-01-18
---

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:

![Production Tree](/img/Energy_tree.png)


### 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

```fortran
    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:

```fortran
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

```fortran
E_nucl = compute_E_nucl()
```

has to be positioned before

```fortran
E      = compute_E(E_nucl,E_elec)
```

otherwise the program is wrong.

If the code is written in this way:

```fortran
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

```fortran
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

```irpf90
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

```bash
$ irpf90 --init
```

Two temporary directories are created:

```bash
$ ls
IRPF90_man  IRPF90_temp  Makefile
```

and a standard `Makefile` is build, using the gfortran compiler:

```make
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:

```irpf90
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.

```irpf90
BEGIN_PROVIDER [ real, surface ]

 BEGIN_DOC
! Surface of the cube
 END_DOC

 surface = 6. * edge**2
END_PROVIDER
```

```irpf90
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``:

```irpf90
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:

```bash
$ 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:

```bash
$ 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:

```irpf90
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:

```irpf90
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``:

```irpf90
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:

```irpf90
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:

```irpf90
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:

```irpf90
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

```bash
$ 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

```irpf90
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.

```irpf90
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:

```irpf90
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 :

```bash
$ 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:

```irpf90
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``:

```irpf90
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

```bash
$ 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.

```irpf90
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:

```irpf90
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:

```python
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 integers
* ``dsort`` for double precision
* ``sort`` 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:

```python
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:

```irpf90
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:

```bash
$ ./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

```irpf90
 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.

```irpf90
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.