All temporary arrays are allocated dynamically. In addition, arrays allocated dynamically are always deallocated by the same routine that allocates them.
Fortran 90 constructs (allocatable and automatic) are used by default with compilers that allow such constructs. Malloc and free (or comparable routines) are used with the other compilers, and can be used with those compilers that allow F90 constructs if desired by modifying cpp defines in the makefiles.
The actual lines of code that declare, allocate, and deallocate the dynamically-allocated arrays are generated automatically by a gawk script that parses a file containing information about the arrays that need to be allocated dynamically for each routine. This makes porting to new machines and compilers easier, and minimizes chances for mistakes like memory leaks (forgetting to add a call to free for an array you added), but complicates the makefiles somewhat due to the complicated dependencies.
[Need to describe this in detail, with code examples.]
c***********************************************************************
c*
c* Computes y=Ax where A is a matrix stored in ELL format.
c*
c* Note that the dynamically-allocated vector "xtmp" is dimensioned
c* to start at zero. This is because in ELL format, the index map
c* matrix "ja" contains 0's where elements of "a" are zero. Thus the
c* 0th element of "xtmp" can be referenced.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idim - leading dimension of a and ja
c* n - number of unknowns
c* maxnz - maximum number of non-zero elements in any row
c* a - matrix
c* ja - map array
c* x - x-vector
c*
c* Output:
c* y - y-vector
c* status - return status
c* -2 ==> memory allocation failure
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
c* <SUBROUTINES REQUIRED>
c*
c* JT_FillVectorFloat
c* JT_y_eq_x
c*
#include "arrays-y_eq_Ax_ELL.inc"
c*
c***********************************************************************
subroutine JT_y_eq_Ax_ELL (idim, n, maxnz, a, ja, x, y, status)
implicit none
c
c ... Input:
integer idim, maxnz, n
integer ja(idim,maxnz)
real a(idim,maxnz), x(n)
c
c ... Output:
integer status
real y(n)
c
c ... Local:
integer i, j
real zero
#include "declare-y_eq_Ax_ELL.inc"
c
parameter (zero=0.0d0)
c
c ... Initialize return status.
status = 0
c
#include "allocate-y_eq_Ax_ELL.inc"
c
xtmp(0) = zero
call JT_y_eq_x (n, x, xtmp(1), status)
c
call JT_FillVectorFloat (n, zero, y, status)
do j=1,maxnz
do i=1,n
y(i) = y(i) + a(i,j)*xtmp(ja(i,j))
enddo
enddo
c
9999 continue
#include "deallocate-y_eq_Ax_ELL.inc"
c
return
end