Disclaimer
this is an informal talk by a poor Fortraner, do not aspect to learn a lot,
probably you already know more Fortran than me...
Fortran was originally developed by IBM (the FORmula TRANslating system) in the mid of 50’
Three main reasons to learn Fortran:
PROGRAM
is the same as ProGram
; program My_bAD
integer :: NabCD
parameter :: naBcd = 10
character(nabcD) :: abc = 'hello'
print*, ABC
program MY_bad
program my_good
integer, parameter :: chars_number=10
character(chars_number) :: greetings='hello'
print '(A)', trim(greetings)
program my_good
program
...
end program
module
...
end module
subroutine
...
end subroutine
function
...
end function
PROGRAM
unit is the main, there should be only one main unit;PROGRAM
unit can contain one or more subprogram units such as SUBROUTINE
or FUNCTION
;MODULE
: program the_greeter
print '(A)', 'Hello World!'
end program the_greeter
module the_greeter_worker
contains
subroutine say_something(what)
character(*), intent(in) :: what !< What I have to say?
print '(A)', what
end subroutine say_something
endmodule the_greeter_worker
program the_modular_greeter
use the_greeter_worker
call say_something(what='Hello modular World!')
end program the_modular_greeter
The standard Fortran file extension is .f90
(for standard 90+), e.g. the_greeter.f90
Building a Fortran application means to pass through 3 steps:
# generate the_greeter.o compiled object
gfortran -c the_greeter.f90
# generate the_greeter executable
gfortran the_greeter.o -o the_greeter
# generate the_greeter executable in one command
gfortran the_greeter.f90 -o the_greeter
The best Fortran compilers you will ever see are also FREE!
INTEGER
s are integer values, which can be represented exactly by the computerREAL
s are floating point numbers, representing real numbers; they are not exact representations of real numbers, having only a finite number of decimal places, the FINITE PRECISION issueCHARACTER
s are text characters, usually encoded as ASCIILOGICAL
s are boolean variables, which can store a value of either .TRUE.
or .FALSE.
COMPLEX
variables are effectively just two real variables, one storing the real component of a number, the other storing the imaginary componentFortran is capable of implicit typing of variables, i.e. variables obtain a type on the basis of ITS FIRST CHARCTER… a nightmare!
|ABCDEFGH|IJKLMN |OPQRSTUVWXYZ|
|________|_______|____________|
| real |integer| real |
implicit none
to disable the default implicity program explicit
implicit none
end program explicit
module explicit
implicit none
end module explicit
subroutine explicit
! if not contained into
! a program/module with
! "implicit none"
implicit none
end subroutine explicit
TYPE[, attributes] [::] variable
where TYPE
is an intrinsic type (or a derived one not yet seen), attributes
are optional attributes of the variable declaration and variables
is the, obviously, the variable name. The double semicolon ::
are optional, but highly recommended!
It is possible and highly recommended to assign and initial value to all variables. It could be very helpful, to initialize them into their definition
integer, parameter :: cars_number=100
The assignment on definition is NOT ALWAYS POSSIBLE, neither DESIRABLE.
Constants are classified into two categories:
100
3.14159265
'Hello World!'
-3.134324390984353454e23 ! subtle mistakes...
Avoid to use literal constants: MAGIC NUMBERS are not clear!
integer, parameter :: cars_number=100
real, parameter :: pi=3.14159265
character(12), parameter :: greetings='Hello World!'
real, parameter :: wrong=-3.134324390984353454e23 ! subtle mistakes...
The PARAMETER
attribute is your CONSTANT friend!
Arrays are declared or by means of dimension attribute, either by directly define a dimension along-side the array name, e.g.
integer, dimension(2:15) :: array_with_explicit_bounds
integer, dimension(11) :: array_with_implicit_bounds
integer :: array_with_explicit_dimensions(3)
integer, dimension(5, 3) :: array ! avoid magic numbers
integer :: c
integer :: r
do r=1, 5
do c=1, 3
print*, array(r, c) ! wrong memory access
enddo
enddo
integer, parameter :: rows_number=5
integer, parameter :: cols_number=3
integer, dimension(rows_number, cols_number) :: array
integer :: c
integer :: r
do c=1, cols_number
do r=1, rows_number
print*, array(r, c) ! right column-major memory access
enddo
enddo
real :: a(0:20), b(3, 0:5, -10:10)
rank
: number of dimensions; `a` has rank 1 and `b` has rank 3
bounds
: upper and lower limits of each dimension of the array; `a` has bounds 0:20 and `b` has bounds 1:3, 0:5 and -10:10
extent
: number of element in each dimension; `a` has extent 21 and `b` has extents 3,6 and 21
size
: total number of elements; `a` has size 21 and `b` has 30
shape
: the shape of an array is its rank and extent; `a` has shape 21 and `b` has shape (3, 6, 21)
Note that
c(4:6) = d(1:3)
The arrays handling is superb in Fortran, computations on them are very efficient
integer, parameter :: elems_number=10
real, parameter :: delta=1.0/elems_number
real, dimension(0:elems_number) :: abscissa
real, dimension(0:elems_number) :: temperature
integer :: e
! array computations
! array constructor
abscissa = [(delta*e, e=0, elems_number)]
! elemental function operating on each array element
temperature = sin(abscissa)/2.0
print "(A)", "Abscissa, Temperature"
do e=0, elems_number
print "(E15.6, 1X, E15.6)", abscissa(e), temperature(e)
enddo
Arrays can have static defined sized (at compile time) or a dynamic size (at run time). Dynamic allocation allow efficient handling of memory with a very little (negligible or null) overhead: dynamic arrays allocation is not the evil!
integer :: elems_number
real :: delta
real, allocatable, dimension(:) :: abscissa
real, allocatable, dimension(:) :: temperature
integer :: e
! input from user
print "(A)", "Enter the number of elements"
read *, elems_number
delta=1.0/elems_number
! allocate arrays
allocate(abscissa(0:elems_number))
allocate(temperature(0:elems_number))
! array computations
! array constructor
abscissa = [(delta*e, e=0, elems_number)]
! elemental function operating on each array element
temperature = sin(abscissa)/2.0
print "(A)", "Abscissa, Temperature"
do e=0, elems_number
print "(E15.6, 1X, E15.6)", abscissa(e), temperature(e)
enddo
real, parameter :: pi_greek=3.14159265358979323 ! you will not get what you expect...
real*8, parameter :: pi_greek=3.14159265358979323 ! still wrong...
real*16, parameter :: pi_greek=3.14159265358979323 ! no improvements, still wrong...
real*16, parameter :: pi_greek=3.14159265358979323_16 ! ok, sometimes works, other fails...
USE KIND SPECIFICATION!
integer(kind=...) :: a_parametric_integer
real(kind=...) :: a_parametric_real
character(kind=..., len=...) :: a_parametric_character
logical(kind=...) :: a_parametric_logical
complex(kind=...) :: a_parametric_complex
The default value of the kind of each intrinsic type is PROCESSOR DEPENDENT: do not expect to obtain the same default kind on all OS/Compiler/CPU combinations (hereafter architecture).
Which Parameter I have to use?
integer(kind=4) :: a_bad_parametric_integer
real(kind=4) :: a_bad_parametric_real
...
The definition KIND=4
is not standard portable definition, it does not mean anything on which you can rely on: on some (the most part) architectures it means 4-bytes, but it is not ensured to happen.
By means of the built-in functions selected_xxx_kind
! parametric portable kinds
integer, parameter :: I4 = selected_int_kind(9) ! 4-bytes
integer, parameter :: R8 = selected_real_kind(15,307) ! 8-bytes
integer, parameter :: CK = selected_char_kind(name='ASCII') ! ascii kind
! parametrized portable definition of variables
integer(kind=I4) :: a_good_parametric_integer
real(kind=R8) :: a_good_parametric_real
character(kind=CK, len=2) :: a_good_parametric_character
By means of the intrinsic module (F03+) iso_fortran_env
! parametric portable kinds
use, intrinsic :: iso_fortran_env, only : I4 => int32, R8 => real64, character_kinds
integer, parameter :: CK = character_kinds(1)
! parametrized portable definition of variables
integer(kind=I4) :: a_good_parametric_integer
real(kind=R8) :: a_good_parametric_real
character(kind=CK, len=2) :: a_good_parametric_character
use, intrinsic :: iso_fortran_env, only : R8 => real64
implicit none
real, parameter :: pi_greek_1=3.14159265358979323
real*8, parameter :: pi_greek_2=3.14159265358979323
real(kind=R8), parameter :: pi_greek_3=atan(1.)*4.
real(kind=R8), parameter :: pi_greek_4=atan(1._R8)*4._R8
print '(E23.15E3)', pi_greek_1
print '(E23.15E3)', pi_greek_2
print '(E23.15E3)', pi_greek_3
print '(E23.15E3)', pi_greek_4
What you expect? The output could be surprising…
0.314159274101257E+001
0.314159274101257E+001
0.314159274101257E+001
0.314159265358979E+001
+
: addition-
: subtraction*
: multiplication/
: division**
: exponentiation==
: equal to/=
: not equal to<
: less than<=
: less than or equal to>
: greater than>=
: greater than or equal to.and.
: intersection.or.
: union.not.
: negation.eqv.
: logical equivalence.neqv.
: exclusive or//
: concatenationThe King of all conditionals, the if
construct. The minimal usage
if (expression) statement
where expression
is a valid logical expression. If the expression is .true.
(when evaluated), the statement
is executed. A slightly more complex usage
[name_of_if:] if (expression1) then
statement1
else if (expression2) then
statement2
else
statement3
end if [name_of_if]
if-then-else
construct allows branching of algorithms;expression1, epression2,...
) are executed sub-sequentially;statement
is executed the if-then-else
is left (i.e. the control goes to end if
);else
is executed (if present) if all the previous expression are .false.
.Often if-then-else
construct is not the optimal branching control, select case
construct could be more efficient
[name_of_case:] select case(expression)
case(selectors1)
statement1
case(selectors2)
statement2
case default
statement3
end select [name_of_case]
where expression
is a generic expression (not only of logical type), selectors1, selectors2,...
are list of (constant) values or range of values.
case(2:11)
means for all cases where expression
is in between [2,11];expression
is evaluated only one time! (efficient branching);case default
is executed (if present) if all the previous cases are not matched.The minimal definite loop is of the form
[name_of_do:] do counter=lower_bound, upper_bound [, stride]
end do [name_of_do]
where counter
is an integer variable, lower/upper_bound
are the loop bounds and the optional stride
is the stride used (default 1). If the stride is negative a count-down is performed, e.g.
[name_of_do:] do counter=upper_bound, lower_bound, negative_stride
end do [name_of_do]
[name_of_do:] do while (expression)
end do [name_of_do]
where expression
is evaluated at the beginning of each loop step and if .true.
the statements are executed.
It is possible, and often desirable, to loop indefinitely at compile time and check for the occurrences of some conditions at run time
[name_of_do:] do
end do [name_of_do]
The loop execution is controlled by means of the intrinsic function exit
(to break the loop) and cycle
to continue the loop with the next step
indefinite do
if (expression1) exit indefinite
if (expression2) cycle indefinite
end do indefinite
If the name of the loop is omitted, the loop exited or cycled is the nearest.
read
, Output => print, write
;open
, close
;inquire
;read
has a long list of arguments…
Formatted
READ (eunit, format [, advance] [, asynchronous]
[, blank] [, decimal] [, id] [, pad] [, pos]
[, round] [, size] [, iostat] [, err] [, end]
[, eor] [, iomsg]) [io-list]
Formatted, List-Directed
READ (eunit, *[, asynchronous] [, blank]
[, decimal] [, id] [, pad] [, pos]
[, round] [, size] [, iostat] [, err] [, end]
[, iomsg]) [io-list]
Formatted, Namelist
READ (eunit, nml-group[, asynchronous] [, blank]
[, decimal] [, id] [, pad] [, pos] [, round]
[, size] [, iostat] [, err] [, end] [, iomsg])
Unformatted
READ (eunit [, asynchronous] [, id] [, pos]
[, iostat] [, err][, end] [, iomsg]) [io-list]
Formatted
READ (eunit, format, rec [, asynchronous] [, blank]
[, decimal] [, id] [, pad] [, pos] [, round]
[, size] [, iostat] [, err] [, iomsg])
[io-list]
Unformatted
READ (eunit, rec [, asynchronous] [, id] [, pos]
[, iostat] [, err] [, iomsg]) [io-list]
READ (iunit, format [, nml-group] [, iostat] [, err]
[, end] [, iomsg]) [io-list]
Namelist
READ (iunit, nml-group [, iostat] [, err] [, end]
[, iomsg]) [io-list]
write
has a long list of arguments too…
Formatted
WRITE (eunit, format [, advance] [, asynchronous]
[, decimal] [, id] [, pos] [, round]
[, sign] [, iostat] [ , err] [, [iomsg])
[io-list]
Formatted, List-Directed
WRITE (eunit, * [, asynchronous] [, decimal]
[, delim] [, id] [, pos] [, round] [, sign]
[, iostat] [, err] [, iomsg])[io-list]
Formatted, Namelist
WRITE (eunit, nml-group [, asynchronous] [, decimal]
[, delim] [, id] [, pos] [, round] [, sign]
[, iostat] [, err] [, iomsg])
Unformatted
WRITE (eunit [, asynchronous] [, id] [, pos]
[, iostat] [, err] [, iomsg])[io-list]
Formatted
WRITE (eunit, format, rec [, asynchronous]
[, decimal] [, delim] [, id] [, pos]
[, round] [, sign] [, iostat] [, err]
[, iomsg])[io-list]
Unformatted
WRITE (eunit, rec [, asynchronous] [, id] [, pos]
[, iostat] [ , err] [, iomsg])[io-list]
WRITE (iunit, format [, nml-group] [, iostat]
[ , err] [, iomsg])[io-list]
Namelist
WRITE (iunit, nml-group [, iostat] [ , err]
[, iomsg])[io-list]
iso_fortran_env
read(5, *) user_input
if (error_occur) then
write(*, '(A)') error_message
else
write(6, '(A)') output_message
endif
use, intrinsic :: iso_fortran_env, only : error_unit, input_unit, output_unit
read(input_unit, *) user_input
if (error_occur) then
write(error_unit, '(A)') error_message
else
write(output_unit, '(A)') output_message
endif
open
, connect Fortran logical unit to a file or device; declare attributes for read and write operations;close
, terminate connection between logical unit and file or devicerewind
, position sequential or direct access file at the beginning of the file (the initial point);inquire
, request information on the status of specified properties of a file or logical unit. OPEN ([UNIT=] io-unit[, FILE= name] [, ERR= label]
[, IOMSG=msg-var] [, IOSTAT=i-var], slist)
where slist
is a list of one or more specifier, e.g. access='stream', form='unformatted'
etc…
CLOSE ([UNIT=] io-unit[, STATUS = p] [, ERR= label]
[, IOMSG=msg-var] [, IOSTAT=i-var])
REWIND ([UNIT=] io-unit[, ERR= label]
[, IOMSG=msg-var] [, IOSTAT=i-var])
Inquiring by File
INQUIRE (FILE=name[, ERR=label] [, ID=id-var]
[, IOMSG=msg-var] [, SIZE=sz]
[, IOSTAT=i-var] slist)
Inquiring by Unit
INQUIRE ([UNIT=]io-unit [, ERR=label] [, ID=id-var]
[, IOMSG=msg-var] [, SIZE=sz]
[, IOSTAT=i-var] slist)
Inquiring by Output List
INQUIRE (IOLENGTH=len) out-item-list
This allows to
do i=1, N
a(i) = cos(b(i-1)) + sin(c(i+1))
end do
do i=1, N
d(i) = cos(e(i-1)) + sin(f(i+1))
end do
do i=1, N
g(i) = cos(h(i-1)) + sin(p(i+1))
end do
subroutine do_something(a, b, d)
...
end subroutine do_something
call do_something(a=a, b=b, c=c)
call do_something(a=d, b=e, c=f)
call do_something(a=g, b=h, c=i)
[prefix [prefix]] subroutine subroutine_name [([d-arg-list]) [lang-binding]]
! modules access, if any
use module_1
...
implicit none
! variables (dummy arguments and local ones) definition
...
! executable statements
[contains]
! internal subprogram units
end subroutine subroutine_name
where
elemental, impure, module, pure, recursive
;BIND (C[, name=ext-name])
.Exchange data through its dummy arguments list or by IO (impure)
call subroutine_name [([d-arg-list])
[prefix [prefix]] FUNCTION function_name [([d-arg-list]) [suffix] [lang-binding]]
! modules access, if any
use module_1
...
implicit none
! variables (dummy arguments and local ones) definition
...
! executable statements
[contains]
! internal subprogram units
end function function_name
where
data type-specifier, elemental, impure, module, pure, recursive
;result(result_name)
;BIND (C[, name=ext-name])
.Exchange data through its returned value or by dummy arguments list (impure) or by IO (impure)
my_resutl = function_name [([d-arg-list])
intent(in)
intent(out)
intent(inout)
intent(in)
in a procedure cannot be changed during the execution of the procedure.present
intrinsic can be used to check for missing arguments subroutine get_temperature(vel, boltz, temperature)
use parameters
real(dp), dimension(:,:), intent(in) :: vel
real(dp), optional, intent(in) :: boltz
real(dp), intent(out) :: temperature
integer(ip) :: i
real(dp) :: ke
if (present(boltz)) then
kb = boltz
else
kb = boltz_parameter
endif
ke = 0._dp
do i = 1, natom
ke = ke + dot_product(vel(i,:), vel(i,:))
end do
temperature = mass * ke / (3._dp * kb * real(natom - 1, dp))
end subroutine get_temperature
! equivalent callings
call get_temperature(velocity, boltzman, temp)
call get_temperature(vel=velocity, boltz=boltzman, temperature=temp)
call get_temperature(vel=velocity, temperature=temp, boltz=boltzman)
call get_temperature(velocity, temperature=temp, boltz=boltzman)
Always prefer KEYWORD ARGUMENT style!
! a lot of this kind of callings
call baseline_sub(arg1=foo, arg2=bar, arg4=baz)
...
! modify the baseline_sub adding a new optional argument do not
! comport to modify all previous callings
! new callings
call baseline_sub(arg1=foo, arg2=bar, arg4=baz, arg3=fofo)
newunit=
in open statementg0
edit descriptorexecute_command_line
iso_fortran_env
: compiler_version
and compiler_options
module
is a program unit whose functionality can be exploited by other programs which attaches to it via the use
statement.COMMON
and INCLUDE
statementsmodule
, functions and subroutines are called module proceduresmodule
procedures can contain internal proceduresmodule
objects that retain their values should be given a save
attributemodule
s can be used by procedures and other modulesmodule
s can be compiled separatelyprivate
statementOld-fashioned Fortran library
subroutine foo(arg1, arg2, arg3)
real :: arg1, agr2(*), arg3
integer :: i
! GOD is real...
god = sum(arg2)
i = bar(arg2)
arg3 = god/i
end subroutine foo
integer function bar(arg)
real :: arg(*)
bar = size(arg, dim=1)
end function bar
Modern, robust Fortran module library
module library
implicit none
private
public :: foo
contains
pure subroutine foo(arg1, arg2, arg3)
real, intent(in) :: arg1
real, intent(in) :: agr2(*)
real, intent(out) :: arg3
real :: god
integer :: i
god = sum(arg2)
i = bar(arg2)
arg3 = god/i
end subroutine foo
pure function bar(arg) result(bar_size)
real, intent(in) :: arg(*)
integer :: bar_size
bar_size = size(arg, dim=1)
end function bar
end module library
Always start with Implicit none
statement
module library
! optional use-association of other modules
implicit none
end module library
Hide ALL! Expose only the minimum
module library
private ! all PRIVATE!
public :: what_really_concern ! expose only what really concern
protected :: read_only_objects
end module library
Rename entities on use-association for names-collision handling
module library
use another_library, only : external_concern => what_really_concern
private
public :: what_really_concern
end module library
module globals
implicit none
private
public :: foo, bar
protected :: baz
integer :: foo = 0
real :: bar = 0.0
logical :: baz = .false.
contains
subroutine initialize
foo = 1
bar = 2.0
baz = .true.
end subroutine initialize
end module globals
type :: line
real :: x(2) ! abscissa
real :: y(2) ! ordinate
end type line
type(line) :: a_line
integer :: Ni
integer :: Nj
integer :: Nk
real, allocatable :: centers(:,:,:)
real, allocatable :: nodes(:,:,:)
real, allocatable :: normals(:,:,:)
real, allocatable :: faces(:,:,:)
real, allocatable :: volumes(:,:,:)
call compute_metrics(Ni, Nj, Nk, &
centers,nodes, &
normals, fases &
volumes)
type :: mesh
integer :: Ni = 0
integer :: Nj = 0
integer :: Nk = 0
real, allocatable :: centers(:,:,:)
real, allocatable :: nodes(:,:,:)
real, allocatable :: normals(:,:,:)
real, allocatable :: faces(:,:,:)
real, allocatable :: volumes(:,:,:)
end type mesh
type(mesh) :: a_mesh
call compute_metrics(grid=a_mesh)
module globals
implicit none
integer :: Ni
integer :: Nj
integer :: Nk
real, allocatable :: centers(:,:,:)
real, allocatable :: nodes(:,:,:)
real, allocatable :: normals(:,:,:)
real, allocatable :: faces(:,:,:)
real, allocatable :: volumes(:,:,:)
end module globals
program bad
use globals
implicit none
call compute_metrics(Ni, Nj, Nk, &
centers,nodes, &
normals, fases &
volumes)
end program bad
module mesh_t
implicit none
private
public :: mesh
type :: mesh
integer :: Ni = 0
integer :: Nj = 0
integer :: Nk = 0
real, allocatable :: centers(:,:,:)
real, allocatable :: nodes(:,:,:)
real, allocatable :: normals(:,:,:)
real, allocatable :: faces(:,:,:)
real, allocatable :: volumes(:,:,:)
end type mesh
end module mesh_t
program good
use mesh_t, only: mesh
implicit none
type(mesh) :: a_mesh
call compute_metrics(grid=a_mesh)
end program good
Apply the same PROCEDURE to different arguments (for type, kind and or rank), e.g. convert number to string
integer :: step_number
real :: time
real :: delta
integer :: unit
integer :: t
delta = 0.1
time = 0.0
do t=1, 100
time = time + delta
open(newunit=unit, file=’clock-‘//string(t)//’.dat’)
write(unit, "(A)")’current time: ‘//string(time)
close(unit)
enddo
Say, you want the STRING
function to accept both the integer t
and the real time
, in the same of the many built-ins like cos, sin, etc...
You can develop a STRING_INTEGER, STRING_REAL
, but if you need to accept also different kinds? You then need a STRING_INTEGER_8, STRING_INTEGER_16, STRING_INTEGER_32...
and you need a different call for each argument type/kind/rank,
... file='clock-'//string_integer_32(t)//'.dat')
... 'current time: '//string_real_32(time)
Too much ERROR PRONE!
Generic Names avoid the plethora of different names exploiting dynamic dispatching at compile-time (no penalty!)
module generic
implicit none
private
public :: string
interface string
module procedure string_integer, string_real
endinterface
contains
elemental function string_integer(n) result(str)
integer, intent(in) :: n
character(11) :: str
write(str, ‘(I11)’) n
end function string_integer
elemental function string_real(n) result(str)
real, intent(in) :: n
character(13) :: str
write(str, ‘(E13.6E2)’) n
end function string_real
end module generic
Very CONCISE, CLEAR and ROBUST
use generic
integer :: step_number
real :: time
real :: delta
integer :: unit
integer :: t
delta = 0.1
time = 0.0
do t=1, 100
time = time + delta
open(newunit=unit, &
file='clock-'//string(t)//'.dat')
write(unit, "(A)") &
'current time: '//string(time)
close(unit)
enddo
A KISS library for exploiting codes portability for modern (2003+) Fortran projects
Create an algebra for your new shining type, e.g. for your vector class
type, public :: Vector
real :: x = 0.0 ! Cartesian component in x direction.
real :: y = 0.0 ! Cartesian component in y direction.
real :: z = 0.0 ! Cartesian component in z direction.
end type Vector
It would be fantastic to be able to do VECTORIAL COMPUTATION… isn’t it?
type(vector) :: vector1
type(vector) :: vector2
type(vector) :: vector3
...
vector3 = vector1 + vector2
vector3 = vector1 - vector2
vector3 = 2.0*vector1 + vector2
vector3 = vector1 - 1.4
vector3 = vector1.cross.vector2
vector3 = vector1.dot.vector2
Well, you can! Exploit OPERATORS OVERLOADING (without penalties!)
module vector_t
implicit none
private
public :: vector, operator (+), operator(.cross.)
type, public :: Vector
real :: x = 0.0 ! Cartesian component in x direction.
real :: y = 0.0 ! Cartesian component in y direction.
real :: z = 0.0 ! Cartesian component in z direction.
end type Vector
interface operator (+)
module procedure add
end interface
interface operator (.cross.)
module procedure crossproduct
end interface
contains
elemental function add(left, rigth) result(summ)
type(Vector), intent(in) :: left
type(Vector), intent(in) :: rigth
type(Vector) :: summ
summ%x = left%x + rigth%x
summ%y = left%y + rigth%y
summ%z = left%z + rigth%z
end function add
elemental function crossproduct(vec1, vec2) result(cross)
type(Vector), intent(in) :: vec1
type(Vector), intent(in) :: vec2
type(Vector) :: cross
cross%x = (vec1%y * vec2%z) - (vec1%z * vec2%y)
cross%y = (vec1%z * vec2%x) - (vec1%x * vec2%z)
cross%z = (vec1%x * vec2%y) - (vec1%y * vec2%x)
end function crossproduct
end module vector_t
A KISS pure Fortran OOD class for computing Vectorial (3D) algebra
subroutine print_char(this, header)
character(len=*), intent(in) :: this
logical, optional, intent(in) :: header
if (present(header).and.header) then
print *, 'This is the header'
endif
print *, this
end subroutine print_char
subroutine print_char(this, header)
character(len=*), intent(in) :: this
logical, optional, intent(in) :: header
if (present(header)) then
if (header) print *, 'This is the header'
endif
print *, this
end subroutine print_char
program intent_gotcha
type mytype
integer :: x
real :: y
end type mytype
type (mytype) :: a
a%x = 1 ; a%y = 2.
call assign(a)
! a%y COULD BE UNDEFINED HERE
print *, a
contains
subroutine assign(this)
type (mytype), intent (out) :: this
this%x = 2
end subroutine assign
end program intent_gotcha
program intent_gotcha
type mytype
integer :: x
real :: y
end type mytype
type (mytype) :: a
a%x = 1 ; a%y = 2.
call assign(a)
print *, a
contains
subroutine assign(this)
type (mytype), intent (out) :: this
this%x = 2 ; this%y = 2.
end subroutine assign
end program intent_gotcha
real function kinetic_energy(v)
real, dimension(:), intent(in) :: v
integer :: i
real :: ke = 0.0
do i = 1, size(v)
ke = ke + v(i)**2
enddo
kinetic_energy = .5*ke
end function kinetic_energy
real function kinetic_energy(v)
real, dimension(:), intent(in) :: v
integer :: i
real :: ke
ke = 0.
do i = 1, size(v)
ke = ke + v(i)**2
enddo
kinetic_energy = .5*ke
end function kinetic_energy
As poor Fortran people we are used to simply translate formula, we do not care about paradigms… What Paradigm is the following?
subroutine solr(toll,Qm,gamma1,a1,u1,p1,rho1,gamma4,a4,u4,p4,rho4,u,p,a,fF)
real*8 a2,p2,rho2,a3,p3,rho3,ps,us
real*8 pn,pn_1
real*8 a,u,p
real*8 w1,w4,toll,Qm
real*8 delta,gamma,gm1
...
gm1=gamma-1.d0
delta=gm1/2.d0
alfa=(p1/p4)**(gm1/(2.d0*gamma))
beta=alfa*delta/a1+delta/a4
...
Aside to be very ugly, the above example shows that data is HIGHLY DECOUPLED. This is typical of PROCEDURAL PROGRAMMING paradigm
The real world can be accurately described as a collection of OBJECTS that interact
IMPLEMENTATION is hidden from the rest of the program, aka Application Program Interface.
Do you need to know exactly how every aspect of a car works (engine, carburettor, alternator…)? No - you need to know how to use the steering wheel, brakes, accelerator…
Define What is your Object and How it interacts
module gas_t
use vecfor, only : vector
implicit none
private
public :: gas ! Only the class GAS is public!
type :: gas
private ! hide all members!
! the suffix '_' means internal implementation
real :: Cp_=0.0
real :: Cv_=0.0
real :: density_=0.0
real :: energy_=0.0
type(vector) :: momentum_
contains
private
! the setter (constructor)
procedure, public :: initialize
! return the members value (the getters)
procedure, public :: density
procedure, public :: momentum
procedure, public :: energy
! interact with the rest of the World...
procedure, public :: temperature
procedure, public :: speed_of_sound
...
end type gas
contains
! follow the implementation
...
end module gas_t
Define your World using the Contract
use gas_t, only : gas
implicit none
type(gas) :: air ! a 'gas' object (instance)
real :: density ! an alert dentisyt value
! respect the contract, use only public methods
call air%initialize(Cp=1004.8, &
Cv=716.0, &
density=0.125)
! well designed API is auto-explicative...
if (air%speed()>air%speed_of_sound()) then
! do your supersonic stuffs...
end if
So, where is encapsulation? The user of GAS class
gas_t
module:Extend a Previous Contract adding new behaviors
module gas_multifluid_t
use gas_t, only : gas
implicit none
private
public :: gas_multifluid
! inherit from gas
type, extends(gas) :: gas_multifluid
private
integer :: species_number=0
real, allocatable :: concentrations(:)
contains
private
! overridded methods
procedure, public :: initialize
...
end type gas_multifluid
contains
! follow the implementation
...
end module gas_multifluid_t
Define your new World using the new Contract
use gas_t, only : gas
use gas_multifluid_t, only : gas_multifluid
implicit none
type(gas) :: air
type(gas_multifluid) :: mixture
! respect the contract, use only public methods
call air%initialize(Cp=1004.8, &
Cv=716.0, &
density=0.125)
call mixture%initialize(Cp=1004.8, &
Cv=716.0, &
species_number=3, &
concentrations=[1.0, &
0.0, &
0.0], &
density=0.125)
! well designed API is auto-explicative...
if (air%speed()>air%speed_of_sound()) then
! do your supersonic stuffs...
end if
So, where is inheritance? The developer of MULTIFLUID GAS class
While the user of MULTIFLUID GAS class
If you needed to write a message on a piece of paper, you could use a pen, pencil, marker or even a crayon. You only require that the item you use can fit in your hand and can make a mark when pressed against the paper.
In other words it means, one method with multiple implementation, for a certain class of action. And which implementation to be used is decided at both runtime/compiletime (dynamic/static). Overloading is static polymorphism while, overriding is dynamic polymorphism.
Develop your algorithms for family of classes… as you do for your daily-physics
module gas_fluxes_t
use gas_t, only : gas
implicit none
private
public :: fluxes
type :: fluxes
real :: mass_=0.0
real :: energy_=0.0
type(vector) :: momentum_
contains
...
procedure, public :: compute
end type fluxes
contains
subroutine compute(self, gas_cell)
! fluxes at interfaces to be computed
class(fluxes), intent(inout) :: self
! gas into the cell from which compute fluxes
class(gas), intent(in) :: gas_cell
...
select type(gas_cell)
class is (gas_multifluid)
...
class is (gas)
...
end select
end subroutine compute
end module gas_fluxes_t
use gas_t, only : gas
use gas_multifluid_t, only : gas_multifluid
use gas_library, only : compute_fluxes
implicit none
type(gas) :: air
type(gas_multifluid) :: mixture
type(fluxes) :: air_fluxes
type(fluxes) :: mixture_fluxes
! respect the contract, use only public methods
call air%initialize(Cp=1004.8, &
Cv=716.0, &
density=0.125)
call mixture%initialize(Cp=1004.8, &
Cv=716.0, &
species_number=3, &
concentrations=[1.0, &
0.0, &
0.0], &
density=0.125)
! well designed API is auto-explicative…
if (air%speed()>air%speed_of_sound()) then
call air_fluxes%compute(gas_cell=air)
call mixture_fluxes%compute(gas_cell=mixture)
end if
Such a mathematical formulation is ubiquitous in the mathematical modelling of physical problems (fluid dynamics, chemistry, biology, evolutionary-anthropology, etc…)
U, R and F have different forms for each problem: in a procedural paradigm programming YOU NEED TO DEVELOP A DIFFERENT IMPLEMENTATION FOR EACH PROBLEM!… no, we can be more smart…
The mathematical formulation of ODE is already ABSTRACT:
Develop an ODE solver for each system of equations:
One Ring to rule them all, One Ring to find them, One Ring to bring them all and in the darkness bind them
We want to code a formula like the following (Euler solver)
u = u + Dt * u%t(t=t)
The code is very, very similar to the Math… this is our aim!
FOODIE, Fortran Object-Oriented Differential-equations Integration Environment
https://github.com/Fortran-FOSS-Programmers/FOODIE
Authors:
How to define an Abstract Class in (modern) Fortran?
type, abstract :: integrand
contains
procedure(time_derivative), deferred, public :: t
...
end type integrand
The ABSTRACT keyword implies
integrand
, but only of its (concrete) extensions:integrand
, you can ignore the implementation:The DEFERRED keyword implies
integrand
must provide the implementation of each deferred method.How look our (abstract) Integrand?
type, abstract :: integrand
contains
! deferred methods
procedure(time_derivative), deferred, public :: t
procedure(symmetric_operator), deferred :: integrand_multiply_integrand
procedure(integrand_op_real), deferred :: integrand_multiply_real
procedure(real_op_integrand), deferred :: real_multiply_integrand
procedure(symmetric_operator), deferred :: add
procedure(symmetric_operator), deferred :: sub
procedure(assignment_integrand), deferred :: assign_integrand
! public interface
generic, public :: operator(+) => add
generic, public :: operator(-) => sub
generic, public :: operator(*) => integrand_multiply_integrand, &
real_multiply_integrand, &
integrand_multiply_real
generic, public :: assignment(=) => assign_integrand
end type integrand
The Public Interface is very, very simple, and matches the Math definition of U
+-x
); generic, public :: operator(+)
generic, public :: operator(-)
generic, public :: operator(*)
=
); generic, public :: assignment(=)
%t
). procedure(time_derivative), deferred, public :: t
With only the knowledge of the INTERFACE of these 4 ingredients we can build ALL ODE solvers… just remember the Euler formula…
u = u + Dt * u%t(t=t)
abstract interface
function real_op_integrand(lhs, rhs) &
result(operator_result)
import :: integrand
real, intent(in) :: lhs
class(integrand), intent(in) :: rhs
class(integrand), allocatable :: operator_result
end function real_op_integrand
function time_derivative(self, t) &
result(dState_dt)
import :: integrand
class(integrand), intent(in) :: self
real, optional, intent(in) :: t
class(integrand), allocatable :: dState_dt
end function time_derivative
function integrand_op_real(lhs, rhs) &
result(operator_result)
import :: integrand
class(integrand), intent(in) :: lhs
real, intent(in) :: rhs
class(integrand), allocatable :: operator_result
end function integrand_op_real
end interface
type, abstract :: integrand
contains
procedure(symmetric_operator), deferred :: &
integrand_multiply_integrand
procedure(integrand_op_real), deferred :: &
integrand_multiply_real
procedure(real_op_integrand), deferred :: &
real_multiply_integrand
procedure(symmetric_operator), deferred :: add
…
end type integrand
For all DEFERRED methods you must provide an
EXPLICIT ABSTRACT INTERFACE
abstract interface
function symmetric_operator(lhs, rhs) &
result(operator_result)
import :: integrand
class(integrand), intent(in) :: lhs
class(integrand), intent(in) :: rhs
class(integrand), allocatable :: operator_result
endfunction symmetric_operator
subroutine assignment_integrand(lhs, rhs)
import :: integrand
class(integrand), intent(inout) :: lhs
class(integrand), intent(in) :: rhs
endsubroutine assignment_integrand
end interface
module foodie_adt_integrand
implicit none
private
public :: integrand
type, abstract :: integrand
contains
procedure(time_derivative), deferred :: t
procedure(symmetric_operator), deferred :: &
integrand_multiply_integrand
procedure(integrand_op_real), deferred :: &
integrand_multiply_real
procedure(real_op_integrand), deferred :: &
real_multiply_integrand
procedure(symmetric_operator), deferred :: add
procedure(symmetric_operator), deferred :: sub
procedure(assignment_integrand), deferred :: &
assign_integrand
generic :: operator(+) => add
generic :: operator(-) => sub
generic :: operator(*) => integrand_multiply_integrand,&
real_multiply_integrand, &
integrand_multiply_real
generic :: assignment(=) => assign_integrand
endtype integrand
abstract interface
function time_derivative(self, t) result(dState_dt)
import :: integrand, R_P
class(integrand), intent(IN) :: self
real(R_P), optional, intent(IN) :: t
class(integrand), allocatable :: dState_dt
end function time_derivative
...
…
function integrand_op_real(lhs, rhs) &
result(operator_result)
import :: integrand, R_P
class(integrand), intent(IN) :: lhs
real(R_P), intent(IN) :: rhs
class(integrand), allocatable :: operator_result
end function integrand_op_real
!
function real_op_integrand(lhs, rhs) &
result(operator_result)
import :: integrand, R_P
real(R_P), intent(IN) :: lhs
class(integrand), intent(IN) :: rhs
class(integrand), allocatable :: operator_result
end function real_op_integrand
!
function symmetric_operator(lhs, rhs) &
result(operator_result)
import :: integrand
class(integrand), intent(IN) :: lhs
class(integrand), intent(IN) :: rhs
class(integrand), allocatable :: operator_result
end function symmetric_operator
!
subroutine assignment_integrand(lhs, rhs)
import :: integrand
class(integrand), intent(INOUT) :: lhs
class(integrand), intent(IN) :: rhs
end subroutine assignment_integrand
end interface
end module foodie_adt_integrand
module foodie_integrator_euler_explicit
use foodie_adt_integrand, only : integrand
implicit none
private
public :: euler_explicit_integrator
!
type :: euler_explicit_integrator
private
contains
private
procedure, nopass, public :: integrate
end type euler_explicit_integrator
contains
subroutine integrate(U, Dt, t)
class(integrand), intent(inout) :: U
real(R_P), intent(in) :: Dt
real(R_P), optional, intent(in) :: t
U = U + U%t(t=t) * Dt
return
end subroutine integrate
end module foodie_integrator_euler_explicit
Compare with the math
U%t, U+U, U*Dt...
;euler
integrator accepts ALL integrand extensions!use foodie_integrator_euler_explicit
use my_integrand
type(euler_explicit_integrator) :: solver
type(my_integrand) :: U
...
call solver%integrate(U=U, Dt=0.1)
module foodie_integrator_low_storage_runge_kutta
use foodie_adt_integrand, only : integrand
implicit none
private
public :: ls_runge_kutta_integrator
!
type :: ls_runge_kutta_integrator
private
integer :: stages=0 !< Number of stages.
real, allocatable :: A(:) !< Low storage *A* coefficients.
real, allocatable :: B(:) !< Low storage *B* coefficients.
real, allocatable :: C(:) !< Low storage *C* coefficients.
contains
private
...
procedure, pass(self), public :: integrate
endtype ls_runge_kutta_integrator
contains
...
subroutine integrate(self, U, stage, Dt, t)
class(ls_runge_kutta_integrator), intent(in) :: self
class(integrand), intent(inout) :: U
class(integrand), intent(inout) :: stage(1:)
real(R_P), intent(in) :: Dt
real(R_P), intent(in) :: t
integer(I_P) :: s
select type(stage)
class is(integrand)
stage(1) = U
stage(2) = U*0._R_P
do s=1, self%stages
stage(2) = stage(2) * self%A(s) + &
stage(1)%t(t=t + self%C(s) * Dt) * Dt
stage(1) = stage(1) + stage(2) * self%B(s)
enddo
U = stage(1)
endselect
end subroutine integrate
end module foodie_integrator_low_storage_runge_kutta
Compare with the math
U%t, U+U, U*Dt...
;euler
integrator accepts ALL integrand extensions!use foodie_integrator_low_storage_runge_kutta
use my_integrand
type(ls_runge_kutta_integrator) :: solver
type(my_integrand) :: U
...
call solver%integrate(U=U, Dt=0.1)
module oscillation_ode
use foodie, only : integrand
implicit none
private
public :: oscillation
type, extends(integrand) :: oscillation
private
real :: f=0.0
real :: U(1:2)=[0.0, 0.0]
contains
procedure, pass(self), public :: t => &
dOscillation_dt
procedure, pass(lhs), public :: &
integrand_multiply_integrand => &
oscillation_multiply_oscillation
procedure, pass(lhs), public :: &
integrand_multiply_real => &
oscillation_multiply_real
procedure, pass(rhs), public :: &
real_multiply_integrand => &
real_multiply_oscillation
procedure, pass(lhs), public :: add => &
add_oscillation
procedure, pass(lhs), public :: sub => &
sub_oscillation
procedure, pass(lhs), public :: &
assign_integrand => &
oscillation_assign_oscillation
end type oscillation
contains
...
end module oscillation_ode
...
function dOscillation_dt(self, t) result(dState_dt)
class(oscillation), intent(in) :: self
real(R_P), optional, intent(in) :: t
class(integrand), allocatable :: dState_dt
allocate(oscillation :: dState_dt)
select type(dState_dt)
class is(oscillation)
dState_dt = self
dState_dt%U(1) = -self%f * self%U(2)
dState_dt%U(2) = self%f * self%U(1)
endselect
end function dOscillation_dt
function oscillation_multiply_oscillation(lhs, rhs) &
result(opr)
class(oscillation), intent(in) :: lhs
class(integrand), intent(in) :: rhs
class(integrand), allocatable :: opr
allocate(oscillation :: opr)
select type(opr)
class is(oscillation)
opr = lhs
select type(rhs)
class is (oscillation)
opr%U = lhs%U * rhs%U
endselect
endselect
end function oscillation_multiply_oscillation
So, to integrate the Oscillation ODE we have passed through the EXTENSION of ABSTRACT INTEGRAND, this is an overhead. why bothering with it? you are thinking. Because:
OpenMP benchmarks
CAF benchmarks
Sometimes called CC-UMA - Cache Coherent UMA. Cache coherent means if one processor updates a location in shared memory, all the other processors know about the update. Cache coherency is accomplished at the hardware level.
The number of parallel regions and the threads that comprise them are arbitrary.
It is not Fortran Standard!
!$OMP some_directives
! fortran statements to be done in parallel
do i=1...
directive [clause[ clause]...]
ATTENTION: they may alter code semantic: the code can be corrected in serial but not in parallel or viceversa.
OMP_NUM_THREADS
: sets number of threadsOMP_STACKSIZE size [B|K|M|G]
: size of the stack for threadsOMP_DYNAMIC {TRUE|FALSE}
: dynamic thread adjustmentOMP_SCHEDULE "schedule[,chunk]
: iteration scheduling schemeOMP_PROC_BIND {TRUE|FALSE}
: bound threads to processorsOMP_NESTED {TRUE|FALSE}
: nested parallelismTo set them
setenv OMP_NUM_THREADS 4
export OMP_NUM_THREADS=4
omp_get_thread_num()
: get thread ID (0 for master thread)omp_get_num_threads()
: get number of threads in the teamomp_set_num_threads(int n)
: set number of threadsomp_init_lock(lock_var)
: initializes the OpenMP lock variable lock_var
of type omp_lock_t
omp_get_wtime()
: returns elapsed wallclock timeomp_get_wtick()
: returns timer precisionFunctions interface
use omp_lib
! or old-style include
include 'omp_lib.h'
!$omp parallel
! some code to execute in parallel
print "(A)", 'Hello (parallel) World!'
!$omp end parallel
Variables are shared by default but with the clause default(none)
no implicit default, you have to scope all variables explicitly, sane way!
sum = 0
!$omp parallel private(i, MyThreadID, psum)
MyThreadID = omp_get_thread_num()
NumThreads = omp_get_num_threads()
psum =0
do i=MyThreadID*N/NumThreads+1, min((MyThreadID+1)*N/NumThreads, N)
psum = psum + x(i)
end do
!$omp critical
sum = sum + psum;
!$omp end critical
!$omp end parallel
parallel
constructcritical
constructomp_get_thread_num()
functionomp_get_num_threads()
functiondo
loop constructsections
constructsingle
constructworkshare
construct program doexample
integer, parameter:: n=10
integer:: i, a(n), b(n), c(n)
!$omp parallel
!$omp do
do i=1, n
a(i) = i
b(i) = i
c(i) = 0
end do
!$omp end do
!$omp do
do i=1, n
c(i) = a(i) + b(i);
end do
!$omp end do
!$omp end parallel
end program doexample
program vec_add_sections
integer, parameter :: n=100
integer :: i
real :: a(n), b(n), c(n), d(n)
do i = 1, n
a(i) = i * 1.5
b(i) = i + 22.35
enddo
!$omp parallel shared(a,b,c,d), private(i)
!$omp sections
!$omp section
do i = 1, n
c(i) = a(i) + b(i)
enddo
!$omp section
do i = 1, n
d(i) = a(i) * b(i)
enddo
!$omp end sections nowait
!$omp end parallel
end
integer, parameter :: n=100
integer :: i
real :: a(n), b(n), c(n), d(n)
...
!$omp parallel
!$omp single
a = b / c
print*, a
!$omp end single
!$omp do
do i = 1, n
d(i) = a(i) * b(i)
enddo
!$omp end parallel
program workshare
integer, parameter :: n=100
integer :: i, j
real :: aa(n,n), bb(n,n)
real :: cc(n,n), dd(n,n)
real :: first, last
do i = 1, n
do j = 1, n
aa(j,i) = i * 1.0
bb(j,i) = j + 1.0
enddo
enddo
!$omp parallel shared(aa,bb,cc,dd,first,last)
!$omp workshare
cc = aa * bb
dd = aa + bb
first = cc(1,1) + dd(1,1)
last = cc(n,n) + dd(n,n)
!$omp end workshare nowait
!$omp end parallel
end
shared(list)
: there is only one instance of the objects in the list accessible by all threads in the teamprivate(list)
: each thread has a copy of the variables in the listfirstprivate(list)
: same as private but all variables in the list are initialized with the value that the original object had before entering the parallel constructlastprivate(list)
: same as private but the thread that executes the sequentially last iteration or section updates the value of the objects in the listdefault(none|shared|private|firstprivate)
!$omp parallel do reduction(+:sum)
do i=1, n
sum = sum + x(i)
end do
!$omp end parallel do
reduction(op:list)
+
they are initialized to 0
, if op is *
they are initialized to 1
)+, *, -, .and., .or., .eqv., .neqv., max, min, iand, ior, ieor
parallel
, do
loop and sections
constructsbarrier
construct;atomic
construct;Not covered today!
subroutine prime_number(n, total)
! Return the number of primes between 1 and N.
implicit none
integer, intent(in) :: n
integer, intent(out) :: total
integer :: i
integer :: j
integer :: prime
total = 0
!$omp parallel &
!$omp shared(n) &
!$omp private(i, j, prime)
!$omp do reduction (+: total)
do i=2, n
prime = 1
do j=2, i - 1
if (mod(i, j) == 0) then
prime = 0
exit
end if
end do
total = total + prime
end do
!$omp end do
!$omp end parallel
return
end
use omp_lib
real, parameter :: pi = 3.141592653589793D+00
integer, parameter :: n=2000
real :: a(n,n), b(n,n), c(n,n)
integer :: i, j, k
real :: s, angle
integer ::
s = 1.0D+00 / sqrt(real(n, kind = 8))
!$omp parallel shared (a, b, c, n, s) &
!$omp private (angle, i, j, k)
!$omp do
do i = 1, n
do j = 1, n
angle = 2.0D+00 * pi * (i - 1) * (j - 1) / real(n)
a(i,j) = s * (sin(angle) + cos(angle))
end do
end do
!$omp end do
!$omp do
do i = 1, n
do j = 1, n
b(i,j) = a(i,j)
end do
end do
!$omp end do
!$omp do
do i = 1, n
do j = 1, n
c(i,j) = 0.0D+00
do k = 1, n
c(i,j) = c(i,j) + a(i,k) * b(k,j)
end do
end do
end do
!$omp end do
!$omp end parallel