Difference between revisions of "Modern Fortran"
m |
|||
Line 13: | Line 13: | ||
* [[#A Brief History of Fortran|Intro, History]] (10 min) | * [[#A Brief History of Fortran|Intro, History]] (10 min) | ||
* [[#New Format, New Syntax|New syntax]] (30 min) | * [[#New Format, New Syntax|New syntax]] (30 min) | ||
− | * [[#Hands | + | * [[#Hands On 1|Hands On 1]](60 min) |
* [[#Procedures and Modules|Functions, Modules]] (45 min) | * [[#Procedures and Modules|Functions, Modules]] (45 min) | ||
− | * Hands | + | * [[#Hands On 2|#Hands On 2]] (30 min) |
* Lunch (1 hr) | * Lunch (1 hr) | ||
* New Array Features (15 min) | * New Array Features (15 min) | ||
Line 632: | Line 632: | ||
|} | |} | ||
− | ==Hands | + | ==Hands On 1== |
* In workedexample/f77 is a simplified, F77ized version of a fluid-dynamics code from Ue-Li Pen, CITA, U of Toronto ([http://www.cita.utoronto.ca/~pen/MHD/ http://www.cita.utoronto.ca/~pen/MHD/]) | * In workedexample/f77 is a simplified, F77ized version of a fluid-dynamics code from Ue-Li Pen, CITA, U of Toronto ([http://www.cita.utoronto.ca/~pen/MHD/ http://www.cita.utoronto.ca/~pen/MHD/]) | ||
* For the purposes of this class, we've turned it from a perfectly good f90 code to something that looks more like something your supervisor would dust off and give to you. | * For the purposes of this class, we've turned it from a perfectly good f90 code to something that looks more like something your supervisor would dust off and give to you. | ||
Line 642: | Line 642: | ||
* ~1 hr (for getting logged in and everything working) | * ~1 hr (for getting logged in and everything working) | ||
− | =Procedures and | + | =Procedures and Modules= |
{| | {| | ||
| | | | ||
Line 736: | Line 736: | ||
==Private and public== | ==Private and public== | ||
− | * Not all of a module’s | + | {| |
− | content need be public | + | |- |
− | * Can give individual | + | | |
− | items public or private | + | * Not all of a module’s content need be public |
− | attribute | + | * Can give individual items public or private attribute |
− | * “private” makes | + | * “private” makes everything private by default |
− | everything private by | + | * Allows hiding implementation specific routines |
− | default | + | | |
− | * Allows hiding | + | ( from samples/procedures/multifilemod/gravity.f90 ) |
− | + | |} | |
− | samples/procedures/multifilemod/gravity.f90 | ||
==Procedures== | ==Procedures== | ||
− | * We’ve already seen | + | {| |
− | procedures defined in | + | |- |
− | new style; let’s look | + | | |
− | more closely. | + | * We’ve already seen procedures defined in new style; let’s look more closely. |
− | * Biggest change: intent | + | * Biggest change: intent attribute to “dummy variables” (eg, parameters passed in/out). |
− | attribute to “dummy | + | * Again, make expectations more explicit, compiler can catch errors, optimize. |
− | variables” (eg, | + | * Intent(in) - read only. Error to change. |
− | parameters passed in/ | + | * Intent(out) - write only. Value undefined before set. |
− | out). | + | * Intent(inout) - read/write. (eg, modify region of an array) |
− | + | * Also documentation of a sort. | |
− | + | | | |
− | + | ( from samples/procedures/funcsub/procedures.f90 ) | |
− | * Again, make expectations | + | |} |
− | more explicit, compiler can | ||
− | catch errors, optimize. | ||
− | * Intent(in) - read only. Error | ||
− | to change. | ||
− | * Intent(out) - write only. | ||
− | Value undefined before set. | ||
− | * Intent(inout) - read/write. | ||
− | (eg, modify region of an | ||
− | array) | ||
− | * Also documentation of a | ||
− | sort. | ||
− | samples/procedures/funcsub/procedures.f90 | ||
==Functions== | ==Functions== | ||
+ | {| | ||
+ | |- | ||
+ | | | ||
* Can be typed a couple of ways. | * Can be typed a couple of ways. | ||
− | * Old-style still works (real | + | * Old-style still works (real function square..) |
− | function square..) | + | * Give a result variable different from function name; set that, type that result (xsquared) |
− | * Give a result variable different | + | * Explicitly type the function name, set that as return value (cube) |
− | from function name; set that, | + | * Function return values don’t take intent - always out |
− | type that | + | | |
− | result (xsquared) | + | ( from samples/procedures/funcsub/procedures.f90 ) |
− | + | |} | |
− | * Explicitly type the function | ||
− | name, set that as return value | ||
− | |||
− | * Function values don’t take | ||
− | |||
− | samples/procedures/funcsub/procedures.f90 | ||
==Procedure interfaces== | ==Procedure interfaces== | ||
− | * The interface to a procedure | + | {| |
− | consists of | + | |- |
− | * A procedure’s name | + | | |
− | * The arguments, their names, | + | * The interface to a procedure consists of |
− | types and all attributes | + | ** A procedure’s name |
− | * For functions, the return value | + | ** The arguments, their names, types and all attributes |
− | name and type | + | ** For functions, the return value name and type |
− | * Eg, the procedure, with all the | + | * Eg, the procedure, with all the real code stripped out. |
− | real code stripped out. | + | * Like a C prototype, but more detailed info |
− | * Like a C prototype, but more | + | * .mod files contain explicit interfaces to all public module procedures. |
− | detailed info | + | |} |
− | * .mod files contain explicit | ||
− | interfaces to all public module | ||
− | procedures. | ||
==Procedure interfaces== | ==Procedure interfaces== | ||
− | * To see where interfaces | + | {| |
− | become necessary, consider | + | |- valign="top" |
− | this sketch of a routine to do | + | | |
− | trapezoid-rule integration | + | * To see where interfaces become necessary, consider this sketch of a routine to do trapezoid-rule integration |
− | * We want to | + | * We want to integrate a passed-in function f, but we don’t know anything about it - type, # of arguments, etc. |
− | function f, but we don’t know | + | * Need to “type” f the same way you do with xlo, xhi, n. |
− | anything about it - type, # of | + | * You do that for procedures with interfaces |
− | arguments, etc. | + | | |
− | * Need to “type” f the same way | + | http://en.wikipedia.org/wiki/File:Trapezoidal_rule_illustration_small.svg |
− | you do with xlo, xhi, n. | + | |} |
− | * You do that for procedures | ||
− | with interfaces | ||
− | http://en.wikipedia.org/wiki/ | ||
− | File:Trapezoidal_rule_illustration_small.svg | ||
==Procedure interfaces== | ==Procedure interfaces== | ||
− | * Define f as a parameter, give its | + | {| |
− | type via an interface. | + | |- valign="top" |
− | * Can then use it, and at compile | + | | |
− | time compiler ensures function | + | * Define f as a parameter, give its type via an interface. |
− | passed in matches | + | * Can then use it, and at compile time compiler ensures function passed in matches thisinterface. |
− | + | | | |
− | + | (from samples/procedures/interface/integrate.f90 ) | |
− | integrate.f90 | + | |} |
==Recursive procedures== | ==Recursive procedures== | ||
− | * By default, Fortran procedures | + | {| |
− | cannot call themselves | + | |- valign="top" |
− | (recursion) | + | | |
− | * Can be enabled by giving the | + | * By default, Fortran procedures cannot call themselves (recursion) |
− | procedure the recursive | + | * Can be enabled by giving the procedure the recursive attribute |
− | attribute | ||
* Subroutines, functions | * Subroutines, functions | ||
− | * Recursive functions must use | + | * Recursive functions '''must''' use “result” keyword to return value. |
− | “result” keyword to return | + | | |
− | value. | + | ( from samples/procedures/recursive/integrate.f90) |
− | + | |} | |
− | samples/procedures/recursive/integrate.f90 | ||
==Pure procedures== | ==Pure procedures== | ||
− | * Procedures are pure or | + | {| |
− | impure depending on | + | |- valign="top" |
− | whether or not they have | + | | |
− | “side effects”: | + | * Procedures are pure or impure depending on whether or not they have “side effects”: |
− | * Changing things other | + | ** Changing things other than their dummy arguments |
− | than their dummy | ||
− | arguments | ||
* Modifying save variables | * Modifying save variables | ||
* Modifying module data | * Modifying module data | ||
* Printing, etc. | * Printing, etc. | ||
− | + | * Optimizations can be made for pure routines which can’t for impure | |
− | + | * Label known-pure routines with the pure attribute. | |
− | + | * Almost all the procedures we’ve seen so far are pure. | |
− | * Optimizations can be made | + | | |
− | for pure routines which | + | (from samples/procedures/purity/purity.f90) |
− | can’t for impure | + | |} |
− | * Label known-pure routines | ||
− | with the pure attribute. | ||
− | * Almost all the procedures | ||
− | we’ve seen so far are pure. | ||
− | samples/procedures/purity/purity.f90 | ||
==Optional Arguments== | ==Optional Arguments== | ||
− | * Can make | + | {| |
− | arguments optional | + | |- valign="top" |
− | by using the optional | + | | |
− | attribute. | + | * Can make arguments optional by using the optional attribute. |
* Use present to test. | * Use present to test. | ||
− | * Can’t use tol if not | + | * Can’t use tol if not present; have to use another variable. |
− | present; have to use | + | | |
− | another variable. | + | (from samples/procedures/optional/integrate.f90) |
− | samples/procedures/optional/integrate.f90 | + | |} |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
==Keyword Arguments== | ==Keyword Arguments== | ||
− | * To avoid ambiguity with | + | {| |
− | omitted arguments - or | + | |- valign="top" |
− | really whenever you | + | | |
− | want - you can specify | + | * When calling the procedure, can use the optional argument or not. |
− | which value is which | + | * Makes sense to leave optional arguments at end - easier to figure out what’s what when it’s omitted. |
− | explicitly. | + | * To avoid ambiguity with omitted arguments - or really whenever you want - you can specify which value is which explicitly. |
− | * Don’t have to be in | + | * Don’t have to be in order. |
− | order. | + | * Can clarify calls of routines with many arguments. |
− | * Can clarify calls of | + | | |
− | routines with many | + | (from samples/procedures/optional/optional.f90) |
− | + | |} | |
− | + | ==Procedures and Modules Summary== | |
− | + | {| | |
− | ==Procedures | + | |- valign="top" |
− | + | | | |
− | * Modules let you bundle procedures, constants | + | * Modules let you bundle procedures, constants in useful packages. |
− | in useful packages. | ||
* Can have public, private components | * Can have public, private components | ||
− | * Compiling them generates a .mod file | + | * Compiling them generates a .mod file (needed for compiling anything that does a “use modulename”) and an .o file (where the code goes, needed to link together the program). |
− | (needed for compiling anything that does a | + | * New syntax for functions/subroutines: intent (IN/OUT/INOUT) |
− | “use modulename”) and an .o file (where the | + | * New syntax for function return values; result or explicit typing of function in argument list. |
− | code goes, needed to link together the | + | * Procedures have interfaces, which are needed for (eg) passing functions |
− | program). | ||
− | |||
− | |||
− | |||
− | * New syntax for functions/subroutines: intent | ||
− | (IN/OUT/INOUT) | ||
− | * New syntax for function return values; result | ||
− | or explicit typing of function in argument list. | ||
− | * Procedures have interfaces, which are needed | ||
− | for (eg) passing functions | ||
* Optional/keyword arguments | * Optional/keyword arguments | ||
* Pure/recursive procedures | * Pure/recursive procedures | ||
+ | |} | ||
− | ==Hands | + | ==Hands On 2== |
− | * In workedexamples/modules, have have | + | {| |
− | pulled the PBM stuff out into a module. | + | |- valign="top" |
− | * Do the same with the hydro routines. What | + | | |
− | needs to be private? Public? | + | * In workedexamples/modules, have have pulled the PBM stuff out into a module. |
− | * The common block (thankfully) only contains | + | * Do the same with the hydro routines. What needs to be private? Public? |
− | constants, can make those module parameters | + | * The common block (thankfully) only contains constants, can make those module parameters |
* ~30 min | * ~30 min | ||
+ | |} | ||
==Fortran arrays== | ==Fortran arrays== |
Revision as of 22:38, 13 August 2012
Modern Fortran for FORTRAN77 users
Jonathan Dursi
This is the Wiki-fied version of the slides used for the Modern Fortran course at SciNet. To follow along, with both the examples and the hands on, you will want to download the source code and make sure you have a working Fortran 2003 compiler.
Course Overview
|
A Brief History of Fortran
- Only major compiled programming language designed specifically for scientific programming.
- Powerful array operations; many mathematical functions (Bessel functions!) built in; designed to enable compiler optimizations for fast code
- Oldest (54-57 yrs) still-used programming language.
- Most people come to Fortran via being given old code by someone.
- Can’t understand the old code, or quirks of modern language, without understanding it’s history
- 1957 - J.W. Backus et al. In Proceedings Western Joint Computer Conference, Los Angeles, California, February 1957.
- IBM 704
- (Arguably) first modern compiled programming language.
- Idea of compilers at all was controversial at time.
FORTRAN (1957)
|
Incremental changes
|
FORTRAN66
|
FORTRAN77
- The most common to see “in the wild” of old code today
- if/else/endif, better do loops, control of implicit typing
- Character strings, saved variables, IO improvements
- Approved in 1978, beginning long tradition of “optimistic” naming of standards by year.
The interregnum
- Programming languages and techniques were moving quite quickly
- Several attempts were made to make new version, but standardization process very slow, failed repeatedly.
- Absent new real standard, implementations began to grow in many different directions
- Some extensions became quasi-standard, many were peculiar to individual compilers.
Fortran90
- Enormous changes; the basis of modern Fortran (not FORTRAN!)
- Free form, array slices, modules, dynamic memory allocation, derived types...
- Changes so major that took several years for compilers to catch up.
- Modern fortran
And since...
- Fortran95 - modest changes to Fortran90, killed off some deprecated F77 constructs.
- Fortran 2003 - bigger change; object-oriented, C interoperability. Most compilers have pretty good F2003 support.
- Fortran 2008 - mostly minor changes, with one big addition (coarray), other parallel stuff. Compiler-writers getting there.
F90, F95, F2003, F2008..
- We won’t distinguish between versions; we’ll just show you a lot of useful features of modern fortran.
- Will only show widely-implemented features from 2003 and 8, with exception of coarrays; these are being implemented and are very important.
New Format, New Syntax
Free Format: some highlights
|
<source lang="fortran"> program example implicit none integer, parameter :: npts = 10000 real, parameter :: startx=0., endx=1. real, parameter :: dx = (endx-startx)/npts real :: integral, xleft, xright, xmid integer :: i if (npts < 2) then print *,'Too few points!' else integral = 0. ! Simpson's Rule xleft = 0. int: do i=0,npts-1 xright = (i+1)*dx xmid = (xleft+xright)/2. integral = integral + (dx/6.)*(f(xleft) + 4.*f(xmid) + & f(xright)) xleft = xright end do int print *,'Numerical integral is ', integral print *,'Exact soln is ', (endx-startx)/2. - & (sin(2*endx)-sin(2*startx))/4. endif
function f(x) implicit none real :: f real, intent(in) :: x f = sin(x)**2 end function f end program example </source> |
Variable Declarations
<source lang="fortran"> integer i parameter (i=5) </source>
|
<source lang="fortran"> implicit none integer, parameter :: npts = 10000 real, parameter :: startx=0., endx=1. real, parameter :: dx = (endx-startx)/npts real :: integral, xleft, xright, xmid integer :: i </source> |
Variable Initialization
|
<source lang="fortran"> |
...
subroutine testvarinit implicit none integer :: i = 5 print '(A,I3)', 'On entry; i = ', i i = 7 print '(A,I3)', 'Now set; i = ', i end subroutine testvarinit |
...
use inittest |
...
call testvarinit call testvarinit |
...
end program initialization </source> ( From samples/variables/initialization/initialization.f90) <source lang="bash"> $ gfortran initialization.f90 -o initialization $ ./initialization On entry; i = 5 Now set; i = 7 On entry; i = 7 Now set; i = 7 </source> |
---|
Real Kinds
|
<source lang="fortran"> program realkinds use iso_fortran_env implicit none real :: x real(kind=real32) :: x32 real(kind=real64) :: x64 real(kind=real128):: x128 real(kind=selected_real_kind(6)) :: y6 real(kind=selected_real_kind(15)):: y15 print *,'Default:' print *, precision(x), range(x) print *,'Real32:' print *, precision(x32), range(x32) print *,'Real64:' print *, precision(x64), range(x64) print *,'Real128:' print *, precision(x128), range(x128) print *, print *,'Selected Real Kind 6:' print *, precision(y6), range(y6) print *,'Selected Real Kind 15:' print *, precision(y15), range(y15) end program realkinds </source> (from samples/variables/kinds/realkinds.f90) <source lang="bash"> $ ./realkinds Default: 6 37 Real32: 6 37 Real64: 15 307 Real128: 18 4931 Selected Real Kind 6: 6 37 Selected Real Kind 15: 15 307 </source> |
Integer Kinds
|
<source lang="fortran"> program integerkinds use iso_fortran_env implicit none integer :: i integer(kind=int8) :: i8 integer(kind=int16) :: i16 integer(kind=int32) :: i32 integer(kind=int64) :: i64 integer(kind=selected_int_kind(6)) :: j6 integer(kind=selected_int_kind(15)):: j15 print *,'Default:' print *, huge(i) print *,'Int8:' print *, huge(i8) print *,'Int16:' print *, huge(i16) print *,'Int32:' print *, huge(i32) print *,'Int64:' print *, huge(i64) print *, print *,'Selected Integer Kind 6:' print *, huge(j6) print *,'Selected Integer Kind 15:' print *, huge(j15) end program integerkinds </source> (from samples/variables/kinds/intkinds.f90) <source lang="bash"> $ ./intkinds Default: 2147483647 Int8: 127 Int16: 32767 Int32: 2147483647 Int64: 9223372036854775807 Selected Integer Kind 6: 2147483647 Selected Integer Kind 15: 9223372036854775807 </source> |
Strings
|
<source lang="fortran"> program strings implicit none character(len=20) :: hello character(len=20) :: world character(len=30) :: helloworld hello = "Hello" world = "World!" helloworld = trim(hello) // " " // trim(world) print *, helloworld if (hello < world) then print *, '<', hello, '> is smaller.' else print *, '<', world, '> is larger.' endif end program strings </source> (from samples/variables/strings/strings.f90) <source lang="bash"> $ ./strings Hello World! <Hello > is smaller. </source> |
Array declarations
|
<source lang="fortran"> program arrays implicit none real, dimension(3) :: x, y x = [1,2,3] y = 2*x print *, x print *, y end program arrays </source> ( from samples/variables/arrays/arrays.f90) <source lang="bash"> $ ./arrays 1.0000000 2.0000000 3.0000000 2.0000000 4.0000000 6.0000000 </source> |
Do loops
|
<source lang="fortran"> program doi implicit none integer :: i do i=1,10 print *, i, i**2, i**3 enddo end program doi </source> ( from samples/variables/doloops/doi.f90 ) <source lang="bash"> $ ./doi 1 1 1 2 4 8 3 9 27 4 16 64 5 25 125 6 36 216 7 49 343 8 64 512 9 81 729 10 100 1000 </source> |
Named loops
|
<source lang="fortran"> program nameddo implicit none integer :: i, j outer: do i=1,3 inner: do j=1,3 print *, i, j, i*i+j*j enddo inner end do outer end program nameddo </source> ( from samples/variables/doloops/nameddo.f90 ) <source lang="bash"> $ ./nameddo 1 1 2 1 2 5 1 3 10 2 1 5 2 2 8 2 3 13 3 1 10 3 2 13 3 3 18 </source> |
Cycle/exit
|
<source lang="fortran"> program cycleexit implicit none integer :: i do print *, 'Enter a number between 1-13' read *, i if (i>=1 .and. i<=13) exit print *, 'Wrong; try again.' enddo print *, 'Good; you entered ', i end program cycleexit </source> ( from samples/variables/doloops/cycleexit.f90 ) <source lang="bash"> $ more cycleexit-out.txt $ ./cycleexit Enter a number between 1-13 23 Wrong; try again. Enter a number between 1-13 -1 Wrong; try again. Enter a number between 1-13 12 Good; you entered 12 </source> |
Do while
|
<source lang="fortran"> program dowhile implicit none integer :: i i = -1 do while (i < 1 .or. i > 13) print *, 'Enter a number between 1-13' read *, i if (i<1 .or. i>13) print *, 'Wrong; try again.' enddo print *, 'Good; you entered ', i end program dowhile </source> ( from samples/variables/doloops/dowhile.f90 ) |
Hands On 1
- In workedexample/f77 is a simplified, F77ized version of a fluid-dynamics code from Ue-Li Pen, CITA, U of Toronto (http://www.cita.utoronto.ca/~pen/MHD/)
- For the purposes of this class, we've turned it from a perfectly good f90 code to something that looks more like something your supervisor would dust off and give to you.
- Today we’ll be translating this version into a very modern Fortran
- Compile (using make) and run (./hydro)
- Outputs a .pbm file; use “display dens.pbm” to see the result of dense blob of fluid moving through a light medium.
- In workedexamples/freeform, have partly converted the program to new freeform format, with enddos, ending procedures, implicit nones, and new variable declaration syntax.
- Finish doing so - just need to do program hydro, subroutine color, subroutine outputpbm, function cfl. Fix indenting (Don't need to start at col 7 anymore).
- ~1 hr (for getting logged in and everything working)
Procedures and Modules
|
Modules
|
samples/procedures/simplemod/simplemod.f90 |
Compiling & Running
|
Modules
|
samples/procedures/simplemod/simplemod.f90 |
use module, only :
- Best practice is to only pull in from the module what you need
- Otherwise, everything.
- Adds some clarity and documentation, good for maintainability
- (Note syntax for continuation of a string...)
samples/procedures/simplemod/simplemod2.f90 |
Modules usually get their own files
|
samples/procedures/multifilemod/gravity.f90 |
Modules usually get their own files
- Compiling gravity.f90 now gives both an .o file (containing the code) and the .mod file as before.
- Compiling the main program (multifilemod.f90) requires the .mod file.
samples/procedures/multifilemod/Makefile |
.mod needed for compilation
|
samples/procedures/multifilemod/multifilemod.f90 |
.o needed for linking
- Linking, however, doesn’t require the .mod file
- Only requires the .o file from the module code.
- .mod file analogous (but better than) .h files for C code.
samples/procedures/multifilemod/Makefile |
Compiling and running
- So compile files with modules first, so those that need them have the .mod files
- Link the .o files
Private and public
|
( from samples/procedures/multifilemod/gravity.f90 ) |
Procedures
|
( from samples/procedures/funcsub/procedures.f90 ) |
Functions
|
( from samples/procedures/funcsub/procedures.f90 ) |
Procedure interfaces
|
Procedure interfaces
|
http://en.wikipedia.org/wiki/File:Trapezoidal_rule_illustration_small.svg |
Procedure interfaces
|
(from samples/procedures/interface/integrate.f90 ) |
Recursive procedures
|
( from samples/procedures/recursive/integrate.f90) |
Pure procedures
|
(from samples/procedures/purity/purity.f90) |
Optional Arguments
|
(from samples/procedures/optional/integrate.f90) |
Keyword Arguments
|
(from samples/procedures/optional/optional.f90) |
Procedures and Modules Summary
|
Hands On 2
|
Fortran arrays
- Fortran made for dealing with scientific data
- Arrays built into language
- The type information associated with an
array includes rank (# of dimension), size, element type, stride..
- Enables powerful optimizations,
programmer-friendly features.
Fortran arrays
- Can be manipulated
like simple scalar variables
- Elementwise
addition, multiplication.. samples/arrays/basic.f90
Array constructors
- Can have array
constants like numerical constants
- use [] or (/ /), then
comma-separated list of values.
- Implied do loops
can be used in constructors
- (Variables have to
be defined)
[1,2,3,4,5] or (/1,2,3,4,5/) [ (i,i=1,5)] [ ((i*j,j=1,3),i=1,5)]
Elementwise operations
- Elementwise operations
can be */+-, or application of an elemental function.
- Math intrinsics are all
elemental - applied to array, applies to every element.
- Order of execution
undefined - allows vectorization, parallelization.
samples/arrays/elementwise.f90
Elemental Functions
- User can create their
own elemental functions
- Label any scalar function
with “elemental” should (until recently, must) be pure, so can be applied everywhere at same time.
- Faster than in loop.
- Can also take multiple
arguments: eg z = addsquare(x,y)
samples/arrays/elemental.f90
Array comparisons
- Array comparisons
return an array of logicals of the same size of the arrays.
- Can use any and all to
see if any or all of those logicals are true. samples/arrays/compare.f90
Array masks
- These logical arrays can
be used to mask several operations
- Only do sums, mins, etc
where the mask is true
- Eg, only pick out positive
values.
- Many array intrinsics
have this mask option
samples/arrays/mask.f90
Where construct
- The where construct
can be used to easily manipulate sections of array based on arbitrary comparisons.
- Where construct => for
whatever indices the comparison is true, set values as follow; otherwise, set other values.
samples/arrays/where.f90
Forall construct
- Forall is an array
assignment statement
- Each line in forall has to be
independent. All done “at once” - no guarantees as to order
- If (say) 2 lines in the forall,
all of the first line is done, then all of the second.
- Any functions called must
be pure
- Can be vectorized or
parallelized by compiler
samples/arrays/forall.f90
Array Sections
- Generalization of array
indexing
- Familiar to users of
Matlab, IDL, Python..
- Can use “slices” of an
array using “index triplet”
- [start]:[end][:step]
- Default start=1, default
end=size, default step=1.
- Can be used for each
index of multid array
a([start]:[end][:step]) a = [1,2,3,4,5,6,7,8,9,10] a(7:) == [7,8,9,10] a(:3) == [1,2,3] a(2:4) == [2,3,4] a(::3) == [1,4,7,10] a(2:4:2) == [2,4] a(2) == 2 a(:) == [1,2,3,4,5,6,7,8,9,10]
Array Sections
- This sort of thing is very
handy in numerical computation
- Replace do-loops with
clearer, shorter, possibly vectorized array operations
- Bigger advantage for
multidimensional arrays
samples/arrays/derivative.f90
Array Sections
- The previous sorts of array
sections - shifting things leftward and rightward - are so common there are intrinsics for them
- +ve shift shifts elements
leftwards (or array bounds rightwards).
- cshift does circular shift - shifting
off the end of the array “wraps around”.
- eoshift fills with zeros, or
optional filling.
- Can work on given dimension
a = [1,2,3,4,5] cshift(a,1) == [2,3,4,5,1] cshift(a,-1) == [5,1,2,3,4] eoshift(a,1) ==[2,3,4,5,0] eoshift(a,-1)==[0,1,2,3,4]
Other important array
intrinsics
- minval/maxval - finds min, max element in an array.
- minloc/maxloc - finds location of min/max element
- product/sum - returns product/sum of array elements
- reshape - Adjusts shape of array data. Eg:
1,4 reshape([1,2,3,4,5,6],[3,2]) == 2,5 3,6
Linear algebra in Fortran
- Comes built in with transpose, matmul,
dot_product for dealing with arrays.
- matmul also does matrix-vector multiplication
- Either use these or system-provided BLAS
libraries - never, ever write yourself.
Matmul times
$ ./matmul 2500 Experiment with matrix size 2500 Triple-loop time: 149.63400 matmul intrinsic time: 10.370000 SGEMM lapack call time: 1.4809999
(gfortran 4.6, compiled -O3 -march=native using Intel MKL 10.3 for sgemm)
samples/arrays/matmul.f90
Linear algebra in Fortran
samples/arrays/matrix.f90
Array sizes and Assumed Shape
- Printmat routine here is
interesting - don’t pass (a,rows,cols), just a.
- Can assume a rank-2 array,
and get size at runtime.
- Simplifies call, and eliminates
possible inconsistency: what if rows, cols is wrong?
- size(array,dim) gets the size of
array in the dim dimension.
samples/arrays/matrix.f90
Array sizes and Assumed Shape
- Assumed shape arrays (eg,
dimension(:,:)) much better than older ways of passing arrays: integer nx, ny integer a(nx,ny) or worse, integer a(*,ny)
- Information is thrown away,
possibility of inconsistency.
- Here, (:,:) means we know the
rank, but don’t know the size yet.
samples/arrays/matrix.f90
Allocatable Arrays
- So far, all our programs have had fixed-size
arrays, set at compile time.
- To change problem size, have to edit code,
recompile.
- Has some advantages (optimization,
determinism) but very inflexible.
- Would like to be able to request memory
at run time, make array of desired size.
- Allocatable arrays are arguably most
important addition to Fortran.
Allocate(), Deallocate()
- Give array a deferred size (eg,
dimension(:)) and the attribute allocatable.
- When time to allocate it, use
allocate(a(n)).
- Deallocate with deallocate(a).
samples/arrays/allocatable.f90
- In between, arrays can be used
as any other array.
Allocate(), Deallocate()
- If allocation fails (not enough memory available for
request), program will exit.
- Can control this by checking for an optional error code,
allocate(a(n),stat=ierr)
- Can then test if ierr>0 (failure condition) and handle
gracefully.
- In scientific programming, the default behaviour is often
fine, if abrupt - you either have enough memory to run the problem, or you don’t.
get_command_argument()
- Previous version still
depended on a compiled-in number.
- Can read from file or from
console, but Fortran now has standard way to get command-line arguments
- Get the count of arguments,
and if there’s at least one argument there, get it, read it as integer, and allocate array. samples/arrays/allocatable2.f90
get_command_argument()
samples/arrays/allocatable2.f90
Hands on #3
- Use array functionality to simplify hydro code
-- don't need to pass, array size, and can simplify mathematics using array operations.
- In workedexamples/arrays, have modified
hydro to allocate u, and pbm to just take array.
- Do the same with the fluid dynamic routines in
solver.f90
- ~30 min
Fortran Pointers
- Pointers, or references,
refer to another variable.
- Eg, p does not contain a
real value, but a reference to another real variable.
- Once associated with
another variable, can read/write to it as if it were stored “in” p.
real, target :: x = 3.2 real, pointer:: p p => x p
x 3.2
Fortran Pointers
samples/pointers/ptr1.f90
Fortran Pointers
- Pointers are either
associated, null, or undefined; start out life undefined.
- Can associate them to
a variable with => , or mark them as not associated with any valid variable by pointing it to null().
real, target :: x = 3.2 real, pointer:: p p => null() p
x
Fortran Pointers
real, target :: x = 3.2 real, pointer:: p
- Reading value from or
writing value to a null pointer will cause errors, probably crash.
p => null() p
x
Fortran Pointers
- Fortran pointers can’t
point just anywhere.
- Must reference a
variable with the same type, that has the target attribute.
real, target :: x = 3.2 real, pointer:: p p => x p
x 3.2
Fortran Pointers
real, target :: x = 3.2 real, pointer:: p1, p2
- Pointers can reference
other pointers.
- Actually references
what they’re pointing to.
p1 => x p2 => p1 p1 p2
x 3.2
Allocating a pointer
- Pointer doesn’t
necessarily have to have another variable to target
- Can allocate memory
for p to point to that does not belong to any other pointer.
- Must deallocate it when
done
real, pointer:: p allocate(p) p = 7.9 p 7.9
Allocating a Pointer
samples/pointers/ptr2.f90
What are they good
for? (1)
- Pointers are essential for
creating, maintaining dynamic data structures
- Linked lists, trees, heaps..
- Some of these can be
sort-of implemented in arrays, but very awkward
- Adaptive meshes, treebased particle solvers
need these structures.
http://en.wikipedia.org/wiki/File:Singly-linked-list.svg
http://en.wikipedia.org/wiki/File:Max-Heap.svg
What are they good
for? (2)
- A pointer can be of
array type, not just scalar
- Fortran pointers +
fortran arrays are quite interesting; can create “views” of subarrays
real, target, dimension(7) :: x real, pointer:: p(:) p => x(2:6)
p x 1
2
3
4
5
6
7
Array Views
samples/pointers/views.f90
Hands on #4
- Use pointers to provide views into subsets of
the arrays in solver.f90 to clarify the functions.
- In workedexamples/pointers, have started
the process with cfl, hydroflux; try tackling tvd1d, others.
- ~30 min
Derived Types and Objects
- Often, groups of
variables naturally go together to represent a larger structure
- Whenever you find
yourself passing the same group of variables to several routines, a good candidate for a derived type.
type griddomain real :: xmin, xmax real :: ymin, ymax real :: nx, ny real, dimension(:,:) :: u endtype griddomain type(griddomain) :: g g % xmin = -1 g % xmax = +1
Derived Types and Objects
- Consider interval
arithmetic (good for quantification of uncertainties, etc).
- An interval inherently
has two values associated with it - the end points.
- Can make this a type.
samples/derivedtypes/simple/intervalmath.f90
Derived Types and Objects
- Note can access the
fields in the type with “%”
- typename
(field1val,field2val..) initializes a value of that type.
- Can pass values of this
type to functions, etc., just like a built-in type. samples/derivedtypes/simple/intervalmath.f90
Procedures using types
- Can start creating
library of routines that operate on these new interval types.
- Procedures can take
the new type as arguments, functions can return the new type. samples/derivedtypes/intervalfunctions/intervalmath.f90
Procedures using types
- Can start creating
library of routines that operate on these new interval types.
- Procedures can take
the new type as arguments, functions can return the new type.
samples/derivedtypes/intervalfunctions/interval1.f90
Procedures using types
- Would prefer not to
have to treat integer and real intervals so differently in main program
- Different types, but
adding should be similar. samples/derivedtypes/intervalfunctions/intervalmath.f90
Procedures using types
|
samples/derivedtypes/genericintervals/interval2.f90 |
Generic Interfaces
- Generic Interfaces
- addintintervals and addrealintervals share the same interface, (two input parameters, one function return), but different types.
- Put them behind the same interface.
- Now, a call to addintervals is resolved at compile time to one or the other.
samples/derivedtypes/genericintervals/intervalmath.f90 |
Generic Interfaces
- Note that everything is
private except what is explicitly made public.
- Types are public.
- Generic interfaces are
public.
- Type specific routines
are not.
- Program using interval
math sees only the generic interfaces. samples/derivedtypes/genericintervals/intervalmath.f90
Generic interfaces
- Call to addintervals or
subtract intervals goes to the correct typespecific routine.
- As does print interval.
- Could create routines
to add real to int interval, etc and add to the same interface. samples/derivedtypes/genericintervals/interval2.f90
Operator overloading
- An infix operator is
really just “syntactic sugar” for a function which takes two operands and returns a third.
a = b (op) c => function op(b,c) returns a
Operator overloading
- An assignment
operator is really just “syntactic sugar” for a subroutine which takes two operands and sets the first from the second.
a=b => subroutine assign(a,b)
Operator overloading
- Here, we’ve defined
two subroutines which set intervals based on an array - 2 ints for an integer interval, or 2 reals for a real interval
samples/derivedtypes/intervaloperators/intervalmath.f90
Generic interfaces
- Once this is done, can
use assignment operator,
- Or add, subtract
multiply intervals.
- Can even compose
them in complex expressions! Functions automatically composed. samples/derivedtypes/intervaloperators/interval3.f90
Type bound procedures
- Types can have not only
variables, but procedures.
- Takes us from a type to
what is usually called a class.
samples/derivedtypes/intervaltypebound/ intervalmath.f90
Type bound procedures
- Called like one accesses
a field - %
- Operates on data of
the particular variable it is invoked from
samples/derivedtypes/intervaltypebound/interval3.f90
Type bound procedures
- It is implicitly passed as
it’s first parameter the variable itself.
- Can take other
arguments as well.
samples/derivedtypes/intervaltypebound/ intervalmath.f90
Object oriented
programming
- F2003 onwards can do full object oriented
programing.
- Types can be derived from other types, inheriting
their fields and type-bound procedures, and extending them.
- Goes beyond scope of today, but this is the
starting-off point.
Interoperability with
other languages
- Large scientific software now frequently uses
multiple languages, either within a single code or between codes.
- Right tool for the job!
- Need to know how to make software
interact.
- Here we’ll look at C/Fortran code calling
each other, and calling Fortran code from python.
C-interoperability
- iso_c_binding module contains definitions for
interacting with C
- Types, ways to bind procedures to C, etc.
- Allows you to call C routines, or bind Fortran
routines in a way that they can be called by C.
Calling a C routine
from Fortran
- As with the case of
calling a passed-in function, need an explicit prototype.
- Tells compiler what
to do with “sqrtc” function when called. samples/C-interop/call-c-fn/main.f90
Calling a C routine
from Fortran
- BIND(C) - tells compiler/
linker will have a C-style, rather than fortran-style name in object file.
- Can optionally supply
another name; otherwise, will default to procedure name.
- Here we’re calling it sqrtc
to avoid Fortran sqrt() function.
Calling a C routine
from Fortran
- The value the function
takes and returns are C doubles; that is, reals of kind(c_double).
- Also defined: c_float,
integers of kind c_int, c_long, etc.
Calling a C routine
from Fortran
- C function arguments by
default are passed by value; Fortran by default are passed by reference.
- Passed by value - values
copied into function
- Passed by reference pointers to values copied
in.
- value attribute for C-style
arguments.
Calling a C routine
from Fortran C math lib.
- Compiling - just make
sure to link in C library you’re calling
- And that’s it.
$ make gfortran -c main.f90 gfortran -o main main.o -lm
$ ./main x = 2.0000000000000000 C sqrt(x) = 1.4142135623730951 Fortran sqrt(x) = 1.4142135623730951
C strings
- When using C
strings, you have to remember that C terminates strings with a special character
- C_NULL_CHAR
- Dealing with functions
that return strings is hairier, as they return a pointer, not an array.
$ make gfortran-mp-4.4 -c main.f90 gfortran-mp-4.4 -o main main.o -lc $ ./main 1234 = 1234
samples/C-interop/c-strings/main.f90
Calling Fortran from C
- To write Fortran
routines callable from C, bind the subroutine to a C name
- Can optionally give it a
different name, as above
- And make sure
arguments are of C types.
samples/C-interop/c-call-fortran/froutine.f90
Calling Fortran from C
- Here, we’ll do
something a little trickier and pass C dynamic arrays
samples/C-interop/c-call-fortran/froutine.f90
Calling Fortran from C
- In Fortran, we accept
the arguments as C pointers
- We then associate
them with fortran pointers to arrays of shape [nx] (1d arrays here)
- Then we can do the
usual Fortran array operations with them.
samples/C-interop/c-call-fortran/froutine.f90
Calling Fortran from C
- The single scalar
argument passed back we just treat as an intent(out) variable
- Of type c_float.
samples/C-interop/c-call-fortran/froutine.f90
More advanced
- Handling arbitrary C code is possible
- Passing C structs -- create a Fortran derived
type with the same fields, use BIND(C) in defining the type.
- C “multidimensional arrays” - have to be
careful, they may not be contiguous! Follow pointers.
- Even taking passed function pointers to C
functions is possible (samples/C-interop/ functionptr)
Example of
Fortran calling C, which calls Fortran
samples/C-interop/valueref/driver.f90
samples/C-interop/valueref/croutine.c
samples/C-interop/valueref/froutine.f90
$ ./main
(1) In the FORTRAN routine before call to C (1) i = 15 x = 2.0000000000000000 (2) In the C routine, got i=15, *x=2.000000 (2) Changed x (3) In the FORTRAN routine called by C (3) Got product p = 30.000000000000000 (1) In the FORTRAN routine after call to C (1) i = 15 x = 3.0000000000000000 prod = 30.000000000000000
samples/C-interop/valueref/Makefile
F2py
- Comes with scipy,
a widely-installed (by scientists) python package.
- Wraps fortran in
python, allows you to easily interface.
- fwrap is another
option
* Will only use solver module.
- Unfortunately, f2py isn’t quite smart enough to understand
using parameters to size arrays, so global replace ‘nvars’=4
- module load use.experimental gcc python/2.7.1 intel/12
- “f2py -m hydro_solver -h hydro_solver.pyf solver.f90”
* generates the following header file (hydro_solver.pyf)
- Comment out stuff we don’t need (lower-level routines)
- f2py -c --fcompiler=gfortran hydro_solver.pyf solver.f90
$ ipython
In [1]: import hydro_solver In [2]: hydro_solver.solver. hydro_solver.solver.cfl hydro_solver.solver.gamma hydro_solver.solver.hydroflux hydro_solver.solver.idens hydro_solver.solver.iener hydro_solver.solver.imomx
hydro_solver.solver.imomy hydro_solver.solver.initialconditions hydro_solver.solver.timestep hydro_solver.solver.tvd1d hydro_solver.solver.xsweep hydro_solver.solver.xytranspose
In [2]: hydro_solver.solver.idens Out[2]: array(1, dtype=int32) In [3]: import numpy In [4]: u = hydro_solver.solver.initialconditions(100,100) In [5]: import pylab In [6]: pylab.imshow(u[1,:,:]) In [7]: for i in range(100) ...dt = hydro_solver.solver.timestep(u) In [8]: pylab.imshow(u[1,:,:])
Coarrays
- In F2008, objects
can also have a “codimension”.
- An object with a
codimension can be “seen” across processes
- Right now, only intel
v 12 supports this samples/coarrays/simple.f90
Coarrays
1
A B
2
A B
N
...
A B
num_images independant processes
Coarrays
1
A B
2
A B
N
...
A B
num_images independant processes
A[1] B[N]
Independent, parallel tasks can “see” each other’s data as easily as an array index.
Sychronization
- When accessing other
processor’s data, must ensure that tasks are synchronized
- Don’t want to read
task 1’s data before read is complete!
- sync all
samples/coarrays/broadcast.f90
Sychronization
- Other synchronization
tools:
- sync images([..]) syncs
only images in the list provided
- critical - only one
image can enter block at a time
- lock - finer-grained
control than critical.
- atomic operations.
samples/coarrays/broadcast.f90
1d Diffusion Eqn
- Calculate heat
diffusion along a bar
- Finite difference;
calculate second derivative using neighbours
- Get new
temperature from old temperatures
1 -2 1
1d Diffusion Eqn
- Initial conditions:
- Use pointers to
point to new, old temperatures
- 1..totpoints+2 (pts 1,
totpoints+1 “ghost cells)
- Setup x, orig
temperature, expected values
samples/coarrays/diffusion/diffusion.f90
1d Diffusion Eqn
- Evolution:
- Apply BCs
- Apply finite
difference stencil to all real data points samples/coarrays/diffusion/diffusion.f90
1d Diffusion Eqn
- Output calculated
values and theory to file output.txt
- Note: Fortran2008,
can use “newunit” to find, use an unused IO unit
samples/coarrays/diffusion/diffusion.f90
$ ./diffusion
[..] $ gnuplot [..] gnuplot> plot 'output.txt' using 1:2 with points title 'Simulated', 'output.txt' using 1:3 with lines title 'Theory'
Coarray Diffusion
- Now we’ll try this in
parallel
- Each image runs it’s part
of the problem (totpoints/num_images())
- Communications is like
boundary condition handling - except boundary data must be obtained from neighbouring image.
Coarray Diffusion
- Everything’s the same
except we are only responsible for locpoints points, starting at start.
- Once calculated, we
never need totpoints again. All arrays are of size (locpoints), etc.
- For simplicity, we
assume here everything divides evenly.
samples/coarrays/diffusion/diffusion-coarray.f90
Coarray Diffusion
- Boundary exchange; if we
have interior neighbours, get updated data from them so we can calculate our new data
- Note: can’t use pointers
here, coarrays can’t be targets.
- Sync all around boundary
exchange a little heavy handed; could just sync neighbours. samples/coarrays/diffusion/diffusion-coarray.f90
Coarray Diffusion
temperature()[leftneigh]
temperature() temperature()[rightneigh]
1 2 3 ...
lnpts+2
Coarray Diffusion
- Everything else exactly
the same.
- For I/O, we each output
our own file, prefixed by our image number
- (eg, 1-output.txt, 2output.txt...)
samples/coarrays/diffusion/diffusion-coarray.f90
$ export FOR_COARRAY_NUM_IMAGES=3
- 3 images
$ ./diffusion-coarray [..] $ gnuplot [..] gnuplot> plot '1-ics.txt' using 1:2, '2-ics.txt' using 1:2, '3-ics.txt' using 1:2
gnuplot> plot '1-output.txt' using 1:2, '2-output.txt' using 1:2, '3output.txt' using 1:2, '1-output.txt' using 1:3 with lines title 'Theory',
'2-output.txt' using 1:3 with lines title 'Theory', '3-output.txt' using 1:3 with lines title 'Theory'
Parallel Matrix Multiplication
- Consider matrix
operations, where matrix is decomposed onto images in blocks
- A given image has
corresponding blocks in A, B, responsible for that block in C.
A
B X C =
Parallel Matrix Multiplication
- Block matrix multiply exactly like matrix multiply, except each operation is a submatrix operation
- matmul instead of “*”.
- For simplicity, we’ll assume square decomposition.
A
B X C =
Parallel Matrix Multiplication
- Each image gets its own block of a, b, c.
- Note [nrows,*]. Coarrays can have any corank(eg, number of codimensions)
- Here, we make decomposition nrows x ncols.
- Can set decomposition array-by-array.
samples/coarrays/blockmatrix.f90
Parallel Matrix Multiplication
- Allocation, deallocation of coarrays is like a sync all
- Forces synchronization across all images
- Any particular coarray must have same size, co-size across all images
samples/coarrays/blockmatrix.f90
Parallel Matrix Multiplication
- This is the entire parallel multiplication operation.
- a[myrow,k] gets the entire a block from image at myrow,k.
- matmul works with those blocks, updates c with new term.
samples/coarrays/blockmatrix.f90
Coarray Summary
- Coarrays are powerful, simple-to-use addition to parallel programming models available
- Will be widely available soon
- Basic model implemented in Cray fortran for ~10 yrs; well tested and understood.
- Typically implemented with MPI under the hood; can give performance similar to MPI with fewer lines of code.
- Other languages have similar extensions (UPC based on C, Titanium based on Java) but are not part of languages’ standard.
Closing Hints
- Always give the compiler info it needs to help you by being as explicit as possible
- implicit none, end [construct] [name], parameters for constants, intent in/out, use only, etc.
- Always get as much info from compiler as possible - always use -Wall (gfortran) or -warn all (ifort).
Useful Resources
- http://fortranwiki.org/
- Reference source; has all standards; Fortran2003/2008 status of major compilers
- http://en.wikipedia.org/wiki/Fortran_language_features
- Succinct summary of new features (spotty past F95)
- http://stackoverflow.com/questions/tagged/fortran
- ( Programmers Questions & Answers