Fortran basics
Program template
program program_name
implicit none
! program body
end program program_name
Comment: !
Case insensitive
Line continuation
Use &
at the end and at the beginning to break a long line into 2.
E.g.
real, dimension (1:16) :: pad1 = (/ 1, 2, 3, 4, &
& 5, 6, 7, 8 /)
Declare variable
real :: var1,x,y
integer type
huge()
gives the largest number that can be held by the specified integer type.
E.g.
integer :: largeval
print *, huge(largeval)
To specify the number of bytes, use kind
specifier:
integer(kind=2) :: shortvar
integer(kind=8) :: longvar
logical type
.TURE.
, .FALSE.
Complex type
cx = 1.0 + 2.0i
cx = cmplx(1.0, 2.0)
Custom type
To define a custom type:
type :: type_name
integer :: id
integer :: age
character(len=50) :: name
character(len=2) :: sex
end type
The double ::
in type :: type_name
is optional.
To declare a custom type variable:
type(type_name) :: var1
To assign values:
var1%id=1
var1%age=20
var1="mark"
var1%sex="M"
or
var1=type_name(1,20,"mark","M")
To access a field:
var1%field_name
To create an array of custom types:
type(type_name), dimension(10) :: type_array
To access an element:
type_array(2)%fieldname
Print string
print *, 'This is a string.'
print *, "The answer is ", answer
or
write(*,*) "new string."
where *
is for the default output option (screen).
Format output
Integer and real number specification
syntax:
write(output_device, label) vars label format_spec
E.g.
integer :: num1, num2
write(*,1) num1 1 format(2i10)
write(*,2) num2 2 format(f6.2)
Another general form:
read (unit=u, fmt=fmt) variable_list
write (unit=u, fmt=fmt) variable_list
where fmt
is formatting string. E.g. (i3,f6.4)
.
u
is a unit number, a number associated with a file.
E.g.
write(*,"(i10.2)") number
*
: output device is screen.2i10
: right justified integer, <= 10 digits, 2 digits per line.f6.2
: no more than 6 digits in total (including sign if there is one and decimal point), 2 decimal places.i10.2
: <= 10 digits, >= 2 digits.
Exponential specification
nEm.d
:
d
: number of decimal places.m
: total width (including sign if any, E and its sign, decimal point and all numbers).n
: number of exponential numbers per line.
E.g.
format(2E14.5)
Character specification
nAm
m
: max number of characters.n
: number of strings per line.
E.g. format(2a10)
Print new line in formatted output
Use //
in the format string:
write(*,"(i3//i3)") 10, 10
Print spaces in formatted output:
Use nx
in format string, n
is the number of spaces. E.g.
write(*,"(A,2x,A,2x,A,i3)") 'name', trim(tt), 'length', tt%length
The output looks like:
name mark length 5
More on formatted I/O
Input
read *, x
read *, y
Character variables
To declare a character variable with a max length:
character :: var*10
or
character (len=10) :: var
character(10) :: var
character(len=*) :: var ! unknown length
To get substring, use indexing:
var(1:4)
Both ends included.
To concatenate characters, use //
:
name=title//firstname//surname
Some string functions
len(str)
: length.len_trim(str)
: length after removing trailing blanks.index(str,substr)
: find index location of a substring. Return 0 if not found.achar(int)
: convert an integer to a character. E.g. 65->A, 90->Z.ichar(c)
: convert a character to an integer.trim(str)
: remove trailing blanks (similar asrstrip()
in python).scan(str,chars)
: find the index location of the 1st occurrence ofchars
. Return 0 if not found.verify(str,chars)
: find the 1st occurrence of any character not contained inchars
. Return 0 if only the characters inchars
have been found.repeat(str,n)
: repeatstr
byn
times.adjusl(str)
: remove leading blanks, and append them to the end.adjustr(str)
: remove trailing blanks, and insert at the front.
Named constants
Define constants use the parameter
keyword:
real, parameter :: pi = 3.1415
Power (same as in python)
=10**2=
Common math functions
NOTE that input args should be real or complex, e.g. sin(10.0) rather than sin(10).
sin(x), cos(x), tan(x), atan(x), abs(x), sqrt(x), exp(x), log10(x) floor(x), int(x).
nint(x)
: round to nearest integer.
mod(a,p)
: remainder of a/p.
if block
if (expression) then
...
else if (expression) then
...
else
...
end if
Simpler if statement
if (num <0) stop
if (num <10) print *, "less than 10"
if (num >10) print *, "> 10"
Logical operators
-
== or .eq.
-
/= or .ne.
-
or .gt.
-
= or .ge.
-
< or .lt.
-
<= .le.
-
- or .not.
-
.and.
-
.or.
if (x>y .and. y>z) then … end if
stop
stop
: stops the program.
pause
Pause the program, until user hit Enter.
Integer division
Same as Python2, need to involve one real to get real division result. E.g.
integer :: x
real :: y
x=1
y=x/3 ! gives 0.00000
y=x/3.0 ! gives 0.33333
do loop
do i=start,stop[,step]
...
end do
E.g.
do i=0,20
print *, i
end do
do i=5,-5,-2
print *, i
end do
exit: same as break
cycle: same as continue
Named do loop
When need to nest do-loops, it is usually helpful to make the code clearer by naming the outer and inner loops, E.g.
outer: do i=1,10
inner: do j=1,20
if (i==3) exit outer
if (j==2) cycle inner
array(i,j)=2*i+j
end do inner
end do outer
do while loop
do while (expression)
...
end do
forall loop
forall(i=1:5, j=1:10, a(i,j)<10)
a(i,j)=1
end forall
Read from file
Use open()
for file handling. The general form:
open(unit=u, file=filename, status=st, action=act, iostat=ios)
where
u
: a positive integer, unit number. Often the range 1-99 is allowed (processor-specific). 0 is often special. 5 is for console input. 6 for output to screen. Therefore, don’t use 0, 5, 6.st
: str,- "new": create new, the file should not yet exists.
- "old": open existing file.
- "replace": create new if not exits, otherwise overwrite.
- "scratch": a file is used only during the execution of the program.
filename
: must not use whenst
is "scratch".act
: "read", "write", "readwrite".ios
: an integer variable. It’s set to 0 if open is successful, a positive integer otherwise.
E.g.
real :: x,y,z
open(10, file='mydata.txt')
read(10,*) x,y,z
where:
mydata.txt
: has 3 values which are read and saved to x,y,z.10
: named input device.
Write to file
E.g.
open(12, file="myoutput")
write(12,*) var
E.g.
open(8, file="myoutput")
write(8, "(a2, 3(f10.6))") dtype,x,y,z
E.g.
program file_example
implicit none
integer :: ierror
open(unit=13, file="test.dat", status="new", action="write", iostat=ierror)
if (ierror /= 0) then
print *, "Failed to open test.data"
stop
end if
write (unit=13, fmt=*) "Hello world"
close (unit=13)
end program file_example
It is good programming practice to test if the open statement was successful. If failed, the program would have crashed when trying to write into it.
Close opened file
When a file is no longer needed, use close()
to close it:
close (unit=u, iostat=ios, status=st)
where:
- st:
- "keep": keep the file. Default, except for "scratch" files.
- "delete": delete the file. Default for "scratch" files.
Internal files
read()
and write()
can also read or write into internal files. An internal file is a character variable.
E.g.
character (len=8) :: string
write (string, "a", iostat=ios) "abcdfefgh"
Extend single precision
In Fortran 77, the prevision of real
type can be increased by using the double precision
attribute:
double precision :: x
In Fortran 90, the types of integer, real, complex and logical have a "default kind" and a number of other kinds. How many other kinds there are for a certain type depends on the particular processor. Each kind has its own kind type parameter, which is an integer of positive value. For example, if, for a certain processor, a kind value of 8 yields a precision equivalent to the old double precision type of Fortran 77, then the following statement
real (kind=8) :: x1
is equivalent to the Fortran 77 statement
double precision x1
However, this is not very portable, because the required kind value may be different on another computer. Although many computers use kind values that indicate the number of bytes used for storage of the variable, you cannot rely on this.
In Fortran 90, there are 2 intrinsic functions to obtain the kind value for the required precision:
selected_real_kind()
: for precision of real.selected_int_kind()
: for precision of int.
The selected_real_kind()
function returns an integer that is the kind type parameter value necessary for a given decimal precision p
and decimal exponent range r
.
E.g.
integer, parameter :: ikind=selected_real_kind(p=7, r=10)
The integer ikind
now contains the kind value needed for a precision
of 7 decimal places, and a range of at least (10^{-10} to 10^{10}).
If the kind value for the required precision or range is not available, a negative integer is returned.
To use the ikind value in a type declaration statement, it has to be a constant. E.g.
integer, parameter :: ikind=selected_real_kind(p=15)
real (kind=ikind) :: sum,x
Valid values for p
are 6, 15 and 18. Unlike variables, parameters can’t change once declared.
To use constants in a program that uses extended precision, append _ikind
after the constant, where ikind
is the precision parameter.
E.g.
integer, parameter :: ikind=selected_real_kind(p=15)
real (kind=ikind) :: x
x=10.0_ikind/3
The selected_int_kind
function returns the lowest kind value needed
for the integers with the specified range:
E.g.
integer, parameter :: ikind=selected_int_kind(r=10)
integer (kind=ikind) :: big_number
Now the integer big_number
can represent integers in the range (-10^{10} – 10^{10}).
Random number
random_seed()
: start random number generator, must call this before calling other random number functions.random_number(t)
: [0,1) uniform random.
subroutines and functions
Define subroutine
subroutine NAME ()
implicit none
...
end subroutine
Can use return
in subroutine to return at given point. If not give, return at the last line.
Subroutine with arguments
Main program and subroutines don’t share namespaces.
To make subroutines take arguments:
subroutine NAME(arg1,arg2...)
implicit none
real :: arg1, arg2,...
end subroutine
Argument passing are achieved by passing their ram addresses. NOTE that still need to declare args inside subroutine. NOTE that arguments changed inside subroutine will be changed in the main program.
Call subroutine
call NAME()
call NAME(arg1,arg2,...)
When it reaches end subroutine
, the values of the input args are
returned to the main program.
Custom functions
Syntax
function FUNC(arg1, arg2,...)
implicit none
real :: FUNC
...
end function
Call a function
program NAME
implicit none
real, external :: FUNC
...
FUN(arg1,arg2,...)
end program
NOTE:
- must declare the data type of the function both in the main program and in the function itself as if it were a variable.
external
is optional, but helpful in distinguishing variables and functions.- functions return 1 value, which is assigned to the function name.
- variables inside functions, except input arguments, have local scope (not accessible from outside).
If you want the returned value to be stored in some other name than
(can’t be the same as function name) the function name, use the result
option:
function funcname(arg1,arg2,...) result (return_name)
...
end function
Pass functions as function arguments
Can pass a function name as argument to another function. E.g.
program main
implicit none
real, external :: func ! declare a func
real, intrinsic :: sin ! build-in sin func
call Exefunc(func)
call Exefunc(sin)
end program
subroutine Exefunc(f)
implicit none
real, external :: f
write (*,*) f(1.0)
end subroutine
real function func(num)
implicit none
real :: num
func=num*2
end function
NOTE that external
cannot be omitted here. Subroutines names can also
be passed to other subroutines.
One-liner function
real, external :: func
func(a,b) = a+b
NOTE that such functions cannot be called from outside the main program.
Intent attribute
Use the intent
attribute to specify the intention of the arguments in
subroutine or function.
The possible attribute values:
- in: used as input values, not changed in function (read only)
- out: use as output, they are overwritten
- inout: arguments are used and overwritten, same as not specifying intent.
E.g.
subroutine intent_example (a,b,c,d)
implicit none
real, intent (in) :: a
real, intent (in) :: b
real, intent (in) :: c
real, intent (out) :: d
d= b * b - 4.0 * a * c
end subroutine intent_example
It’s good practice to use intent
attributes, which make the code
clearer and more transparent. It also makes the compiler catch
mistakes, for instance, when modifying an argument which is denoted as
intent(in)
.
Optional arguments
Use optional
attribute to declare an argument as optional, and use
present()
function to check whether the optional argument is passed
in.
E.g.
program main
implicit none
interface
subroutine sub(a,b)
implicit none
interger :: a
interger, optional :: b
end subroutine
end interface
call sub(1,2)
call sub(2)
subroutine sub(a,b)
implicit none
integer :: a
integer, optional :: b
if (present(b)) then
...
else
...
end if
end subroutine
NOTE that variable-argument type function/subroutine has to be inside
interface
block, unless inside a module. This can be used to create
functions with default arguments, and function call can use keyword
argument syntax, where the order is insignificant:
call func(a=1,b=2,c=3)
call func(b=2,a=1,c=3)
call func(b=2,a=1) ! make c an optional arg, and give it a default value if present(c) is .false.
NOTE that you can’t use a var that is not present, have to define one for receiving the optional argument, and another for actual computation:
E.g.
subroutine sub1(arr,opt_x)
implicit none
real, dimension(:,:), intent(in) :: arr
integer, optional :: opt_x
integer :: x
if (present(opt_x)) then
x=opt_x
else
x=10
end if
...
! use x so on
...
end subroutine
Recursive procedures
Declare a recursive function/subroutine using the recursive
keyword. When a function is defined as recursive, the result
option
must be used.
E.g.
recursive function factorial(n) result(fac)
implicit none
integer :: fac
integer, intent(in) :: n
select case(n)
case (0:1)
fac=1
case default
fac=n*factorial(n-1)
end select
end function factorial
Local functions/subroutines
Put function/subroutine definitions inside a contains
– end contains
block to make it private (cannot access from other
program/function/subroutine).
Elemental functions/subroutines
Put elemental
keyword before function/subroutine definition to make it elmental E.g.
interface
elemental function func(num)
implicit none
integer, intent(in) :: num
end function
end interface
elemental function func(num)
implicit none
integer, intent(in):: num
func=...
end function
NOTE that elemental func/sub cannot take arrays as argument, and has
to be inside interface
block.
To call an elemental function:
a = func(a)
where a
is an array, and func()
is defined based on scalar
computations, but its elemental
feature broadcasts the computation to
an array, equivalent to do
a do loop over all elements.
Interface
Interface blocks are used to specify the types of functionn arguments and return values. When a procedure is external (defined outside of the main programe), the compiler will (in most cases) know about the procedure’s interface, and can’t check if the precedure call is consistent with the procedure declaration. It’s such good practice to provide interfaces for external procedures.
These following scenarios require an interface:
- function returns an array.
- function has variable number of arguments
- function takes pointer as arguments or returns pointers.
Syntax:
interface
function func_name(arg1,arg2,...)
implicit none
!------ declare arg types --------
real :: arg1
integer :: arg2
...
!------- declare return value type ------
real :: func_name(10) ! an array
end function
subroutine sub_name
implicit none
interger :: ....
...
end subroutine
end interface
The interface block is placed at the beginning of the program.
If a procedure proc1
calls another procedure proc2
, then the interface
block of proc2
should be placed at the beginning of the procedure
proc1
.
Date and time
The date_and_time()
function returns the date and time.
E.g.
charater(len=8) :: dateinfo ! ccyymmdd
charater(len=4) :: year, month*2, day*2
charater(len=10) :: timeinfo ! hhmmss.sss
charater(len=2) :: hour, minute, second*6
call date_and_time(dateinfo, timeinfo)
year=dateinfo(1:4)
month=dateinfo(5:6)
day=dateinfo(7:8)
Array
Declare array
type, dimension(n) :: x
E.g.
real, dimension(10) :: x
Use parameter in defining array size, e.g.
integer, parameter :: imax = 10
real, dimension(imax) :: x
real :: sum, average
integer :: i
sum=0.0
do i=1,imax
sum=sum+x(i)
end do
average=sum/imax
Declare dynamic array when size is unknown
The size of a dynamic array is unknown until execution time.
To declare:
integer, dimension(:), allocatable :: array1
real, dimension(:,:), allocatable :: array2
NOTE that the number of dimensions has to be specified.
To allocate memory, use allocate()
.
E.g.
integer, allocatable, dimension(:) :: vector
integer :: size
read *, size
allocate(vector(size))
To prevent memory shortage error, get the return value of allocate()
, if it is 0, then memory allocation is successful.
E.g.
integer :: size, return=0
integer, allocatable, dimension(:) :: array
allocate(array(size),stat=return)
if (error/=0) exit
To clear up the memory of an array, use deallocate()
deallocate(vecotr)
Declare array with explicit bounds:
real, dimension(2:6) :: numbers
integer, dimension(-3:2, 0:4) :: matrix
Multi-dimensional arrays
integer, dimension(m,n) :: a
Fortran allows the dimension size up to 7. To declare an extended precision 3-d array:
real (kind=ikind), dimension(3,4,5) :: a
Assign values
do i=1,5
number(i) = i*20
end do
1-D array can be assigned using a shorthand:
numbers = (/1,2,3,4,5/)
numbers2 = [ 1,2,3,4,5 ]
numbers3 = (/ (i, i=1,10,2) /) ! (1,3,5,7,9)
numbers4 = (/ (i*2. i=1,5) /) ! (2.0,4.0,6.0,8.0,10.0)
Use array constructors to define array constants:
integer, dimension(8), parameter :: p = (/1,2,3,5,7,11,13,17,19/)
Indexing
Array indexing starts from 1, using parentheses: e.g. x(2)
.
Both ends are included.
Syntax: array(start:stop:step)
E.g.
a(1:7)
b(1:10:2)
c(5:)
d(1, 1:4:2)
Array functions
size(a[,dim])
: total size.dot_product(vector_a, vector_b)
: dot product.matmul(matrix_a, matrix_b)
: matrix multiplication.transpose()
: transpose.all(condition[,dim])
: similar tonumpy.all()
. Ifdim
is given, only
along given dimension.any(condition[,dim])
: similar tonumpy.any()
.count(condition[,dim])
: count the number of.true.
values returned by condition.maxval(array[,dim,condition])
: return max value of array, optionally along dimensiondim
, satisfying requirementcondition
.minval(array[,dim,condition])
.product(array[,dim,condition])
: return product of all elements, optionally along given dimension and filtered by mask.sum(array[,dim,condition])
.allocated(array)
: check if the array if allocated.lbound(array[,dim])
: get the lower dimension limit for the array. Ifdim
not given, get an integer vector. Ifdim
given, return lower bound of given dimension.ubound(array[,dim])
: upper bound.shape(array)
: same as python.merge(x,y,condition)
: same asnp.where(condition,x,y)
.
reshape(array,shape[,pad,order])
E.g.
real, dimension(3,2) :: a
a=reshape((/1,2,3,4,5,6/), (/3, 2/))
- if
pad
is given, will padarray
with numbers frompad
to fill up the new shape. order
is like the order innumpy.ndarray
, and specifies the order in which numbers are filling. Default is column-wise.
E.g.
integer, dimension(9) :: a = (/1,2,3,4,5,6,7,8,9/)
integer, dimension(10) :: pad = (/-1,-2,-3,-4,-5,-6,-7,-8,-9,-10/)
integer, dimension(4,4) :: c
c=reshape(a, (/4,4/), pad=pad, order=(/2,1/))
c will be:
1 2 3 4 5 6 7 8 9 -1 -2 -3 -4 -5 -6 -7
spread(array,dim,n):
array/slab grow. Grow array
along the dim
dimension by n
copies.
E.g.
integer, dimension(10) :: a = (/1,2,3,4,5,6,7,8,9,10/)
integer, dimension(2,5) :: b
integer, dimension(2,2,5) :: c
b=reshape(a,(/2,5/))
c=spread(b,1,2)
b will be:
1 3 5 7 9 2 4 6 8 10
c(1) and c(2) will both be:
1 2 3 4 5 6 7 8 9 10
cshift(array,shift,dim)
Similar to np.roll(array,shift,axis)
. If shift
is positive, shift to the left, otherwise to the right.
E.g.
integer, dimension(5) :: a = (/1,2,3,4,5/)
integer, dimension(5) :: b
b=cshift(a,1,1)
b will be 2 3 4 5 1
.
eoshift(array,shift,boundary,dim)
Similar to cshift
except that new elements are taken from boundary
which is a scalar. If boundary is not given, pad with 0s.
E.g.
integer, dimension(5) :: a = (/1,2,3,4,5/)
integer, dimension(5) :: b
b=eoshift(a,boundary=-10,shift=2)
b will be 3 4 5 -10 -10
.
maxloc(array[,condition])
Find the max element indices in array
, optionally satisfying the condition
. Equivalent to np.where(a==np.max(a))
.
Return a vector containing the indices.
E.g.
integer, dimension(6) :: a=(/1,2,3,4,5,6/)
integer, dimension(2,3) :: b
integer, dimension(2) :: c
b=reshape(a,shape(b))
c=maxloc(b,b<=5)
b is:
1 3 5 2 4 6
c will be 1 3
.
minloc(array[,condition])
Similar to maxloc()
.
pack(array,mask)
Fancy indexing an array using a logical mask. E.g.
integer, dimension(5) :: data = (/ 1,2,3,4,5 /)
logical, dimension(5) :: mask
mask=.TRUE.
mask(1)=.FALSE.
mask(3)=.FALSE.
write(*,*) pack(data,mask)
will print 2,4,5
.
NOTE that if both array
and mask
are n-dimensional (e.g. both being (3×3)), pack(array,mask)
always returns an 1D array.
where statement
Syntax:
where (condition)
...
[elsewhere (condition)
...]
[elsewhere
...]
end where
E.g.
integer, dimension(6) :: a = (/1,2,3,4,5,6/)
where (a<4)
a=1
elsewhere
a=10
end where
a will be 1 1 1 10 10 10
.
Compiling
File name convention
File names ending in .f90
and .f95
are assumed to be free form, suitable for Fortran 90/95 compilation.
File names ending in .f
and .for
are assumed to be fixed form, compatible with Fortran 77.
Compile a program
Basic compilation
The simplest command:
gfortran myprog.f90
This creates an executable a.out
in the current directory.
To specify the name of the executable:
gfortran myprog.f90 -o myprog
Compile with include files
The path of include files can be given with the -I
option:
gfortran myprog.f90 -o myprog -I/home/user/fortran/inc
or
gfortran myprog.f90 -o myprog -I$MYINC
where the environment variable MYINC
is set with
MYINC=/home/user/fortran/inc/
Speed optimization
The -Olevel
option performs some optimization of the executable and can lead to significant increases in execution speed.
E.g.
gfortran myprog.f90 -o myprog -O2
Warning options
The -Wlevel
option enables most warning messages. Such messages are generated at compile-time warning of, for example, unused or unset variables.
E.g.
gfortran myprog.f90 -o myprog -O2 -Wall
Runtime options
Various run-time options can be selected, these options cause extra code to be added to the executable and so can cause significant decreases in execution speed. However these options can be very useful during program development and debugging.
E.g.
gfortran myprog.f90 -o myprog -O2 -fbounds-check
This causes the executable to check for "array index out of bounds conditions".
Recommended options
For better safety:
gfortran myprog.f90 -o myprog -Wall -fbounds-check -ftrace=full -O2
If speed of execution is important then use the following options:
gfortran myprog.f90 -o myprog -Wuninitialized -Wimplicit-none -Wunused-vars -Wunset-vars -O2
Compile subprograms
If a Fortran 90/95 project contains more than one program source files, then to compile all source files to an executable, use when:
subprograms are external
gfortran main.f90 sub1.f90 sub2.f90 sub3.f90 -o myprog
NOTE main file in front of subprogram files.
subprograms are defined in module
gfortran mod.f90 main.f90 -o myprog
NOTE module file in front of main file.
Create and link object files
When working on large projects, especially projects that are shared between a number of programmers/users, it is useful to create object files or libraries of object files that can be shared between programmers and users.
If a subprogram has been compiled, it’s not necessary to recompile it again during the development of other parts of the project. Completed subprograms can be compiled into object files that can be linked during the compilation of the whole program.
E.g.
gfortran main.f90 sub1.o sub2.o sub3.o -o myprog
where the .o
object files are created by
gfortran -c sub1.f90
gfortran -c sub2.f90
gfortran -c sub3.f90
Create and link library files
If many object files are to be linked then it is convenient to place them in a library.
Libraries can be created using the ar
command. To place object files in a library:
ar rcvf libsubs.a sub1.o sub2.o sub3.o
The library created is named libsubs.a
, containing 3 objects: sub1.o
, sub2.o
, sub3.o
.
NOTE that the name of the library must begin with lib
and end with .a
.
To add more objects into the library file:
ar rcvf libsubs.a sub4.o
To list the contents of a library:
ar tv libsubs.a
To list subroutines and functions in the library:
nm libsubs.a
To delete object from the library:
ar dv libsubs.a sub.o
To link a library while compiling a program:
gfortran main.f90 -o myprog -L. -lsubs
(The -L.
option tells the linker that the library is in the current directory.)
If the library is not in the current directory, need to give the path:
gfortran main.f90 -o myprog -L/home/user/fortran/lib/ -lsubs
Linux library types
There are 2 types of library files:
1. static library .a
Library of object code which is linked at compilation and becomes part of the application. If there is any change in the lib, you need to recompile the application.
Two possible ways of usage:
-
Dynamically linked at run time but statically aware.
The libraries must be available during compile/link phase. The shared objects are not included into the executable component but are tied to the execution.
-
Dynamically loaded/unloaded and linked during execution.
Use the dynamic linking loader system function. If there is any change in the
.so
file, no need to recompile the main program, but do need to make sure the main program is linked to the updated lib file.
At compile time, the linker insures that all the necessary symbols are either linked into the executable, or can be linked at runtime from the shared library. Executables compiled from shared libraries are smaller, but the shared libraries must be included with the executable to function correctly. When multiple programs use the same shared library, only one copy of the library is required in memory.
Library naming conventions
Libraries are typically named with a lib
prefix. When linking, the command line reference to the library will NOT contain the lib
prefix. E.g.
gcc src-file.c -lm -lpthread
The 2 libraries linked in the above example are:
m
->libm.a
pthread
->libpthread.a
.
Both files can be found in /usr/lib/
by default.
- Use the
-shared
option. - On Linux systems using IA-32 architecture and Intel@64 architecture, you must also specify
-fpic
for the compilation of each object file you want to include in the shared library. - Use the
-o
option to specify the output file name.
E.g.
gfortran -shared -fpic -o octagon.so octagon.f90
This will create octagon.so
.
Once the shared library is created, it must be installed for private or system-wide use before you run a program that refers to it:
-
To install a private shared library (when you are testing, for example), set the environment variable
LD_LIBRARY_PATH
, as described inld(1)
. -
To install a system-wide shared library, place the shared library file in one of the standard directory paths used by
ld
orlibtool
.
- Create the shared lib
.so
. - Compile the main program.
- Link using
gfortran
.
E.g.
gfortran -shared -fpic -o mylib.so mylib.f90
gfortran -c main.f90
gfortran -o main.exe main.o mylib.so
NOTE that in the last command, the main goes in front of the linked lib. This also holds for linking to LAPACK or BLAS: write the name of the source file first. E.g.
gfortran mylapack.f90 -llapack -lblas -o mylapack
NOTE that if the shared lib files are not in the system’s dynamic linker path, need to either modify the LD_LIBRARY_PATH
variable to include it, or give the path to the gfortran
command:
gfortran -o main.exe main.o -L/path/to/lib/ mylib.so
Debugging
gdb – GNU debugger
Generally comes with Linux. xxgdb
is a GUI to gdb for X window.
In order to use a debugger like gdb
to track the execution of your Fortran program, it is necessary to compile the program with the -g
option.
E.g.
suppose your source file is called foo.f
. To compile it:
gfortran -g foo.f -o foo
To start execution of a program named foo
under gdb
:
gdb foo
Then it will enter the gdb
session. Then type:
break main
run
This will start execution of the program, but execution will pause just before the first executable statement.
Commands in gdb
and xxgdb
:
break
: set a breakpoint.run
: begin.cont
(orc
) : continue.next
(orn
) : next line.step
(ors
) : step in.print
(orp
) expr : display the value of the expression.list
(orl
) : same asl
in pdb.
dbx debugger
run
: begin.cont
(orc
) : continue.next
(orn
) : next line.step
(ors
) : step in.stop [var]
: breakpoint when the variable [var] changes.stop in [proc]
: stop when procedure [proc] is entered.stop at [line]
: breakpoint at line [line].
Fortran libraries
Some free libraries:
- RANDLIB: random number and statistical distribution generators
- BLAS
- EISPACK
- GAMS–NIST: Guide to Available Math Software
- Some statistical and other routines from NIST
- LAPACK
- LINPACK
- MINPACK
- MUDPACK
- NCAR Mathematical Library
- The Netlib collection of mathematical software, papers, and databases.
- ODEPACK
- ODERPACK: a set of routines for rankning and ording
- Expokit for computing matrix exponentials
- SLATEC
- SPECFUN: Common Mathematical Library: a comprehensive software library containing over 1400 general purpose mathematical and statistical routines written in Fortran 77 (http://www.netlib.org/slatec/index.html) and Fotran 90 (http://people.sc.fsu.edu/~jburkardt/f_src/slatec/slatec.html).
- STARPAC
- StatLib statistical library
- TOMS
- Sorting and merging strings
Modules
A module serves as a packaging means for subprogrames, data and
interface blocks, a new feature of Fortran 90. They make common
blocks
(routinely used in Fortran 77) and include
statements obsolete.
module structure
A typical module has 2 parts:
-
declaration (including
interface
blocks,type
andparameter
declarations). -
a
contains
part for subroutine and function definitions.module mod_name implicit none [declarations] contains ! function definitions ! subroutine definitions end contains end module
A module can have just the declaration or subprogram part, or both.
Procedures defined inside a program are internal procedures.
An internal procedure is local to its host, and have access to host’s variables.
It is often better to use external (define outside main program) rather than internal, because they can be called from more than one program, and they are safer that variables from the calling program are hidden from the procedure.
If implicit none
is used at the beginning of module definition, it’s
not necessary to put implicit none
in the subprograms in the module.
When a subroutine is defined in a module, then there is no need to
provide an explicit interface in the calling program (as long as the
module’s contents are made available to the program via the use
statement).
use module
To use a module in a program or a subroutine
use mod_name1
use mod_name2
...
NOTE use
statement has to be before implicit none
in the main program (only comments are allowed before use
).
accessibility of module variables and procedures
By default, all variables, functions and subroutines are available to
the program that uses the module, by the use
statement.
- Variables declared in a module are global to the module.
- Variables declared in a module are global in the importing program/routine.
Can use the private
and public
attributes to control accessibility.
E.g.
module my_module
implicit none
real, parameter, private :: pi = 3.1415
contains
subroutine show_const()
write (*,*) pi
end subroutine
end module
The default accessibility of the module can be set by the public
or
private
statements.
If private
is specified, then all module contents are private, except those that are explicitly defined as public:
E.g.
module my_module
implicit none
private
integer :: status, count
! public variables
public :: ndim, nvar, max_number, data_array
! public procedures
public :: getData, getList
contains
subroutine getData()
...
end subroutine
subroutine getList()
...
end subroutine
end my_module
Selecting module elements
Generally, the use
statement makes available all (public) elements of a module.
However, when only a subset of the module is needed, the accessibility
can be restricted with only
.
E.g.
program convert_temperature
use convertT, only: CtoF, FtoC
...
end program convert_temperature
Aliasing element names in use
use A, only: Afoo=>foo
is equivalent to
from A import foo as Afoo
Data encapsulation
Modules allow data and operations to be hidden from the rest of the program. Data encapsulation functions as a security tool, because the data in the object is only available through the methods.
E.g.
module student_class
implicit none
private
public :: create_student, get_mark
type student_data
character (len=50) :: name
real :: mark
end type student_data
type (student_data), dimension (100) :: student
contains
subroutine create_student (student_n, name, mark)
! here some code to set the name and mark of a student
end subroutine create_student
subroutine get_mark (name, mark)
! dummy arguments
character (len=*), intent (in) :: name
real, intent (out) :: mark
! local variables
integer :: i
do i = 1,100
if (student(i) == name) then
mark = student(i)%mark
end if
end do
end subroutine get_mark
end module students
The student_class
module defines a data type (student_data
) to hold
information of a student (name and a mark). Only the subroutines,
create_student
and get_mark
, are accessible from outside the module,
all other module contents are private. Thus, one cannot obtain the
mark of a student by writing:
mark1 = student(1)%mark
because the array student
is private.
Global data management – no more common blocks
Procedures can communicate with each other via their argument lists. However, a program may consist of many procedures that require access to the same data. It would be convenient if this data were globally accessible to the whole program.
In FORTRAN 77, this was accomplished by common
blocks. However,
modules can replace all uses of common
blocks. Global data can be
packed in a module, and all procedures requiring this data can simply
use the module.
Modules are much safer and cleaner than common
blocks. Common
blocks
have no mechanisms to check errors, variables can be renamed
implicitly, and there are no access restrictions. So, don’t use common
blocks, use modules instead!
Pointers
A pointer in Fortran is a data object with the pointer attribute, containing information about a particular object, such as its type, rank, extent, and memory address.
A pointer variable is declared with the pointer
attribute. A pointer
variable that is an array must be a deferred-shape
array (only the
rank of the array is specified).
E.g.
integer, pointer :: p ! pointer to integer
real, pointer, dimension(:) :: rp ! pointer to 1d array
real, pointer, dimension(:,:) :: rp2 ! pointer to 2d array
In contrast to a normal data object, a pointer has initially no space set aside for its contents. It can only be used after space has been associated with it.
A target is the space that becomes associated with the pointer.
A pointer can point to:
- an area of dynamically allocated memory, as illustrated in the next section.
- a data object of the same type as the pointer, with the target attribute (see section on targets below)
allocating space for a pointer
E.g.
integer, pointer :: p
allocate(p)
p=1
The allocated storage space can be deallocated by the deallocate
statement. It is a good idea to deallocate storage space that is not
any more needed, to avoid accumulation of unused and unusable memory
space.
target
To act as a target for a pointer it must be declared with the target attribute. This is to allow code optimization by the compiler.
integer, pointer :: p1
integer, target :: t1
p1=>t1
p1=1
After the statement `p1=>t1`, the pointer p1
acts as an alias for t1
. Changing p1
has the effect of changing t1
.
Association
The association status of a pointer can be:
- undefined
- associated
- disassociated
If the association status is not undefined, it can be tested using the associated function:
associated(ptr)
: check if pointer ptr is associated with any target.associated(ptr,target)
: check if pointer is associated with target.
A pointer can be explicitly disassociated with a target using the
nullify
statement:
nullify(ptr)
It is a good idea to nullify pointers instead of leaving their status undefined, because they can then be tested with the associated function.
Nullify does not deallocate the targets (because there can be more than one pointer pointing to the same target). Deallocate implies nullification as well.
Linked lists
A linked list is a special kind of data storage structure. It is a series of derived type that are linked together by pointers.
There are several kinds of linked lists (single-linked lists, double-linked lists, binary trees).
Each element (also called node or link) of a linked list is an object of derived type that consists of a part with data and a pointer to the next object of the same list:
The pointer is of the same type as the other elements of the list.
The derived type can for example be something like:
type node
integer :: i
real :: value
type(node), pointer :: next
end type node
The “next” pointer of the last link in the list should be nullified. You also need a pointer (often referred to as head pointer) that refers to the first item in the list.
Differences between linked lists and arrays
- Linked lists can be allocated dynamically, so you don’t need to know before the program is executed how many elements are needed (this also saves memory space).
- The size of the list can change during execution (links can be added and removed).
- Links can be added at any position in the list.
- The links are not necessarily stored contiguously in memory.
Example
program linkedlist
implicit none
type :: link
integer :: i
type (link), pointer :: next
end type link
type (link), pointer :: first, current
integer :: number
nullify (first)
nullify (current)
! read in a number, until 0 is entered
do
read*, number
if (number == 0) then
exit
end if
allocate (current) ! create new link
current%i = number
current%next => first ! point to previous link
first => current ! update head pointer
end do
! print the contents of the list
current => first
! point to beginning of the list
do
if (.not. associated (current)) then ! end of list reached
exit
end if
print*, current%i
current => current%next ! go the next link in the list
end do
end program linkedlist
In this program a link is defined that can hold an integer. The pointer “first” is the head pointer. In the first do loop, numbers are read in until a 0 is entered. After each number is read in, a new link is created and added before the previous link.
INCLUDE vs USE module
See this explanation from stackoverflow.
The short conclusion is: prefer modules to include
in most cases.
One possible valid case to use include is: when a single module file
gets too long, then separate it into different files, and use include
in the module file to package them together.
Fortran namelist
namelist input format
! comments
&GROUP_NAME
fields = 'cin' ! use single/double quotes for strings
mercator_defs = .true. ! no need to use quotes on boolean
interp_method = 1 ! no need to use qutoes on numbers
interp_levels = 1000.,900.,500. ! use comma to separate
...
/
&GROUP_NAME2
...
/
- Order of key-value pairs in group is not important.
- Order of groups is important.
Usage of namelist
In the main program:
namelist /namelist_name/ key1, key2, ...
open(1,file='namelist_file')
read(1, namelist_name)
Code snippets
linspace
subroutine linspace(x1, x2, n, arr)
! linspace, real
implicit none
integer :: n, i
real :: x1, x2, dx
real, dimension(n), intent(out) :: arr
dx=(x2-x1)/n
arr=(/ (i*dx, i=1, n) /)
end subroutine linspace
Measure execution time
Use CPU_TIME()
to get cpu time.
E.g.
real :: t1, t2
call CPU_TIME(t1)
! do stuff
call CPU_TIME(t2)
write (*,*) 'cpu time:', t2-t1
pretty print of 2d matrix
NOTE that this subroutine has optional arguments and inferred input dimension (arr), need to put it in a module, or add an interface block.
subroutine print2D_r(arr,pre,space)
! Pretty print of 2D matrix
implicit none
real, dimension(:,:), intent(in) :: arr
integer, optional :: pre, space ! precision, spacing
character(len=20) :: f_str
character(len=1) :: s1,s2
integer :: nx, ny, i, j
integer :: pre2, space2 ! max len, spacing
if (present(pre) .eqv. .false.) then
pre2=4
else
pre2=pre
end if
if (present(space) .eqv. .false.) then
space2=4
else
space2=space
end if
ny=size(arr,1)
nx=size(arr,2)
!--------------Compose format string--------------
! f_str = '(f8.n,mx)'
! n = decimal places (pre)
! m = spaces between numbers (space)
write (s1, "(I1)") pre2
write (s2, "(I1)") space2
f_str='(f8.'//s1//','//s2//'x)'
do i = 1,ny
do j = 1,nx
write(*,trim(f_str),advance='no') arr(i,j)
end do
write(*,*)
end do
write(*,*)
end subroutine print2D_r
Trouble shooting
use with netcdf fortran
Install netcdf c and fortran libraries from the package manager.
To compile a script with use netcdf
in it:
f95 -o outfile -I/usr/include fortran_code.f90 -lnetcdff -lnetcdf
or try this:
gfortran -o outfile -I/usr/include/ -L/usr/lib/ -lnetcdf simple_xy_wr.f90 /usr/lib/libnetcdff.a
See the explanations here: https://stackoverflow.com/questions/13941549/compiling-fortran-netcdf-programs-on-ubuntu
Specifically:
Ordering of object files and archives on the linker command line is very important on Unix systems since the default linker behaviour is to search for symbol definitions only in archives that follow the object file or archive, where an unresolved reference was found, referred to single pass linking.
This means that if your code references _<sub>netcdf</sub><sub>MOD</sub><sub>nf90</sub><sub>strerror</sub>, then the archive that contains the definition of this symbol (libnetcdff.a) must appear after the list of object files from your program. libnetcdff.a itself references symbols from the C library libnetcdf.a, hence it must be linked after libnetcdff.a. So the correct link order is:
f95 -o xy -I/usr/include/ simple<sub>xy</sub><sub>wr.f90</sub> -lnetcdff -lnetcdf