Difference between revisions of "Modern Fortran"
(23 intermediate revisions by the same user not shown) | |||
Line 12: | Line 12: | ||
| | | | ||
* [[#A Brief History of Fortran|Intro, History]] (10 min) | * [[#A Brief History of Fortran|Intro, History]] (10 min) | ||
− | * New syntax (30 min) | + | * [[#New Format, New Syntax|New syntax]] (30 min) |
− | * [[#Hands | + | * [[#Hands On 1|Hands On 1]](60 min) |
− | * 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 29: | Line 29: | ||
=A Brief History of Fortran= | =A Brief History of Fortran= | ||
{| | {| | ||
− | |||
− | |||
* Only major compiled programming language designed specifically for scientific programming. | * 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 | * Powerful array operations; many mathematical functions (Bessel functions!) built in; designed to enable compiler optimizations for fast code | ||
Line 40: | Line 38: | ||
* (Arguably) first modern compiled programming language. | * (Arguably) first modern compiled programming language. | ||
* Idea of compilers at all was controversial at time. | * Idea of compilers at all was controversial at time. | ||
+ | |} | ||
+ | {| | ||
+ | |- | ||
| | | | ||
[[File:ModernFortran-010.jpg|thumb|(from http://www.softwarepreservation.org/projects/FORTRAN/ )]] | [[File:ModernFortran-010.jpg|thumb|(from http://www.softwarepreservation.org/projects/FORTRAN/ )]] | ||
+ | | | ||
[[File:ModernFortran-013.jpg|thumb|(from http://en.wikipedia.org/wiki/File:Fortran_acs_cover.jpeg )]] | [[File:ModernFortran-013.jpg|thumb|(from http://en.wikipedia.org/wiki/File:Fortran_acs_cover.jpeg )]] | ||
|} | |} | ||
Line 126: | Line 128: | ||
==Free Format: some highlights== | ==Free Format: some highlights== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Columns no longer significant; can start at left margin | * Columns no longer significant; can start at left margin | ||
Line 146: | Line 149: | ||
* (From samples/freeform/freeform.f90) | * (From samples/freeform/freeform.f90) | ||
| | | | ||
− | + | <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 | ||
+ | |||
+ | |||
+ | contains | ||
+ | |||
+ | function f(x) | ||
+ | implicit none | ||
+ | real :: f | ||
+ | real, intent(in) :: x | ||
+ | f = sin(x)**2 | ||
+ | end function f | ||
+ | |||
+ | end program example | ||
+ | </source> | ||
|} | |} | ||
==Variable Declarations== | ==Variable Declarations== | ||
− | {| | + | {| |
+ | |- valign="top" | ||
| | | | ||
* Implicit none turns off all implicit typing. | * Implicit none turns off all implicit typing. | ||
Line 166: | Line 211: | ||
* Useful for things which shouldn’t change. | * Useful for things which shouldn’t change. | ||
* F77 equivalent: | * F77 equivalent: | ||
− | < | + | <source lang="fortran"> |
− | parameter (i=5)</ | + | integer i |
+ | parameter (i=5) | ||
+ | </source> | ||
* Initialization of variables at declaration time | * Initialization of variables at declaration time | ||
− | * Required for parameters (because can’t change them later), can be done for other variables. | + | ** Required for parameters (because can’t change them later), can be done for other variables. |
− | * Can do anything that compiler can figure out at compile time, including math with other parameters. | + | ** Can do anything that compiler can figure out at compile time, including math with other parameters. |
+ | | | ||
+ | <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== | ==Variable Initialization== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Initializing variables this way gives unexpected behaviour in functions/subroutines “for historical reasons”. | * Initializing variables this way gives unexpected behaviour in functions/subroutines “for historical reasons”. | ||
Line 186: | Line 241: | ||
* Not a problem for main program, parameters. | * Not a problem for main program, parameters. | ||
| | | | ||
− | + | <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 | ||
+ | |||
+ | ! ... | ||
+ | |||
+ | |||
+ | program initialization | ||
+ | use inittest | ||
+ | |||
+ | !... | ||
+ | call testvarinit | ||
+ | call testvarinit | ||
+ | !... | ||
− | + | end program initialization | |
− | + | </source> | |
( From samples/variables/initialization/initialization.f90) | ( 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== | ==Real Kinds== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Reals, Double precisions, really different “kinds” of same “type” - floating point real numbers. | * Reals, Double precisions, really different “kinds” of same “type” - floating point real numbers. | ||
Line 209: | Line 295: | ||
* several others | * several others | ||
| | | | ||
− | + | <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) | (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== | ==Integer Kinds== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Similar constructs for integers | * Similar constructs for integers | ||
Line 224: | Line 355: | ||
* huge(N): largest positive number of that type samples/variables/kinds/intkinds.f90 | * huge(N): largest positive number of that type samples/variables/kinds/intkinds.f90 | ||
| | | | ||
− | + | <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) | (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== | ==Strings== | ||
− | {| | + | {| |
+ | |- valign="top" | ||
| | | | ||
* Character types are usually used for strings | * Character types are usually used for strings | ||
Line 244: | Line 426: | ||
* gfortran has partial support for selected_char_kind(“ISO_10646”) for unicode strings. | * gfortran has partial support for selected_char_kind(“ISO_10646”) for unicode 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== | ==Array declarations== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Array declarations have changed, too: | * Array declarations have changed, too: | ||
Line 253: | Line 465: | ||
* Can easily declare several arrays with same dimension | * Can easily declare several arrays with same dimension | ||
| | | | ||
+ | <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> | ||
|} | |} | ||
Line 258: | Line 488: | ||
==Do loops== | ==Do loops== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Do loops syntax has had some changes | * Do loops syntax has had some changes | ||
Line 265: | Line 496: | ||
* Helps catch various simple errors (mismatched ends, etc.) | * Helps catch various simple errors (mismatched ends, etc.) | ||
| | | | ||
− | samples/variables/doloops/doi.f90 | + | <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> | ||
+ | |||
|} | |} | ||
Line 271: | Line 527: | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Can name control structures like do, if statements now, too. | * Can name control structures like do, if statements now, too. | ||
Line 277: | Line 534: | ||
* enddo or end do; fortran isn’t picky about spaces. | * enddo or end do; fortran isn’t picky about spaces. | ||
| | | | ||
− | samples/variables/doloops/nameddo.f90 | + | <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== | ==Cycle/exit== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Do loops don’t even need a i=1,n | * Do loops don’t even need a i=1,n | ||
Line 290: | Line 574: | ||
* cycle - jumps back to do | * cycle - jumps back to do | ||
| | | | ||
− | samples/variables/doloops/cycleexit.f90 | + | <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> | ||
|} | |} | ||
Line 296: | Line 609: | ||
==Do while== | ==Do while== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* do while - repeats as long as precondition is true. | * do while - repeats as long as precondition is true. | ||
* Seems like it should be useful, but in practice, just do/enddo with exit condition is usually cleaner. | * Seems like it should be useful, but in practice, just do/enddo with exit condition is usually cleaner. | ||
| | | | ||
− | samples/variables/doloops/dowhile.f90 | + | <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 | + | ==Hands On 1== |
− | * In workedexample/f77 is a simplified, F77ized version of a fluid-dynamics code from | + | * 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/]) |
− | * Today we’ll be translating | + | * 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) | * 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. | * Outputs a .pbm file; use “display dens.pbm” to see the result of dense blob of fluid moving through a light medium. | ||
Line 312: | Line 643: | ||
* ~1 hr (for getting logged in and everything working) | * ~1 hr (for getting logged in and everything working) | ||
− | =Procedures and | + | =Procedures and Modules= |
{| | {| | ||
| | | | ||
Line 322: | Line 653: | ||
==Modules== | ==Modules== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* Easiest to show by example | * Easiest to show by example | ||
Line 328: | Line 660: | ||
* “Use” goes before “implicit none” | * “Use” goes before “implicit none” | ||
| | | | ||
− | samples/procedures/simplemod/simplemod.f90 | + | <source lang="fortran"> |
+ | module gravity | ||
+ | implicit none | ||
+ | real, parameter :: G = 6.67e-11 ! MKS units | ||
+ | |||
+ | contains | ||
+ | real function gravforce(x1,x2,m1,m2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1,x2 | ||
+ | real, intent(in) :: m1, m2 | ||
+ | real :: dist | ||
+ | |||
+ | dist = sqrt(sum((x1-x2)**2)) | ||
+ | gravforce = G * m1 * m2 / dist**2 | ||
+ | end function gravforce | ||
+ | end module gravity | ||
+ | |||
+ | program simplemod | ||
+ | use gravity | ||
+ | implicit none | ||
+ | |||
+ | print *, 'Gravitational constant = ', G | ||
+ | print *, 'Force between 2 1kg masses at [1,0,0] & | ||
+ | &and [0,0,1] is' | ||
+ | |||
+ | print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.) | ||
+ | |||
+ | end program simplemod | ||
+ | </source> | ||
+ | (from samples/procedures/simplemod/simplemod.f90) | ||
| | | | ||
|} | |} | ||
Line 334: | Line 695: | ||
==Compiling & Running== | ==Compiling & Running== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* When compiling the code a gravity.mod file is created | * When compiling the code a gravity.mod file is created | ||
Line 339: | Line 701: | ||
* Not compatible between different compilers, versions. | * Not compatible between different compilers, versions. | ||
| | | | ||
+ | <source lang="bash"> | ||
+ | $ ls | ||
+ | simplemod.f90 | ||
+ | |||
+ | $ gfortran -o simplemod simplemod.f90 -Wall | ||
+ | |||
+ | $ ls | ||
+ | gravity.mod simplemod simplemod.f90 | ||
+ | |||
+ | $ ./simplemod | ||
+ | Gravitational constant = 6.6700000E-11 | ||
+ | Force between 2 1kg masses at [1,0,0] and [0,0,1] is | ||
+ | 3.3350003E-11 | ||
+ | </source> | ||
|} | |} | ||
==Modules== | ==Modules== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* function gravforce can “see” the modulewide parameter defined above. | * function gravforce can “see” the modulewide parameter defined above. | ||
* So can main program, through use statement. | * So can main program, through use statement. | ||
| | | | ||
− | samples/procedures/simplemod/simplemod.f90 | + | <source lang="fortran"> |
+ | <source lang="fortran"> | ||
+ | module gravity | ||
+ | implicit none | ||
+ | real, parameter :: G = 6.67e-11 ! MKS units | ||
+ | |||
+ | contains | ||
+ | real function gravforce(x1,x2,m1,m2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1,x2 | ||
+ | real, intent(in) :: m1, m2 | ||
+ | real :: dist | ||
+ | |||
+ | dist = sqrt(sum((x1-x2)**2)) | ||
+ | gravforce = G * m1 * m2 / dist**2 | ||
+ | end function gravforce | ||
+ | end module gravity | ||
+ | |||
+ | program simplemod | ||
+ | use gravity | ||
+ | implicit none | ||
+ | |||
+ | print *, 'Gravitational constant = ', G | ||
+ | print *, 'Force between 2 1kg masses at [1,0,0] & | ||
+ | &and [0,0,1] is' | ||
+ | |||
+ | print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.) | ||
+ | |||
+ | end program simplemod | ||
+ | </source> | ||
+ | (from samples/procedures/simplemod/simplemod.f90) | ||
|} | |} | ||
==use module, only :== | ==use module, only :== | ||
{| | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Best practice is to only pull in from the module what you need | * Best practice is to only pull in from the module what you need | ||
* Otherwise, everything. | * Otherwise, everything. | ||
Line 357: | Line 766: | ||
* (Note syntax for continuation of a string...) | * (Note syntax for continuation of a string...) | ||
| | | | ||
+ | <source lang="fortran"> | ||
+ | module gravity | ||
+ | implicit none | ||
+ | real, parameter :: G = 6.67e-11 ! MKS units | ||
+ | |||
+ | contains | ||
+ | real function gravforce(x1,x2,m1,m2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1,x2 | ||
+ | real, intent(in) :: m1, m2 | ||
+ | real :: dist | ||
+ | |||
+ | dist = sqrt(sum((x1-x2)**2)) | ||
+ | gravforce = G * m1 * m2 / dist**2 | ||
+ | end function gravforce | ||
+ | end module gravity | ||
+ | |||
+ | program simplemod2 | ||
+ | use gravity, only : G, gravforce | ||
+ | implicit none | ||
+ | |||
+ | print *, 'Gravitational constant = ', G | ||
+ | print *, 'Force between 2 1kg masses at [1,0,0] & | ||
+ | &and [0,0,1] is' | ||
+ | |||
+ | print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.) | ||
+ | |||
+ | end program simplemod2 | ||
+ | </source> | ||
samples/procedures/simplemod/simplemod2.f90 | samples/procedures/simplemod/simplemod2.f90 | ||
|} | |} | ||
Line 362: | Line 800: | ||
==Modules usually get their own files== | ==Modules usually get their own files== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* For encapsulation | * For encapsulation | ||
Line 369: | Line 808: | ||
* (Main program hasn’t changed much). | * (Main program hasn’t changed much). | ||
| | | | ||
− | samples/procedures/multifilemod/gravity.f90 | + | <source lang="fortran"> |
+ | module gravity | ||
+ | implicit none | ||
+ | private | ||
+ | |||
+ | character (len=8), parameter, public :: massunit="kilogram" | ||
+ | character (len=8), parameter, public :: forceunit="Newton" | ||
+ | public :: gravforce | ||
+ | |||
+ | real, parameter :: G = 6.67e-11 ! MKS units | ||
+ | |||
+ | contains | ||
+ | real function distance(x1,x2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1, x2 | ||
+ | |||
+ | distance = sqrt(sum((x1-x2)**2)) | ||
+ | end function distance | ||
+ | |||
+ | real function gravforce(x1,x2,m1,m2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1,x2 | ||
+ | real, intent(in) :: m1, m2 | ||
+ | real :: dist | ||
+ | |||
+ | dist = distance(x1,x2) | ||
+ | gravforce = G * m1 * m2 / dist**2 | ||
+ | end function gravforce | ||
+ | end module gravity | ||
+ | </source> | ||
+ | (from samples/procedures/multifilemod/gravity.f90 ) | ||
|} | |} | ||
==Modules usually get their own files== | ==Modules usually get their own files== | ||
{| | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Compiling gravity.f90 now gives both an .o file (containing the code) and the .mod file as before. | * 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. | * Compiling the main program (multifilemod.f90) requires the .mod file. | ||
| | | | ||
− | samples/procedures/multifilemod/Makefile | + | <source lang="make"> |
+ | FC=gfortran | ||
+ | FFLAGS=-O3 -Wall | ||
+ | |||
+ | multifilemod: multifilemod.o gravity.o | ||
+ | $(FC) -o $@ multifilemod.o gravity.o | ||
+ | |||
+ | %.mod: %.f90 | ||
+ | $(FC) $(FFLAGS) -c $< | ||
+ | |||
+ | multifilemod.o: multifilemod.f90 gravity.mod | ||
+ | $(FC) $(FFLAGS) -c $< | ||
+ | |||
+ | clean: | ||
+ | rm -f *.o *~ *.mod multifilemod | ||
+ | </source> | ||
+ | (from samples/procedures/multifilemod/Makefile) | ||
|} | |} | ||
==.mod needed for compilation== | ==.mod needed for compilation== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
* ...because needs the type information of the constants, | * ...because needs the type information of the constants, | ||
Line 387: | Line 875: | ||
* Can’t compile without these | * Can’t compile without these | ||
| | | | ||
− | samples/procedures/multifilemod/multifilemod.f90 | + | <source lang="fortran"> |
+ | program simplemod2 | ||
+ | use gravity, only : gravforce, massunit, forceunit | ||
+ | implicit none | ||
+ | |||
+ | print *, 'Force between 2 1 ', massunit ,' masses ', & | ||
+ | ' at [1,0,0] and [0,0,1] is' | ||
+ | |||
+ | print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.), forceunit | ||
+ | |||
+ | end program simplemod2 | ||
+ | </source> | ||
+ | (from samples/procedures/multifilemod/multifilemod.f90) | ||
|} | |} | ||
==.o needed for linking== | ==.o needed for linking== | ||
{| | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Linking, however, doesn’t require the .mod file | * Linking, however, doesn’t require the .mod file | ||
* Only requires the .o file from the module code. | * Only requires the .o file from the module code. | ||
* .mod file analogous (but better than) .h files for C code. | * .mod file analogous (but better than) .h files for C code. | ||
| | | | ||
− | samples/procedures/multifilemod/Makefile | + | <source lang="make"> |
+ | FC=gfortran | ||
+ | FFLAGS=-O3 -Wall | ||
+ | |||
+ | multifilemod: multifilemod.o gravity.o | ||
+ | $(FC) -o $@ multifilemod.o gravity.o | ||
+ | |||
+ | %.mod: %.f90 | ||
+ | $(FC) $(FFLAGS) -c $< | ||
+ | |||
+ | multifilemod.o: multifilemod.f90 gravity.mod | ||
+ | $(FC) $(FFLAGS) -c $< | ||
+ | |||
+ | clean: | ||
+ | rm -f *.o *~ *.mod multifilemod | ||
+ | </source> | ||
+ | (from samples/procedures/multifilemod/Makefile) | ||
|} | |} | ||
==Compiling and running== | ==Compiling and running== | ||
{| | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* So compile files with modules first, so those that need them have the .mod files | * So compile files with modules first, so those that need them have the .mod files | ||
* Link the .o files | * Link the .o files | ||
+ | | | ||
+ | <source lang="bash"> | ||
+ | $ make | ||
+ | gfortran -O3 -Wall -c gravity.f90 | ||
+ | gfortran -O3 -Wall -c multifilemod.f90 | ||
+ | gfortran -o multifilemod multifilemod.o gravity.o | ||
+ | reposado-$ ./multifilemod | ||
+ | Force between 2 1 kilogram masses at [1,0,0] and [0,0,1] is | ||
+ | 3.33500033E-11 Newton | ||
+ | </source> | ||
|} | |} | ||
==Private and public== | ==Private and public== | ||
− | * Not all of a module’s | + | {| |
− | content need be public | + | |- valign="top" |
− | * 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 | + | <source lang="fortran"> |
− | + | module gravity | |
− | samples/procedures/multifilemod/gravity.f90 | + | implicit none |
+ | private | ||
+ | |||
+ | character (len=8), parameter, public :: massunit="kilogram" | ||
+ | character (len=8), parameter, public :: forceunit="Newton" | ||
+ | public :: gravforce | ||
+ | |||
+ | real, parameter :: G = 6.67e-11 ! MKS units | ||
+ | |||
+ | contains | ||
+ | real function distance(x1,x2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1, x2 | ||
+ | |||
+ | distance = sqrt(sum((x1-x2)**2)) | ||
+ | end function distance | ||
+ | |||
+ | real function gravforce(x1,x2,m1,m2) | ||
+ | implicit none | ||
+ | real, dimension(3), intent(in) :: x1,x2 | ||
+ | real, intent(in) :: m1, m2 | ||
+ | real :: dist | ||
+ | |||
+ | dist = distance(x1,x2) | ||
+ | gravforce = G * m1 * m2 / dist**2 | ||
+ | end function gravforce | ||
+ | end module gravity | ||
+ | </source> | ||
+ | ( from samples/procedures/multifilemod/gravity.f90 ) | ||
+ | |} | ||
==Procedures== | ==Procedures== | ||
− | * We’ve already seen | + | {| |
− | procedures defined in | + | |- valign="top" |
− | 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. | |
+ | | | ||
+ | <source lang="fortran"> | ||
+ | module procedures | ||
+ | |||
+ | contains | ||
+ | |||
+ | function square(x) result(xsquared) | ||
+ | implicit none | ||
+ | real :: xsquared | ||
+ | real, intent(IN) :: x | ||
+ | |||
+ | xsquared = x*x | ||
+ | end function square | ||
+ | |||
+ | function cube(x) | ||
+ | implicit none | ||
+ | real :: cube | ||
+ | real, intent(IN) :: x | ||
+ | |||
+ | cube = x*x*x | ||
+ | end function cube | ||
+ | |||
+ | subroutine squareAndCube(x, squarex, cubex) | ||
+ | implicit none | ||
+ | real, intent(in) :: x | ||
+ | real, intent(out) :: squarex | ||
+ | real, intent(out) :: cubex | ||
+ | |||
+ | squarex = square(x) | ||
+ | cubex = cube(x) | ||
+ | end subroutine squareAndCube | ||
− | + | end module procedures | |
− | + | </source> | |
− | + | ( from samples/procedures/funcsub/procedures.f90 ) | |
− | + | |} | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | ( | ||
− | |||
− | |||
− | |||
− | samples/procedures/funcsub/procedures.f90 | ||
==Functions== | ==Functions== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* 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) | + | <source lang="fortran"> |
− | + | function square(x) result(xsquared) | |
− | * Explicitly type the function | + | implicit none |
− | name, set that as return value | + | real :: xsquared |
− | real :: cube | + | real, intent(IN) :: x |
− | * | + | |
− | intent | + | xsquared = x*x |
− | samples/procedures/funcsub/procedures.f90 | + | end function square |
+ | |||
+ | function cube(x) | ||
+ | implicit none | ||
+ | real :: cube | ||
+ | real, intent(IN) :: x | ||
+ | |||
+ | cube = x*x*x | ||
+ | end function cube | ||
+ | |||
+ | subroutine squareAndCube(x, squarex, cubex) | ||
+ | implicit none | ||
+ | real, intent(in) :: x | ||
+ | real, intent(out) :: squarex | ||
+ | real, intent(out) :: cubex | ||
+ | |||
+ | squarex = square(x) | ||
+ | cubex = cube(x) | ||
+ | end subroutine squareAndCube | ||
+ | </source> | ||
+ | ( from samples/procedures/funcsub/procedures.f90 ) | ||
+ | |} | ||
==Procedure interfaces== | ==Procedure interfaces== | ||
− | * The interface to a procedure | + | {| |
− | consists of | + | |- valign="top" |
− | * 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 | + | <source lang="fortran"> |
− | interfaces to all public module | + | function square(x) result(xsquared) |
− | procedures. | + | implicit none |
+ | real :: xsquared | ||
+ | real, intent(IN) :: x | ||
+ | |||
+ | xsquared = x*x | ||
+ | end function square | ||
+ | |||
+ | function cube(x) | ||
+ | implicit none | ||
+ | real :: cube | ||
+ | real, intent(IN) :: x | ||
+ | |||
+ | cube = x*x*x | ||
+ | end function cube | ||
+ | |||
+ | subroutine squareAndCube(x, squarex, cubex) | ||
+ | implicit none | ||
+ | real, intent(in) :: x | ||
+ | real, intent(out) :: squarex | ||
+ | real, intent(out) :: cubex | ||
+ | |||
+ | squarex = square(x) | ||
+ | cubex = cube(x) | ||
+ | end subroutine squareAndCube | ||
+ | </source> | ||
+ | ( from samples/procedures/funcsub/procedures.f90 ) | ||
+ | |} | ||
==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 | + | <source lang="fortran"> |
− | you do with xlo, xhi, n. | + | function integratefx(xlo, xhi, f, n) |
− | * You do that for procedures | + | ! integrate with trapezoid rule |
− | with | + | ! .... |
− | http://en.wikipedia.org/wiki/ | + | |
− | File:Trapezoidal_rule_illustration_small.svg | + | integer :: i |
+ | real :: dx, xleft, xright | ||
+ | |||
+ | integratefx = 0. | ||
+ | dx = (xhi-xlo)/n | ||
+ | xleft = xlo | ||
+ | do i=0, n-1 | ||
+ | xright = xleft + dx | ||
+ | integratefx = integratefx + dx*(f(xright)+f(xleft))/2. | ||
+ | xleft = xright | ||
+ | enddo | ||
+ | end function integratefx | ||
+ | </source> | ||
+ | (from samples/procedures/interface/integrate.f90 ) | ||
+ | |||
+ | [[File:Wiki_trapezoidal.png|300px]] | ||
+ | |||
+ | (from 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. |
− | interface. | + | | |
− | * samples/procedures/interface/ | + | <source lang="fortran"> |
− | integrate.f90 | + | function integratefx(xlo, xhi, f, n) |
+ | ! integrate with trapezoid rule | ||
+ | implicit none | ||
+ | real, intent(in) :: xlo, xhi | ||
+ | interface | ||
+ | function f(x) | ||
+ | implicit none | ||
+ | real :: f | ||
+ | real, intent(in) :: x | ||
+ | end function f | ||
+ | end interface | ||
+ | integer, intent(in) :: n | ||
+ | real :: integratefx | ||
+ | |||
+ | integer :: i | ||
+ | real :: dx, xleft, xright | ||
+ | |||
+ | integratefx = 0. | ||
+ | dx = (xhi-xlo)/n | ||
+ | xleft = xlo | ||
+ | do i=0, n-1 | ||
+ | xright = xleft + dx | ||
+ | integratefx = integratefx + dx*(f(xright)+f(xleft))/2. | ||
+ | xleft = xright | ||
+ | enddo | ||
+ | end function integratefx | ||
+ | |||
+ | </source> | ||
+ | (from samples/procedures/interface/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 | + | | |
− | + | <source lang="fortran"> | |
+ | recursive function integratefx(xlo, xhi, f, tol) result(integral) | ||
+ | ! integrate with trapezoid rule, simpsons rule; | ||
+ | ! if difference between two is larger than | ||
+ | ! relevant tolerance, subdivide region. | ||
− | samples/procedures/recursive/integrate.f90 | + | ! ...variable declarations as before... |
+ | |||
+ | dx = xhi-xlo | ||
+ | xmid = (xlo+xhi)/2. | ||
+ | trapezoid = dx*(f(xlo)+f(xhi))/2. | ||
+ | simpsons = dx/6.*(f(xlo)+4.*f(xmid)+f(xhi)) | ||
+ | error = abs(trapezoid-simpsons)/(0.5*(trapezoid+simpsons)) | ||
+ | |||
+ | if (error > tol) then | ||
+ | ! too coarse; subdivide | ||
+ | integral = integratefx(xlo,xmid,f,tol) + & | ||
+ | integratefx(xmid,xhi,f,tol) | ||
+ | else | ||
+ | integral = trapezoid | ||
+ | endif | ||
+ | end function integratefx | ||
+ | </source> | ||
+ | ( from 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. | ||
+ | | | ||
+ | <source lang="fortran"> | ||
+ | pure subroutine axpy(a, x, y) | ||
+ | ! y = y + a*x | ||
+ | implicit none | ||
+ | real, intent(IN) :: a, x | ||
+ | real, intent(INOUT) :: y | ||
− | == | + | y = y + a*x |
− | + | end subroutine axpy | |
− | + | ||
− | + | subroutine printaxpy(a, x, y) | |
− | * | + | ! y = y + a*x |
− | + | implicit none | |
− | * | + | real, intent(IN) :: a, x |
− | + | real, intent(INOUT) :: y | |
− | samples/procedures/purity/purity.f90 | + | |
+ | print *, a, '*', x, ' + ', y, & | ||
+ | ' = ', a*x+y | ||
+ | y = a*x + y | ||
+ | end subroutine printaxpy | ||
+ | </source> | ||
+ | (from 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 | + | * Often useful to clarify calling of procedure if many parameters have sensible defaults. |
− | present; have to use | + | * Avoid code duplication, wrappers of having one version of routine with default values, one with user-supplied |
− | another variable. | + | * Can’t use tol if not present; have to use another variable. |
− | samples/procedures/optional/integrate.f90 | + | | |
+ | <source lang="fortran"> | ||
+ | recursive function integratefx(xlo, xhi, f, tol) result(integral) | ||
+ | ! .... | ||
+ | real, intent(in), optional :: tol | ||
+ | ! .... | ||
+ | |||
+ | real :: errtol | ||
+ | |||
+ | ! use parameter if passed, | ||
+ | ! else use default | ||
+ | if (present(tol)) then | ||
+ | errtol = tol | ||
+ | else | ||
+ | errtol = 1.e-6 | ||
+ | endif | ||
+ | |||
+ | ! .... | ||
+ | |||
+ | if (error > errtol) then | ||
+ | ! too coarse; subdivide | ||
+ | integral = integratefx(xlo,xmid,f,errtol) + & | ||
+ | integratefx(xmid,xhi,f,errtol) | ||
+ | else | ||
+ | integral = trapezoid | ||
+ | endif | ||
+ | end function integratefx | ||
+ | |||
+ | </source> | ||
+ | (from samples/procedures/optional/integrate.f90) | ||
+ | |} | ||
==Optional Arguments== | ==Optional Arguments== | ||
− | * When calling the | + | {| |
− | procedure, can use | + | |- valign="top" |
− | the optional | + | | |
− | argument or not. | + | * When calling the procedure, can use the optional argument or not. |
− | * Makes sense to leave | + | * Makes sense to leave optional arguments at end - easier to figure out what’s what when it’s omitted. |
− | optional arguments at | + | |
− | end - easier to figure | + | | |
− | out what’s what | + | <source lang="fortran"> |
− | when it’s omitted. | + | print *, 'Integrating using default tol' |
+ | approx = integratefx(0., 2*pi, sinesquared) | ||
+ | print *, 'Approximate integral = ', approx | ||
+ | print *, 'Exact integral = ', exact | ||
− | samples/procedures/optional/optional.f90 | + | print *, '' |
+ | print *, 'Integrating using coarser tol' | ||
+ | approx = integratefx(0., 2*pi, sinesquared, 0.01) | ||
+ | print *, 'Approximate integral = ', approx | ||
+ | |||
+ | ! .... | ||
+ | </source> | ||
+ | (from samples/procedures/optional/optional.f90) | ||
+ | |} | ||
==Keyword Arguments== | ==Keyword Arguments== | ||
− | * To avoid ambiguity with | + | {| |
− | omitted arguments - or | + | |- valign="top" |
− | really whenever you | + | | |
− | want - you can specify | + | * To avoid ambiguity with omitted arguments - or really whenever you want - you can specify which value is which explicitly. |
− | which value is which | + | * Don’t have to be in order. |
− | explicitly. | + | * Can clarify calls of routines with many arguments |
− | * Don’t have to be in | + | | |
− | order. | + | <source lang="fortran"> |
− | * Can clarify calls of | + | ! .... |
− | routines with many | + | print *, 'Integrating using still coarser tol' |
− | + | approx = integratefx(xhi=2*pi, xlo=0., tol=0.5, & | |
+ | f=sinesquared) | ||
+ | print *, 'Approximate integral = ', approx | ||
+ | </source> | ||
+ | (from samples/procedures/optional/optional.f90) | ||
+ | |} | ||
− | |||
− | ==Procedures | + | ==Procedures and Modules Summary== |
− | + | {| | |
− | * Modules let you bundle procedures, constants | + | |- valign="top" |
− | in useful packages. | + | | |
+ | * Modules let you bundle procedures, constants 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= | |
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Fortran made for dealing with scientific data | * Fortran made for dealing with scientific data | ||
* Arrays built into language | * Arrays built into language | ||
− | * The type information associated with an | + | * The type information associated with an array includes rank (# of dimension), size, element type, stride.. |
− | array includes rank (# of dimension), size, | + | * Enables powerful optimizations, programmer-friendly features. |
− | element type, stride.. | + | * Can be manipulated like simple scalar variables |
− | * Enables powerful optimizations, | + | * Elementwise addition, multiplication.. |
− | programmer-friendly features. | + | | |
+ | <source lang="fortran"> | ||
+ | program basicarrays | ||
+ | implicit none | ||
+ | integer, dimension(5) :: a, b, c | ||
+ | integer :: i | ||
+ | |||
+ | a = [1,2,3,4,5] | ||
+ | b = [(2*i+1, i=1,5)] | ||
+ | |||
+ | |||
+ | print *, 'a = ', a | ||
+ | print *, 'b = ', b | ||
+ | |||
+ | c = a+b | ||
+ | print *, 'c = ', c | ||
− | = | + | c = a*b + 1 |
− | * | + | print *, 'a*b+1=', c |
− | + | end program basicarrays | |
− | + | </source> | |
− | * | + | (from samples/arrays/basic.f90 ) |
− | + | |} | |
− | |||
− | samples/arrays/basic.f90 | ||
==Array constructors== | ==Array constructors== | ||
− | * Can have array | + | {| |
− | constants like | + | |- valign="top" |
− | numerical constants | + | | |
− | * use [] or (/ /), then | + | * Can have array constants like numerical constants |
− | comma-separated | + | * use [] or (/ /), then comma-separated list of values. |
− | list of values. | + | * Implied do loops can be used in constructors |
− | * Implied do loops | + | * (Implied do-loop variables have to be defined) |
− | can be used in | + | | |
− | constructors | + | <source lang="fortran"> |
− | * ( | + | x = [1,2,3,4,5] |
− | be defined) | + | x = (/1,2,3,4,5/) |
+ | x = [ (i,i=1,5)] | ||
+ | a = [ ((i*j,j=1,3),i=1,5)] | ||
+ | </source> | ||
+ | |} | ||
− | + | ==Elementwise operations== | |
− | + | {| | |
− | + | |- valign="top" | |
+ | | | ||
+ | * 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. | ||
+ | | | ||
+ | <source lang="fortran"> | ||
+ | program elementwise | ||
+ | implicit none | ||
+ | real, dimension(10) :: x,y,z | ||
+ | integer :: i | ||
+ | real, parameter:: pi = 4.*atan(1.) | ||
− | = | + | x = [(2*pi*(i-1)/9.,i=1,10)] |
− | * | + | |
− | + | y = sin(x) | |
− | + | z = x*x | |
− | |||
− | * | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | samples/arrays/elementwise.f90 | + | print *, x |
+ | print *, y | ||
+ | print *, z | ||
+ | end program elementwise | ||
+ | </source> | ||
+ | (from samples/arrays/elementwise.f90 ) | ||
+ | |} | ||
==Elemental Functions== | ==Elemental Functions== | ||
− | * User can create their | + | {| |
− | own elemental functions | + | |- valign="top" |
− | * Label any scalar function | + | | |
− | with “elemental” should (until recently, | + | * User can create their own elemental functions |
− | must) be pure, so can be | + | * Label any scalar function with “elemental” should (until recently, must) be pure, so can be applied everywhere at same time. |
− | applied everywhere at | + | * Can be faster than loop. |
− | same time. | + | * Can also take multiple arguments: eg z = addsquare(x,y) |
− | * | + | | |
− | * Can also take multiple | + | <source lang="fortran"> |
− | arguments: eg | + | program elementalfn |
− | z = | + | implicit none |
+ | real, dimension(10) :: x,y,z | ||
+ | integer :: i | ||
+ | real, parameter:: pi = 4.*atan(1.) | ||
+ | |||
+ | x = [(2*pi*(i-1)/9.,i=1,10)] | ||
+ | |||
+ | y = sinesquared(x) | ||
+ | z = sin(x)*sin(x) | ||
− | samples/arrays/elemental.f90 | + | print *, x |
+ | print *, y | ||
+ | print *, z | ||
+ | print *,z(::3) | ||
+ | contains | ||
+ | elemental function sinesquared(x) | ||
+ | implicit none | ||
+ | real :: sinesquared | ||
+ | real, intent(in) :: x | ||
+ | |||
+ | sinesquared = sin(x)**2 | ||
+ | end function sinesquared | ||
+ | end program elementalfn | ||
+ | </source> | ||
+ | (from samples/arrays/elemental.f90 ) | ||
+ | |} | ||
==Array comparisons== | ==Array comparisons== | ||
− | * Array comparisons | + | {| |
− | return an array of | + | |- valign="top" |
− | logicals of the same size | + | | |
− | of the arrays. | + | * Array comparisons return an array of logicals of the same size of the arrays. |
− | * Can use any and all to | + | * Can use any and all to see if any or all of those logicals are true. |
− | see if any or all of those | + | | |
− | logicals are true. | + | <source lang="fortran"> |
− | samples/arrays/compare.f90 | + | program comparearrays |
+ | implicit none | ||
+ | integer, dimension(5) :: a, b | ||
+ | integer :: i | ||
+ | |||
+ | a = [1,2,3,4,5] | ||
+ | b = [(2*i-3, i=1,5)] | ||
+ | |||
+ | print *, 'A = ', a | ||
+ | print *, 'B = ', b | ||
+ | |||
+ | if (any(a > b)) then | ||
+ | print *, 'An A is larger than a B' | ||
+ | endif | ||
+ | if (all(a > b)) then | ||
+ | print *, 'All As ares larger than Bs' | ||
+ | else if (all(b > a)) then | ||
+ | print *, 'All Bs are larger than As' | ||
+ | else | ||
+ | print *, 'A, B values overlap' | ||
+ | endif | ||
+ | |||
+ | end program comparearrays | ||
+ | </source> | ||
+ | (from samples/arrays/compare.f90) | ||
+ | |} | ||
==Array masks== | ==Array masks== | ||
− | * These logical arrays can | + | {| |
− | be used to mask several | + | |- valign="top" |
− | operations | + | | |
− | * Only do sums, mins, etc | + | * These logical arrays can be used to mask several operations |
− | where the mask is true | + | * Only do sums, mins, etc where the mask is true |
− | * Eg, only pick out positive | + | * Eg, only pick out positive values. |
− | values. | + | * Many array intrinsics have this mask option |
− | * Many array intrinsics | + | | |
− | have this mask option | + | <source lang="fortran"> |
+ | program mask | ||
+ | implicit none | ||
+ | integer, dimension(10) :: a | ||
+ | logical, dimension(10) :: pos | ||
+ | integer :: i | ||
− | samples/arrays/mask.f90 | + | a = [(2*i-7, i=1,10)] |
+ | pos = (a > 0) | ||
+ | |||
+ | print '(A,10(I4,1X))','A = ', a | ||
+ | print *,'# of positive values: ', count(pos) | ||
+ | print *,'Sum of positive values: ', sum(a,pos) | ||
+ | print *,'Minimum positive value: ', minval(a,pos) | ||
+ | |||
+ | end program mask | ||
+ | </source> | ||
+ | ( from samples/arrays/mask.f90 ) | ||
+ | <source lang="bash"> | ||
+ | $ ./mask | ||
+ | A = -5 -3 -1 1 3 5 7 9 11 13 | ||
+ | # of positive values: 7 | ||
+ | Sum of positive values: 49 | ||
+ | Minimum positive value: 1 | ||
+ | </source> | ||
+ | |} | ||
==Where construct== | ==Where construct== | ||
− | * The where construct | + | {| |
− | can be used to easily | + | |- valign="top" |
− | manipulate sections of | + | | |
− | array based on arbitrary | + | * The where construct can be used to easily manipulate sections of array based on arbitrary comparisons. |
− | comparisons. | + | * Where construct => for whatever indices the comparison is true, set values as follow; otherwise, set other values. |
− | * Where construct => for | + | | |
− | whatever indices the | + | <source lang="fortran"> |
− | comparison is true, set | + | program wherearray |
− | values as follow; | + | implicit none |
− | otherwise, set other | + | real, dimension(6) :: a, diva |
− | + | integer :: i | |
+ | |||
+ | a = [(2*i-6, i=1,6)] | ||
+ | where (a /= 0) | ||
+ | diva = 1/a | ||
+ | elsewhere | ||
+ | diva = -999 | ||
+ | endwhere | ||
+ | |||
+ | print *,'a = ' | ||
+ | print '(6(F8.3,1X))',a | ||
+ | print *,'1/a = ' | ||
+ | print '(6(F8.3,1X))',diva | ||
− | samples/arrays/where.f90 | + | end program wherearray |
+ | </source> | ||
+ | (from samples/arrays/where.f90) | ||
+ | <source lang="bash"> | ||
+ | $ ./where | ||
+ | a = | ||
+ | -4.000 -2.000 0.000 2.000 4.000 6.000 | ||
+ | 1/a = | ||
+ | -0.250 -0.500 -999.000 0.500 0.250 0.167 | ||
+ | </source> | ||
+ | |} | ||
==Forall construct== | ==Forall construct== | ||
− | * Forall is an array | + | {| |
− | assignment statement | + | |- valign="top" |
− | * Each line in forall has to be | + | | |
− | independent. All done “at | + | * Forall is an array assignment statement |
− | once” - no guarantees as | + | * Each line in forall has to be independent. All done “at once” - no guarantees as to order |
− | to order | + | * If (say) 2 lines in the forall, all of the first line is done, then all of the second. |
− | * If (say) 2 lines in the forall, | + | * Any functions called must be pure |
− | all of the first line is done, | + | * Can be vectorized or parallelized by compiler |
− | then all of the second. | + | | |
− | * Any functions called must | + | <source lang="fortran"> |
− | be pure | + | program forallarray |
− | * Can be vectorized or | + | implicit none |
− | parallelized by compiler | + | integer, dimension(6,6) :: a |
+ | integer :: i,j | ||
− | samples/arrays/forall.f90 | + | a = -999 |
+ | forall (i=1:6, j=1:6, i/=j) | ||
+ | a(i,j) = i-j | ||
+ | endforall | ||
+ | |||
+ | do i=1,6 | ||
+ | print '(6(I5,1X))',(a(i,j),j=1,6) | ||
+ | enddo | ||
+ | end program forallarray | ||
+ | </source> | ||
+ | (from samples/arrays/forall.f90) | ||
+ | <source lang="bash"> | ||
+ | $ ./forall | ||
+ | -999 -1 -2 -3 -4 -5 | ||
+ | 1 -999 -1 -2 -3 -4 | ||
+ | 2 1 -999 -1 -2 -3 | ||
+ | 3 2 1 -999 -1 -2 | ||
+ | 4 3 2 1 -999 -1 | ||
+ | 5 4 3 2 1 -999 | ||
+ | </source> | ||
+ | |} | ||
==Array Sections== | ==Array Sections== | ||
− | * Generalization of array | + | {| |
− | indexing | + | |- valign="top" |
− | * Familiar to users of | + | | |
− | Matlab, IDL, Python.. | + | * Generalization of array indexing |
− | * Can use “slices” of an | + | * Familiar to users of Matlab, IDL, Python.. |
− | array using “index | + | * Can use “slices” of an array using “index triplet” |
− | triplet” | ||
* [start]:[end][:step] | * [start]:[end][:step] | ||
− | * Default start=1, default | + | * Default start=1, default end=size, default step=1. |
− | end=size, default step=1. | + | * Can be used for each index of multid array |
− | * Can be used for each | + | | |
− | index of multid array | + | <source lang="fortran"> |
a([start]:[end][:step]) | a([start]:[end][:step]) | ||
+ | |||
a = [1,2,3,4,5,6,7,8,9,10] | a = [1,2,3,4,5,6,7,8,9,10] | ||
+ | |||
a(7:) == [7,8,9,10] | a(7:) == [7,8,9,10] | ||
a(:3) == [1,2,3] | a(:3) == [1,2,3] | ||
Line 754: | Line 1,647: | ||
a(2) == 2 | a(2) == 2 | ||
a(:) == [1,2,3,4,5,6,7,8,9,10] | a(:) == [1,2,3,4,5,6,7,8,9,10] | ||
+ | </source> | ||
+ | |} | ||
==Array Sections== | ==Array Sections== | ||
− | * This sort of thing is very | + | {| |
− | handy in numerical | + | |- valign="top" |
− | computation | + | | |
− | * Replace do-loops with | + | * This sort of thing is very handy in numerical computation |
− | clearer, shorter, possibly | + | * Replace do-loops with clearer, shorter, possibly vectorized array operations |
− | vectorized array | + | * Bigger advantage for multidimensional arrays |
− | operations | + | | |
− | * Bigger advantage for | + | <source lang="fortran"> |
− | multidimensional arrays | + | program derivative |
+ | implicit none | ||
+ | real, dimension(10) :: x | ||
+ | real, dimension(9) :: derivx | ||
+ | integer :: i | ||
+ | real, parameter:: pi = 4.*atan(1.), h=1. | ||
+ | |||
+ | x = [(2*pi*(i-1)/9.,i=1,10)] | ||
+ | |||
+ | derivx = ((x(2:10)-x(1:9))/h) | ||
+ | print *, derivx | ||
− | samples/arrays/derivative.f90 | + | do i=1,9 |
+ | derivx(i) = (x(i+1)-x(i))/h | ||
+ | enddo | ||
+ | print *, derivx | ||
+ | end program derivative | ||
+ | </source> | ||
+ | (from samples/arrays/derivative.f90) | ||
+ | |} | ||
==Array Sections== | ==Array Sections== | ||
− | * The previous sorts of array | + | {| |
− | sections - shifting things leftward | + | |- valign="top" |
− | and rightward - are so common | + | | |
− | there are intrinsics for them | + | * The previous sorts of array sections - shifting things leftward and rightward - are so common there are intrinsics for them |
− | * | + | * positive shift shifts elements leftwards (or array bounds rightwards). |
− | leftwards (or array bounds | + | * cshift does circular shift - shifting off the end of the array “wraps around”. |
− | rightwards). | + | * eoshift fills with zeros, or optional filling. |
− | * cshift does circular shift - shifting | ||
− | off the end of the array “wraps | ||
− | around”. | ||
− | * eoshift | ||
− | optional filling. | ||
* Can work on given dimension | * Can work on given dimension | ||
+ | | | ||
+ | <source lang="fortran"> | ||
+ | |||
a = [1,2,3,4,5] | a = [1,2,3,4,5] | ||
+ | |||
cshift(a,1) == [2,3,4,5,1] | cshift(a,1) == [2,3,4,5,1] | ||
cshift(a,-1) == [5,1,2,3,4] | cshift(a,-1) == [5,1,2,3,4] | ||
eoshift(a,1) ==[2,3,4,5,0] | eoshift(a,1) ==[2,3,4,5,0] | ||
eoshift(a,-1)==[0,1,2,3,4] | eoshift(a,-1)==[0,1,2,3,4] | ||
+ | </source> | ||
+ | |} | ||
− | ==Other important array== | + | ==Other important array intrinsics== |
− | + | {| | |
+ | |- valign="top" | ||
+ | | | ||
* minval/maxval - finds min, max element in an array. | * minval/maxval - finds min, max element in an array. | ||
* minloc/maxloc - finds location of min/max element | * minloc/maxloc - finds location of min/max element | ||
* product/sum - returns product/sum of array elements | * product/sum - returns product/sum of array elements | ||
* reshape - Adjusts shape of array data. Eg: | * reshape - Adjusts shape of array data. Eg: | ||
− | 1,4 | + | |
+ | <source lang="fortran"> | ||
+ | |||
+ | |||
+ | 1,4 | ||
reshape([1,2,3,4,5,6],[3,2]) == 2,5 | reshape([1,2,3,4,5,6],[3,2]) == 2,5 | ||
− | 3,6 | + | 3,6 |
+ | </source> | ||
+ | |} | ||
==Linear algebra in Fortran== | ==Linear algebra in Fortran== | ||
− | * Comes built in with transpose, matmul, | + | {| |
− | dot_product for dealing with arrays. | + | |- valign="top" |
+ | | | ||
+ | * Comes built in with transpose, matmul, dot_product for dealing with arrays. | ||
* matmul also does matrix-vector multiplication | * matmul also does matrix-vector multiplication | ||
− | * Either use these or system-provided BLAS | + | * Either use these or system-provided BLAS libraries - never, ever write yourself. |
− | libraries - never, ever write yourself. | + | * Consider the following timings for a matrix multiplication - a naive triple loop, a clal to that matmul intrinsic, or an explicit call to a BLAS package. |
+ | * factor of 100 difference! | ||
+ | * Note that you can build gfortran to use fast (BLAS) routines for intrinsics... | ||
+ | | | ||
+ | |} | ||
+ | |||
+ | ==Matrix Multiplication Times== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | <source lang="fortran"> | ||
+ | !... | ||
+ | |||
+ | print *, 'Experiment with matrix size ', n | ||
+ | print *, 'Times in seconds.' | ||
+ | |||
+ | allocate(a(n,n)) | ||
+ | allocate(b(n,n)) | ||
+ | allocate(c(n,n)) | ||
+ | call random_number(a) | ||
+ | call random_number(b) | ||
+ | |||
+ | call tick(starttime) | ||
+ | do j=1,n | ||
+ | do i=1,n | ||
+ | c(i,j) = 0. | ||
+ | do k=1,n | ||
+ | c(i,j) = c(i,j)+a(i,k)*b(k,j) | ||
+ | enddo | ||
+ | enddo | ||
+ | enddo | ||
+ | looptime = tock(starttime) | ||
+ | |||
+ | call tick(starttime) | ||
+ | c = matmul(a,b) | ||
+ | matmultime = tock(starttime) | ||
+ | |||
+ | call tick(starttime) | ||
+ | call sgemm('N','N',n,n,n,1.,a,n,b,n,0.,c,n) | ||
+ | sgemmtime = tock(starttime) | ||
− | = | + | |
+ | print *, 'Triple-loop time: ', looptime | ||
+ | print *, 'matmul intrinsic time: ', matmultime | ||
+ | print *, 'SGEMM lapack call time:', sgemmtime | ||
+ | !... | ||
+ | </source> | ||
+ | | | ||
+ | <source lang="bash"> | ||
$ ./matmul 2500 | $ ./matmul 2500 | ||
Experiment with matrix size | Experiment with matrix size | ||
− | 2500 | + | 2500 |
Triple-loop time: | Triple-loop time: | ||
− | 149.63400 | + | 149.63400 |
matmul intrinsic time: | matmul intrinsic time: | ||
− | 10.370000 | + | 10.370000 |
SGEMM lapack call time: | SGEMM lapack call time: | ||
− | 1.4809999 | + | 1.4809999 |
+ | </source> | ||
+ | (gfortran 4.6, compiled -O3 -march=native using Intel MKL 10.3 for sgemm) | ||
+ | (program from samples/arrays/matmul.f90) | ||
+ | |} | ||
− | ( | + | ==Linear algebra in Fortran== |
− | + | {| | |
+ | |- valign="top" | ||
+ | | | ||
+ | * Things like transposes work, too: | ||
+ | * see samples/arrays/matrix.f90 | ||
+ | <source lang="fortran"> | ||
+ | program matvec | ||
+ | implicit none | ||
+ | integer, dimension(4,5) :: a | ||
+ | integer, dimension(5,4) :: at | ||
+ | integer, dimension(4,4) :: aat | ||
+ | integer :: i | ||
− | + | a = reshape([(i,i=1,4*5)],[4,5]) | |
+ | at = transpose(a) | ||
+ | print *,'A = ' | ||
+ | call printmat(a) | ||
+ | print *,'A^T = ' | ||
+ | call printmat(at) | ||
− | == | + | aat = matmul(a,at) |
− | + | print *,'A . A^T = ' | |
− | + | call printmat(aat) | |
+ | !... | ||
+ | </source> | ||
+ | | | ||
+ | <source lang="bash"> | ||
+ | $ ./matrix | ||
+ | A = | ||
+ | 1 5 9 13 17 | ||
+ | 2 6 10 14 18 | ||
+ | 3 7 11 15 19 | ||
+ | 4 8 12 16 20 | ||
+ | A^T = | ||
+ | 1 2 3 4 | ||
+ | 5 6 7 8 | ||
+ | 9 10 11 12 | ||
+ | 13 14 15 16 | ||
+ | 17 18 19 20 | ||
+ | A . A^T = | ||
+ | 565 610 655 700 | ||
+ | 610 660 710 760 | ||
+ | 655 710 765 820 | ||
+ | 700 760 820 880 | ||
+ | A . A^T subarray = | ||
+ | 610 660 710 | ||
+ | 655 710 765 | ||
+ | 700 760 820 | ||
+ | </source> | ||
+ | |} | ||
==Array sizes and Assumed Shape== | ==Array sizes and Assumed Shape== | ||
− | * Printmat routine here is | + | {| |
− | interesting - don’t pass | + | |- valign="top" |
− | (a,rows,cols), just a. | + | | |
− | * Can assume a rank-2 array, | + | * Printmat routine here is interesting - don’t pass (a,rows,cols), just a. |
− | and get size at runtime. | + | * Can assume a rank-2 array, and get size at runtime. |
− | * Simplifies call, and eliminates | + | * Simplifies call, and eliminates possible inconsistency: what if rows, cols is wrong? |
− | possible inconsistency: what if | + | * size(array,dim) gets the size of array in the dim dimension. |
− | rows, cols is wrong? | + | * Assumed shape arrays (eg, dimension(:,:)) much better than older ways of passing arrays: |
− | * size(array,dim) gets the size of | ||
− | array in the dim dimension. | ||
− | + | integer nx, ny | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
integer a(nx,ny) | integer a(nx,ny) | ||
+ | |||
or worse, | or worse, | ||
+ | |||
integer a(*,ny) | integer a(*,ny) | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | samples/arrays/matrix.f90 | + | * Information is thrown away, possibility of inconsistency. |
+ | * Here, (:,:) means we know the rank, but don’t know the size yet. | ||
+ | | | ||
+ | <source lang="fortran"> | ||
+ | subroutine printmat(a) | ||
+ | implicit none | ||
+ | integer, dimension(:,:) :: a | ||
+ | integer :: nr, nc, i, j | ||
+ | |||
+ | nr = size(a,1) | ||
+ | nc = size(a,2) | ||
+ | do i=1,nr | ||
+ | print '(99(I4,1X))', (a(i,j), j=1,nc) | ||
+ | enddo | ||
+ | end subroutine printmat | ||
+ | </source> | ||
+ | (from samples/arrays/matrix.f90) | ||
+ | |} | ||
+ | |||
==Allocatable Arrays== | ==Allocatable Arrays== | ||
− | * So far, all our programs have had fixed-size | + | {| |
− | arrays, set at compile time. | + | |- valign="top" |
− | * To change problem size, have to edit code, | + | | |
− | recompile. | + | * So far, all our programs have had fixed-size arrays, set at compile time. |
− | * Has some advantages (optimization, | + | * To change problem size, have to edit code, recompile. |
− | determinism) but very inflexible. | + | * Has some advantages (optimization, determinism) but very inflexible. |
− | * Would like to be able to request memory | + | * Would like to be able to request memory at run time, make array of desired size. |
− | at run time, make array of desired size. | + | * Allocatable arrays are arguably most important addition to Fortran. |
− | * Allocatable arrays are arguably most | + | |} |
− | important addition to Fortran. | ||
==Allocate(), Deallocate()== | ==Allocate(), Deallocate()== | ||
− | * Give array a deferred size (eg, | + | {| |
− | dimension(:)) and the attribute | + | |- valign="top" |
− | allocatable. | + | | |
− | * When time to allocate it, use | + | * Give array a deferred size (eg, dimension(:)) and the attribute allocatable. |
− | allocate(a(n)). | + | * When time to allocate it, use allocate(a(n)). |
* Deallocate with deallocate(a). | * Deallocate with deallocate(a). | ||
− | + | * In between, arrays can be used as any other array. | |
− | * In between, arrays can be used | + | * If allocation fails (not enough memory available for request), program will exit. |
− | as any other array. | + | * 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. | ||
+ | | | ||
+ | <source lang="fortran"> | ||
+ | program allocarray | ||
+ | implicit none | ||
+ | integer :: i, n | ||
+ | integer, dimension(:), allocatable :: a | ||
+ | |||
+ | n = 10 | ||
+ | allocate(a(n)) | ||
+ | a = [(i, i=2,20,2)] | ||
+ | |||
+ | print *,'A = ' | ||
+ | print *,a | ||
− | + | deallocate(a) | |
− | + | end program allocarray | |
− | + | </source> | |
− | + | (from samples/arrays/allocatable.f90 ) | |
− | + | |} | |
− | |||
− | |||
− | |||
− | |||
− | |||
==get_command_argument()== | ==get_command_argument()== | ||
− | * Previous version still | + | {| |
− | depended on a compiled-in | + | |- valign="top" |
− | number. | + | | |
− | * Can read from file or from | + | * Previous version still depended on a compiled-in number. |
− | console, but Fortran now has | + | * Can read from file or from console, but Fortran now has standard way to get command-line arguments |
− | standard way to get | + | * Get the count of arguments, and if there’s at least one argument there, get it, read it as integer, and allocate array. |
− | command-line arguments | + | | |
− | * Get the count of arguments, | + | <source lang="fortran"> |
− | and if there’s at least one | + | program allocarray2 |
− | argument there, get it, read it | + | implicit none |
− | as integer, and allocate array. | + | integer :: i, n |
− | + | integer, dimension(:), allocatable :: a | |
+ | character(len=30) :: arg | ||
− | == | + | if (command_argument_count() < 1) then |
+ | print *,'Use: allocatable N, '//& | ||
+ | ' where N is array size.' | ||
+ | stop | ||
+ | endif | ||
+ | |||
+ | call get_command_argument(1, arg) | ||
+ | read( arg,'(I30)'), n | ||
+ | |||
+ | print *,'Allocating array of size ', n | ||
+ | allocate(a(n)) | ||
+ | |||
+ | a = [(i,i=1,n)] | ||
+ | print *, a | ||
+ | |||
+ | deallocate(a) | ||
+ | end program allocarray2 | ||
+ | </source> | ||
+ | (from samples/arrays/allocatable2.f90) | ||
+ | <source lang="bash"> | ||
+ | $ ./allocatable2 | ||
+ | Use: allocatable N, where N is array size. | ||
+ | |||
+ | $ ./allocatable2 3 | ||
+ | Allocating array of size 3 | ||
+ | 1 2 3 | ||
− | + | $ ./allocatable2 5 | |
+ | Allocating array of size 5 | ||
+ | 1 2 3 4 5 | ||
+ | </source> | ||
+ | |} | ||
==Hands on #3== | ==Hands on #3== | ||
− | * Use array functionality to simplify hydro code | + | {| |
− | -- don't need to pass, array size, and can | + | |- valign="top" |
− | simplify mathematics using array operations. | + | | |
− | * In workedexamples/arrays, have modified | + | * Use array functionality to simplify hydro code -- don't need to pass, array size, and can simplify mathematics using array operations. |
− | hydro to allocate u, and pbm to just take array. | + | * In workedexamples/arrays, have modified hydro to allocate u, and pbm to just take array. |
− | * Do the same with the fluid dynamic routines in | + | * Do the same with the fluid dynamic routines in solver.f90 |
− | solver.f90 | ||
* ~30 min | * ~30 min | ||
+ | | | ||
+ | |} | ||
− | + | =Fortran Pointers= | |
− | * Pointers, or references, | + | {| |
− | refer to another | + | |- valign="top" |
− | variable. | + | | |
− | * Eg, p does not contain a | + | * Pointers, or references, refer to another variable. |
− | real value, but a | + | * Eg, p does not contain a real value, but a reference to another real variable. |
− | reference to another | + | * Once associated with another variable, can read/write to it as if it were stored “in” p. |
− | real variable. | + | | |
− | * Once associated with | + | <source lang="fortran"> |
− | another variable, can | ||
− | read/write to it as if it | ||
− | were stored “in” p. | ||
− | |||
real, target :: x = 3.2 | real, target :: x = 3.2 | ||
real, pointer:: p | real, pointer:: p | ||
p => x | p => x | ||
+ | |||
p | p | ||
x | x | ||
3.2 | 3.2 | ||
− | + | </source> | |
− | |||
− | |||
samples/pointers/ptr1.f90 | samples/pointers/ptr1.f90 | ||
+ | |} | ||
==Fortran Pointers== | ==Fortran Pointers== | ||
− | * Pointers are either | + | {| |
− | associated, null, or | + | |- valign="top" |
− | undefined; start out life | + | | |
− | undefined. | + | * Pointers are either associated, null, or undefined; start out life undefined. |
− | * Can associate them to | + | * Can associate them to a variable with => , or mark them as not associated with any valid variable by pointing it to null(). |
− | a variable with => , or | + | * Reading value from or writing value to a null pointer will cause errors, probably crash. |
− | mark them as not | + | * Fortran pointers can’t point just anywhere. |
− | associated with any | + | * Must reference a variable with the same type, that has the target attribute. |
− | valid variable by | + | | |
− | pointing it to null(). | + | <source lang="fortran"> |
− | |||
real, target :: x = 3.2 | real, target :: x = 3.2 | ||
real, pointer:: p | real, pointer:: p | ||
Line 962: | Line 2,020: | ||
x | x | ||
+ | </source> | ||
+ | | | ||
+ | |} | ||
==Fortran Pointers== | ==Fortran Pointers== | ||
− | + | {| | |
− | + | |- valign="top" | |
− | * | + | | |
− | + | * Pointers can reference other pointers. | |
− | + | * Actually references what they’re pointing to. | |
− | + | | | |
− | + | <source lang="fortran"> | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | * | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
real, target :: x = 3.2 | real, target :: x = 3.2 | ||
real, pointer:: p1, p2 | real, pointer:: p1, p2 | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
p1 => x | p1 => x | ||
p2 => p1 | p2 => p1 | ||
Line 1,008: | Line 2,041: | ||
x | x | ||
3.2 | 3.2 | ||
+ | </source> | ||
+ | |} | ||
==Allocating a pointer== | ==Allocating a pointer== | ||
− | * Pointer doesn’t | + | {| |
− | necessarily have to have | + | |- valign="top" |
− | another variable to | + | | |
− | target | + | * Pointer doesn’t necessarily have to have another variable to target |
− | * Can allocate memory | + | * Can allocate memory for p to point to that does not belong to any other pointer. |
− | for p to point to that | + | * Must deallocate it when done |
− | does not belong to any | + | | |
− | other pointer. | + | <source lang="fortran"> |
− | * Must deallocate it when | ||
− | |||
− | |||
real, pointer:: p | real, pointer:: p | ||
allocate(p) | allocate(p) | ||
Line 1,026: | Line 2,058: | ||
p | p | ||
7.9 | 7.9 | ||
− | + | | | |
− | + | </source> | |
− | |||
samples/pointers/ptr2.f90 | samples/pointers/ptr2.f90 | ||
+ | |} | ||
− | ==What are they good | + | ==What are they good for? (1)== |
− | for? (1) | + | {| |
− | * Pointers are essential for | + | |- valign="top" |
− | creating, maintaining | + | | |
− | dynamic data structures | + | * Pointers are essential for creating, maintaining dynamic data structures |
* Linked lists, trees, heaps.. | * Linked lists, trees, heaps.. | ||
− | * Some of these can be | + | * Some of these can be sort-of implemented in arrays, but very awkward |
− | sort-of implemented in | + | * Adaptive meshes, treebased particle solvers need these structures. |
− | 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:Singly-linked-list.svg | ||
− | |||
http://en.wikipedia.org/wiki/File:Max-Heap.svg | http://en.wikipedia.org/wiki/File:Max-Heap.svg | ||
+ | |} | ||
− | ==What are they good | + | ==What are they good for? (2)== |
− | for? (2) | + | {| |
− | * A pointer can be of | + | |- valign="top" |
− | array type, not just | + | | |
− | scalar | + | * A pointer can be of array type, not just scalar |
− | * Fortran pointers + | + | * Fortran pointers + fortran arrays are quite interesting; can create “views” of subarrays |
− | fortran arrays are quite | + | | |
− | interesting; can create | + | <source lang="fortran"> |
− | “views” of subarrays | ||
− | |||
real, target, dimension(7) :: x | real, target, dimension(7) :: x | ||
real, pointer:: p(:) | real, pointer:: p(:) | ||
Line 1,077: | Line 2,103: | ||
7 | 7 | ||
− | + | </source> | |
− | |||
− | |||
samples/pointers/views.f90 | samples/pointers/views.f90 | ||
+ | |} | ||
==Hands on #4== | ==Hands on #4== | ||
− | * Use pointers to provide views into subsets of | + | {| |
− | the arrays in solver.f90 to clarify the functions. | + | |- valign="top" |
− | * In workedexamples/pointers, have started | + | | |
− | the process with cfl, hydroflux; try tackling | + | * Use pointers to provide views into subsets of the arrays in solver.f90 to clarify the functions. |
− | tvd1d, others. | + | * In workedexamples/pointers, have started the process with cfl, hydroflux; try tackling tvd1d, others. |
* ~30 min | * ~30 min | ||
+ | | | ||
+ | |} | ||
− | + | =Derived Types and Objects= | |
− | * Often, groups of | + | {| |
− | variables naturally go | + | |- valign="top" |
− | together to represent a | + | | |
− | larger structure | + | * Often, groups of variables naturally go together to represent a larger structure |
− | * Whenever you find | + | * Whenever you find yourself passing the same group of variables to several routines, a good candidate for a derived type. |
− | yourself passing the | + | <source lang="fortran"> |
− | same group of variables | ||
− | to several routines, a | ||
− | good candidate for a | ||
− | derived type. | ||
− | |||
type griddomain | type griddomain | ||
− | real :: xmin, xmax | + | real :: xmin, xmax |
− | real :: ymin, ymax | + | real :: ymin, ymax |
− | real :: nx, ny | + | real :: nx, ny |
− | real, dimension(:,:) :: u | + | real, dimension(:,:) :: u |
endtype griddomain | endtype griddomain | ||
+ | |||
type(griddomain) :: g | type(griddomain) :: g | ||
g % xmin = -1 | g % xmin = -1 | ||
g % xmax = +1 | g % xmax = +1 | ||
+ | </source> | ||
+ | |} | ||
==Derived Types and Objects== | ==Derived Types and Objects== | ||
− | * Consider interval | + | {| |
− | arithmetic (good for | + | |- valign="top" |
− | quantification of | + | | |
− | uncertainties, etc). | + | * Consider interval arithmetic (good for quantification of uncertainties, etc). |
− | * An interval inherently | + | * An interval inherently has two values associated with it - the end points. |
− | has two values | ||
− | associated with it - the | ||
− | end points. | ||
* Can make this a type. | * Can make this a type. | ||
− | + | * 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. | |
− | * Note can access the | + | | |
− | fields in the type with | + | (from samples/derivedtypes/simple/intervalmath.f90 ) |
− | “%” | + | | |
− | * 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== | ==Procedures using types== | ||
− | + | {| | |
− | + | |- valign="top" | |
− | + | | | |
− | + | * 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. | |
− | + | | | |
− | + | (from samples/derivedtypes/intervalfunctions/intervalmath.f90 ) | |
− | + | |} | |
− | |||
− | |||
− | |||
− | |||
− | * 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/ | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
==Procedures using types== | ==Procedures using types== | ||
{| | {| | ||
+ | |- valign="top" | ||
| | | | ||
+ | * Would prefer not to have to treat integer and real intervals so differently in main program | ||
+ | * Different types, but adding should be similar. | ||
* Would like to be able to call “addintervals” and have language call the right subroutine, do the right thing. | * Would like to be able to call “addintervals” and have language call the right subroutine, do the right thing. | ||
* Similar to how intrinsics work - sin() works on any kind of real, matmult() works on integer, real, or complex matricies. | * Similar to how intrinsics work - sin() works on any kind of real, matmult() works on integer, real, or complex matricies. | ||
Line 1,183: | Line 2,175: | ||
==Generic Interfaces== | ==Generic Interfaces== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
{| | {| | ||
* Generic Interfaces | * Generic Interfaces | ||
Line 1,193: | Line 2,188: | ||
==Generic Interfaces== | ==Generic Interfaces== | ||
− | * Note that everything is | + | {| |
− | private except what is | + | |- valign="top" |
− | explicitly made public. | + | | |
+ | * Note that everything is private except what is explicitly made public. | ||
* Types are public. | * Types are public. | ||
− | * Generic interfaces are | + | * Generic interfaces are public. |
− | public. | + | * Type specific routines are not. |
− | * Type specific routines | + | * Program using interval math sees only the generic interfaces. |
− | are not. | + | | |
− | * Program using interval | ||
− | math sees only the | ||
− | generic interfaces. | ||
samples/derivedtypes/genericintervals/intervalmath.f90 | samples/derivedtypes/genericintervals/intervalmath.f90 | ||
+ | |} | ||
==Generic interfaces== | ==Generic interfaces== | ||
− | * Call to addintervals or | + | {| |
− | subtract intervals goes | + | |- valign="top" |
− | to the correct typespecific routine. | + | | |
+ | * Call to addintervals or subtract intervals goes to the correct typespecific routine. | ||
* As does print interval. | * As does print interval. | ||
− | * Could create routines | + | * Could create routines to add real to int interval, etc and add to the same interface. |
− | to add real to int | + | | |
− | interval, etc and add to | ||
− | the same interface. | ||
samples/derivedtypes/genericintervals/interval2.f90 | samples/derivedtypes/genericintervals/interval2.f90 | ||
+ | |} | ||
==Operator overloading== | ==Operator overloading== | ||
− | * An infix operator is | + | {| |
− | really just “syntactic | + | |- valign="top" |
− | sugar” for a function | + | | |
− | which takes two | + | * An infix operator is really just “syntactic sugar” for a function which takes two operands and returns a third. |
− | operands and returns a | + | | |
− | third. | + | <source lang="fortran" |
− | |||
a = b (op) c | a = b (op) c | ||
=> | => | ||
function op(b,c) | function op(b,c) | ||
returns a | returns a | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
a=b | a=b | ||
=> | => | ||
subroutine assign(a,b) | subroutine assign(a,b) | ||
+ | </source> | ||
+ | |} | ||
==Operator overloading== | ==Operator overloading== | ||
− | * Here, we’ve defined | + | {| |
− | two subroutines which | + | |- valign="top" |
− | set intervals based on | + | | |
− | an array - 2 ints for an | + | * 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 |
− | integer interval, or 2 | + | | |
− | reals for a real interval | ||
− | |||
samples/derivedtypes/intervaloperators/intervalmath.f90 | samples/derivedtypes/intervaloperators/intervalmath.f90 | ||
+ | |} | ||
==Generic interfaces== | ==Generic interfaces== | ||
− | * Once this is done, can | + | {| |
− | use assignment | + | |- valign="top" |
− | operator, | + | | |
− | * Or add, subtract | + | * Once this is done, can use assignment operator, |
− | multiply intervals. | + | * Or add, subtract multiply intervals. |
− | * Can even compose | + | * Can even compose them in complex expressions! Functions automatically composed. |
− | them in complex | + | | |
− | expressions! Functions | ||
− | automatically | ||
− | composed. | ||
samples/derivedtypes/intervaloperators/interval3.f90 | samples/derivedtypes/intervaloperators/interval3.f90 | ||
+ | |} | ||
==Type bound procedures== | ==Type bound procedures== | ||
− | * Types can have not only | + | {| |
− | variables, but | + | |- valign="top" |
− | procedures. | + | | |
− | * Takes us from a type to | + | * Types can have not only variables, but procedures. |
− | what is usually called a | + | * Takes us from a type to what is usually called a class. |
− | class. | + | | |
− | + | (from samples/derivedtypes/intervaltypebound/intervalmath.f90 ) | |
− | samples/derivedtypes/intervaltypebound/ | + | |} |
− | intervalmath.f90 | ||
==Type bound procedures== | ==Type bound procedures== | ||
− | * Called like one accesses | + | {| |
− | a field - % | + | |- valign="top" |
− | * Operates on data of | + | | |
− | the particular variable it | + | * Called like one accesses a field - % |
− | is invoked from | + | * Operates on data of the particular variable it is invoked from |
− | + | | | |
samples/derivedtypes/intervaltypebound/interval3.f90 | samples/derivedtypes/intervaltypebound/interval3.f90 | ||
+ | |} | ||
==Type bound procedures== | ==Type bound procedures== | ||
− | * It is implicitly passed as | + | {| |
− | it’s first parameter the | + | |- valign="top" |
− | variable itself. | + | | |
− | * Can take other | + | * It is implicitly passed as it’s first parameter the variable itself. |
− | arguments as well. | + | * Can take other arguments as well. |
+ | | | ||
+ | samples/derivedtypes/intervaltypebound/intervalmath.f90 | ||
+ | |} | ||
− | + | ==Object oriented programming== | |
− | + | {| | |
+ | |- valign="top" | ||
+ | | | ||
+ | * 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= |
− | + | {| | |
− | + | |- valign="top" | |
− | + | | | |
− | + | * Large scientific software now frequently uses multiple languages, either within a single code or between codes. | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | * Large scientific software now frequently uses | ||
− | multiple languages, either within a single code | ||
− | or between codes. | ||
* Right tool for the job! | * Right tool for the job! | ||
− | * Need to know how to make software | + | * Need to know how to make software interact. |
− | interact. | + | * Here we’ll look at C/Fortran code calling each other, and calling Fortran code from python. |
− | * Here we’ll look at C/Fortran code calling | + | | |
− | each other, and calling Fortran code from | + | |} |
− | python. | ||
==C-interoperability== | ==C-interoperability== | ||
− | * iso_c_binding module contains definitions for | + | {| |
− | interacting with C | + | |- valign="top" |
+ | | | ||
+ | * iso_c_binding module contains definitions for interacting with C | ||
* Types, ways to bind procedures to C, etc. | * Types, ways to bind procedures to C, etc. | ||
− | * Allows you to call C routines, or bind Fortran | + | * Allows you to call C routines, or bind Fortran routines in a way that they can be called by C. |
− | routines in a way that they can be called by C. | + | | |
+ | |} | ||
− | ==Calling a C routine== | + | ==Calling a C routine from Fortran== |
− | + | {| | |
− | * As with the case of | + | |- valign="top" |
− | calling a passed-in | + | | |
− | function, need an | + | * As with the case of calling a passed-in function, need an explicit prototype. |
− | explicit prototype. | + | * Tells compiler what to do with “sqrtc” function when called. |
− | * Tells compiler what | + | | |
− | to do with “sqrtc” | ||
− | function when called. | ||
samples/C-interop/call-c-fn/main.f90 | samples/C-interop/call-c-fn/main.f90 | ||
+ | |} | ||
− | ==Calling a C routine== | + | ==Calling a C routine from Fortran== |
− | + | {| | |
− | * BIND(C) - tells compiler/ | + | |- valign="top" |
− | linker will have a C-style, | + | | |
− | rather than fortran-style | + | * BIND(C) - tells compiler/ linker will have a C-style, rather than fortran-style name in object file. |
− | name in object file. | + | * Can optionally supply another name; otherwise, will default to procedure name. |
− | * Can optionally supply | + | * Here we’re calling it sqrtc to avoid Fortran sqrt() function. |
− | another name; otherwise, | + | * The value the function takes and returns are C doubles; that is, reals of kind(c_double). |
− | will default to procedure | + | * Also defined: c_float, integers of kind c_int, c_long, etc. |
− | name. | + | | |
− | * Here we’re calling it sqrtc | + | |} |
− | to avoid Fortran sqrt() | ||
− | function. | ||
− | ==Calling a C routine== | + | ==Calling a C routine from Fortran== |
− | + | {| | |
− | + | |- valign="top" | |
− | + | | | |
− | + | * 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. | |
− | + | * Compiling - just make sure to link in C library you’re calling | |
− | |||
− | |||
− | |||
− | * 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. | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | * Compiling - just make | ||
− | sure to link in C library | ||
− | you’re calling | ||
* And that’s it. | * And that’s it. | ||
− | + | | | |
+ | <source lang="bash"> | ||
$ make | $ make | ||
gfortran -c main.f90 | gfortran -c main.f90 | ||
Line 1,390: | Line 2,352: | ||
2.0000000000000000 | 2.0000000000000000 | ||
C | C | ||
− | sqrt(x) = | + | sqrt(x) = 1.4142135623730951 |
− | 1.4142135623730951 | + | Fortran sqrt(x) = 1.4142135623730951 |
− | Fortran sqrt(x) = | + | </source> |
− | 1.4142135623730951 | + | |} |
==C strings== | ==C strings== | ||
− | * When using C | + | {| |
− | strings, you have to | + | |- valign="top" |
− | remember that C | + | | |
− | terminates strings | + | * When using C strings, you have to remember that C terminates strings with a special character |
− | with a special | ||
− | character | ||
* C_NULL_CHAR | * C_NULL_CHAR | ||
− | * Dealing with functions | + | * Dealing with functions that return strings is hairier, as they return a pointer, not an array. |
− | that return strings is | + | | |
− | hairier, as they return | + | <source lang="bash"> |
− | a pointer, not an array. | ||
− | |||
$ make | $ make | ||
gfortran-mp-4.4 -c main.f90 | gfortran-mp-4.4 -c main.f90 | ||
Line 1,414: | Line 2,372: | ||
1234 = | 1234 = | ||
1234 | 1234 | ||
− | + | </source> | |
− | samples/C-interop/c-strings/main.f90 | + | (from samples/C-interop/c-strings/main.f90 ) |
+ | |} | ||
==Calling Fortran from C== | ==Calling Fortran from C== | ||
− | * To write Fortran | + | {| |
− | routines callable from | + | |- valign="top" |
− | C, bind the subroutine | + | | |
− | to a C name | + | * To write Fortran routines callable from C, bind the subroutine to a C name |
− | * Can optionally give it a | + | * Can optionally give it a different name, as above |
− | different name, as above | + | * And make sure arguments are of C types. |
− | * And make sure | + | | |
− | arguments are of C | ||
− | types. | ||
− | |||
samples/C-interop/c-call-fortran/froutine.f90 | samples/C-interop/c-call-fortran/froutine.f90 | ||
+ | |} | ||
==Calling Fortran from C== | ==Calling Fortran from C== | ||
− | * Here, we’ll do | + | {| |
− | something a little | + | |- valign="top" |
− | trickier and pass C | + | | |
− | dynamic arrays | + | * Here, we’ll do something a little trickier and pass C dynamic arrays |
− | + | * 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. | |
− | + | * The single scalar argument passed back we just treat as an intent(out) variable | |
− | * 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. | ||
− | |||
− | |||
− | |||
− | |||
− | * The single scalar | ||
− | argument passed back | ||
− | we just treat as an | ||
− | intent(out) variable | ||
* Of type c_float. | * Of type c_float. | ||
− | samples/C-interop/c-call-fortran/froutine.f90 | + | | |
+ | (from samples/C-interop/c-call-fortran/froutine.f90 ) | ||
+ | |} | ||
==More advanced== | ==More advanced== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Handling arbitrary C code is possible | * Handling arbitrary C code is possible | ||
− | * Passing C structs -- create a Fortran derived | + | * Passing C structs -- create a Fortran derived type with the same fields, use BIND(C) in defining the type. |
− | type with the same fields, use BIND(C) in | + | * C “multidimensional arrays” - have to be careful, they may not be contiguous! Follow pointers. |
− | defining the type. | ||
− | * C “multidimensional arrays” - have to be | ||
− | careful, they may not be contiguous! Follow | ||
− | pointers. | ||
* Even taking passed function pointers to C | * Even taking passed function pointers to C | ||
− | functions is possible (samples/C-interop/ | + | functions is possible |
− | functionptr) | + | | |
− | + | (from samples/C-interop/functionptr ) | |
− | + | |} | |
− | |||
− | |||
− | |||
+ | ==Fortran calling C, which calls Fortran== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
samples/C-interop/valueref/driver.f90 | samples/C-interop/valueref/driver.f90 | ||
+ | | | ||
+ | |} | ||
==samples/C-interop/valueref/croutine.c== | ==samples/C-interop/valueref/croutine.c== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | | | ||
+ | |} | ||
==samples/C-interop/valueref/froutine.f90== | ==samples/C-interop/valueref/froutine.f90== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | | | ||
+ | |} | ||
==$ ./main== | ==$ ./main== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | <source lang="bash"> | ||
+ | $ ./main | ||
(1) In the FORTRAN routine before call to C | (1) In the FORTRAN routine before call to C | ||
(1) i = | (1) i = | ||
Line 1,495: | Line 2,452: | ||
30.000000000000000 | 30.000000000000000 | ||
(1) In the FORTRAN routine after call to C | (1) In the FORTRAN routine after call to C | ||
− | (1) i = | + | (1) i = 15 x = 3.0000000000000000 prod = 30.000000000000000 |
− | 15 x = | + | </source> |
− | 3.0000000000000000 | + | | |
− | prod = | + | |} |
− | 30.000000000000000 | ||
==samples/C-interop/valueref/Makefile== | ==samples/C-interop/valueref/Makefile== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | | | ||
+ | |} | ||
==F2py== | ==F2py== | ||
− | * Comes with scipy, | + | {| |
− | a widely-installed | + | |- valign="top" |
− | (by scientists) | + | | |
− | python package. | + | * Comes with scipy, a widely-installed (by scientists) python package. |
− | * Wraps fortran in | + | * Wraps fortran in python, allows you to easily interface. |
− | python, allows you | + | * fwrap is another option |
− | to easily interface. | + | | |
− | * fwrap is another | ||
− | option | ||
− | |||
http://www.scipy.org/F2py | http://www.scipy.org/F2py | ||
+ | |} | ||
==* Will only use solver module.== | ==* Will only use solver module.== | ||
− | * Unfortunately, f2py isn’t quite smart enough to understand | + | {| |
− | using parameters to size arrays, so global replace ‘nvars’=4 | + | |- valign="top" |
+ | | | ||
+ | * 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 | * module load use.experimental gcc python/2.7.1 intel/12 | ||
* “f2py -m hydro_solver -h hydro_solver.pyf solver.f90” | * “f2py -m hydro_solver -h hydro_solver.pyf solver.f90” | ||
+ | | | ||
+ | |} | ||
==* generates the following header file (hydro_solver.pyf)== | ==* generates the following header file (hydro_solver.pyf)== | ||
− | + | {| | |
+ | |- valign="top" | ||
+ | | | ||
* Comment out stuff we don’t need (lower-level routines) | * Comment out stuff we don’t need (lower-level routines) | ||
* f2py -c --fcompiler=gfortran hydro_solver.pyf solver.f90 | * f2py -c --fcompiler=gfortran hydro_solver.pyf solver.f90 | ||
+ | | | ||
+ | |} | ||
==$ ipython== | ==$ ipython== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | <source lang="bash"> | ||
In [1]: import hydro_solver | In [1]: import hydro_solver | ||
In [2]: hydro_solver.solver. | In [2]: hydro_solver.solver. | ||
Line 1,553: | Line 2,524: | ||
...dt = hydro_solver.solver.timestep(u) | ...dt = hydro_solver.solver.timestep(u) | ||
In [8]: pylab.imshow(u[1,:,:]) | In [8]: pylab.imshow(u[1,:,:]) | ||
+ | </source> | ||
+ | |} | ||
− | + | =Coarrays= | |
− | * In F2008, objects | + | {| |
− | can also have a | + | |- valign="top" |
− | “codimension”. | + | | |
− | * An object with a | + | * In F2008, objects can also have a “codimension”. |
− | codimension can be | + | * An object with a codimension can be “seen” across processes |
− | “seen” across | + | * Right now, only intel v 12 supports this |
− | processes | + | | |
− | * Right now, only intel | ||
− | v 12 supports this | ||
samples/coarrays/simple.f90 | samples/coarrays/simple.f90 | ||
+ | |} | ||
==Coarrays== | ==Coarrays== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
1 | 1 | ||
Line 1,587: | Line 2,562: | ||
independant | independant | ||
processes | processes | ||
+ | | | ||
+ | |} | ||
==Coarrays== | ==Coarrays== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
1 | 1 | ||
Line 1,615: | Line 2,595: | ||
Independent, parallel tasks can | Independent, parallel tasks can | ||
“see” each other’s data as easily as an array index. | “see” each other’s data as easily as an array index. | ||
+ | | | ||
+ | |} | ||
==Sychronization== | ==Sychronization== | ||
− | * When accessing other | + | {| |
− | processor’s data, must | + | |- valign="top" |
− | ensure that tasks are | + | | |
− | synchronized | + | * When accessing other processor’s data, must ensure that tasks are synchronized |
− | * Don’t want to read | + | * Don’t want to read task 1’s data before read is complete! |
− | task 1’s data before | ||
− | read is complete! | ||
* sync all | * sync all | ||
+ | | | ||
samples/coarrays/broadcast.f90 | samples/coarrays/broadcast.f90 | ||
+ | |} | ||
==Sychronization== | ==Sychronization== | ||
− | * Other synchronization | + | {| |
− | tools: | + | |- valign="top" |
− | * sync images([..]) syncs | + | | |
− | only images in the list | + | * Other synchronization tools: |
− | provided | + | * sync images([..]) syncs only images in the list provided |
− | * critical - only one | + | * critical - only one image can enter block at a time |
− | image can enter block | + | * lock - finer-grained control than critical. |
− | at a time | ||
− | * lock - finer-grained | ||
− | control than critical. | ||
* atomic operations. | * atomic operations. | ||
− | + | | | |
samples/coarrays/broadcast.f90 | samples/coarrays/broadcast.f90 | ||
+ | |} | ||
==1d Diffusion Eqn== | ==1d Diffusion Eqn== | ||
− | * Calculate heat | + | {| |
− | diffusion along a bar | + | |- valign="top" |
− | * Finite difference; | + | | |
− | calculate second | + | * Calculate heat diffusion along a bar |
− | derivative using | + | * Finite difference; calculate second derivative using neighbours |
− | neighbours | + | * Get new temperature from old temperatures |
− | * Get new | ||
− | temperature from | ||
− | old temperatures | ||
− | |||
1 -2 1 | 1 -2 1 | ||
+ | | | ||
+ | |} | ||
==1d Diffusion Eqn== | ==1d Diffusion Eqn== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Initial conditions: | * Initial conditions: | ||
− | * Use pointers to | + | * Use pointers to point to new, old temperatures |
− | point to new, old | + | * 1..totpoints+2 (pts 1, totpoints+1 “ghost cells) |
− | temperatures | + | * Setup x, orig temperature, expected values |
− | * 1..totpoints+2 (pts 1, | + | | |
− | totpoints+1 “ghost | ||
− | cells) | ||
− | * Setup x, orig | ||
− | temperature, expected | ||
− | values | ||
− | |||
samples/coarrays/diffusion/diffusion.f90 | samples/coarrays/diffusion/diffusion.f90 | ||
+ | |} | ||
==1d Diffusion Eqn== | ==1d Diffusion Eqn== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Evolution: | * Evolution: | ||
* Apply BCs | * Apply BCs | ||
− | * Apply finite | + | * Apply finite difference stencil to all real data points |
− | difference stencil to | + | | |
− | all real data points | ||
samples/coarrays/diffusion/diffusion.f90 | samples/coarrays/diffusion/diffusion.f90 | ||
+ | |} | ||
==1d Diffusion Eqn== | ==1d Diffusion Eqn== | ||
− | * Output calculated | + | {| |
− | values and theory to | + | |- valign="top" |
− | file output.txt | + | | |
− | * Note: Fortran2008, | + | * Output calculated values and theory to file output.txt |
− | can use “newunit” | + | * Note: Fortran2008, can use “newunit” to find, use an unused IO unit |
− | to find, use an | + | | |
− | unused IO unit | ||
− | |||
samples/coarrays/diffusion/diffusion.f90 | samples/coarrays/diffusion/diffusion.f90 | ||
+ | |} | ||
==$ ./diffusion== | ==$ ./diffusion== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
[..] | [..] | ||
$ gnuplot | $ gnuplot | ||
Line 1,694: | Line 2,675: | ||
gnuplot> plot 'output.txt' using 1:2 with points title 'Simulated', | gnuplot> plot 'output.txt' using 1:2 with points title 'Simulated', | ||
'output.txt' using 1:3 with lines title 'Theory' | 'output.txt' using 1:3 with lines title 'Theory' | ||
+ | | | ||
+ | |} | ||
==Coarray Diffusion== | ==Coarray Diffusion== | ||
− | * Now we’ll try this in | + | {| |
− | parallel | + | |- valign="top" |
− | * Each image runs it’s part | + | | |
− | of the problem | + | * Now we’ll try this in parallel |
− | (totpoints/num_images()) | + | * Each image runs it’s part of the problem (totpoints/num_images()) |
− | * Communications is like | + | * Communications is like boundary condition handling - except boundary data must be obtained from neighbouring image. |
− | boundary condition | + | | |
− | handling - except | + | |} |
− | boundary data must be | ||
− | obtained from | ||
− | neighbouring image. | ||
==Coarray Diffusion== | ==Coarray Diffusion== | ||
− | * Everything’s the same | + | {| |
− | except we are only | + | |- valign="top" |
− | responsible for locpoints | + | | |
− | points, starting at start. | + | * Everything’s the same except we are only responsible for locpoints points, starting at start. |
− | * Once calculated, we | + | * Once calculated, we never need totpoints again. All arrays are of size (locpoints), etc. |
− | never need totpoints | + | * For simplicity, we assume here everything divides evenly. |
− | again. All arrays are of | + | | |
− | size (locpoints), etc. | ||
− | * For simplicity, we | ||
− | assume here everything | ||
− | divides evenly. | ||
− | |||
samples/coarrays/diffusion/diffusion-coarray.f90 | samples/coarrays/diffusion/diffusion-coarray.f90 | ||
+ | |} | ||
==Coarray Diffusion== | ==Coarray Diffusion== | ||
− | * Boundary exchange; if we | + | {| |
− | have interior neighbours, | + | |- valign="top" |
− | get updated data from | + | | |
− | them so we can calculate | + | * Boundary exchange; if we have interior neighbours, get updated data from them so we can calculate our new data |
− | our new data | + | * Note: can’t use pointers here, coarrays can’t be targets. |
− | * Note: can’t use pointers | + | * Sync all around boundary exchange a little heavy handed; could just sync neighbours. |
− | 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 | samples/coarrays/diffusion/diffusion-coarray.f90 | ||
+ | |} | ||
==Coarray Diffusion== | ==Coarray Diffusion== | ||
− | + | {| | |
+ | |- valign="top" | ||
+ | | | ||
+ | <source lang="fortran"> | ||
temperature()[leftneigh] | temperature()[leftneigh] | ||
Line 1,748: | Line 2,723: | ||
lnpts+2 | lnpts+2 | ||
+ | </source> | ||
+ | |} | ||
==Coarray Diffusion== | ==Coarray Diffusion== | ||
− | * Everything else exactly | + | {| |
− | the same. | + | |- valign="top" |
− | * For I/O, we each output | + | | |
− | our own file, prefixed by | + | * Everything else exactly the same. |
− | our image number | + | * For I/O, we each output our own file, prefixed by our image number |
* (eg, 1-output.txt, 2output.txt...) | * (eg, 1-output.txt, 2output.txt...) | ||
− | + | | | |
samples/coarrays/diffusion/diffusion-coarray.f90 | samples/coarrays/diffusion/diffusion-coarray.f90 | ||
+ | |} | ||
==$ export FOR_COARRAY_NUM_IMAGES=3== | ==$ export FOR_COARRAY_NUM_IMAGES=3== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
+ | <source lang="fortran"> | ||
# 3 images | # 3 images | ||
$ ./diffusion-coarray | $ ./diffusion-coarray | ||
Line 1,767: | Line 2,749: | ||
gnuplot> plot '1-ics.txt' using 1:2, '2-ics.txt' using 1:2, '3-ics.txt' | gnuplot> plot '1-ics.txt' using 1:2, '2-ics.txt' using 1:2, '3-ics.txt' | ||
using 1:2 | using 1:2 | ||
+ | </source> | ||
+ | |} | ||
==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',== | ==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',== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
'2-output.txt' using 1:3 with lines title 'Theory', '3-output.txt' using 1:3 | '2-output.txt' using 1:3 with lines title 'Theory', '3-output.txt' using 1:3 | ||
with lines title 'Theory' | with lines title 'Theory' | ||
+ | | | ||
+ | |} | ||
==Parallel Matrix Multiplication== | ==Parallel Matrix Multiplication== | ||
− | * Consider matrix | + | {| |
− | operations, where | + | |- valign="top" |
− | matrix is | + | | |
− | decomposed onto | + | * Consider matrix operations, where matrix is decomposed onto images in blocks |
− | images in blocks | + | * A given image has corresponding blocks in A, B, responsible for that block in C. |
− | * A given image has | ||
− | corresponding blocks | ||
− | in A, B, responsible | ||
− | for that block in C. | ||
A | A | ||
Line 1,789: | Line 2,774: | ||
C | C | ||
= | = | ||
+ | | | ||
+ | |} | ||
==Parallel Matrix Multiplication== | ==Parallel Matrix Multiplication== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Block matrix multiply exactly like matrix multiply, except each operation is a submatrix operation | * Block matrix multiply exactly like matrix multiply, except each operation is a submatrix operation | ||
* matmul instead of “*”. | * matmul instead of “*”. | ||
Line 1,801: | Line 2,791: | ||
C | C | ||
= | = | ||
+ | | | ||
+ | |} | ||
==Parallel Matrix Multiplication== | ==Parallel Matrix Multiplication== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Each image gets its own block of a, b, c. | * Each image gets its own block of a, b, c. | ||
* Note [nrows,*]. Coarrays can have any corank(eg, number of codimensions) | * Note [nrows,*]. Coarrays can have any corank(eg, number of codimensions) | ||
Line 1,809: | Line 2,804: | ||
samples/coarrays/blockmatrix.f90 | samples/coarrays/blockmatrix.f90 | ||
+ | | | ||
+ | |} | ||
==Parallel Matrix Multiplication== | ==Parallel Matrix Multiplication== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Allocation, deallocation of coarrays is like a sync all | * Allocation, deallocation of coarrays is like a sync all | ||
* Forces synchronization across all images | * Forces synchronization across all images | ||
Line 1,816: | Line 2,816: | ||
samples/coarrays/blockmatrix.f90 | samples/coarrays/blockmatrix.f90 | ||
+ | | | ||
+ | |} | ||
==Parallel Matrix Multiplication== | ==Parallel Matrix Multiplication== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* This is the entire parallel multiplication operation. | * This is the entire parallel multiplication operation. | ||
* a[myrow,k] gets the entire a block from image at myrow,k. | * a[myrow,k] gets the entire a block from image at myrow,k. | ||
Line 1,823: | Line 2,828: | ||
samples/coarrays/blockmatrix.f90 | samples/coarrays/blockmatrix.f90 | ||
+ | | | ||
+ | |} | ||
==Coarray Summary== | ==Coarray Summary== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Coarrays are powerful, simple-to-use addition to parallel programming models available | * Coarrays are powerful, simple-to-use addition to parallel programming models available | ||
* Will be widely available soon | * Will be widely available soon | ||
Line 1,830: | Line 2,840: | ||
* Typically implemented with MPI under the hood; can give performance similar to MPI with fewer lines of code. | * 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. | * Other languages have similar extensions (UPC based on C, Titanium based on Java) but are not part of languages’ standard. | ||
+ | | | ||
+ | |} | ||
==Closing Hints== | ==Closing Hints== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* Always give the compiler info it needs to help you by being as explicit as possible | * 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. | * 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). | * Always get as much info from compiler as possible - always use -Wall (gfortran) or -warn all (ifort). | ||
+ | | | ||
+ | |} | ||
==Useful Resources== | ==Useful Resources== | ||
+ | {| | ||
+ | |- valign="top" | ||
+ | | | ||
* http://fortranwiki.org/ | * http://fortranwiki.org/ | ||
− | * Reference source; has all standards; Fortran2003/2008 status of major compilers | + | ** Reference source; has all standards; Fortran2003/2008 status of major compilers |
* http://en.wikipedia.org/wiki/Fortran_language_features | * http://en.wikipedia.org/wiki/Fortran_language_features | ||
− | * Succinct summary of new features (spotty past F95) | + | ** Succinct summary of new features (spotty past F95) |
* http://stackoverflow.com/questions/tagged/fortran | * http://stackoverflow.com/questions/tagged/fortran | ||
− | * Programmers Questions & Answers | + | ** ( Programmers Questions & Answers |
+ | |} |
Latest revision as of 14:59, 23 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
|
<source lang="fortran"> module gravity implicit none real, parameter :: G = 6.67e-11 ! MKS units contains real function gravforce(x1,x2,m1,m2) implicit none real, dimension(3), intent(in) :: x1,x2 real, intent(in) :: m1, m2 real :: dist dist = sqrt(sum((x1-x2)**2)) gravforce = G * m1 * m2 / dist**2 end function gravforce end module gravity program simplemod use gravity implicit none print *, 'Gravitational constant = ', G print *, 'Force between 2 1kg masses at [1,0,0] & &and [0,0,1] is' print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.) end program simplemod </source> (from samples/procedures/simplemod/simplemod.f90) |
Compiling & Running
|
<source lang="bash"> $ ls simplemod.f90 $ gfortran -o simplemod simplemod.f90 -Wall $ ls gravity.mod simplemod simplemod.f90 $ ./simplemod Gravitational constant = 6.6700000E-11 Force between 2 1kg masses at [1,0,0] and [0,0,1] is 3.3350003E-11 </source> |
Modules
|
<source lang="fortran"> <source lang="fortran"> module gravity implicit none real, parameter :: G = 6.67e-11 ! MKS units contains real function gravforce(x1,x2,m1,m2) implicit none real, dimension(3), intent(in) :: x1,x2 real, intent(in) :: m1, m2 real :: dist dist = sqrt(sum((x1-x2)**2)) gravforce = G * m1 * m2 / dist**2 end function gravforce end module gravity program simplemod use gravity implicit none print *, 'Gravitational constant = ', G print *, 'Force between 2 1kg masses at [1,0,0] & &and [0,0,1] is' print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.) end program simplemod </source> (from samples/procedures/simplemod/simplemod.f90) |
use module, only :
|
<source lang="fortran"> module gravity implicit none real, parameter :: G = 6.67e-11 ! MKS units contains real function gravforce(x1,x2,m1,m2) implicit none real, dimension(3), intent(in) :: x1,x2 real, intent(in) :: m1, m2 real :: dist dist = sqrt(sum((x1-x2)**2)) gravforce = G * m1 * m2 / dist**2 end function gravforce end module gravity program simplemod2 use gravity, only : G, gravforce implicit none print *, 'Gravitational constant = ', G print *, 'Force between 2 1kg masses at [1,0,0] & &and [0,0,1] is' print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.) end program simplemod2 </source> samples/procedures/simplemod/simplemod2.f90 |
Modules usually get their own files
|
<source lang="fortran"> module gravity implicit none private character (len=8), parameter, public :: massunit="kilogram" character (len=8), parameter, public :: forceunit="Newton" public :: gravforce real, parameter :: G = 6.67e-11 ! MKS units contains real function distance(x1,x2) implicit none real, dimension(3), intent(in) :: x1, x2 distance = sqrt(sum((x1-x2)**2)) end function distance real function gravforce(x1,x2,m1,m2) implicit none real, dimension(3), intent(in) :: x1,x2 real, intent(in) :: m1, m2 real :: dist dist = distance(x1,x2) gravforce = G * m1 * m2 / dist**2 end function gravforce end module gravity </source> (from samples/procedures/multifilemod/gravity.f90 ) |
Modules usually get their own files
|
<source lang="make"> FC=gfortran FFLAGS=-O3 -Wall multifilemod: multifilemod.o gravity.o $(FC) -o $@ multifilemod.o gravity.o %.mod: %.f90 $(FC) $(FFLAGS) -c $< multifilemod.o: multifilemod.f90 gravity.mod $(FC) $(FFLAGS) -c $< clean: rm -f *.o *~ *.mod multifilemod </source> (from samples/procedures/multifilemod/Makefile) |
.mod needed for compilation
|
<source lang="fortran"> program simplemod2 use gravity, only : gravforce, massunit, forceunit implicit none print *, 'Force between 2 1 ', massunit ,' masses ', & ' at [1,0,0] and [0,0,1] is' print *, gravforce([1.,0.,0.],[0.,0.,1.],1.,1.), forceunit end program simplemod2 </source> (from samples/procedures/multifilemod/multifilemod.f90) |
.o needed for linking
|
<source lang="make"> FC=gfortran FFLAGS=-O3 -Wall multifilemod: multifilemod.o gravity.o $(FC) -o $@ multifilemod.o gravity.o %.mod: %.f90 $(FC) $(FFLAGS) -c $< multifilemod.o: multifilemod.f90 gravity.mod $(FC) $(FFLAGS) -c $< clean: rm -f *.o *~ *.mod multifilemod </source> (from samples/procedures/multifilemod/Makefile) |
Compiling and running
|
<source lang="bash"> $ make gfortran -O3 -Wall -c gravity.f90 gfortran -O3 -Wall -c multifilemod.f90 gfortran -o multifilemod multifilemod.o gravity.o reposado-$ ./multifilemod Force between 2 1 kilogram masses at [1,0,0] and [0,0,1] is 3.33500033E-11 Newton </source> |
Private and public
|
<source lang="fortran"> module gravity implicit none private character (len=8), parameter, public :: massunit="kilogram" character (len=8), parameter, public :: forceunit="Newton" public :: gravforce real, parameter :: G = 6.67e-11 ! MKS units contains real function distance(x1,x2) implicit none real, dimension(3), intent(in) :: x1, x2 distance = sqrt(sum((x1-x2)**2)) end function distance real function gravforce(x1,x2,m1,m2) implicit none real, dimension(3), intent(in) :: x1,x2 real, intent(in) :: m1, m2 real :: dist dist = distance(x1,x2) gravforce = G * m1 * m2 / dist**2 end function gravforce end module gravity </source> ( from samples/procedures/multifilemod/gravity.f90 ) |
Procedures
|
<source lang="fortran"> module procedures contains function square(x) result(xsquared) implicit none real :: xsquared real, intent(IN) :: x xsquared = x*x end function square function cube(x) implicit none real :: cube real, intent(IN) :: x cube = x*x*x end function cube subroutine squareAndCube(x, squarex, cubex) implicit none real, intent(in) :: x real, intent(out) :: squarex real, intent(out) :: cubex squarex = square(x) cubex = cube(x) end subroutine squareAndCube end module procedures </source> ( from samples/procedures/funcsub/procedures.f90 ) |
Functions
|
<source lang="fortran"> function square(x) result(xsquared) implicit none real :: xsquared real, intent(IN) :: x xsquared = x*x end function square function cube(x) implicit none real :: cube real, intent(IN) :: x cube = x*x*x end function cube subroutine squareAndCube(x, squarex, cubex) implicit none real, intent(in) :: x real, intent(out) :: squarex real, intent(out) :: cubex squarex = square(x) cubex = cube(x) end subroutine squareAndCube </source> ( from samples/procedures/funcsub/procedures.f90 ) |
Procedure interfaces
|
<source lang="fortran"> function square(x) result(xsquared) implicit none real :: xsquared real, intent(IN) :: x xsquared = x*x end function square function cube(x) implicit none real :: cube real, intent(IN) :: x cube = x*x*x end function cube subroutine squareAndCube(x, squarex, cubex) implicit none real, intent(in) :: x real, intent(out) :: squarex real, intent(out) :: cubex squarex = square(x) cubex = cube(x) end subroutine squareAndCube </source> ( from samples/procedures/funcsub/procedures.f90 ) |
Procedure interfaces
|
<source lang="fortran"> function integratefx(xlo, xhi, f, n) |
integrate with trapezoid rule | ....
integer :: i real :: dx, xleft, xright integratefx = 0. dx = (xhi-xlo)/n xleft = xlo do i=0, n-1 xright = xleft + dx integratefx = integratefx + dx*(f(xright)+f(xleft))/2. xleft = xright enddo end function integratefx </source> (from samples/procedures/interface/integrate.f90 ) (from http://en.wikipedia.org/wiki/File:Trapezoidal_rule_illustration_small.svg ) |
---|
Procedure interfaces
|
<source lang="fortran"> function integratefx(xlo, xhi, f, n) |
integrate with trapezoid rule
implicit none real, intent(in) :: xlo, xhi interface function f(x) implicit none real :: f real, intent(in) :: x end function f end interface integer, intent(in) :: n real :: integratefx integer :: i real :: dx, xleft, xright integratefx = 0. dx = (xhi-xlo)/n xleft = xlo do i=0, n-1 xright = xleft + dx integratefx = integratefx + dx*(f(xright)+f(xleft))/2. xleft = xright enddo end function integratefx </source> (from samples/procedures/interface/integrate.f90 ) |
---|
Recursive procedures
|
<source lang="fortran"> recursive function integratefx(xlo, xhi, f, tol) result(integral) |
integrate with trapezoid rule, simpsons rule; | if difference between two is larger than | relevant tolerance, subdivide region. | ...variable declarations as before...
dx = xhi-xlo xmid = (xlo+xhi)/2. trapezoid = dx*(f(xlo)+f(xhi))/2. simpsons = dx/6.*(f(xlo)+4.*f(xmid)+f(xhi)) error = abs(trapezoid-simpsons)/(0.5*(trapezoid+simpsons)) if (error > tol) then |
too coarse; subdivide
integral = integratefx(xlo,xmid,f,tol) + & integratefx(xmid,xhi,f,tol) else integral = trapezoid endif end function integratefx </source> ( from samples/procedures/recursive/integrate.f90) |
---|
Pure procedures
|
<source lang="fortran"> pure subroutine axpy(a, x, y) |
y = y + a*x
implicit none real, intent(IN) :: a, x real, intent(INOUT) :: y y = y + a*x end subroutine axpy subroutine printaxpy(a, x, y) |
y = y + a*x
implicit none real, intent(IN) :: a, x real, intent(INOUT) :: y print *, a, '*', x, ' + ', y, & ' = ', a*x+y y = a*x + y end subroutine printaxpy </source> (from samples/procedures/purity/purity.f90) |
---|
Optional Arguments
|
<source lang="fortran"> recursive function integratefx(xlo, xhi, f, tol) result(integral) |
....
real, intent(in), optional :: tol |
....
real :: errtol |
use parameter if passed, | else use default
if (present(tol)) then errtol = tol else errtol = 1.e-6 endif |
....
if (error > errtol) then |
too coarse; subdivide
integral = integratefx(xlo,xmid,f,errtol) + & integratefx(xmid,xhi,f,errtol) else integral = trapezoid endif end function integratefx </source> (from samples/procedures/optional/integrate.f90) |
---|
Optional Arguments
|
<source lang="fortran"> print *, 'Integrating using default tol' approx = integratefx(0., 2*pi, sinesquared) print *, 'Approximate integral = ', approx print *, 'Exact integral = ', exact print *, print *, 'Integrating using coarser tol' approx = integratefx(0., 2*pi, sinesquared, 0.01) print *, 'Approximate integral = ', approx |
....
</source> (from samples/procedures/optional/optional.f90) |
---|
Keyword Arguments
|
<source lang="fortran"> |
....
print *, 'Integrating using still coarser tol' approx = integratefx(xhi=2*pi, xlo=0., tol=0.5, & f=sinesquared) print *, 'Approximate integral = ', approx </source> (from samples/procedures/optional/optional.f90) |
---|
Procedures and Modules Summary
|
Hands On 2
|
Fortran arrays
|
<source lang="fortran"> program basicarrays implicit none integer, dimension(5) :: a, b, c integer :: i a = [1,2,3,4,5] b = [(2*i+1, i=1,5)]
print *, 'a = ', a print *, 'b = ', b c = a+b print *, 'c = ', c c = a*b + 1 print *, 'a*b+1=', c end program basicarrays </source> (from samples/arrays/basic.f90 ) |
Array constructors
|
<source lang="fortran"> x = [1,2,3,4,5] x = (/1,2,3,4,5/) x = [ (i,i=1,5)] a = [ ((i*j,j=1,3),i=1,5)] </source> |
Elementwise operations
|
<source lang="fortran"> program elementwise implicit none real, dimension(10) :: x,y,z integer :: i real, parameter:: pi = 4.*atan(1.) x = [(2*pi*(i-1)/9.,i=1,10)] y = sin(x) z = x*x print *, x print *, y print *, z end program elementwise </source> (from samples/arrays/elementwise.f90 ) |
Elemental Functions
|
<source lang="fortran"> program elementalfn implicit none real, dimension(10) :: x,y,z integer :: i real, parameter:: pi = 4.*atan(1.) x = [(2*pi*(i-1)/9.,i=1,10)] y = sinesquared(x) z = sin(x)*sin(x) print *, x print *, y print *, z print *,z(::3) contains elemental function sinesquared(x) implicit none real :: sinesquared real, intent(in) :: x sinesquared = sin(x)**2 end function sinesquared end program elementalfn </source> (from samples/arrays/elemental.f90 ) |
Array comparisons
|
<source lang="fortran"> program comparearrays implicit none integer, dimension(5) :: a, b integer :: i a = [1,2,3,4,5] b = [(2*i-3, i=1,5)] print *, 'A = ', a print *, 'B = ', b if (any(a > b)) then print *, 'An A is larger than a B' endif if (all(a > b)) then print *, 'All As ares larger than Bs' else if (all(b > a)) then print *, 'All Bs are larger than As' else print *, 'A, B values overlap' endif end program comparearrays </source> (from samples/arrays/compare.f90) |
Array masks
|
<source lang="fortran"> program mask implicit none integer, dimension(10) :: a logical, dimension(10) :: pos integer :: i a = [(2*i-7, i=1,10)] pos = (a > 0) print '(A,10(I4,1X))','A = ', a print *,'# of positive values: ', count(pos) print *,'Sum of positive values: ', sum(a,pos) print *,'Minimum positive value: ', minval(a,pos) end program mask </source> ( from samples/arrays/mask.f90 ) <source lang="bash"> $ ./mask A = -5 -3 -1 1 3 5 7 9 11 13 # of positive values: 7 Sum of positive values: 49 Minimum positive value: 1 </source> |
Where construct
|
<source lang="fortran"> program wherearray implicit none real, dimension(6) :: a, diva integer :: i a = [(2*i-6, i=1,6)] where (a /= 0) diva = 1/a elsewhere diva = -999 endwhere print *,'a = ' print '(6(F8.3,1X))',a print *,'1/a = ' print '(6(F8.3,1X))',diva end program wherearray </source> (from samples/arrays/where.f90) <source lang="bash"> $ ./where a = -4.000 -2.000 0.000 2.000 4.000 6.000 1/a = -0.250 -0.500 -999.000 0.500 0.250 0.167 </source> |
Forall construct
|
<source lang="fortran"> program forallarray implicit none integer, dimension(6,6) :: a integer :: i,j a = -999 forall (i=1:6, j=1:6, i/=j) a(i,j) = i-j endforall do i=1,6 print '(6(I5,1X))',(a(i,j),j=1,6) enddo end program forallarray </source> (from samples/arrays/forall.f90) <source lang="bash"> $ ./forall -999 -1 -2 -3 -4 -5 1 -999 -1 -2 -3 -4 2 1 -999 -1 -2 -3 3 2 1 -999 -1 -2 4 3 2 1 -999 -1 5 4 3 2 1 -999 </source> |
Array Sections
|
<source lang="fortran"> 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] </source> |
Array Sections
|
<source lang="fortran"> program derivative implicit none real, dimension(10) :: x real, dimension(9) :: derivx integer :: i real, parameter:: pi = 4.*atan(1.), h=1. x = [(2*pi*(i-1)/9.,i=1,10)] derivx = ((x(2:10)-x(1:9))/h) print *, derivx do i=1,9 derivx(i) = (x(i+1)-x(i))/h enddo print *, derivx end program derivative </source> (from samples/arrays/derivative.f90) |
Array Sections
|
<source lang="fortran">
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] </source> |
Other important array intrinsics
<source lang="fortran">
1,4 reshape([1,2,3,4,5,6],[3,2]) == 2,5 3,6 </source> |
Linear algebra in Fortran
|
Matrix Multiplication Times
<source lang="fortran"> |
...
print *, 'Experiment with matrix size ', n print *, 'Times in seconds.' allocate(a(n,n)) allocate(b(n,n)) allocate(c(n,n)) call random_number(a) call random_number(b) call tick(starttime) do j=1,n do i=1,n c(i,j) = 0. do k=1,n c(i,j) = c(i,j)+a(i,k)*b(k,j) enddo enddo enddo looptime = tock(starttime) call tick(starttime) c = matmul(a,b) matmultime = tock(starttime) call tick(starttime) call sgemm('N','N',n,n,n,1.,a,n,b,n,0.,c,n) sgemmtime = tock(starttime)
print *, 'Triple-loop time: ', looptime print *, 'matmul intrinsic time: ', matmultime print *, 'SGEMM lapack call time:', sgemmtime |
...
</source> |
<source lang="bash"> $ ./matmul 2500 Experiment with matrix size 2500 Triple-loop time: 149.63400 matmul intrinsic time: 10.370000 SGEMM lapack call time: 1.4809999 </source> (gfortran 4.6, compiled -O3 -march=native using Intel MKL 10.3 for sgemm) (program from samples/arrays/matmul.f90) |
---|
Linear algebra in Fortran
<source lang="fortran"> program matvec implicit none integer, dimension(4,5) :: a integer, dimension(5,4) :: at integer, dimension(4,4) :: aat integer :: i a = reshape([(i,i=1,4*5)],[4,5]) at = transpose(a) print *,'A = ' call printmat(a) print *,'A^T = ' call printmat(at) aat = matmul(a,at) print *,'A . A^T = ' call printmat(aat) |
...
</source> |
<source lang="bash"> $ ./matrix A = 1 5 9 13 17 2 6 10 14 18 3 7 11 15 19 4 8 12 16 20 A^T = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 A . A^T = 565 610 655 700 610 660 710 760 655 710 765 820 700 760 820 880 A . A^T subarray = 610 660 710 655 710 765 700 760 820 </source> |
---|
Array sizes and Assumed Shape
integer nx, ny integer a(nx,ny) or worse, integer a(*,ny)
|
<source lang="fortran"> subroutine printmat(a) implicit none integer, dimension(:,:) :: a integer :: nr, nc, i, j nr = size(a,1) nc = size(a,2) do i=1,nr print '(99(I4,1X))', (a(i,j), j=1,nc) enddo end subroutine printmat </source> (from samples/arrays/matrix.f90) |
Allocatable Arrays
|
Allocate(), Deallocate()
|
<source lang="fortran"> program allocarray implicit none integer :: i, n integer, dimension(:), allocatable :: a n = 10 allocate(a(n)) a = [(i, i=2,20,2)] print *,'A = ' print *,a deallocate(a) end program allocarray </source> (from samples/arrays/allocatable.f90 ) |
get_command_argument()
|
<source lang="fortran"> program allocarray2 implicit none integer :: i, n integer, dimension(:), allocatable :: a character(len=30) :: arg if (command_argument_count() < 1) then print *,'Use: allocatable N, '//& ' where N is array size.' stop endif call get_command_argument(1, arg) read( arg,'(I30)'), n print *,'Allocating array of size ', n allocate(a(n)) a = [(i,i=1,n)] print *, a deallocate(a) end program allocarray2 </source> (from samples/arrays/allocatable2.f90) <source lang="bash"> $ ./allocatable2 Use: allocatable N, where N is array size. $ ./allocatable2 3 Allocating array of size 3 1 2 3 $ ./allocatable2 5 Allocating array of size 5 1 2 3 4 5 </source> |
Hands on #3
|
Fortran Pointers
|
<source lang="fortran"> real, target :: x = 3.2 real, pointer:: p p => x p x 3.2 </source> samples/pointers/ptr1.f90 |
Fortran Pointers
|
<source lang="fortran"> real, target :: x = 3.2 real, pointer:: p p => null() p x </source> |
Fortran Pointers
|
<source lang="fortran"> real, target :: x = 3.2 real, pointer:: p1, p2 p1 => x p2 => p1 p1 p2 x 3.2 </source> |
Allocating a pointer
|
<source lang="fortran"> real, pointer:: p allocate(p) p = 7.9 p 7.9 |
</source> samples/pointers/ptr2.f90 |
What are they good for? (1)
|
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)
|
<source lang="fortran"> real, target, dimension(7) :: x real, pointer:: p(:) p => x(2:6) p x 1 2 3 4 5 6 7 </source> samples/pointers/views.f90 |
Hands on #4
|
Derived Types and Objects
<source lang="fortran"> 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 </source> |
Derived Types and Objects
|
(from samples/derivedtypes/simple/intervalmath.f90 ) |
Procedures using types
|
(from samples/derivedtypes/intervalfunctions/intervalmath.f90 ) |
Procedures using types
|
samples/derivedtypes/genericintervals/interval2.f90 |
Generic Interfaces
Generic Interfaces
Generic interfaces
Operator overloading
Operator overloading
Generic interfaces
Type bound procedures
Type bound procedures
Type bound procedures
Object oriented programming
Interoperability with other languages
C-interoperability
Calling a C routine from Fortran
Calling a C routine from Fortran
Calling a C routine from Fortran
C strings
Calling Fortran from C
Calling Fortran from C
More advanced
Fortran calling C, which calls Fortran
samples/C-interop/valueref/croutine.csamples/C-interop/valueref/froutine.f90$ ./main
samples/C-interop/valueref/MakefileF2py
* Will only use solver module.
* generates the following header file (hydro_solver.pyf)
$ ipython
Coarrays
Coarrays
Coarrays
Sychronization
Sychronization
1d Diffusion Eqn
1d Diffusion Eqn
1d Diffusion Eqn
1d Diffusion Eqn
$ ./diffusion
Coarray Diffusion
Coarray Diffusion
Coarray Diffusion
Coarray Diffusion
Coarray Diffusion
$ export FOR_COARRAY_NUM_IMAGES=3
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',
Parallel Matrix Multiplication
Parallel Matrix Multiplication
Parallel Matrix Multiplication
Parallel Matrix Multiplication
Parallel Matrix Multiplication
Coarray Summary
Closing Hints
Useful Resources
|