Fortran code examples
From Wikipedia, the free encyclopedia
The following sample programs can be compiled and run with any standard Fortran compiler (see the end of the main Fortran article for lists of compilers). Most modern Fortran compilers expect a file with a .f or .for extension (for FORTRAN 66 or FORTRAN 77 source, although the FORTRAN 66 dialect may have to be selected specifically with a command-line option) or .f90/.f95 extension (for Fortran 90/95 source, respectively).
[edit] "Retro" FORTRAN IV
A retro example of a FORTRAN IV (later evolved into FORTRAN 66) program deck is available on the IBM 1130 page, including the IBM 1130 DM2 JCL required for compilation and execution. An IBM 1130 emulator is available at IBM 1130.org that will allow the FORTRAN IV program to be compiled and run on a PC.
The Snoopy Calendar program is the classic Fortran program often referenced in computer history folklore (e.g. Real Programmers Don't Use Pascal). This version can also be run using the IBM 1130 emulator above.
[edit] Hello, world
In keeping with computing tradition, the first example presented is a simple program to display the words "Hello, world" on the screen (or printer).
FORTRAN 66 (also FORTRAN IV)
C FORTRAN IV WAS ONE OF THE FIRST PROGRAMMING
C LANGUAGES TO SUPPORT SOURCE COMMENTS
WRITE (6,7)
7 FORMAT(13H HELLO, WORLD)
STOP
END
This program prints "HELLO, WORLD" to Fortran unit number 6, which on most machines was the line printer or terminal. (The card reader or keyboard was usually connected as unit 5). The number 7 in the WRITE statement refers to the statement number of the corresponding FORMAT statement. FORMAT statements may be placed anywhere in the same program or function/subroutine block as the WRITE statements which reference them. Typically a FORMAT statement is placed immediately following the WRITE statement which invokes it; alternatively, FORMAT statements are grouped together at the end of the program or subprogram block. If execution flows into a FORMAT statement, it is a no-op; thus, the example above has only two executable statements, WRITE and STOP.
The initial 13H in the FORMAT statement in the above example defines a Hollerith constant, here meaning that the 13 characters immediately following are to be taken as a character constant (note that the Hollerith constant is not surrounded by delimiters). (Some compilers also supported character literals enclosed in single quotes, a practice that came to be standard with FORTRAN 77.)
The space immediately following the 13H is a carriage control character, telling the I/O system to advance to a new line on the output. A zero in this position advances two lines (double space), a 1 advances to the top of a new page and + character will not advance to a new line, allowing overprinting.
FORTRAN 77
As of FORTRAN 77, single quotes are used to delimit character literals, and inline character strings may be used instead of references to FORMAT statements. Comment lines may be indicated with either a C or an asterisk (*) in column 1.
PROGRAM HELLOW
* The PRINT statement is like WRITE,
* but prints to the standard output unit
PRINT '(A)', 'Hello, world'
STOP
END
Fortran 90
As of Fortran 90, double quotes are allowed in addition to single quotes. An updated version of the Hello, world example (which here makes use of list-directed I/O, supported as of FORTRAN 77) could be written in Fortran 90 as follows:
program HelloWorld write (*,*) 'Hello, world!' ! This is an inline comment end program HelloWorld
[edit] Greatest common divisor (FORTRAN 77)
The following introductory example in FORTRAN 77 finds the greatest common divisor for two numbers A and B using a verbatim implementation of Euclid's algorithm.
* euclid.f (FORTRAN 77)
* Find greatest common divisor using the Euclidean algorithm
PROGRAM EUCLID
PRINT *, 'A?'
READ *, NA
IF (NA.LE.0) THEN
PRINT *, 'A must be a positive integer.'
STOP
END IF
PRINT *, 'B?'
READ *, NB
IF (NB.LE.0) THEN
PRINT *, 'B must be a positive integer.'
STOP
END IF
PRINT *, 'The GCD of', NA, ' and', NB, ' is', NGCD(NA, NB), '.'
STOP
END
FUNCTION NGCD(NA, NB)
IA = NA
IB = NB
1 IF (IB.NE.0) THEN
ITEMP = IA
IA = IB
IB = MOD(ITEMP, IB)
GOTO 1
END IF
NGCD = IA
RETURN
END
The above example is intended to illustrate the following:
- The
PRINTandREADstatements in the above use '*' as a format, specifying list-directed formatting. List-directed formatting instructs the compiler to make an educated guess about the required input or output format based on the following arguments. - As the earliest machines running Fortran had restricted character sets, FORTRAN 77 uses abbreviations such as
.EQ.,.NE.,.LT.,.GT.,.LE., and.GE.to represent the relational operators =, ≠, <, >, ≤, and ≥, respectively. - This example relies on the implicit typing mechanism described previously to specify the INTEGER types of
NA,NB,IA,IB, andITEMP. - In the function
NGCD(NA, NB), the values of the function argumentsNAandNBare copied into the local variablesIAandIBrespectively. This is necessary as the values ofIAandIBare altered within the function. Because argument passing in Fortran functions and subroutines utilize call by reference by default (rather than call by value, as is the default in languages such as C), modifyingNAandNBfrom within the function would effectively have modified the corresponding actual arguments in the mainPROGRAMunit which called the function.
The following shows the results of compiling and running the program.
mac:~/Desktop $ g77 -o euclid euclid.f mac:~/Desktop $ euclid A? 24 B? 36 The GCD of 24 and 36 is 12.
[edit] Complex numbers (FORTRAN 77)
The following FORTRAN 77 example prints out the values of ejiπ / 4 (where 
) for values of 
.
* cmplxd.f (FORTRAN 77)
* Demonstration of COMPLEX numbers
*
* Prints the values of e ** (j * i * pi / 4) for i = 0, 1, 2, ..., 7
* where j is the imaginary number sqrt(-1)
PROGRAM CMPLXD
IMPLICIT COMPLEX(X)
PARAMETER (PI = 3.141592653589793, XJ = (0, 1))
DO 1, I = 0, 7
X = EXP(XJ * I * PI / 4)
IF (AIMAG(X).LT.0) THEN
PRINT 2, 'e**(j*', I, '*pi/4) = ', REAL(X), ' - j',-AIMAG(X)
ELSE
PRINT 2, 'e**(j*', I, '*pi/4) = ', REAL(X), ' + j', AIMAG(X)
END IF
2 FORMAT (A, I1, A, F10.7, A, F9.7)
1 CONTINUE
STOP
END
The above example is intended to illustrate the following:
- The
IMPLICITstatement can be used to specify the implicit type of variables based on their initial letter if different from the default implicit typing scheme described above. In this example, this statement specifies that the implicit type of variables beginning with the letterXshall beCOMPLEX. - The
PARAMETERstatement may be used to specify constants. The second constant in this example (XJ) is given the complex-valued value 0 + j1, where j is the imaginary unit
. - The first number in the
DOstatement specifies the number of the last statement considered to be within the body of theDOloop. In this example, as neither theEND IFnor theFORMATis a single executable statement, theCONTINUEstatement (which does nothing) is used simply in order for there to be some statement to denote as the final statement of the loop. EXP()corresponds to the exponential function ex. In FORTRAN 77, this is a generic function, meaning that it accepts arguments of multiple types (such asREALand, in this example,COMPLEX). In FORTRAN 66, a specific function would have to be called by name depending on the type of the function arguments (for this example,CEXP()for aCOMPLEX-valued argument).- When applied to a
COMPLEX-valued argument,REAL()andAIMAG()return the values of the argument's real and imaginary components, respectively.
Incidentally, the output of the above program is as follows (see the article on Euler's formula for the geometric interpretation of these values as eight points spaced evenly about a unit circle in the complex plane).
mac:~/Desktop $ cmplxd e**(j*0*pi/4) = 1.0000000 + j0.0000000 e**(j*1*pi/4) = 0.7071068 + j0.7071068 e**(j*2*pi/4) = 0.0000000 + j1.0000000 e**(j*3*pi/4) = -0.7071068 + j0.7071068 e**(j*4*pi/4) = -1.0000000 - j0.0000001 e**(j*5*pi/4) = -0.7071066 - j0.7071069 e**(j*6*pi/4) = 0.0000000 - j1.0000000 e**(j*7*pi/4) = 0.7071070 - j0.7071065
Error can be seen occurring in the last decimal place in some of the numbers above, a result of the COMPLEX data type representing its real and imaginary components in single precision. Incidentally, Fortran 90 also made standard a double-precision complex-number data type (although several compilers provided such a type even earlier).
[edit] Summations with a DO loop (FORTRAN 77)
In this example of FORTRAN 77 code, the programmer has written the bulk of the code inside of a DO loop. Upon execution, instructions are printed to the screen and a SUM variable is initialized to zero outside the loop. Once the loop begins, it asks the user to input any number. This number is added to the variable SUM every time the loop repeats. If the user inputs 0, the EXIT statement terminates the loop, and the value of SUM is displayed on screen.
Also apparent in this program is a data file. Before the loop begins, the program creates (or opens, if it has already been run before) a text file called "SumData.DAT". During the loop, the WRITE statement stores any user-inputted number in this file, and upon termination of the loop, also saves the answer.
c sum.f
c Performs summations using in a loop using EXIT statement
c Saves inputted information and the summation in a data file
PROGRAM SUMMATION
IMPLICIT INTEGER (A-Z)
PRINT*,'This program performs summations. Enter 0 to stop.'
OPEN(UNIT=2, FILE='SumData.DAT')
SUM = 0
DO
PRINT*,'Add:'
READ*,A
IF(A.EQ.0) EXIT
IF(A.GT.0) THEN
SUM = SUM + A
ENDIF
WRITE(2,*)A
ENDDO
PRINT*,'Summation =', SUM
WRITE(2,*)'Summation =', SUM
CLOSE(2)
END
When executed, the console would display the following:
This program performs summations. Enter 0 to stop. Add: 1 Add: 2 Add: 3 Add: 0 Summation = 6
And the file SumData.DAT would contain:
1 2 3 Summation = 6
[edit] Calculating cylinder area (Fortran 90)
The following program, which calculates the surface area of a cylinder, illustrates free-form source input and other features introduced by Fortran 90.
program cylinder
! Calculate the surface area of a cylinder.
!
! Declare variables and constants.
! constants=pi
! variables=radius squared and height
implicit none ! Require all variables to be explicitly declared
integer :: ierr
character :: yn
real :: radius, height, area
real, parameter :: pi = 3.141592653589793
interactive_loop: do
! Prompt the user for radius and height
! and read them.
write (*,*) 'Enter radius and height.'
read (*,*,iostat=ierr) radius,height
! If radius and height could not be read from input,
! then cycle through the loop.
if (ierr /= 0) then
write(*,*) 'Error, invalid input.'
cycle interactive_loop
end if
! Compute area. The ** means "raise to a power."
area = 2 * pi * (radius**2 + radius*height)
! Write the input variables (radius, height)
! and output (area) to the screen.
write (*,'(1x,a7,f6.2,5x,a7,f6.2,5x,a5,f6.2)') &
'radius=',radius,'height=',height,'area=',area
yn = ' '
yn_loop: do
write(*,*) 'Perform another calculation? y[n]'
read(*,'(a1)') yn
if (yn=='y' .or. yn=='Y') exit yn_loop
if (yn=='n' .or. yn=='N' .or. yn==' ') exit interactive_loop
end do yn_loop
end do interactive_loop
end program cylinder
[edit] Dynamic memory allocation and arrays (Fortran 90)
The following program illustrates dynamic memory allocation and array-based operations, two features introduced with Fortran 90. Particularly noteworthy is the absence of DO loops and IF/THEN statements in manipulating the array; mathematical operations are applied to the array as a whole. Also apparent is the use of descriptive variable names and general code formatting that comport with contemporary programming style. This example computes an average over data entered interactively.
program average
! read in some numbers and take the average
! As written, if there are no data points (or no positive/negative points)
! an average of zero is returned.
! While this may not be desired behavior, it keeps this example simple.
implicit none
integer :: NumberOfPoints
real, dimension(:), allocatable :: Points
real :: AveragePoints=0., PositiveAverage=0., NegativeAverage=0.
write (*,*) 'Input number of points to average:'
read (*,*) NumberOfPoints
allocate (Points(NumberOfPoints))
write (*,*) 'Enter the Points to average:'
read (*,*) Points
! take the average by summing Points and dividing by NumberOfPoints
if (NumberOfPoints>0) AveragePoints = sum(Points)/NumberOfPoints
! now form average over positive and negative points only
if (count(Points>0.)>0) PositiveAverage = sum(Points, Points>0.)/count(Points>0.)
if (count(Points<0.)>0) NegativeAverage = sum(Points, Points<0.)/count(Points<0.)
deallocate (Points)
! print result to terminal
write (*,'(''Average = '', 1g12.4)') AveragePoints
write (*,'(''Average of positive points = '', 1g12.4)') PositiveAverage
write (*,'(''Average of negative points = '', 1g12.4)') NegativeAverage
end program average
[edit] Writing functions (Fortran 90)
Modern Fortran features available for use with procedures, including deferred-shape, protected, and optional arguments, are illustrated in the following example, a function to solve a system of linear equations.
function GaussSparse(NumIter, Tol, b, A, X, ActualIter)
! This function solves a system of equations (Ax = b) by using the Gauss-Jordan Method
implicit none
real GaussSparse
! Input: its value cannot be modified from within the function
integer, intent(in) :: NumIter
real, intent(in) :: Tol
real, intent(in), dimension(1:) :: b
real, intent(in), dimension(1:,1:) :: A
! Input/Output: its input value is used within the function, and can be modified
real, intent(inout), dimension(1:) :: X
! Output: its value is modified from within the function, only if the argument is required
integer, optional, intent(out) :: ActualIter
! Locals
integer i, n, Iter
real TolMax, Xk
! Initialize values
n = ubound(b, dim = 1) ! Size of array, obtained using ubound intrinsic function
TolMax = 2. * Tol
Iter = 0
! Compute solution until convergence
convergence_loop: do while (TolMax >= Tol.AND.Iter < NumIter); Iter = Iter + 1
TolMax = -1. ! Reset the tolerance value
! Compute solution for the k-th iteration
iteration_loop: do i = 1, n
! Compute the current x-value
Xk = (b(i) - sum(A(i,1:i-1) * X(1:i-1)) - sum(A(i,i+1:n) * X(i+1:n))) / A(i, i)
! Compute the error of the solution
TolMax = max((abs(X(i) - Xk)/(1. + abs(Xk))) ** 2, abs(A(i, i) * (X(i) - Xk)), TolMax)
X(i) = Xk
enddo iteration_loop
enddo convergence_loop
if (present(ActualIter)) ActualIter = Iter
GaussSparse = TolMax
end function GaussSparse
Note that an explicit interface to this routine must be available to its caller so that the type signature is known. This is preferably done by placing the function in a MODULE and then USEing the module in the calling routine. An alternative is to use an INTERFACE block.
[edit] Writing subroutines (Fortran 90)
In those cases where it is desired to return values via a procedure's arguments, a subroutine is preferred over a function; this is illustrated by the following subroutine to swap the contents of two arrays (this example assumes that the arrays are of equal size):
subroutine Swap_Real(a1, a2)
implicit none
! Input/Output
real, intent(inout) :: a1(:), a2(:)
! Locals
integer :: lb(1), & !Lower bound
ub(1) !Upper bound
integer i
real a
! Get bounds
lb = lbound(a1)
ub = ubound(a1)
! Swap
do i = lb(1), ub(1)
a = a1(i)
a1(i) = a2(i)
a2(i) = a
end do
end subroutine Swap_Real
As in the previous example, an explicit interface to this routine must be available to its caller so that the type signature is known. As before, this is preferably done by placing the function in a MODULE and then USEing the module in the calling routine. An alternative is to use a INTERFACE block.
[edit] Pointers and targets (Fortran 90)
In Fortran, the concept of pointers differs from that in C-like languages. A Fortran 90 pointer does not merely store the memory address of a target variable; it also contains additional descriptive information such as the target's rank, the upper and lower bounds of each dimension, and even strides through memory. This allows a Fortran 90 pointer to point at submatrices.
Fortran 90 pointers are "associated" with well-defined "target" variables, via either the pointer assignment operator (=>) or an ALLOCATE statement. When appearing in expressions, pointers are always dereferenced; no "pointer arithmetic" is possible.
The following example illustrates the concept:
program Test
! NOTE: Variable expressions in format statements (e.g., <m> or <n>)
! are compiler-dependent.
! Also, array notation (e.g., [1,2,3]) is a Fortran 2003 feature;
! for Fortran 95, use (/1,2,3/) instead.
use FunctionsModule, only: DoSomething => A !This function performs some
!operation on the integer
!input and returns its result
implicit none
integer, parameter :: m = 3, n = 3
integer, pointer :: p(:)=>null(), q(:,:)=>null()
integer, allocatable, target :: A(:,:)
integer ios = 0
allocate(A(1:m, 1:n), q(1:m, 1:n), stat = ios)
if (ios /= 0) stop 'Error during allocation of A and q'
! Assign the matrix
! A = [[1 4 7]
! [2 5 8]
! [3 6 9]]
A = reshape([(i, i = 1, m*n)], [m, n])
q = A
! p will be associated with the first column of A
p => A(:, 1)
! This operation on p has a direct effect on matrix A
p = p ** 2
! This will end the association between p and the first column of A
nullify(p)
! Matrix A becomes:
! A = [[1 4 7 ]
! [4 5 8 ]
! [9 6 9 ]]
write(*, '(/, "Matrix A becomes:",/,"A = [",<m>("[",<n>(i1,2x),"]",/,5x)"]")') &
((A(i, j), j = 1, n), i = 1, m)
! Perform some array operation
q = q + A
! Matrix q becomes:
! q = [[ 2 8 14 ]
! [ 6 10 16 ]
! [12 12 18 ]]
write(*, '(/, "Matrix q becomes:",/,"q = [",<m>("[",<n>(i2,2x),"]",/,5x)"]")') &
((q(i, j), j = 1, n), i = 1, m)
! Use p as an ordinary array
allocate (p(1:m*n), stat = ios)
if (ios /= 0) stop 'Error during allocation of p'
! Perform some array operation
p = [((DoSomething(A(i, j) + A(i, j) ** 2), i = 1, m), j = 1, n)]
write(*, '(<m*n>(i1,4x,"p[",i1,"] = ",i5))') (i, p(i), i = 1, m * n)
deallocate(A, p, q, stat = ios)
if (ios /= 0) stop 'Error during deallocation'
end program Test
[edit] Modules (Fortran 90)
A module is a program unit which contains data definitions, global data, and CONTAINed procedures. Unlike a simple INCLUDE file, a module is an independent program unit that can be compiled separately and linked in its binary form. Once compiled, a module's public contents can be made visible to a calling routine via the USE statement.
The module mechanism makes the explicit interface of procedures easily available to calling routines. In fact, modern Fortran encourages every SUBROUTINE and FUNCTION to be CONTAINed in a MODULE. This allows the programmer to use the newer argument passing options and allows the compiler to perform full type checking on the interface.
The following example also illustrates derived types and overloading of procedures and operators.
module GlobalModule
! Reference to a pair of procedures included in a previously compiled
! module named PortabilityLibrary
use PortabilityLibrary, only: GetLastError, & ! Generic procedure
Date ! Specific procedure
! Constants
integer, parameter :: dp_k = kind (1.0d0) ! Double precision kind
real, parameter :: zero = (0.)
real(dp_k), parameter :: pi = 3.141592653589793_dp_k
! Variables
integer :: n, m, retint
logical :: status, retlog
character(50) :: AppName
! Arrays
real, allocatable, dimension(:,:,:) :: a, b, c, d
complex(dp_k), allocatable, dimension(:) :: z
! Derived type definitions
type ijk
integer i
integer j
integer k
end type ijk
type matrix
integer m, n
real, allocatable :: a(:,:) ! Fortran 2003 feature. For Fortran 95,
! use the pointer attribute instead
end type matrix
! All the variables and procedures from this module can be accessed
! by other program units, except for AppName
public
private AppName
! Procedure overloading
interface swap
module procedure swap_Integer, swap_Real
end interface swap
interface GetLastError ! This adds a new, additional procedure to the
! generic procedure GetLastError
module procedure GetLastError_GlobalModule
end interface GetLastError
! Operator overloading
interface operator(+)
module procedure add_ijk
end interface
! Prototype for external procedure
interface
real function GaussSparse(NumIter, Tol, b, A, X)
integer, intent(in) :: NumIter
real, intent(in) :: Tol
real, intent(in), dimension(1:) :: b
real, intent(in), dimension(1:,1:) :: A
real, intent(inout), dimension(1:) :: X
end function GaussSparse
end interface
! Procedures included in the module
contains
! Internal function
function add_ijk(ijk_1, ijk_2)
type(ijk) add_ijk, ijk_1, ijk_2
intent(in) :: ijk_1, ijk_2
add_ijk = ijk(ijk_1%i + ijk_2%i, ijk_1%j + ijk_2%j, ijk_1%k + ijk_2%k)
end function add_ijk
! Include external files
include 'Swap_Integer.f90' ! Comments SHOULDN'T be added on include lines
include 'Swap_Real.f90'
end module GlobalModule
