Physics Department, University of California at Los Angeles
Los Angeles, CA 90095-1547
and
Jet Propulsion Laboratory, California Institute of Technology
Pasadena, California 91109
email: decyk@physics.ucla.edu
Jet Propulsion Laboratory, California Institute of Technology
MS 168-522, 4800 Oak Grove Drive
Pasadena, California 91109
email: Charles.D.Norton@jpl.nasa.gov
Department of Computer Science
and Scientific Computation Research Center (SCOREC)
110 Eighth Street, Rensselaer Polytechnic Institute, Troy, NY 12180-3590
email: szymansk@cs.rpi.edu
Abstract
Fortran90 is a modern, powerful language with features that support important new programming concepts, including those used in object-oriented programming. This paper explains the concepts of data encapsulation, function overloading, classes, objects, inheritance, and dynamic dispatching, and how to implement them in Fortran90. As a result, a methodology can be developed to do object-oriented programming in the language.
A postscript version of this paper is also available.
There are several properties of program design that are useful in the programming and maintenance of large computational projects. First, the code which modifies a given data structure should be localized or confined to one location in the program and not spread, uncontrollably, throughout the program. This property is called an encapsulation. In a sense, it is a generalization of the familiar notion of a function or subroutine. However, modern encapsulation allows all related operations to be grouped together around a data type. For example, given a data type representing a logical, encapsulation may group together a set of operations that perform logical operations. A related notion is that of information hiding. Once all the operations are encapsulated, one can hide the details by which the results are computed, allowing only well defined operations to be applied to the data. For example, one can use the .not. operation on the intrinsic type logical, without worrying about the way the compiler writers implemented the operation or the internal representation of logicals. The divide operation on logicals, on the other hand, is not defined in Fortran and therefore cannot be used.
Using encapsulation and information hiding, the users can define their own data types, called abstract data types. For example, one can define a data type representing an angle. Abstract data types are further defined by the operations that can be applied to them. A class definition encapsulates the code for the allowable operations, or methods, written by the programmer along with the abstract data type, hiding the implementation details there. For instance, users of the data type angle can apply the sine and cosine functions to angle variables (objects) if such methods have been defined for angles, but cannot multiply two angles if such a method has not been defined. Encapsulated code can be safely changed (for example, a different method of computing sines could be used on a new architecture) without the need to change the code using these functions. The other advantage of abstract data types is that the software designer can describe the program in terms resembling the application domain. For example, if a data type electron is defined, one could define what it means to push an electron object a given distance in space. (The electron is pushed in this view, whereas in Fortran77 the values of arrays representing the electron coordinates are updated).
Some operations are defined for more than one standard type. For example multiplication can be defined for integer arguments, real arguments, vectors and arrays. This property is called overloading and the concept extends to abstract data types. Consider an electron type and an ion type, both of which permit a push operation to be applied, but using a different procedure in each case. Overloading helps to write shorter and clearer programs by providing a general operation in which the specific action taken is defined by the type of the arguments at compile- time.
Often in science, properties of some entities are abstracted and grouped into classes. For example, electron and ion types can be considered different kinds of species types. If separate data types were created for electron and ion, many operations would be shared by these two types. Overloading can help, but it is even more convenient to introduce the abstract data type species for the common properties. Then, like in everyday life, everything that applies to the species type can also be applicable to electron and ion types (but not the opposite). Once the data type species is defined one can simplify the description of ion and electron types by permitting inheritance of species properties by those two more specialized data types.
All of the ideas presented so far are part of the modern programming paradigm called object-oriented programming (OOP). Many users with large investments in Fortran77 have reason to be reluctant to shift to a very different programming style and language. The demands on development time and the training required to shift to a new approach, combined with the many years invested in existing Fortran77 based applications and the need to continually produce new science are good reasons to be skeptical of experimenting with new ideas, as promising as they might appear.
Fortran90 supports these new important programming concepts, including those used in object-oriented programming. Since Fortran90 is backward compatible with Fortran77, it is possible to incorporate these new ideas into old programs in an incremental fashion, enabling the scientist to continue his or her scientific activities. Some of these ideas are useful for the typical kinds of programs written by individual authors now. The usefulness of other ideas only becomes apparent for more ambitious programs written by multiple authors. These are programs that might never have been written in Fortran77 because the complexity involved would have been unmanageable. These new ideas enable more productive program development, encourage software collaboration and allow the scientist to use the same abstract concepts in a program that have been used so successfully in scientific theory. Scientific productivity will then improve. Additionally, there is also a migration path to parallel computers, since High Performance Fortran (HPF) is also based on Fortran90.
In this paper, we will introduce the concepts of data encapsulation, function overloading, classes and objects, inheritance, and dynamic dispatching in Fortran90. Since these are the fundamental building blocks of object-oriented programming, it is important to understand them before one can effectively use OOP in Fortran90. Many of these ideas are powerful by themselves and are useful even without adopting the object-oriented paradigm. This paper is intended to be introductory in nature. For those who wish to pursue these ideas further, there are a number of references which discuss object- oriented ideas in a language independent manner [1-3]. There are many textbooks available on Fortran90 and C++. Two that we have found useful are those by Ellis [4] and Lippman [5].
subroutine fft1r(f,t,isign,mixup,sct,indx,nx,nxh)
integer isign, indx, nx, nxh, mixup(nxh)
real f(nx)
complex sct(nxh), t(nxh)
c rest of procedure goes here
return
end
Here f is the array to be transformed (and the output), t is a temporary work array, mixup is a bit-reverse table, sct is a sine/cosine table, indx is the power of 2 defining the length of the transform, and nx (>=2**indx) is the size of the f array, and nxh (=nx/2) is the size of the remaining arrays. The variable isign determines the direction of the transform (or if zero, initiates the tables mixup and sct.)
Since the procedure fft1r is designed to work with variable
size data (and might be compiled separately), and Fortran77 cannot
dynamically allocate such data, the work and table arrays t, sct,
and mixup have to be declared in the main program and passed to
the procedure. If the FFT procedure is itself embedded inside
other procedures, then all these arrays have to be passed down the
chain of arguments or else stored in a common block. Thus the
main program might look something like:
program main
integer isign, indx, nx, nxh
parameter(indx=11,nx=2**indx,nxh=nx/2)
integer mixup(nxh)
real f(nx)
complex sct(nxh), t(nxh)
c initialize fft tables
isign = 0
call fft1r(f,t,isign,mixup,sct,indx,nx,nxh)
stop
end
The goal of data encapsulation is make the FFT call look like:
call fft1r(f,isign)
where all the auxiliary arrays and constants which are needed only by the FFT are hidden inside the FFT, and the rest of the program does not have to be concerned about them. This hiding of data greatly simplifies bookkeeping with procedures.
Fortran90 allows dynamic arrays which are only used inside a
procedure to be created and destroyed there, and they are
therefore unknown outside the procedure. One mechanism for such
encapsulation is the automatic array, which is created on entry
and destroyed upon exiting a procedure. It is easy to implement
the work array t as an automatic array. One just omits the name
from the argument list:
subroutine fft1r(f,isign,mixup,sct,indx,nx,nxh)
integer isign, indx, nx, nxh, mixup(nxh)
real f(nx)
complex sct(nxh)
! t is an automatic array, it disappears on exit
complex, dimension(nxh) :: t
! rest of procedure goes here
end subroutine fft1r
Notice that we have begun to use the new Fortran 90 :: syntax for declaring arrays, and a new END SUBROUTINE statement. Automatic arrays differ from local arrays previously available in Fortran77 because their dimensions can now be variables.
Another mechanism for encapsulation of arrays in Fortran90 is the allocatable (or deferred-size) array. They are similar to automatic arrays, except that their creation (allocation) and destruction (deallocation) are entirely under programmer control. Such arrays are created with the ALLOCATE statement. If the SAVE attribute is used in their declaration, then they will not be destroyed on exit from the procedure. They can be explicitly destroyed with the DEALLOCATE statement.
For now, let us assume that the table arrays mixup and sct
do not change between calls to the FFT. We can then remove them
from the argument list in our example and explicitly ALLOCATE them
inside the procedure when the tables are initialized, as follows:
subroutine fft1r(f,isign,indx,nx,nxh)
integer isign, indx, nx, nxh
real f(nx)
! mixup and sct are saved, allocatable arrays
integer, dimension(:), allocatable, save :: mixup
complex, dimension(:), allocatable, save :: sct
! t is an automatic array
complex, dimension(nxh) :: t
! special initialization call
if (isign.eq.0) allocate(mixup(nxh),sct(nxh))
! rest of procedure goes here
end subroutine fft1r
Later, we will add error checking conditions.
A very powerful feature of Fortran90 is that arrays are
actually array objects which contain not only the data itself, but
information about their size. This was previously only available
for character arrays in Fortran77, which supplied the LEN
intrinsic to obtain the character length. In Fortran90, assumed-
shape arrays are available whose dimensions can be obtained from
the SIZE intrinsic. This is a third mechanism useful for data
encapsulation. They are declared like ordinary arrays, except
their dimension lengths are replaced with a colon. We can now
omit the dimensions nx and nxh from the argument list, and obtain
them inside the procedure. The declaration of the automatic array
t also has to be revised to use the SIZE intrinsic. The result
is:
! f is an assumed-shape array
real, dimension(:) :: f
integer :: isign, indx, nx, nxh
integer, dimension(:), allocatable, save :: mixup
complex, dimension(:), allocatable, save :: sct
! t is an automatic array whose size is determined from array f
complex, dimension(size(f)/2) :: t
! size of arrays mixup and sct are determined from array f.
nx = size(f)
nxh = nx/2
! special initialization call
if (isign.eq.0) allocate(mixup(nxh),sct(nxh))
! rest of procedure goes here
end subroutine fft1r
In order to use assumed-shape arrays the compiler must have
knowledge of the actual argument types being used when the
procedure is called. One way this information can be supplied is
with the INTERFACE statement, which declares to the main program
the argument types of the procedure. For example, one can write:
program main
integer :: isign
integer, parameter :: indx=11, nx=2**indx
real, dimension(nx) :: f
! declare interface
interface
subroutine fft1r(a,i,j)
real, dimension(:) :: a
integer :: i, j
end subroutine fft1r
end interface
! initialize fft tables
isign = 0
call fft1r(f,isign,indx)
stop
end
where we have used the new form of the PARAMETER statement. An
additional advantage of explicit INTERFACE blocks is the compiler
will now check and require that the actual arguments passed to the
procedure match in number and type with those declared in the
interface. Thus if one accidentally omits the argument isign in a
procedure call
call fft1r(f,indx)
the compiler will flag this.
One possible source of error is that one can mistakenly
declare the data in an INTERFACE block to be different than the
actual data in the procedure. This source of error is removed if
the procedure is stored in a MODULE which is then "used," because
in this case the compiler creates the INTERFACE automatically.
The USE statement is similar to the INCLUDE extension commonly
found in Fortran77, but it is not a text substitution. Rather, it
makes information "available", and is much more powerful than the
INCLUDE statement. Thus:
module fft1r_module
contains
subroutine fft1r(f,isign,indx)
real, dimension(:) :: f
integer :: isign, indx, nx, nxh
integer, dimension(:), allocatable, save :: mixup
complex, dimension(:), allocatable, save :: sct
complex, dimension(size(f)/2) :: t
nx = size(f)
nxh = nx/2
! special initialization call
if (isign.eq.0) allocate(mixup(nxh),sct(nxh))
! rest of procedure goes here
end subroutine fft1r
end fft1r_module
!
program main
use fft1r_module ! explicit interface not needed now
integer :: isign = 0
integer, parameter :: indx=11, nx=2**indx
real, dimension(nx) :: f
! initialize fft tables
call fft1r(f,isign,indx)
stop
end
where we have used a new way to initialize the integer isign.
Fortran90 supports a number of other statements and
attributes that contribute to programming safety. One is the
IMPLICIT NONE statement that requires all variables to be
explicitly declared. Another is the INTENT attribute for
arguments, that declares whether arguments are intended as input
only, output only, or both. If we declare the arguments in fft1r
as follows:
subroutine fft1r(f,isign,indx)
implicit none
real, dimension(:), intent(inout) :: f
integer, intent(in) :: isign, indx
then the variables isign and indx cannot be modified in this procedure since they were declared with INTENT(IN) attributes. Such features mean that more errors are now caught by the compiler rather than by the operating system when the code is running. Because of this added safety we will use modules for all of the remaining subroutines in this paper.
The encapsulation of data we have illustrated in this example makes it easier for multiple authors to independently develop programs which will be used by others. The FFT program now has a simple interface which is less likely to be changed, so that even if the author of the procedure makes changes internally, users of the procedure do not have to change their code.
Another benefit of such an approach is that one can hide old, ugly code that cannot be changed, perhaps because one does not have access to the source. For example, if one were using a library FFT which was optimized for some specific architecture, one could encapsulate it in a "shell" procedure which allocates any work or table arrays needed and then calls the library FFT. Since the details of the library FFT are hidden from the user, one could replace it with another by making changes only in this "shell" procedure and not impact the rest of the code. This allows a code to remain portable and yet optimized.
Let us now add more error checking when allocating data.
The ALLOCATE statement allows an optional error return code to
check if the data was actually allocated. And the ALLOCATED
statement allows one to check if the data is already allocated
(perhaps we had previously used the FFT with different length
data). Thus the following version of the allocation is safer:
if (isign.eq.0) then
if (allocated(mixup)) deallocate(mixup)
if (allocated(sct)) deallocate(sct)
allocate(mixup(nxh),sct(nxh),stat=ierr)
if (ierr.ne.0) then
print *,'allocation error'
stop
endif
endif
One can simplify the interface even further by noting that the FFT
length parameter indx is only needed during the initialization of
the FFT. In subsequent calls it would be an error to use a
different value without initializing new tables. Fortran90
supports the use of OPTIONAL arguments and the intrinsic PRESENT
to determine if it was passed. With these tools, we can save the
indx parameter passed during initialization and not require it to
be passed subsequently. Furthermore, if we initialize the saved
indx parameter to some nonsense value, we can test it subsequently
to prevent the FFT from being used before the FFT tables were
initialized. The result is:
subroutine fft1r(f,isign,indx)
integer, intent(in), optional :: indx
integer, save :: saved_indx = -1 ! initialize to nonsense
if (isign.eq.0) then
if (.not.present(indx)) then
print *,'indx must be present during initialization!'
stop
endif
saved_indx = indx
else
if (saved_indx.lt.0) then
print *,'fft tables not initialized!'
stop
endif
endif
In the following version of the FFT, the procedure lib_fft1r is
intended to refer to some library FFT where either the source code
is unavailable or one does not desire to change it. This new
version is much safer and easier to use and modify.
subroutine fft1r(f,isign,indx)
implicit none
real, dimension(:), intent(inout) :: f
integer, intent(in) :: isign
integer, intent(in), optional :: indx
integer :: nx, nxh, ierr
integer, save :: saved_indx = -1
integer, dimension(:), allocatable, save :: mixup
complex, dimension(:), allocatable, save :: sct
complex, dimension(size(f)/2) :: t
nx = size(f)
nxh = nx/2
! special initialization call
if (isign.eq.0) then
! indx must be present during initialization
if (.not.present(indx)) then
print *,'indx must be present during initialization'
stop
endif
! indx must be non-negative
if (indx.lt.0) then
print *,'indx must be non-negative'
stop
endif
! save indx for future calls
saved_indx = indx
! deallocate if already allocated
if (allocated(mixup)) deallocate(mixup)
if (allocated(sct)) deallocate(sct)
! allocate table arrays
allocate(mixup(nxh),sct(nxh),stat=ierr)
! check if allocation error
if (ierr.ne.0) then
print *,'allocation error'
stop
endif
! make sure fft tables initialized
else
if (saved_indx.lt.0) then
print *,'fft tables not initialized!'
stop
endif
endif
! call old, ugly but fast fft here
! saved_indx used here instead of indx
call lib_fft1r(f,t,isign,mixup,sct,saved_indx,nx,nxh)
end subroutine fft1r
During initialization one would call:
call fft1r(f,isign,indx)
But subsequently one can use just:
call fft1r(f,isign)
Logically, we have bundled two distinct operations, initialization and performing the FFT, into one procedure. There is a third procedure which would be useful, deallocating the internal table arrays to free up memory if we are done with performing FFTs. This could be done by adding an extra, OPTIONAL argument to the procedure, and adding more code, but this is unattractive. It makes for clearer programming to separate logically distinct operations into distinct procedures. The difficulty with this is how to allow distinct procedures to share access to the internal table arrays mixup and sct. In Fortran77, the only mechanism to do this was common blocks. In Fortran90, there is a new mechanism: modules can contain global data which are shared by all the procedures in the module without explicitly declaring them inside the procedures. This is a new idea in Fortran, although common in other languages such as C++. Furthermore, this global data can be made local to the module and inaccessible to other procedures which USE the module.
In our example, we will make the table arrays mixup and sct,
as well as the integer saved_indx global by moving their
declaration outside the procedure to the declaration section at
the beginning of the module. We will also add the PRIVATE
attribute, to block access to this data from outside the module.
Adding the deallocation procedure, the module looks like:
module fft1r_module
! all module procedures have access to this data
integer, save, private :: saved_indx = -1
integer, dimension(:), allocatable, save, private :: mixup
complex, dimension(:), allocatable, save, private :: sct
contains
subroutine fft1r_end
! this procedure has access to saved_indx, mixup, and sct
! reset saved_indx to nonsense value
saved_indx = -1
! deallocate table arrays
deallocate(mixup,sct)
end subroutine fft1r_end
! other procedures go here
end module fft1r_module
By separating the original fft1r procedure into a new initialization (fft1r_init) and FFT procedure, we no longer need to use optional arguments. In the final version of this module which is shown in Appendix A, we have used the ';' syntax which allows multiple statements on one line.
Allocatable arrays can also be used in the main program,
which allows one to create all arrays at run time rather than at
compile time. In Fortran90, one no longer has to recompile a code
because the dimensions change. (It is the programmer's
responsibility, however, not to use an allocatable array before it
has been allocated or after it has been deallocated.) In the
following main program, we make f an allocatable array, obtain the
value of indx from the input device, allocate f, and use the new
array constructor syntax to initialize it.
program main
use fft1r_module
implicit none
integer :: indx, nx, i
real, dimension(:), allocatable :: f
! write prompt without linefeed
write (6,'(a)',advance='no') 'enter indx: '
! obtain indx from input device
read (5,*) indx
! allocate array f
nx = 2**indx
allocate(f(nx))
! initialize data using array constructor
f = (/(i,i=1,nx)/)
! initialize fft
call fft1r_init(indx)
! call fft
call fft1r(f,-1)
! terminate fft
call fft1r_end
stop
end
Notice that we have not modified the original lib_fft1r procedure, which is a private procedure in the module. Instead we have simplified the user interface to the bare essentials while adding substantial safety to its usage.
In Fortran90, generic functions allow user defined functions to also have this feature. In the case of the FFT example, users of Fortran77 have had to remember to use different function names for every possible type of FFT, such as real to complex, complex to complex, 1 dimensional, 2 dimensional, single precision or double precision FFTs. The generic function facility allows a single name, for example, fft, to be used for all of them, and the compiler will automatically select the correct FFT to use based on the number and types of arguments actually used.
Thus in the case of the 1d real to complex FFT, we had a
procedure with the following interface:
subroutine fft1r(f,isign)
real, dimension(:), intent(inout) :: f
integer, intent(in) :: isign
In a manner similar to what we described in the previous section,
one can construct a 2d real to complex FFT with the following
interface:
subroutine fft2r(f,isign)
real, dimension(:,:), intent(inout) :: f
integer, intent(in) :: isign
where the procedure fft2r "hides" an old, ugly but fast 2d real to
complex FFT which has lots of arguments. In the first case the
argument f is a real, one dimensional array, while in the second
case it is a real, two dimensional array. If both of these
procedures are in the same module, one constructs a generic
function fft by placing the following statements in the
declaration section of the module:
interface fft
module procedure fft1r
module procedure fft2r
end interface
Then in a main program which uses the module, the statement
call fft(f,isign)
will call the procedure fft1r, if f is a real, one dimensional array, or will call fft2r, if f is a real, two dimensional array. If f is any other type of argument, a compile error will be generated.
If the 2 dimensional FFT is contained in a separate module,
then two separate INTERFACE statements are required. In the first
module one includes
interface fft
module procedure fft1r
end interface
and in the second module one includes
interface fft
module procedure fft2r
end interface
The main program then uses both modules, and each module procedure will then be added to the list of procedures that have the generic interface fft.
In a similar manner, one can include in the same interface
all other types of FFTs and the programmer is protected from
making errors in calling the wrong procedure. If desired, one can
even make the specific names fft1r, fft2r inaccessible by adding
the declaration:
private :: fft1r, fft2r
in the module. One advantage of this is that the specific names can be reused in other modules without conflict. Function overloading is also called ad hoc polymorphism.
subroutine push1 (part,fx,qbm,dt,ek,np,idimp,nop,nx)
integer np, idimp, nop, nx
real qbm, dt, ek, part(idimp,nop), fx(nx)
c rest of procedure goes here
return
end
Here part is the array which contains the particle coordinates and
velocities and fx is the electric field array. The integer idimp
is the dimensionality of phase space, nop is the maximum number of
particles allowed, and nx is the size of the electric field array.
As we saw in the FFT example, we do not have to pass these
integers in Fortran90, since they can be determined by the SIZE
intrinsic if part and fx are passed as assumed-shape arrays. The
integer np (np<=nop) is the actual number of valid particles in
the part array, qbm is the charge/mass ratio, ek is the kinetic
energy of the np valid particles, and dt is the time step. All of
the scalars except for the time step describe a group of charged
particles and they usually are passed together whenever the
particles are processed by some procedure. We can use a derived
type to store them together as follows:
type species_descriptor
integer :: number_of_particles
real :: charge, charge_to_mass, kinetic_energy
end type species_descriptor
This is similar to structures and record types which appear in
other programming languages. We have added charge to the list,
since there are some procedures which also require that. Notice
that in this derived type, there are components of both integer
and real type, so that this could not have been implemented with
just an array in Fortran77. To create a variable of this type,
one makes the following declaration:
type (species_descriptor) :: electron_args, ion_args
where we have created two variables of type species_descriptor,
one for electrons and one for ions. The components of this new
type are accessed with the '%' symbol. Thus we can assign values
as follows:
electron_args%number_of_particles = 1000
electron_args%charge = 1.0
It is best to put the definition in the declaration section of a
module along with the new push1 subroutine (to avoid having to declare
an explicit interface for it), as shown in Appendix
B. Then this module is "used" in the main program to give access
to the derived type and new push1 procedure, which can now be called
with a much simpler interface:
call push1(part,fx,electron_args,dt)
We have shown here a simple use of derived types, merely to
reduce bookkeeping when passing arguments to procedures. But
derived types are much more powerful than that. They can be used
to express sophisticated, abstract quantities. In fact, with
derived types it is possible to express in programming the same
high level, abstract quantities that physicists are used to in
their mathematics. To illustrate how one might begin to express
more sophisticated mathematics in programming, let us define a new
private_complex type and the procedures which will operate on that
type. This is, of course, an academic exercise for Fortran
programmers, since the complex type already exists in the
language. Nevertheless, it is a useful example to illustrate the
basic principles involved and will lead to our definition of
classes. This type is defined as follows:
type private_complex
real :: real, imaginary
end type private_complex
To create variables a, b, and c of this new type, and assign
values, one proceeds as before:
type (private_complex) :: a, b, c
! assign values to a
a%real = 1.0
a%imaginary = 2.0
If this private_complex type behaves the same as ordinary complex
numbers, then multiplication of c = a*b can be defined as
follows:
c%real = a%real*b%real - a%imaginary*b%imaginary
c%imaginary = a%real*b%imaginary + a%imaginary*b%real
A new function pc_mult to multiply private_complex numbers could
then be written:
type (private_complex) function pc_mult(a,b)
type (private_complex), intent(in) :: a, b
pc_mult%real = a%real*b%real - a%imaginary*b%imaginary
pc_mult%imaginary = a%real*b%imaginary + a%imaginary*b%real
end function pc_mult
Note that this function returns a variable of type
private_complex. One can thus multiply two numbers of this type
with the following statement:
c = pc_mult(a,b)
It makes sense to place a new derived type together with the
procedures which operate on that type into the same module:
module private_complex_module
! define private_complex type
type private_complex
real :: real, imaginary
end type private_complex
contains
type (private_complex) function pc_mult(a,b)
! multiply private_complex variables
type (private_complex), intent(in) :: a, b
pc_mult%real = a%real*b%real - a%imaginary*b%imaginary
pc_mult%imaginary = a%real*b%imaginary + a%imaginary*b%real
end function pc_mult
end module private_complex_module
A program to illustrate the multiplication of two private_complex
numbers then looks like the following:
program main
! bring in private_complex definition and procedures
use private_complex_module
! define sample variables
type (private_complex):: a, b, c
! initialize sample variables
a%real = 1. ; a%imaginary = -1.
b%real = -1. ; b%imaginary = 2.
! perform multiplication
c = pc_mult(a,b)
print *,'c=', c%real, c%imaginary
stop
end program main
It is also possible to encapsulate the individual components
of a derived type. This is a common and useful practice in
object-oriented programming. It means that when a module is
"used" in another program unit, the private_complex type can be
defined, but the individual components, such as a%real or
a%imaginary are not accessible. In the sample program above, the
individual components were accessed in initializing the data and
in printing the result of the multiplication. If the components
are encapsulated, then additional procedures would have to be
provided in the module to perform this function. The
encapsulation is achieved by adding the PRIVATE attribute to the
derived type definition as follows:
type private_complex
private
real :: real, imaginary
end type private_complex
A procedure to initialize a private_complex number from real
numbers can be written as follows:
subroutine pc_init(a,real,imaginary)
! initialize private_complex variable from reals
type (private_complex), intent(out) :: a
real, intent(in) :: real, imaginary
a%real = real
a%imaginary = imaginary
end subroutine pc_init
while one to display the contents can be written:
subroutine pc_display(a,c)
! display value of private_complex variable with label
type (private_complex), intent(in) :: a
character*(*), intent(in) :: c
print *, c, a%real, a%imaginary
end subroutine pc_display
The main program then looks like the following:
program main
use private_complex_module
type (private_complex) :: a, b, c
! initialize sample variables
call pc_init(a,1.,-1.)
call pc_init(b,-1.,2.)
! perform multiplication
c = pc_mult(a,b)
! display result
call pc_display(c,'c=')
stop
end program main
The advantage of such encapsulation is that procedures in other modules can never impact the internal representation of the private_complex type. Furthermore, any changes made to the internal representation of private_complex type would be confined to this module, and would not impact program units in other modules. This makes it easier to develop software with interchangeable parts.
We have seen earlier how functions can be overloaded. In
Fortran90, operators such as '*' can also be overloaded. This is
also done with the INTERFACE statement, which is placed in the
declaration section of the module, as follows:
interface operator(*)
module procedure pc_mult
end interface
We have now equated the operator '*' with the name pc_mult.
Thus in the main program, one can multiply two private_complex
numbers using the more familiar syntax:
c = a*b
If one adds the declaration:
private :: pc_mult
to the module, one can also make the original name pc_mult no longer accessible.
In the language of object-oriented programming, the module we have just created is known as a class. It consists of a derived type definition, known as a class name, along with the procedures which operate on that class, called class member functions. The components of the derived type are called the class data members, while global data in the module (if any) corresponds to static class data members. The actual variable of type private_complex is known as an object.
To make this appear more familiar to those who already know C++, we will adopt the convention to make the derived type the first argument in all the module procedures, and we will give it the name "this." We will also overload the initialization function pc_init and give it the public name "new," while making the specific name pc_init private. The final version of the private_complex class is listed in Appendix C. In this version we have provided a new public name "display" to the procedure pc_display, but we have not made it private, so both names can still be used.
The main program which uses this class can be written:
program main
use private_complex_class
type (private_complex) :: a, b, c
call new(a,1.,-1.)
call new(b,-1.,2.)
! perform multiplication
c = a*b
call pc_display(c,'c=')
stop
end program main
We will begin with the most general case, where one class makes use of another. Such a form of inheritance is sometimes called class composition, since classes are composed of other classes. To illustrate this, let us create a new class which consists of private_complex arrays. Array syntax for intrinsic data types is a powerful new feature of Fortran90, and it is instructive to see how one can implement it for user defined types. Fortran90 permits one to define arrays of derived types as follows:
program main
use private_complex_class
integer, parameter :: nx = 8
type (private_complex), dimension(nx) :: f
We can initialize the elements of this array type by calling in a
loop the "new" procedure we defined in the private_complex_class:
do i = 1, nx
call new(f(i),real(i),-real(i))
enddo
Note that the do loop without statement numbers is now an official
part of Fortran90. In this way, we can create a procedure called
pcarray_init which will convert two real arrays into an array of
type (private_complex). We place this procedure in a new module
called complex_array_class, as follows:
module complex_array_class
! get private_complex type and its public procedures
use private_complex_class
contains
subroutine pcarray_init(this,real,imaginary)
! initialize private_complex array from real arrays
type (private_complex), dimension(:), intent(out) :: this
real, dimension(:), intent(in) :: real, imaginary
do i = 1, size(this)
! new here uses the pc_init procedure
call new(this(i),real(i),imaginary(i))
enddo
end subroutine pcarray_init
end module complex_array_class
In this example, the private_complex module is "used" (or
inherited) by the complex_array module. This means that all the
public entities in the first module (the base class) will be
available to the second module (the derived class). Since arrays
of derived type are considered a different type than a single
scalar of that type, we can include the following interface
statement in the module to overload the procedure new so that it
can also execute pcarray_init:
interface new
module procedure pcarray_init
end interface
Now, if the argument of new is a scalar of type private_complex,
then pc_init will be executed, but if the argument of new is an
array of type private_complex, then pcarray_init will be executed.
The main program which uses this new module looks like the
following:
program main
! bring in complex_array (and private_complex) procedures
use complex_array_class
integer, parameter :: nx = 8
! define a single scalar of type private_complex
type (private_complex) :: a
! define an array of type private_complex
type (private_complex), dimension(nx) :: f
! initialize a single scalar
call new(a,1.0,-1.0)
! initialize an array
call new(f,(/(real(i),i=1,nx)/),(/(-real(i),i=1,nx)/))
stop
end program main
In a similar fashion we can add a multiply function to the new
module which will multiply arrays. To do this we take advantage
of a new feature of Fortran90, which is that functions can return
entire arrays. This is done by declaring the function name to be
an array, and setting its dimension from another array in the
argument list with the SIZE intrinsic. Since the '*' operator for
scalars of type private_complex has already been defined in the
module private_complex_class, the array function is created by
calling this operator in a loop, as follows:
function pcarray_mult(this,b)
! multiply arrays of private_complex variables
type (private_complex), dimension(:), intent(in) :: this, b
! declare output array to be of same size as "this" array
type (private_complex), dimension(size(this)) ::pcarray_mult
do i = 1, size(this)
! this multiplication is actually performed by function pc_mult
pcarray_mult(i) = this(i)*b(i)
enddo
end function pcarray_mult
Finally, we can overload the '*' operator to call pcarray_mult
when the arguments are arrays of type private_complex, as follows:
interface operator(*)
module procedure pcarray_mult
end interface
In the final version of the derived class complex_array, which is
listed in Appendix D, we have also overloaded the
display function so that we can print out the array elements. In the
following main program, we can now multiply arrays of type
private_complex using array syntax:
program main
use complex_array_class
integer, parameter :: nx = 8
! define sample arrays
type (private_complex), dimension(nx) :: f, g, h
! initialize sample arrays
call new(f,(/(real(i),i=1,nx)/),(/(-real(i),i=1,nx)/))
call new(g,(/(-real(i),i=1,nx)/),(/(real(2*i),i=1,nx)/))
! perform multiplication of arrays
h = f*g
! display the first three elements of the results
call display(h(1:3),'h=')
stop
end program main
Thus we have constructed a derived class built upon the definitions and procedures in a base class using class composition. It was not necessary to construct a new derived type for arrays, since arrays of derived types are automatically supported in the language. Procedures which operate on the arrays, however, had to be constructed. This was done in two stages. First a new procedure for the derived class is written (for example, pcarray_mult) which executes the original base class operator ('*') in a loop. The original base class operator ('*') is then overloaded to include the new derived class procedure. Note that the base class was never modified during this process.
Inheritance is generally used to mean something more restricted than class composition. In this form of the inheritance relation, the base class contains the properties (procedures) which are common to a group of derived classes. Each derived class can modify or extend these procedures for its own needs if necessary.
As an example, consider the private_complex class discussed
in the previous section. Suppose we want to extend this class so
that it keeps track of the last operation performed. Such a
feature could be useful in debugging, for example. Except for the
additional feature of monitoring operations, we would like this
extended class to behave exactly like the private_complex class.
We can accomplish this by creating a new class called the
monitor_complex class. First we create a new derived type, as
follows:
type monitor_complex
type (private_complex) :: pc
character*8 :: last_op
end type monitor_complex
which contains one instance of a private_complex type plus an
additional character component to be used for monitoring. We want
to extend all three procedures of the private_complex class
(new,'*',and display) so that they also work in the new
monitor_complex class. We will accomplish this with methods very
similar to those we used in composition. Initializing the
monitor_complex type can be performed by making use of the new
operator in the private_complex class to initialize the
private_complex component of the monitor_complex type as follows:
type (monitor_complex) :: x
call new(x%pc,1.0,-1.0)
The additional monitor component can be initialized in the usual
way:
x%last_op = 'INIT'
This initialization can be embedded in a procedure called mc_init.
By also adding the interface new to the monitor_complex class, we
can overload the new procedure so that it calls mc_init if the
argument is of type monitor_complex. As a result, the operator
new which previously worked on private_complex type has been
extended to also work on monitor_complex types, as required. The
resulting monitor_complex class looks like:
module monitor_complex_class
! get (inherit) private_complex type and its public procedures
use private_complex_class
! define monitor_complex type
type monitor_complex
private
type (private_complex) :: pc
character*8 :: last_op
end type monitor_complex
interface new
module procedure mc_init
end interface
contains
subroutine mc_init(this,real,imaginary)
! initialize monitor_complex variable from reals
type (monitor_complex), intent(out) :: this
real, intent(in) :: real, imaginary
! initialize private_complex component of monitor_complex
call new(this%pc,real,imaginary)
this%last_op = 'INIT' ! set last operation
end subroutine mc_init
end module monitor_complex_class
In a similar fashion we can extend the multiplication
function to also work on monitor_complex types. Multiplication is
performed by using the '*' operator (defined in the
private_complex class) on the private_complex component of the
monitor_complex type, as follows:
type (monitor_complex) :: x, y, z
z%pc = x%pc*y%pc
This operation can be embedded in a function which returns a
result of type monitor_complex. We also want to add to this
function the operation of setting the monitor component. The
resulting procedure can be written:
type (monitor_complex) function mc_mult(this,b)
type (monitor_complex), intent(in) :: this, b
mc_mult%pc = this%pc*b%pc
mc_mult%last_op = 'MULTIPLY'
end function mc_mult
Finally, we can overload the '*' operator with the interface statement
so that it calls mc_mult when the arguments are of type
monitor_complex. In the final version of the derived class
monitor_complex, which is listed in Appendix E,
we have extended the display function to work with the monitor_complex
class, as well as written an entirely new procedure (last_op) to
display the last operation performed. In order to allow expressions
with mixed types, we have also written a conversion operator mc to
convert from private_complex to monitor_complex types. In the
following main program, all of the procedures in the base class have
been extended to work in the derived class:
program main
! bring in monitor_complex definition and procedures
use monitor_complex_class
! define sample variables
type (private_complex) :: a, b, c
type (monitor_complex) :: x, y, z
! initialize private_complex variables
call new(a,1.,-1.)
call new(b,-1.,2.)
! initialize monitor_complex variables
call new(x,1.,-1.)
call new(y,-1.,2.)
call new(z,0.,0.)
! perform multiplication with private_complex types
c = a*b
! perform multiplication with monitor_complex types
z = x*y
! display results
call display(c,'c=')
call display(z,'z=')
! perform multiplication with mixed types
z = mc(c)*z
! display last operation for monitor_complex type
call last_op(z,'z')
end program main
Thus we have constructed a derived class which has a special form of relationship to its base class. Here, the derived type definition in the derived class contains within it exactly one component of the derived type definition in the base class. In other words, each object of the derived class contains within it exactly one object of the base class. Furthermore all the procedures in the base class have been extended to also work in the derived class, although two have been internally modified, and one new one created. In addition, a conversion operator was supplied. Extending the procedures was accomplished by first writing a derived class procedure which called the base class procedure on the base class component, and then overloading the name to be the same as the base class procedure name. As in class composition, the base class was never modified during this process.
This form of inheritance is sometimes called subtyping because derived class objects can use base class procedures as if they were the same type. Although inheritance by subtyping is rather restrictive, it is also very convenient because special mechanisms exist in some languages (such as C++) to automatically extend unmodified base class procedures to derived types through automatic conversion of types. Fortran90 does not have such mechanisms and therefore inheritance by subtyping must be constructed explicitly, which can be cumbersome.
In many object-oriented languages, composition and subtyping are considered separately because they are implemented by two different language mechanisms, and composition is rarely considered a form of inheritance. However, in Fortran90, subtyping is implemented as a special case of composition and is constructed using similar methods, so that it makes sense here to consider the two ideas as related and describe them jointly. In our experience, composition is a more powerful idea which reflects more closely than subtyping how concepts are constructed in physics.
This example with private_complex types shows one of the reasons object-oriented programming is so powerful: if we decide to change the internal representation of private_complex to use polar coordinates instead of cartesian, we would have to modify only the procedures in the private_complex class to accommodate the change, while the derived classes would still work without modification. This example is perhaps academic. However, one can use the same techniques to create more powerful and interesting classes to represent other kinds of algebras. For example, one can create vector classes with procedures such as gradient operators, or tensor classes and associated procedures. One can then program at the same high level that one can do mathematics, with all the power and safety such abstractions give.
subroutine push1(part,fx,species_args,dt)
The argument species_args was an instance of type species_descriptor that contained certain parameters describing a group of charged particles. The array part contained the particle coordinates for that group. This works well and the subroutine interface is considerably simplified from the original version.
Nevertheless, it is possible to add even more safety and
simplicity. Clearly, the array part needs to have the object
species_args always present. This leads to a possible source of
error, since one can pass a descriptor which is inconsistent with
the particle data. It makes sense therefore to make a new derived
type which unites the two data types together. One might define
such a type as follows:
type species
type (species_descriptor) :: descriptor
real, dimension(idimp,nop) :: coordinates
end type species
If idimp and nop are parameters known at compile time, such a
definition is perfectly valid. However, it is not convenient,
since changing idimp or nop would require much of the code to be
recompiled. One might guess that it would be better to have an
allocatable array in the type definition, such as:
type species
type (species_descriptor) :: descriptor
real, dimension(:,:), allocatable :: coordinates
end type species
It turns out that this is invalid: allocatable arrays cannot be used in derived type definitions. There is an alternative, however, which is valid, namely pointers to arrays. In order to explain how this works, we have to make a digression and explain how Fortran90 pointers work.
Pointers in any language refer to the location in memory of
data, rather than the data itself. Fortran90 pointers are in
reality a special kind of object. Compared to pointers in other
languages, their use is greatly restricted in order to avoid the
kind of performance degradation common with indiscriminant use of
pointers. The programmer does not have access to the value of a
pointer nor is pointer arithmetic permitted. Instead, pointers
are really aliases to other data. Suppose we define two real
arrays, and two pointers to real arrays:
real, dimension(4), target :: a, b
real, dimension(:), pointer :: ptr1, ptr2
One can then associate the first pointer with the array a using
the '=>' operator as follows:
ptr1 => a
Notice that the TARGET attribute must be used for the arrays a and
b, in order to allow them to be pointed to. The kind of data with
which a pointer can be associated is determined in its
declaration. Thus, the pointer ptr1 can only point to real, one
dimensional arrays, and cannot point to a real, two dimensional
array, for example. Once associated with a target, the pointer
can then be used just as if it were the array a, including passing
it as an argument to a procedure. This is called dereferencing.
For example the statement
ptr1(1) = 1.0
will assign the value of 1.0 to a(1). The NULLIFY intrinsic is
used to disassociate a pointer from an array:
nullify(ptr1)
Similarly, if we associate the second pointer with b,
ptr2 => b
then the statement
ptr2(2) = 2.0
will assign the value of 2.0 to b(2). The ASSOCIATED intrinsic
function can be used to determine if a pointer has been associated
with any array:
if (associated(ptr1)) print *,'ptr1 associated!'
It can also be used to determine if two pointers point to the
same location in memory.
if (associated(ptr1,ptr2)) print *,'associated together!'
For our purposes, the most useful feature of pointers is that they
can be associated with an unnamed variable with the ALLOCATE
statement. For example, the statements
nullify(ptr1)
allocate(ptr1(4))
will first disassociate ptr1 from any previous array, then
allocate an unnamed array consisting of 4 real words and associate
the pointer ptr1 to that unnamed array. The pointer ptr1 can then
be used as if it were a normal array. When we are finished with
that array, we can deallocate it with:
deallocate(ptr1)
The ALLOCATED intrinsic does not work with pointers to arrays. However, the ASSOCIATED statement can tell us whether the data had actually been allocated, if we first NULLIFY ptr1 before allocating it. For all practical purposes, pointers to unnamed arrays function just like allocatable arrays, except that they have the advantage they can be used in derived type definitions.
Thus, the correct definition for the new species type is:
type species
type (species_descriptor) :: descriptor
real, dimension(:,:), pointer :: coordinates
end type species
If we define an object of type species called electrons,
type (species) :: electrons
then the pointer array for the particle data can be allocated as
follows:
allocate(electrons%coordinates(idimp,nop),stat=ierr)
The new push1 procedure can now be called with the simple
statement:
call push1(electrons,fx,dt)
This can be organized efficiently by creating two classes. First, we create a species_descriptor class (see Appendix F) which contains the type definition for this class along with procedures to read and write objects of this class. This class is then "inherited" by a derived class called species (see Appendix G), which uses the base class information to define a species type along with a creation procedure and a new push1 subroutine.
There is one subtle issue that occurs when using derived
types containing pointers. When two such types are copied:
type (species) :: electrons, ions
ions = electrons
what actually occurs is
ions%descriptor = electrons%descriptor
ions%coordinates => electrons%coordinates
In the second line, the pointer component is copied, not the data
to which it points. Sometimes this is not the desired behavior. If
we want to copy the data instead of the pointer, we need to
perform the following operation instead:
ions%coordinates = electrons%coordinates
which will dereference the pointer and copy the data. To
accomplish this, we create a procedure:
subroutine copy_species(a,b)
type (species), intent(out) :: a
type (species), intent(in) :: b
a%descriptor = b%descriptor
a%coordinates = b%coordinates
end subroutine copy_species
Fortran90 allows such a procedure to be associated with the '='
operator by means of the interface statement in the module:
interface assignment (=)
module procedure copy_species
end interface
If we implement such a procedure, then the statement:
ions = electrons
will copy the data instead of the pointer.
To illustrate this, let us write a subroutine which does
some kind of work using the methods in our private_complex class
hierarchy. Since this inheritance hierarchy was rather simple,
our work subroutine will also be simple: it will square a number
and then print out the result:
subroutine work(a)
type (private_complex), intent(inout) :: a
a = a*a
call display(a,Õwork:Õ)
end subroutine work
We will define a display procedure so that it works differently in
each class (for example, by modifying the display procedure in the
monitor_complex class listed in Appendix E so that it calls
last_op). The work subroutine is written for the private_complex
type, but we would like it to function correctly even if we pass a
monitor_complex type instead (which would normally be a type
violation). In other words, when we say that this procedure works
on a private_complex type, we actually mean that it is supposed to
work on all types derived from private_complex as well.
Furthermore, we want to decide at run time which we mean. Thus, if
we declare and initialize the class objects as follows:
type (private_complex), pointer :: pc
type (private_complex), target :: a
type (monitor_complex), target :: b
! initialize private_complex variables
call new(a,1.,-1.)
! initialize monitor_complex variables
call new(b,1.,-1.)
we would like to do something like:
! point to private_complex variable
pc => a
call work(pc)
! point to monitor_complex variable
pc => b !error, type violation
call work(pc)
This is not possible in Fortran90, because the pointer pc can only
point to targets of type private_complex, so that the statement:
pc => b
is illegal. The solution is to first define a derived type which
contains a pointer to each of the possible types in the
inheritance hierarchy, such as:
type complex_subtype
type (private_complex), pointer :: pc
type (monitor_complex), pointer :: mc
end type complex_subtype
This subtype has the ability to point to either complex type:
type (complex_subtype) :: generic_pc
! point to private_complex variable
generic_pc%pc => a
! point to monitor_complex variable
generic_pc%mc => b
Dynamic dispatching can then be implemented by defining a
subtype class which inherits all the classes in the inheritance
hierarchy and contains a new version of each class member
function. This new version will test which of the pointers in the
complex_subtype type has been associated and execute the
corresponding version of the function. For example, one can
define a new display procedure as follows:
subroutine display_subtype(a,c)
type (complex_subtype), intent(in) :: a
character*(*), intent(in) :: c
! check if pointer is associated with private_complex type
if (associated(a%pc)) then
! if so, execute private_complex version of display
call display(a%pc,c)
! check if pointer is associated with monitor_complex type
elseif (associated(a%mc)) then
! if so, execute monitor_complex version of display
call display(a%mc,c)
endif
end subroutine display_subtype
The argument to this procedure is a variable of type
complex_subtype. The procedure hides the decision about which
actual function to call. In a similar fashion, one could write a
new multiplication function, which would hide the decisions about
the appropriate multiplication procedure to call. If one
overloads the procedure names to have the same names as the
corresponding class member functions, the resulting work procedure
then has exactly the same form as before, except that the argument
is of type complex_subtype rather than private_complex, to
indicate that this subroutine is intended to be used with the
entire class hierarchy:
subroutine work(a)
type (complex_subtype), intent(inout) :: a
! multiplication operator has been overloaded to cover all types
a = a*a
call display(a,'work:')
end subroutine work
From the above, it is clear that the subtype class will need
an assignment operator. Since we will make the components of the
complex_subtype class private, one has to write a procedure to
dynamically assign one of the possible pointers in the class
hierarchy and nullify the rest. For example, the assignment of
the private_complex pointer would be performed by:
subroutine assign_pc(cs,pc)
type (complex_subtype), intent(out) :: cs
type (private_complex), target, intent(in) :: pc
! assign private_complex to complex_subtype
cs%pc => pc
! nullify monitor_complex pointer
nullify(cs%mc)
end subroutine assign_pc
A similar procedure, which we call assign_mc, must be created to
assign the monitor_complex pointer. Finally, in order for the
following expression:
a = a*a
to produce the expected results, one also needs to define a new copy operator, copy_subtype, which copies data rather than pointers, similar to the one we defined at the end of the section VI. These assignment and copy operators can be overloaded to the assignment operator ('='). The complex_subtype class which combines all these features is shown in Appendix H. It encapsulates all the details of how dynamic dispatching works, so that users of this class can freely write abstract or generic procedures based on the class hierarchy without being aware of how dynamic dispatching is implemented.
The program which calls work looks like:
program main
use complex_subtype_class
type (complex_subtype) :: pc
type (private_complex), target :: a
type (monitor_complex), target :: b
call new(a,1.,-1.)
call new(b,1.,-1.)
call assign(pc,a)
call work(pc)
call assign(pc,b)
call work(pc)
end program main
A distinguishing feature of typed object-oriented languages, is that support for the equivalent of such a subtype class is supported automatically by the compiler. There is always a performance penalty for using dynamic dispatching, however, even in object-oriented languages.
The rules for implementing a subtype class in Fortran90 which supports dynamic dispatching are as follows: create a derived type which contains exactly one pointer to each possible class in the inheritance hierarchy. Then implement a generic method for each class member function which will test which of the possible pointers have been associated and pass the corresponding pointer to the appropriate function. Assignment procedures are also required. The users of the subtype class, however, do not need to concern themselves with the details of how it is constructed.
Clearly, it is more cumbersome to implement dynamic dispatching in Fortran90 than in an object-oriented language. Once implemented, however, the usage is similar. Whether this feature is important depends on the nature of the problem one wants to model. Some studies of object-oriented programming [9] indicate that dynamic dispatching is needed about 15% of the time. There is some debate about what constitutes object-oriented programming. Some would argue that if a program does not use dynamic dispatching, it is not object-oriented. Others, however, argue that object-oriented is a question of design, not a question of what features of an object-oriented language are used. We tend to agree with the latter view.
module fft1r_module
integer, save, private :: saved_indx = -1
integer, dimension(:), allocatable, save, private :: mixup
complex, dimension(:), allocatable, save, private :: sct
contains
subroutine fft1r_init(indx)
! initialization call
integer, intent(in) :: indx
integer :: nx, nxh, ierr, isign=0
! allocate f and t: old, ugly fft requires them as arguments
real, dimension(2**indx) :: f
complex, dimension(2**(indx-1)) :: t
if (indx.lt.0) then
print *,'indx must be non-negative'
stop
endif
nx = 2**indx ; nxh = nx/2 ; saved_indx = indx
if (allocated(mixup)) deallocate(mixup)
if (allocated(sct)) deallocate(sct)
allocate(mixup(nxh),sct(nxh),stat=ierr)
if (ierr.ne.0) then
print *,'allocation error'
stop
endif
! call old, ugly but fast fft here
call lib_fft1r(f,t,isign,mixup,sct,saved_indx,nx,nxh)
end subroutine fft1r_init
!
subroutine fft1r_end
! deallocate internal data
saved_indx = -1
deallocate(mixup,sct)
end subroutine fft1r_end
!
subroutine fft1r(f,isign)
real, dimension(:), intent(inout) :: f
integer, intent(in) :: isign
integer :: nx, nxh
complex, dimension(size(f)/2) :: t
nx = size(f) ; nxh = nx/2
! do nothing if isign is invalid
if (isign.eq.0) return
if (saved_indx.lt.0) then
print *,'fft tables not initialized!'
stop
endif
! call old, ugly but fast fft here
call lib_fft1r(f,t,isign,mixup,sct,saved_indx,nx,nxh)
end subroutine fft1r
end module fft1r_module
module species_module
! define derived type
type species_descriptor
integer :: number_of_particles
real :: charge, charge_to_mass, kinetic_energy
end type species_descriptor
contains
subroutine push1(part,fx,species_args,dt)
! declare assumed-shape arrays
real, dimension(:,:), intent(inout) :: part
real, dimension(:), intent(in) :: fx
! declare argument of derived type
type (species_descriptor), intent(inout) :: species_args
real, intent(in) :: dt
integer :: np, idimp, nop, nx
real :: qbm, ek
! extract array sizes
idimp = size(part,1) ; nop = size(part,2) ; nx = size(fx)
! unpack input scalars from derived type
qbm = species_args%charge_to_mass
np = species_args%number_of_particles
! call old, ugly but fast particle pusher here
call orig_push1(part,fx,qbm,dt,ek,np,idimp,nop,nx)
! pack output scalars into derived type
species_args%kinetic_energy = ek
end subroutine push1
end module species_module
module private_complex_class
private :: pc_init, pc_mult
type private_complex
private
real :: real, imaginary
end type private_complex
interface new
module procedure pc_init
end interface
interface operator(*)
module procedure pc_mult
end interface
interface display
module procedure pc_display
end interface
contains
subroutine pc_init(this,real,imaginary)
! initialize private_complex variable from reals
type (private_complex), intent(out) :: this
real, intent(in) :: real, imaginary
this%real = real
this%imaginary = imaginary
end subroutine pc_init
!
type (private_complex) function pc_mult(this,b)
! multiply private_complex variables
type (private_complex), intent(in) :: this, b
pc_mult%real = this%real*b%real -
&this%imaginary*b%imaginary
pc_mult%imaginary = this%real*b%imaginary +
&this%imaginary*b%real
end function pc_mult
!
subroutine pc_display(this,c)
! display value of private_complex variable with optional label
type (private_complex), intent(in) :: this
character*(*), intent(in), optional :: c
if (present(c)) then
print *, c, this%real, this%imaginary
else
write (6,'(2f14,7)',advance='no') this%real,
&this%imaginary
endif
end subroutine pc_display
end module private_complex_class
module complex_array_class
use private_complex_class
private :: pcarray_init, pcarray_mult
interface new
module procedure pcarray_init
end interface
interface operator(*)
module procedure pcarray_mult
end interface
interface display
module procedure pcarray_display
end interface
contains
subroutine pcarray_init(this,real,imaginary)
! initialize private_complex variable from reals
type (private_complex), dimension(:), intent(out) :: this
real, dimension (:), intent(in) :: real, imaginary
do i = 1, size(this)
call new(this(i),real(i),imaginary(i))
enddo
end subroutine pcarray_init
!
function pcarray_mult(this,b)
! multiply arrays of private_complex variables
type (private_complex), dimension(:),intent(in) :: this,b
type (private_complex), dimension(size(this))::
&pcarray_mult
! this multiplication is actually defined by function pc_mult
do i = 1, size(this)
pcarray_mult(i) = this(i)*b(i)
enddo
end function pcarray_mult
!
subroutine pcarray_display(this,c)
! display value of private_complex array with label
type (private_complex), dimension(:), intent(in) :: this
character*(*), intent(in) :: c
write (6,'(a)',advance='no') c
do i = 1, size(this)
call display(this(i))
enddo
print *
end subroutine pcarray_display
end module complex_array_class
module monitor_complex_class
! get (inherit) private_complex type and its public procedures
use private_complex_class
private :: mc_init, mc_mult, mc_display
! define monitor_complex type
type monitor_complex
private
type (private_complex) :: pc
character*8 :: last_op
end type monitor_complex
interface new
module procedure mc_init
end interface
interface operator(*)
module procedure mc_mult
end interface
interface display
module procedure mc_display
end interface
contains
subroutine mc_init(this,real,imaginary)
! initialize monitor_complex variable from reals
type (monitor_complex), intent(out) :: this
real, intent(in) :: real, imaginary
! initialize private_complex component of monitor_complex
call new(this%pc,real,imaginary)
! set last operation
this%last_op = 'INIT'
end subroutine mc_init
!
type (monitor_complex) function mc_mult(this,b)
! multiply monitor_complex variables
type (monitor_complex), intent(in) :: this, b
! this multiplication is actually defined by function pc_mult
mc_mult%pc = this%pc*b%pc
! set last operation
mc_mult%last_op = 'MULTIPLY'
end function mc_mult
!
subroutine mc_display(this,c)
! display value of monitor_complex variable with label
type (monitor_complex), intent(in) :: this
character*(*), intent(in), optional :: c
call display(this%pc,c)
end subroutine mc_display
!
subroutine last_op(this,c)
! display last operation
type (monitor_complex), intent(in) :: this
character*(*), intent(in) :: c
print *, 'last op for ', c, ' was ', this%last_op
end subroutine last_op
!
type (monitor_complex) function mc(pc)
! convert private_complex object to monitor_complex object
type (private_complex), intent(in) :: pc
mc%pc = pc
mc%last_op = 'INIT'
end function mc
end module monitor_complex_class
module species_descriptor_class
type species_descriptor
private
integer :: number_of_particles
real :: charge, charge_to_mass, kinetic_energy
end type species_descriptor
contains
subroutine get_species(this,np,qm,qbm,ek)
! unpack components of species descriptor
implicit none
type (species_descriptor), intent(in) :: this
integer, intent(out), optional :: np
real, intent(out), optional :: qm, qbm, ek
if (present(np)) np = this%number_of_particles
if (present(qm)) qm = this%charge
if (present(qbm)) qbm = this%charge_to_mass
if (present(ek)) ek = this%kinetic_energy
end subroutine get_species
!
subroutine put_species(this,np,qm,qbm,ek)
! pack components of species descriptor
implicit none
type (species_descriptor), intent(out) :: this
integer, intent(in), optional :: np
real, intent(in), optional :: qm, qbm, ek
if (present(np)) this%number_of_particles = np
if (present(qm)) this%charge = qm
if (present(qbm)) this%charge_to_mass = qbm
if (present(ek)) this%kinetic_energy = ek
end subroutine put_species
end module species_descriptor_class
module species_class
! inherit species descriptor class
use species_descriptor_class
type species
private
type (species_descriptor) :: descriptor
real, dimension(:,:), pointer :: coordinates
end type species
interface new
module procedure species_init
end interface
contains
subroutine species_init(this,species_args,idimp,nop)
! allocate particle coordinate pointer array and store descriptor
type (species), intent(inout) :: this
type (species_descriptor), intent(in) :: species_args
integer, intent(in) :: idimp, nop
integer :: ierr
! allocate pointer array
allocate(this%coordinates(idimp,nop),stat=ierr)
! check for allocation error
if (ierr.ne.0) then
print *,'species allocation error'
! store descriptor
else
this%descriptor = species_args
endif
end subroutine species_init
!
subroutine push1(this,fx,dt)
type (species), intent(inout) :: this
real, dimension (:), intent(in) :: fx
real, intent(in) :: dt
integer :: np, idimp, nop, nx
real :: qbm, ek
! extract array sizes
idimp = size(this%coordinates,1)
nop = size(this%coordinates,2)
nx = size(fx)
! unpack input scalars from derived type with inherited procedure
call get_species(this%descriptor,np=np,qbm=qbm)
! call old, ugly but fast particle pusher here
call orig_push1(this%coordinates,fx,qbm,dt,ek,np,idimp,
&nop,nx)
! pack output scalar into derived type with inherited procedure
call put_species(this%descriptor,ek=ek)
end subroutine push1
end module species_class
module complex_subtype_class
! get (inherit) private_complex type and its public procedures
use private_complex_class
! get (inherit) monitor_complex type and its public procedures
use monitor_complex_class
private
public :: private_complex, monitor_complex, complex_subtype
public :: new, assign, assignment(=), operator(*), display
! define complex_subtype type
type complex_subtype
private
type (private_complex), pointer :: pc
type (monitor_complex), pointer :: mc
end type complex_subtype
interface assign
module procedure assign_pc
module procedure assign_mc
end interface
interface assignment (=)
module procedure copy_subtype
end interface
interface operator(*)
module procedure mult_subtype
end interface
interface display
module procedure display_subtype
end interface
contains
subroutine assign_pc(cs,pc)
! assign private_complex to complex_subtype
type (complex_subtype), intent(out) :: cs
type (private_complex), target, intent(in) :: pc
cs%pc => pc
nullify(cs%mc)
end subroutine assign_pc
!
subroutine assign_mc(cs,mc)
! assign monitor_complex to complex_subtype
type (complex_subtype), intent(out) :: cs
type (monitor_complex), target, intent(in) :: mc
nullify(cs%pc)
cs%mc => mc
end subroutine assign_mc
!
subroutine copy_subtype(this,b)
! assign contents of complex_subtype to complex_subtype
type (complex_subtype), intent(inout) :: this
type (complex_subtype), intent(in) :: b
! check if pointer is associated with private_complex type
if (associated(b%pc)) then
this%pc = b%pc
nullify(this%mc)
! check if pointer is associated with monitor_complex type
elseif (associated(b%mc)) then
this%mc = b%mc
nullify(this%pc)
endif
end subroutine copy_subtype
!
function mult_subtype(this,b) result(output)
! multiply complex_subtype variables
type (complex_subtype), intent(in) :: this, b
type (complex_subtype) :: output
type (private_complex), target, save :: tpc
type (monitor_complex), target, save :: tmc
! check if pointer is associated with private_complex type
if (associated(this%pc)) then
tpc = this%pc*b%pc
output = tpc
! check if pointer is associated with monitor_complex type
elseif (associated(this%mc)) then
tmc = this%mc*b%mc
output = tmc
endif
end function mult_subtype
!
subroutine display_subtype(a,c)
! display value of complex_subtype variable with label
type (complex_subtype), intent(in) :: a
character*(*), intent(in) :: c
! check if pointer is associated with private_complex type
if (associated(a%pc)) then
call display(a%pc,c)
! check if pointer is associated with monitor_complex type
elseif (associated(a%mc)) then
call display(a%mc,c)
endif
end subroutine display_subtype
end module complex_subtype_class