ACM Computing Seminar Fortran Guide

Table of Contents

☰ ×

1 Introduction

This guide is intended to quickly get you up-and-running in scientific computing with Fortran.

1.1 About the language

Fortran was created in the 1950s for mathematical FOR-mula TRAN-slation, and has since gone through a number of revisions (FORTRAN 66, 77, and Fortran 90, 95, 2003, 2008, 2015). The language standards are put forth by the Fortran standards committee J3 in a document (ISO 1539-1:2010) available for purchase. The language syntax and intrinsic procedures make it especially suited for scientific computing. Fortran is a statically-typed and compiled language, like C++. You must declare the type, i.e. integer, real number, etc. of variables in programs you write. Your programs will be translated from human-readable source code into an executable file by software called a compiler. Fortran is not case-sensitive, so matrix and MaTrIx are translated to the same token by the compiler.

2 Getting started

The software that you need to get started comes prepackaged and ready to download on most Linux distributions. There are a few options for emulating a Linux environment in Windows or Mac OS, such as a virtual machine (VirtualBox) or package manager (MinGW or Cygwin on Windows and Brew on Mac OS).

2.1 Text editor

You will write the source code of your programs using a text editor. There are many options that have features designed for programming such as syntax highlighting and auto-completion. If you are an impossible-to-please perfectionist, you might want to check out Emacs. If you are easier to please, you might want to check out Sublime Text.

2.2 Compiler

To translate your source code into an executable, you will need a Fortran compiler. A free option is gfortran, part of the GNU compiler collection (gcc). The features of the Fortran language that are supported by the gfortran compiler are specified in the compiler manual. This is your most complete reference for the procedures intrinsic to Fortran that your programs can use. At the time of this writing, gfortran completely supports Fortran 95 and partially supports more recent standards.

2.3 Writing and compiling a program

A program is delimited by the begin program / end program keywords. A useful construct for keeping code that a program can use is called a module. A module is delimited by the begin module / end module keywords.

2.3.1 Hello world

Let's write a tiny program that prints "hello world" to the terminal screen in hello.f90. 

1: program main 2:  print*, 'hello world' 3: end program main

To compile the program, execute the following command on the command line in the same directory as hello.f90 

gfortran hello.f90

This produces an executable file named a.out by default (On Windows, this is probably named a.exe by default). To run, execute the file. 

./a.out

We could have specified a different name for the executable file during compilation with the -o option of gfortran. 

gfortran hello.f90 -o my_executable_file

On Windows, you should append the .exe extension to my_executable_file.

2.3.2 Template

Now let's write an empty source code template for future projects. Our source code template will consist of two files in the same directory (./source/). In the following files, the contents of a line after a ! symbol is a comment that is ignored by the compiler. One file header.f90 contains a module that defines things to be used in the main program. 

1: module header 2:  implicit none 3:  ! variable declarations and assignments 4: contains 5:  ! function and subroutine definitions 6: end module header

This file should be compiled with the -c option of gfortran. 

gfortran -c header.f90

This outputs the object file named header.o by default. An object file contains machine code that can be linked to an executable. A separate file main.f90 contains the main program. 

1: program main 2:  use header 3:  implicit none 4:  ! variable declarations and assignments 5:  ! function and subroutine calls 6: contains 7:  ! function and subroutine definitions 8: end program main

On line 2 of main.f90, we instruct the main program to use the contents of header.f90, so we must link header.o when compiling main.f90. 

gfortran main.f90 header.o -o main

To run the program, execute the output file main. 

./main

As you get more experience, you may find it cumbersome to repeatedly execute gfortran commands with every modification to your code. A way around this is to use the make command-line utility. Using make, all the of the compilation commands for your project can be coded in a file named makefile in the same directory as your .f90 source files. For example, the template above could use the following makefile. 

 1: COMPILER = gfortran  2: SOURCE = main.f90  3: EXECUTABLE = main  4: OBJECTS = header.o  5:   6: all: $(EXECUTABLE)  7: $(EXECUTABLE): $(OBJECTS)  8: $(COMPILER) $(SOURCE) $(OBJECTS) -o $(EXECUTABLE)  9: %.o: %.f90 10: $(COMPILER) -c $<

Then, to recompile both header.f90 and main.f90 after modifying either file, execute 

make

in the same directory as makefile. The first four lines of the makefile above define the compiler command, file name of the main program, file name of the executable to be created, and file name(s) of linked object file(s), respectively. If you wrote a second module in a separate file my_second_header.f90 that you wanted to use in main.f90, you would modify line 4 of makefile to OBJ = header.o my_second_header.o. The remaining lines of makefile define instructions for compilation.

2.4 Exercises

  1. Compile and run hello.f90.
  2. Execute man gfortran in any directory to bring up the manual for gfortran. Read the description and skim through the options. Do the same for make.

3 Data types

In both programs and modules, variables are declared first before other procedures. A variable is declared by listing its data type followed by :: and the variable name, i.e. integer :: i or real :: x.

We will use the implicit none keyword at the beginning of each program and module as in line 2 of header.f90 and line 3 of main.f90 in Section 2.3.2. The role of this keyword is to suppress implicit rules for interpreting undeclared variables. By including it, we force ourselves to declare each variable we use, which should facilitate debugging when our program fails to compile. Without it, an undeclared variable with a name such as i is assumed to be of the integer data type whereas an undeclared variable with a name such as x is assumed to be of the real data type.

In addition to the most common data types presented below, Fortran has a complex data type and support for data types defined by the programmer (see Section 7.1).

3.1 The logical type

A logical data type can have values .true. or .false.. Logical expressions can be expressed by combining unary or binary operations. 

 1: logical :: a,b,c  2: a = .true.  3: b = .false.  4:   5: ! '.not.' is the logical negation operator  6: c = .not.a ! c is false  7:   8: ! '.and,' is the logical and operator  9: c = a.and.b ! c is false 10:  11: ! '.or.' is the logical or operator 12: c = a.or.b ! c is true 13:  14: ! '==' is the test for equality 15: c = 1 == 2 ! c is false 16:  17: ! '/=' is test for inequality 18: c = 1 /= 2 ! c is true 19: print*, c

Other logical operators include

  • < or .lt.: less than
  • <= or .le.: less than or equal
  • > or .gt.: greater than
  • >= or .ge.: greater than or equal

Logical expressions are often used in control structures.

3.2 The integer type

An integer data type can have integer values. If a real value is assigned to an integer type, the decimal portion is chopped off. 

 1: integer :: a = 6, b = 7 ! initialize a and b to 6 and 7, resp  2: integer :: c  3:   4: c = a + b ! c is 13  5: c = a - b ! c is -1  6: c = a / b ! c is 0  7: c = b / a ! c is 1  8: c = a*b ! c is 42  9: c = a**b ! c is 6^7 10: c = mod(b,a) ! c is (b mod a) = 1 11: c = a > b ! c is 0 (logical gets cast to integer) 12: c = a < b ! c is 1 (logical gets cast to integer)

3.3 Floating point types

The two floating point data types real and double precision correspond to IEEE 32- and 64-bit floating point data types. A constant called machine epsilon is the least positive number in a floating point system that when added to 1 results in a floating point number larger than 1. It is common in numerical analysis error estimates. 

 1: real :: a ! declare a single precision float  2: double precision :: b ! declare a double precision float  3:   4: ! Print the min/max value and machine epsilon  5: ! for the single precision floating point system  6: print*, tiny(a), huge(a), epsilon(a)  7:   8: ! Print the min/max value and machine epsilon  9: ! for the double precision floating point system 10: print*, tiny(b), huge(b), epsilon(b) 
 1.17549435E-38 3.40282347E+38 1.19209290E-07 2.2250738585072014E-308 1.7976931348623157E+308 2.2204460492503131E-016

3.4 The character type

A character data type can have character values, i.e. letters or symbols. A character string is declared with a positive integer specifying its maximum possible length. 

1: ! declare a character variable s at most 32 characters 2: character(32) :: s 3:  4: ! assign value to s 5: s = 'file_name' 6:  7: ! trim trailing spaces from s and 8: ! append a character literal '.txt' 9: print*, trim(s) // '.txt' 
 file_name.txt

3.5 Casting

An integer can be cast to a real and vice versa. 

 1: integer :: a = 1, b  2: real :: c, PI = 3.14159  3:   4: ! explicit cast real to integer  5: b = int(PI) ! b is 3  6:   7: ! explicit cast integer to real then divide  8: c = a/real(b) ! c is .3333...  9:  10: ! divide then implicit cast real to integer 11: c = a/b ! c is 0

3.6 The parameter keyword

The parameter keyword is used to declare constants. A constant must be assigned a value at declaration and cannot be reassigned a value. The following code is not valid because of an attempt to reassign a constant. 

1: ! declare constant variable 2: real, parameter :: PI = 2.*asin(1.) ! 'asin' is arcsine 3:  4: PI = 3 ! not valid

The compiler produces an error like Error: Named constant ‘pi’ in variable definition context (assignment).

3.7 Setting the precision

The kind function returns an integer for each data type. The precision of a floating point number can be specified at declaration by a literal or constant integer of the desired kind. 

 1: ! declare a single precision  2: real :: r  3: ! declare a double precision  4: double precision :: d  5: ! store single precision and double precision kinds  6: integer, parameter :: sp = kind(r), dp = kind(d)  7: ! set current kind  8: integer, parameter :: rp = sp  9:  10: ! declare real b in double precision 11: real(dp) :: b 12:  13: ! declare real a with precision kind rp 14: real(rp) :: a 15:  16: ! cast 1 to real with precision kind rp and assign to a 17: a = 1.0_rp 18:  19: ! cast b to real with precision kind rp and assign to a 20: a = real(b,rp)

To switch the precision of each variable above with kind rp, we would only need to modify the declaration of rp on line 8.

3.8 Pointers

Pointers have the same meaning in Fortran as in C++. A pointer is a variable that holds the memory address of a variable. The implementation of pointers is qualitatively different in Fortran than in C++. In Fortran, the user cannot view the memory address that a pointer stores. A pointer variable is declared with the pointer modifier, and a variable that it points to is declared with the target modifier. The types of a pointer and its target must match. 

 1: ! declare pointer  2: integer, pointer :: p  3: ! declare targets  4: integer, target :: a = 1, b = 2  5:   6: p => a ! p has same memory address as a  7: p = 2 ! modify value at address  8: print*, a==2 ! a is 2  9:  10: p => b ! p has same memory address as b 11: p = 1 ! modify value at address 12: print*, b==1 ! b is 1 13:  14: ! is p associated with a target? 15: print*, associated(p) 16:  17: ! is p associated with the target a? 18: print*, associated(p, a) 19:  20: ! point to nowhere 21: nullify(p) 
 T T T F

3.9 Arrays

The length of an array can be fixed or dynamic. The index of an array starts at 1 by default, but any index range can be specified.

3.9.1 Fixed-length arrays

An array can be declared with a single integer specifying its length in which cast the first index of the array is 1. An array can also be declared with an integer range specifying its first and last index.

Here's a one-dimensional array example. 

 1: ! declare arrray of length 5  2: ! index range is 1 to 5 (inclusive)  3: real :: a(5)  4:   5: ! you can work with each component individually  6: ! set the first component to 1  7: a(1) = 1.0  8:   9: ! or you can work with the whole array 10: ! set the whole array to 2 11: a = 2.0 12:  13: ! or you can with slices of the array 14: ! set elements 2 to 4 (inclusive) to 3 15: a(2:4) = 3.0

And, here's a two-dimensional array example. 

 1: ! declare 5x5 array  2: ! index range is 1 to 5 (inclusive) in both axes  3: real :: a(5,5)  4:   5: ! you can work with each component individually  6: ! set upper left component to 1  7: a(1,1) = 1.0  8:   9: ! or you can work with the whole array 10: ! set the whole array to 2 11: a = 2.0 12:  13: ! or you can with slices of the array 14: ! set a submatrix to 3 15: a(2:4, 1:2) = 3.0

Fortran includes intrinsic functions to operate on an array a such as

  • size(a): number of elements of a
  • minval(a): minimum value of a
  • maxval(a): maximum value of a
  • sum(a): sum of elements in a
  • product(a): product of elements in a

See the gfortran documentation for more.

3.9.2 Dynamic length arrays

Dynamic arrays are declared with the allocatable modifier. Before storing values in such an array, you must allocate memory for the array. After you are finished the array, you ought to deallocate the memory that it occupies.

Here's a one-dimensional array example. 

 1: ! declare a one-dim. dynamic length array  2: real, allocatable :: a(:)  3:   4: ! allocate memory for a  5: allocate(a(5))  6:   7: ! now you can treat a like a normal array  8: a(1) = 1.0  9: ! etc... 10:  11: ! deallocate memory occupied by a 12: deallocate(a) 13:  14: ! we can change the size and index range of a 15: allocate(a(0:10)) 16:  17: a(0) = 1.0 18: ! etc... 19:  20: deallocate(a)

Without the last dellaocate statement on line 20 the code above is valid, but the memory that is allocated for a will not be freed. That memory then cannot be allocated to other resources.

Here's a two-dimensional array example. 

 1: ! declare a two-dim. dynamic length array  2: real, allocatable :: a(:,:)  3:   4: ! allocate memory for a  5: allocate(a(5,5))  6:   7: ! now you can treat a like a normal array  8: a(1,1) = 1.0  9: ! etc... 10:  11: ! deallocate memory occupied by a 12: deallocate(a) 13:  14: ! we can change the size and index range of a 15: allocate(a(0:10,0:10)) 16:  17: a(0,0) = 1.0 18: ! etc... 19:  20: deallocate(a)

4 Control structures

Control structures are used to direct the flow of code execution.

4.1 Conditionals

4.1.1 The if construct

The if construct controls execution of a single block of code. If the block of code is more than one line, it should be delimited by an if / end if pair. If the block of code is one line, it can be written on one line. A common typo is to forget the then keyword following the logical in an if / end if pair. 

1: real :: num = 0.75 2:  3: if (num < .5) then 4:  print*, 'num: ', num 5:  print*, 'num is less than 0.5' 6: end if 7:  8: if (num > .5) print*, 'num is greater than 0.5' 
 num is greater than 0.5

4.1.2 Example: if / else and random number generation

The if / else construct controls with mutually exclusive logic the execution of two blocks of code.

The following code generates a random number between 0 and 1, then prints the number and whether or not the number is greater than 0.5 

 1: real :: num  2:   3: ! seed random number generator  4: call srand(789)  5:   6: ! rand() returns a random number between 0 and 1  7: num = rand()  8:   9: print*, 'num: ', num 10:  11: if (num < 0.5) then 12:  print*, 'num is less than 0.5' 13: else 14:  print*, 'num is greater then 0.5' 15: end if 16:  17: ! do it again 18: num = rand() 19:  20: print*, 'num: ', num 21:  22: if (num < 0.5) then 23:  print*, 'num is less than 0.5' 24: else 25:  print*, 'num is greater then 0.5' 26: end if 
 num: 6.17480278E-03 num is less than 0.5 num: 0.783314705 num is greater then 0.5

Since the random number generator was seeded with a literal integer, the above code will produce the same output each time it is run.

4.1.3 Example: if / else if / else

The if / else if / else construct controls with mutually exclusive logic the execution of three or more blocks of code. The following code generates a random number between 0 and 1, then prints the number and which quarter of the interval \([0,1]\) that the number is in. 

 1: real :: num  2:   3: ! seed random number generator with current time  4: call srand(time())  5:   6: ! rand() returns a random number between 0 and 1  7: num = rand()  8:   9: print*, 'num:', num 10:  11: if (num > 0.75) then 12:  print*, 'num is between 0.75 and 1' 13: else if (num > 0.5) then 14:  print*, 'num is between 0.5 and 0.75' 15: else if (num > 0.25) then 16:  print*, 'num is between 0.25 and 0.5' 17: else 18:  print*, 'num is between 0 and 0.25' 19: end if 
 num: 0.570252180 num is between 0.5 and 0.75

Since the random number generator was seeded with the current time, the above code will produce a different output each time it is run.

4.2 Loops

4.2.1 The do loop

A do loop iterates a block of code over a range of integers. It takes two integer arguments specifying the minimum and maximum (inclusive) of the range and takes an optional third integer argument specifying the iteration stride in the form do i=min,max,stride. If omitted, the stride is 1.

The following code assigns a value to each component of an array then prints it. 

 1: integer :: max = 10, i  2: real, allocatable :: x(:)  3:   4: allocate(x(0:max))  5:   6: do i = 0,max  7:  ! assign to each array component  8:  x(i) = i / real(max)  9:  10:  ! print current component 11:  print "('x(', i0, ') = ', f3.1)", i, x(i) 12: end do 13:  14: deallocate(x) 
 x(0) = 0.0 x(1) = 0.1 x(2) = 0.2 x(3) = 0.3 x(4) = 0.4 x(5) = 0.5 x(6) = 0.6 x(7) = 0.7 x(8) = 0.8 x(9) = 0.9 x(10) = 1.0

An implicit do loop can be used for formulaic array assignments. The following code creates the same array as the last example. 

1: integer :: max = 10 2: real, allocatable :: x(:) 3:  4: allocate(x(0:max)) 5:  6: ! implicit do loop for formulaic array assignment 7: x = [(i / real(max), i=0, max)] 8:  9: deallocate(x)
4.2.1.1 Example: row-major matrix

The following code stores matrix data in a one-dimensional array named matrix in row-major order. This means the first n_cols elements of the array will contain the first row of the matrix, the next n_cols of the array will contain the second row of the matrix, etc. 

 1: integer :: n_rows = 4, n_cols = 3  2: real, allocatable :: matrix(:)  3: ! temporary indices  4: integer :: i,j,k  5:   6: ! index range is 1 to 12 (inclusive)  7: allocate(matrix(1:n_rows*n_cols))  8:   9: ! assign 0 to all elements of matrix 10: matrix = 0.0 11:  12: do i = 1,n_rows 13:  do j = 1,n_cols 14:  ! convert (i,j) matrix index to "flat" row-major index 15:  k = (i-1)*n_cols + j 16:  17:  ! assign 1 to diagonal, 2 to sub/super-diagonal 18:  if (i==j) then 19:  matrix(k) = 1.0 20:  else if ((i==j-1).or.(i==j+1)) then 21:  matrix(k) = 2.0 22:  end if 23:  end do 24: end do 25:  26: ! print matrix row by row 27: do i = 1,n_rows 28:  print "(3(f5.1))", matrix(1+(i-1)*n_cols:i*n_cols) 29: end do 30:  31: deallocate(matrix) 
 1.0 2.0 0.0 2.0 1.0 2.0 0.0 2.0 1.0 0.0 0.0 2.0

4.2.2 The do while loop

A do while loop iterates while a logical condition evaluates to .true..

4.2.2.1 Example: truncated sum

The following code approximates the geometric series

\begin{equation*} \sum_{n=1}^{\infty}\left(\frac12\right)^n=1. \end{equation*}

The do while loop begins with \(n=1\) and exits when the current summand does not increase the current sum. It prints the iteration number, current sum, and absolute error

\begin{equation*} E=1-\sum_{n=1}^{\infty}\left(\frac12\right)^n. \end{equation*} 
 1: real :: sum = 0.0, base = 0.5, tol = 1e-4  2: real :: pow = 0.5  3: integer :: iter = 1  4:   5: do while (sum+pow > sum)  6:  ! add pow to sum  7:  sum = sum+pow  8:  ! update pow by one power of base  9:  pow = pow*base 10:  11:  print "('Iter: ', i3, ', Sum: ', f0.10, ', Abs Err: ', f0.10)", iter, sum, 1-sum 12:  13:  ! update iter by 1 14:  iter = iter+1 15: end do 
 Iter: 1, Sum: .5000000000, Abs Err: .5000000000 Iter: 2, Sum: .7500000000, Abs Err: .2500000000 Iter: 3, Sum: .8750000000, Abs Err: .1250000000 Iter: 4, Sum: .9375000000, Abs Err: .0625000000 Iter: 5, Sum: .9687500000, Abs Err: .0312500000 Iter: 6, Sum: .9843750000, Abs Err: .0156250000 Iter: 7, Sum: .9921875000, Abs Err: .0078125000 Iter: 8, Sum: .9960937500, Abs Err: .0039062500 Iter: 9, Sum: .9980468750, Abs Err: .0019531250 Iter: 10, Sum: .9990234375, Abs Err: .0009765625 Iter: 11, Sum: .9995117188, Abs Err: .0004882812 Iter: 12, Sum: .9997558594, Abs Err: .0002441406 Iter: 13, Sum: .9998779297, Abs Err: .0001220703 Iter: 14, Sum: .9999389648, Abs Err: .0000610352 Iter: 15, Sum: .9999694824, Abs Err: .0000305176 Iter: 16, Sum: .9999847412, Abs Err: .0000152588 Iter: 17, Sum: .9999923706, Abs Err: .0000076294 Iter: 18, Sum: .9999961853, Abs Err: .0000038147 Iter: 19, Sum: .9999980927, Abs Err: .0000019073 Iter: 20, Sum: .9999990463, Abs Err: .0000009537 Iter: 21, Sum: .9999995232, Abs Err: .0000004768 Iter: 22, Sum: .9999997616, Abs Err: .0000002384 Iter: 23, Sum: .9999998808, Abs Err: .0000001192 Iter: 24, Sum: .9999999404, Abs Err: .0000000596 Iter: 25, Sum: 1.0000000000, Abs Err: .0000000000
4.2.2.2 Example: estimating machine epsilon

The following code finds machine epsilon by shifting the rightmost bit of a binary number rightward until it falls off. Think about how it does this. Could you write an algorithm that finds machine epsilon using the function rshift that shifts the bits of float rightward? 

 1: double precision :: eps  2: integer, parameter :: dp = kind(eps)  3: integer :: count = 1  4:   5: eps = 1.0_dp  6: do while (1.0_dp + eps*0.5 > 1.0_dp)  7:  eps = eps*0.5  8:  count = count+1  9: end do 10:  11: print*, eps, epsilon(eps) 12: print*, count, digits(eps) 
 2.2204460492503131E-016 2.2204460492503131E-016 53 53

4.2.3 Example: the exit keyword

The exit keyword stops execution of code within the current scope.

The following code finds the hailstone sequence of \(a_1=6\) defined recursively by

\begin{equation*} a_{n+1} = \begin{cases} a_n/2 & \text{if } a_n \text{ is even}\\ 3a_n+1 & \text{ if } a_n \text{ is odd} \end{cases} \end{equation*}

for \(n\geq1\). It is an open conjecture that the hailstone sequence of any initial value \(a_1\) converges to the periodic sequence \(4, 2, 1, 4, 2, 1\ldots\). Luckily, it does for \(a_1=6\) and the following infinite do loop exits. 

 1: integer :: a = 6, count = 1  2:   3: ! infinite loop  4: do  5:  ! if a is even, divide by 2  6:  ! otherwise multiply by 3 and add 1  7:  if (mod(a,2)==0) then  8:  a = a/2  9:  else 10:  a = 3*a+1 11:  end if 12:  13:  ! if a is 4, exit infinite loop 14:  if (a==4) then 15:  exit 16:  end if 17:  18:  ! print count and a 19:  print "('count: ', i2, ', a: ', i2)", count, a 20:  21:  ! increment count 22:  count = count + 1 23: end do 
 count: 1, a: 3 count: 2, a: 10 count: 3, a: 5 count: 4, a: 16 count: 5, a: 8

5 Input/Output

5.1 File input/output

5.1.1 Reading data from file

The contents of a data file can be read into an array using read. Suppose you have a file ./data/array.txt that contains two columns of data 

 1 1.23 2 2.34 3 3.45 4 4.56 5 5.67

This file can be opened with the open command. The required first argument of open is an integer that specifies a file unit for array.txt. Choose any number that is not in use. The unit numbers 0, 5, and 6 are reserved for system files and should not be used accidentally. Data are read in row-major format, i.e. across the first row, then across the second row, etc.

The following code reads the contents of ./data/array.txt into an array called array. 

 1: ! declare array  2: real :: array(5,2)  3: integer :: row  4:   5: ! open file and assign file unit 10  6: open (10, file='./data/array.txt', action='read')  7:   8: ! read data from file unit 10 into array  9: do row = 1,5 10:  read(10,*) array(row,:) 11: end do 12:  13: ! close file 14: close(10)

5.1.2 Writing data to file

Data can be written to a file with the write command. 

 1: real :: x  2: integer :: i, max = 5  3:   4: ! open file, specify unit 10, overwrite if exists  5: open(10, file='./data/sine.txt', action='write', status='replace')  6:   7: do i = 0,max  8:  x = i / real(max)  9:  10:  ! write to file unit 10 11:  write(10,*) x, sin(x) 12: end do

This produces a file sine.txt in the directory data containing 

 0.00000000 0.00000000 0.200000003 0.198669329 0.400000006 0.389418334 0.600000024 0.564642489 0.800000012 0.717356086 1.00000000 0.841470957

5.2 Formatted input/output

The format of a print, write, or read statement can be specified with a character string. A format character string replaces the * symbol in print* and the second * symbol in read(*,*) or write(*,*). A format string is a list of literal character strings or character descriptors from

  • a: character string
  • iW: integer
  • fW.D: float point
  • esW.DeE: scientific notation
  • Wx: space

where W, D, and E should be replaced by numbers specifying width, number of digits, or number of exponent digits, resp. The width of a formatted integer or float defaults to the width of the number when W is 0. 

 1: character(32) :: fmt, a = 'word'  2: integer :: b = 1  3: real :: c = 2.0, d = 3.0  4:   5: ! character string and 4 space-delimited values  6: print "('four values: ', a, 1x i0, 1x f0.1, 1x, es6.1e1)", trim(a), b, c, d  7:   8: ! character string and 2 space-delimited values  9: fmt = '(a, 2(f0.1, 1x))' 10: print fmt, 'two values: ', c, d 
 four values: word 1 2.0 3.0E+0 two values: 2.0 3.0

5.3 Command line arguments

Arguments can be passed to a program from the command line using get_command_argument. The first argument received by get_command_argument is the program executable file name and the remaining arguments are passed by the user. The following program accepts any number of arguments, each at most 32 characters, and prints them. 

 1: program main  2:  implicit none  3:   4:  character(32) :: arg  5:  integer :: n_arg = 0  6:   7:  do  8:  ! get next command line argument  9:  call get_command_argument(n_arg, arg) 10:  11:  ! if it is empty, exit 12:  if (len_trim(arg) == 0) exit 13:  14:  ! print argument to screen 15:  print"('argument ', i0, ': ', a)", n_arg, trim(arg) 16:  17:  ! increment count 18:  n_arg = n_arg+1 19:  end do 20:  21:  ! print total number of arguments 22:  print "('number of arguments: ', i0)", n_arg 23:  24: end program main

After compiling to a.out, you can pass arguments in the executing command. 

./a.out 1 2 34 
 argument 0: ./a.out argument 1: 1 argument 2: 2 argument 3: 34 number of arguments: 4

6 Functions/Subroutines

Functions and subroutines are callable blocks of code. A function returns a value from a set of arguments. A subroutine executes a block of code from a set of arguments but does not explicitly return a value. Changes to arguments made within a function are not returned whereas changes to arguments made within a subroutine can be returned to the calling program. Both functions and subroutines are defined after the contains keyword in a module or program.

6.1 Writing a function

The definition of a function starts with the name of the function followed by a list of arguments and return variable. The data types of the arguments and return variable are defined within the function body.

6.1.1 Example: linspace: generating a set of equally-space points

The following program defines a function linspace that returns a set of equidistant points on an interval. The main function makes a call to the function. 

 1: program main  2:  implicit none  3:   4:  real :: xs(10)  5:   6:  ! call function linspace to set values in xs  7:  xs = linspace(0.0, 1.0, 10)  8:   9:  ! print returned value of xs 10:  print "(10(f0.1, 1x))" , xs 11:  12: contains 13:  14:  ! linspace: return a set of equidistant points on an interval 15:  ! min: minimum value of interval 16:  ! max: maximum value of interval 17:  ! n_points: number of points in returned set 18:  ! xs: set of points 19:  function linspace(min, max, n_points) result(xs) 20:  real :: min, max, dx 21:  integer :: n_points 22:  integer :: i 23:  real :: xs(n_points) 24:  25:  ! calculate width of subintervals 26:  dx = (max-min) / real(n_points-1) 27:  28:  ! fill xs with points 29:  do i = 1,n_points 30:  xs(i) = min + (i-1)*dx 31:  end do 32:  33:  end function linspace 34:  35: end program main 
 .0 .1 .2 .3 .4 .6 .7 .8 .9 1.0

6.2 Writing a subroutine

The definition of a subroutine begins with the name of the subroutine and list of arguments. Arguments are defined within the subroutine body with one of the following intents

  • intent(in): changes to the argument are not returned
  • intent(inout): changes to the argument are returned
  • intent(out): the initial value of the argument is ignored and changes to the argument are returned.

Subroutines are called using the call keyword followed by the subroutine name.

6.2.1 Example: polar coordinates

The following code defines a subroutine polar_coord that returns the polar coordinates \((r,\theta)\) defined by \(r=\sqrt{x^2+y^2}\) and \(\theta=\arctan(y/x)\) from the rectangular coordinate pair \((x,y)\). 

 1: program main  2:   3:  real :: x = 1.0, y = 1.0, rad, theta  4:   5:  ! call subroutine that returns polar coords  6:  call polar_coord(x, y, rad, theta)  7:  print*, rad, theta  8:   9: contains 10:  11:  ! polar_coord: return the polar coordinates of a rect coord pair 12:  ! x,y: rectangular coord 13:  ! rad,theta: polar coord 14:  subroutine polar_coord(x, y, rad, theta) 15:  real, intent(in) :: x, y 16:  real, intent(out) :: rad, theta 17:  18:  ! compute polar coord 19:  ! hypot = sqrt(x**2+y**2) is an intrinsic function 20:  ! atan2 = arctan with correct sign is an intrinsic function 21:  rad = hypot(x, y) 22:  theta = atan2(y, x) 23:  24:  end subroutine polar_coord 25:  26: end program main 
 1.41421354 0.785398185

6.3 Passing procedures as arguments

An inteface can be used to pass a function or subroutine to another function or a subroutine. For this purpose, an interface is defined in the receiving procedure essentially the same way as the passed procedure itself but with only declarations and not the implementation.

6.3.1 Example: Newton's method for rootfinding

Newton's method for finding the root of a function \(f:\mathbb{R}\rightarrow\mathbb{R}\) refines an initial guess \(x_0\) according to the iteration rule

\begin{equation*} x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} \end{equation*}

for \(n\geq1\) until \(f(x)\) is less than a chosen tolerance or a maximum number of iterations.

The following code defines a subroutine newton_root that returns a root of an input function as well as the number of iterations of Newton's method used to find the root. It is called by the main program to approximate the positive root of \(f(x)=x^2-2\) from an initial guess \(x_0=1\). 

 1: program main  2:  implicit none  3:   4:  character(64) :: fmt  5:  real :: x = 1.0  6:  integer :: iter = 1000  7:   8:  ! call newton rootfinding function  9:  call newton_root(f, df, x, iter, 1e-6, .true.) 10:  11:  ! print found root and number of iterations used 12:  fmt = "('number of iterations: ', i0, ', x: ', f0.7, ', f(x): ', f0.7)" 13:  print fmt, iter, x, f(x) 14:  15: contains 16:  17:  ! function f(x) = x^2 - 2 18:  function f(x) result(y) 19:  real :: x, y 20:  y = x*x - 2 21:  end function f 22:  23:  ! function df(x) = 2x 24:  function df(x) result(dy) 25:  real :: x, dy 26:  dy = 2*x 27:  end function df 28:  29:  ! newton_root: newtons method for rootfinding 30:  ! f: function with root 31:  ! df: derivative of f 32:  ! x: sequence iterate 33:  ! iter: max number of iterations at call, number of iterations at return 34:  ! tol: absolute tolerance 35:  ! print_iters: boolean to toggle verbosity 36:  subroutine newton_root(f, df, x, iter, tol, print_iters) 37:  38:  ! interface to function f 39:  interface 40:  function f(x) result(y) 41:  real :: x, y 42:  end function f 43:  end interface 44:  45:  ! interface to function df 46:  interface 47:  function df(x) result(dy) 48:  real :: x, dy 49:  end function df 50:  end interface 51:  52:  real, intent(inout) :: x 53:  real, intent(in) :: tol 54:  integer, intent(inout) :: iter 55:  logical, intent(in) :: print_iters 56:  integer :: max_iters 57:  58:  max_iters = iter 59:  iter = 0 60:  61:  ! while f(x) greater than absolute tolerance 62:  ! and max number of iterations not exceeded 63:  do while (abs(f(x))>tol.and.iter<max_iters) 64:  ! print current x and f(x) 65:  if (print_iters) print "('f(', f0.7, ') = ', f0.7)", x, f(x) 66:  67:  ! Newton's update rule 68:  x = x - f(x)/df(x) 69:  70:  ! increment number of iterations 71:  iter = iter + 1 72:  end do 73:  74:  end subroutine newton_root 75:  76: end program main 
 f(1.0000000) = -1.0000000 f(1.5000000) = .2500000 f(1.4166666) = .0069444 f(1.4142157) = .0000060 number of iterations: 4, x: 1.4142135, f(x): -.0000001

6.3.2 Example: The midpoint rule for definite integrals

The midpoint rule approximates the definite integral \(\int_a^bf(x)~dx\) with integrand \(f:\mathbb{R}\rightarrow\mathbb{R}\) by

\begin{equation} \Delta x\sum_{i=1}^nf(\bar{x}_i) \end{equation}

where \(\Delta x=(b-a)/n\), \(x_i=a+i\Delta x\) and \(\bar{x}_i=(x_{i-1}+x_i)/2\).

The following code defines a function midpoint that computes the approximation eq. 1 given \(a\), \(b\), and \(n\). The main program calls midpoint to approximate the definite integral of \(f(x)=1/x\) on \([1,e]\) for a range of \(n\). 

 1: program main  2:  implicit none  3:   4:  real, parameter :: E = exp(1.)  5:  integer :: n  6:  real :: integral  7:   8:  ! Approximate the integral of 1/x from 1 to e  9:  ! with the midpoint rule for a range of number of subintervals 10:  do n = 2,20,2 11:  print "('n: ', i0, ', M_n: ', f0.6)", n, midpoint(f, 1.0, E, n) 12:  end do 13:  14: contains 15:  16:  ! function f(x) = 1/x 17:  function f(x) result(y) 18:  real :: x, y 19:  y = 1.0/x 20:  end function f 21:  22:  ! midpoint: midpoint rule for definite integral 23:  ! f: integrand 24:  ! a: left endpoint of interval of integration 25:  ! b: right endpoint of interval of integration 26:  ! n: number of subintervals 27:  ! sum: approximate definite integral 28:  function midpoint(f, a, b, n) result(sum) 29:  30:  ! interface to f 31:  interface 32:  function f(x) 33:  real :: x, y 34:  end function f 35:  end interface 36:  37:  real :: a, b, min, xi, dx, sum 38:  integer :: n, i 39:  40:  ! subinterval increment 41:  dx = (b-a)/real(n) 42:  ! minimum to increment from 43:  min = a - dx/2.0 44:  45:  ! midpoint rule 46:  do i = 1,n 47:  xi = min + i*dx 48:  sum = sum + f(xi) 49:  end do 50:  sum = sum*dx 51:  end function midpoint 52:  53: end program main 
 n: 2, M_n: .976360 n: 4, M_n: .993575 n: 6, M_n: .997091 n: 8, M_n: .998353 n: 10, M_n: .998942 n: 12, M_n: .999264 n: 14, M_n: .999459 n: 16, M_n: .999585 n: 18, M_n: .999672 n: 20, M_n: .999735

6.4 Polymorphism

An interface can be used as an entry into two different implementations of a subroutine or function with the same name so long as the different implementations have different argument signatures. This may be particularly useful for defining both a single precision and double precision version of a function or subroutine.

6.4.1 Example: machine epsilon

The following code implements two versions of a function that computes machine epsilon in either single or double precision. The different implementations are distinguished by their arguments. The single precision version mach_eps_sp accepts one single precision float and the double precision version mach_eps_dp accepts one double precision float. Both functions are listed in the interface and can be called by its name mach_eps. 

 1: program main  2:  implicit none  3:   4:  integer, parameter :: sp = kind(0.0)  5:  integer, parameter :: dp = kind(0.d0)  6:   7:  interface mach_eps  8:  procedure mach_eps_sp, mach_eps_dp  9:  end interface mach_eps 10:  11:  print*, mach_eps(0.0_sp), epsilon(0.0_sp) 12:  print*, mach_eps(0.0_dp), epsilon(0.0_dp) 13:  14: contains 15:  16:  function mach_eps_sp(x) result(eps) 17:  real(sp) :: x, eps 18:  integer :: count = 0 19:  20:  eps = 1.0_sp 21:  do while (1.0_sp + eps*0.5 > 1.0_sp) 22:  eps = eps*0.5 23:  count = count+1 24:  end do 25:  end function mach_eps_sp 26:  27:  function mach_eps_dp(x) result(eps) 28:  real(dp) :: x, eps 29:  integer :: count = 0 30:  31:  eps = 1.0_dp 32:  do while (1.0_dp + eps*0.5 > 1.0_dp) 33:  eps = eps*0.5 34:  count = count+1 35:  end do 36:  end function mach_eps_dp 37:  38: end program main 
 1.19209290E-07 1.19209290E-07 2.2204460492503131E-016 2.2204460492503131E-016

6.5 Recursion

A function or subroutine that calls itself must be defined with the recursive keyword preceding the construct name.

6.5.1 Example: factorial

The following code defines a recursive function factorial that computes \(n!\). If \(n>1\), the function call itself to return \(n(n-1)!\), otherwise the function returns \(1\). The main program calls factorial to compute \(5!\). 

 1: program main  2:  implicit none  3:   4:  ! print 5 factorial  5:  print*, factorial(5)  6:   7: contains  8:   9:  ! factorial(n): product of natural numbers up to n 10:  ! n: integer argument 11:  recursive function factorial(n) result(m) 12:  integer :: n, m 13:  14:  ! if n>1, call factorial recursively 15:  ! otherwise 1 factorial is 1 16:  if (n>1) then 17:  m = n*factorial(n-1) 18:  else 19:  m = 1 20:  end if 21:  22:  end function factorial 23:  24: end program main 
 120

7 Object-oriented programming

7.1 Derived types

Data types can be defined by the programmer. Variables and procedures that belong to a defined data type are declared between a type / end type pair. Type-bound procedures, i.e. functions and subroutines, are defined by the procedure keyword followed by :: and the name of the procedure within the type / end type pair after the contains keyword. A variable with defined type is declared with the type keyword and the name of the type. The variables and procedures of a defined type variable can be accessed by appending a % symbol to the name of the variable. 

 1: ! define a 'matrix' type  2: ! type-bound variables: shape, data  3: ! type-bound procedures: construct, destruct  4: type matrix  5:  integer :: shape(2)  6:  real, allocatable :: data(:,:)  7:  contains  8:  procedure :: construct  9:  procedure :: destruct 10: end type matrix 11:  12: ! declare a matrix variable 13: type(matrix) :: mat 14:  15: ! assign value to type-bound variable 16: mat%shape = [3,3]

7.2 Modules

A type-bound procedure can be defined after the contains keyword in the same program construct, i.e. a module, as the type definition. The first argument in the definition of a type-bound procedure is of the defined type and is declared within the procedure body with the class keyword and the name of the type. 

 1: module matrix_module  2:  implicit none  3:   4:  type matrix  5:  integer :: shape(2)  6:  real, allocatable :: data(:,:)  7:  contains  8:  procedure :: construct  9:  procedure :: destruct 10:  end type matrix 11:  12: contains 13:  14:  ! construct: populate shape and allocate memory for matrix 15:  ! m,n: number of rows,cols of matrix 16:  subroutine construct(this, m, n) 17:  class(matrix) :: this 18:  integer :: m, n 19:  this%shape = [m,n] 20:  allocate(this%data(m,n)) 21:  end subroutine construct 22:  23:  ! destruct: deallocate memory that matrix occupies 24:  subroutine destruct(this) 25:  class(matrix) :: this 26:  deallocate(this%data) 27:  end subroutine destruct 28:  29: end module matrix_module

To define variables of the matrix type in the main program, tell it to use the module defined above with use matrix_module immediately after the program main line. The procedures bound to a defined type can be access through variables of that type by appending the % symbol to the name of the variable. 

 1: program main  2:  use matrix_module  3:  implicit none  4:   5:  type(matrix) :: mat  6:  mat%shape = [3,3]  7:   8:  ! create matrix  9:  call mat%construct(3,3) 10:  11:  ! treat matrix variable 'data' like an array 12:  mat%data(1,1) = 1.0 13:  ! etc... 14:  15:  ! destruct matrix 16:  call matrix%destruct() 17: end program main

7.3 Example: determinant of random matrix

The following module defines a matrix type with two variables: an integer array shape that stores the number of rows and columns of the matrix and a real array data that stores the elements of the matrix. The type has four procedures: a subroutine construct that sets the shape and allocates memory for the data, a subroutine destruct that deallocates memory, a subroutine print that prints a matrix, and a function det that computes the determinant of a matrix. Note det is based on the definition of determinant using cofactors, and is very inefficient. A function random_matrix defined within the module generates a matrix with uniform random entries in \([-1,1]\). 

 1: module matrix_module  2:  implicit none  3:   4:  type matrix  5:  integer :: shape(2)  6:  real, allocatable :: data(:,:)  7:  contains  8:  procedure :: construct  9:  procedure :: destruct  10:  procedure :: print  11:  procedure :: det  12:  end type matrix  13:   14: contains  15:   16:  subroutine construct(this, m, n)  17:  class(matrix) :: this  18:  integer :: m,n  19:  this%shape = [m,n]  20:  allocate(this%data(m,n))  21:  end subroutine construct  22:   23:  subroutine destruct(this)  24:  class(matrix) :: this  25:  deallocate(this%data)  26:  end subroutine destruct  27:   28:  ! print: formatted print of matrix  29:  subroutine print(this)  30:  class(matrix) :: this  31:  ! row_fmt: format character string for row printing  32:  ! fmt: temporary format string  33:  character(32) :: row_fmt, fmt = '(a,i0,a,i0,a,i0,a)'  34:  ! w: width of each entry printed  35:  ! d: number of decimal digits printed  36:  integer :: w, d = 2, row  37:  ! find largest width of element in matrix  38:  w = ceiling(log10(maxval(abs(this%data)))) + d + 2  39:  ! write row formatting to 'row_fmt' variable  40:  write(row_fmt,fmt) '(',this%shape(2),'(f',w,'.',d,',1x))'  41:  ! print matrix row by row  42:  do row = 1,this%shape(1)  43:  print row_fmt, this%data(row,:)  44:  end do  45:  end subroutine print  46:   47:  ! det: compute determinant of matrix  48:  ! using recursive definition based on cofactors  49:  recursive function det(this) result(d)  50:  class(matrix) :: this  51:  type(matrix) :: submatrix  52:  real :: d, sgn, element, minor  53:  integer :: m, n, row, col, i, j  54:   55:  m = this%shape(1)  56:  n = this%shape(2)  57:  d = 0.0  58:   59:  ! compute cofactor  60:  ! if 1x1 matrix, return value  61:  if (m==1.and.n==1) then  62:  d = this%data(1,1)  63:  ! if square and not 1x1  64:  else if (m==n) then  65:  ! cofactor sum down the first column  66:  do row = 1,m  67:  ! sign of term  68:  sgn = (-1.0)**(row+1)  69:  ! matrix element  70:  element = this%data(row,1)  71:  ! construct the cofactor submatrix and compute its determinant  72:  call submatrix%construct(m-1,n-1)  73:  if (row==1) then  74:  submatrix%data = this%data(2:,2:)  75:  else if (row==m) then  76:  submatrix%data = this%data(:m-1,2:)  77:  else  78:  submatrix%data(:row-1,:) = this%data(:row-1,2:)  79:  submatrix%data(row:,:) = this%data(row+1:,2:)  80:  end if  81:  minor = submatrix%det()  82:  call submatrix%destruct()  83:   84:  ! determinant accumulator  85:  d = d + sgn*element*minor  86:  end do  87:  end if  88:  end function det  89:   90:  ! random_matrix: generate matrix with random entries in [-1,1]  91:  ! m,n: number of rows,cols  92:  function random_matrix(m,n) result(mat)  93:  integer :: m,n,i,j  94:  type(matrix) :: mat  95:  ! allocate memory for matrix  96:  call mat%construct(m,n)  97:  ! seed random number generator  98:  call srand(time())  99:  ! populate matrix 100:  do i = 1,m 101:  do j = 1,n 102:  mat%data(i,j) = 2.0*rand() - 1.0 103:  end do 104:  end do 105:  end function random_matrix 106:  107: end module matrix_module

The main program uses the matrix_module defined above to find the determinants of a number of random matrices of increasing size. 

 1: program main  2:  use matrix_module  3:  implicit none  4:   5:  type(matrix) :: mat  6:  integer :: n  7:   8:  ! compute determinants of random matrices  9:  do n = 1,5 10:  ! generate random matrix 11:  mat = random_matrix(n,n) 12:  13:  ! print determinant of matrix 14:  print "('n: ', i0, ', det: ', f0.5)", n, det(mat) 15:  16:  ! destruct matrix 17:  call mat%destruct() 18:  end do 19:  20: end program main 
./main 
 n: 1, det: -.68676 n: 2, det: .45054 n: 3, det: .37319 n: 4, det: -.27328 n: 5, det: .26695

7.4 Example: matrix module

 1: module matrix_module  2:  implicit none  3:   4:  public :: zeros  5:  public :: identity  6:  public :: random  7:   8:  type matrix  9:  integer :: shape(2)  10:  real, allocatable :: data(:,:)  11:  contains  12:  procedure :: construct => matrix_construct  13:  procedure :: destruct => matrix_destruct  14:  procedure :: norm => matrix_norm  15:  end type matrix  16:   17:  type vector  18:  integer :: length  19:  real, allocatable :: data(:)  20:  contains  21:  procedure :: construct => vector_construct  22:  procedure :: destruct => vector_destruct  23:  procedure :: norm => vector_norm  24:  end type vector  25:   26:  ! assignments  27:  interface assignment(=)  28:  procedure vec_num_assign, vec_vec_assign, mat_num_assign, mat_mat_assign  29:  end interface assignment(=)  30:   31:  ! operations  32:  interface operator(+)  33:  procedure vec_vec_sum, mat_mat_sum  34:  end interface operator(+)  35:   36:  interface operator(-)  37:  procedure vec_vec_diff, mat_mat_diff  38:  end interface operator(-)  39:   40:  interface operator(*)  41:  procedure num_vec_prod, num_mat_prod, mat_vec_prod, mat_mat_prod  42:  end interface operator(*)  43:   44:  interface operator(/)  45:  procedure vec_num_quot, mat_num_quot  46:  end interface operator(/)  47:   48:  interface operator(**)  49:  procedure mat_pow  50:  end interface operator(**)  51:   52:  ! functions  53:  interface norm  54:  procedure vector_norm, matrix_norm  55:  end interface norm  56:   57:  ! structured vectors/matrices  58:  interface zeros  59:  procedure zeros_vector, zeros_matrix  60:  end interface zeros  61:   62:  interface random  63:  procedure random_vector, random_matrix  64:  end interface random  65:   66: contains  67:   68:  subroutine matrix_construct(this, m, n)  69:  class(matrix) :: this  70:  integer :: m,n  71:  this%shape = [m,n]  72:  allocate(this%data(m,n))  73:  end subroutine matrix_construct  74:   75:  subroutine vector_construct(this, n)  76:  class(vector) :: this  77:  integer :: n  78:  this%length = n  79:  allocate(this%data(n))  80:  end subroutine vector_construct  81:   82:  subroutine matrix_destruct(this)  83:  class(matrix) :: this  84:  deallocate(this%data)  85:  end subroutine matrix_destruct  86:   87:  subroutine vector_destruct(this)  88:  class(vector) :: this  89:  deallocate(this%data)  90:  end subroutine vector_destruct  91:   92:  ! assignment  93:  subroutine vec_num_assign(vec,num)  94:  type(vector), intent(inout) :: vec  95:  real, intent(in) :: num  96:  vec%data = num  97:  end subroutine vec_num_assign  98:   99:  subroutine vec_vec_assign(vec1,vec2) 100:  type(vector), intent(inout) :: vec1 101:  type(vector), intent(in) :: vec2 102:  vec1%data = vec2%data 103:  end subroutine vec_vec_assign 104:  105:  subroutine mat_num_assign(mat,num) 106:  type(matrix), intent(inout) :: mat 107:  real, intent(in) :: num 108:  mat%data = num 109:  end subroutine mat_num_assign 110:  111:  subroutine mat_mat_assign(mat1,mat2) 112:  type(matrix), intent(inout) :: mat1 113:  type(matrix), intent(in) :: mat2 114:  mat1%data = mat2%data 115:  end subroutine mat_mat_assign 116:  117:  ! operations 118:  function vec_vec_sum(vec1,vec2) result(s) 119:  type(vector), intent(in) :: vec1, vec2 120:  type(vector) :: s 121:  call s%construct(vec1%length) 122:  s%data = vec1%data + vec2%data 123:  end function vec_vec_sum 124:  125:  function mat_mat_sum(mat1,mat2) result(s) 126:  type(matrix), intent(in) :: mat1, mat2 127:  type(matrix) :: s 128:  call s%construct(mat1%shape(1),mat1%shape(2)) 129:  s%data = mat1%data+mat2%data 130:  end function mat_mat_sum 131:  132:  function vec_vec_diff(vec1,vec2) result(diff) 133:  type(vector), intent(in) :: vec1, vec2 134:  type(vector) :: diff 135:  call diff%construct(vec1%length) 136:  diff%data = vec1%data-vec2%data 137:  end function vec_vec_diff 138:  139:  function mat_mat_diff(mat1,mat2) result(diff) 140:  type(matrix), intent(in) :: mat1, mat2 141:  type(matrix) :: diff 142:  call diff%construct(mat1%shape(1),mat1%shape(2)) 143:  diff%data = mat1%data-mat2%data 144:  end function mat_mat_diff 145:  146:  function num_vec_prod(num,vec) result(prod) 147:  real, intent(in) :: num 148:  type(vector), intent(in) :: vec 149:  type(vector) :: prod 150:  call prod%construct(vec%length) 151:  prod%data = num*vec%data 152:  end function num_vec_prod 153:  154:  function num_mat_prod(num,mat) result(prod) 155:  real, intent(in) :: num 156:  type(matrix), intent(in) :: mat 157:  type(matrix) :: prod 158:  call prod%construct(mat%shape(1),mat%shape(2)) 159:  prod%data = num*mat%data 160:  end function num_mat_prod 161:  162:  function mat_vec_prod(mat,vec) result(prod) 163:  type(matrix), intent(in) :: mat 164:  type(vector), intent(in) :: vec 165:  type(vector) :: prod 166:  call prod%construct(mat%shape(1)) 167:  prod%data = matmul(mat%data,vec%data) 168:  end function mat_vec_prod 169:  170:  function mat_mat_prod(mat1,mat2) result(prod) 171:  type(matrix), intent(in) :: mat1, mat2 172:  type(matrix) :: prod 173:  call prod%construct(mat1%shape(1),mat2%shape(2)) 174:  prod%data = matmul(mat1%data,mat2%data) 175:  end function mat_mat_prod 176:  177:  function vec_num_quot(vec,num) result(quot) 178:  type(vector), intent(in) :: vec 179:  real, intent(in) :: num 180:  type(vector) :: quot 181:  call quot%construct(vec%length) 182:  quot%data = vec%data/num 183:  end function vec_num_quot 184:  185:  function mat_num_quot(mat,num) result(quot) 186:  type(matrix), intent(in) :: mat 187:  real, intent(in) :: num 188:  type(matrix) :: quot 189:  call quot%construct(mat%shape(1),mat%shape(2)) 190:  quot%data = mat%data/num 191:  end function mat_num_quot 192:  193:  function mat_pow(mat1,pow) result(mat2) 194:  type(matrix), intent(in) :: mat1 195:  integer, intent(in) :: pow 196:  type(matrix) :: mat2 197:  integer :: i 198:  mat2 = mat1 199:  do i = 2,pow 200:  mat2 = mat1*mat2 201:  end do 202:  end function mat_pow 203:  204:  ! functions 205:  function vector_norm(this,p) result(mag) 206:  class(vector), intent(in) :: this 207:  integer, intent(in) :: p 208:  real :: mag 209:  integer :: i 210:  ! inf-norm 211:  if (p==0) then 212:  mag = 0.0 213:  do i = 1,this%length 214:  mag = max(mag,abs(this%data(i))) 215:  end do 216:  ! p-norm 217:  else if (p>0) then 218:  mag = (sum(abs(this%data)**p))**(1./p) 219:  end if 220:  end function vector_norm 221:  222:  function matrix_norm(this, p) result(mag) 223:  class(matrix), intent(in) :: this 224:  integer, intent(in) :: p 225:  real :: mag, tol = 1e-6 226:  integer :: m, n, row, col, iter, max_iters = 1000 227:  type(vector) :: vec, last_vec 228:  m = size(this%data(:,1)); n = size(this%data(1,:)) 229:  230:  ! entry-wise norms 231:  if (p<0) then 232:  mag = (sum(abs(this%data)**(-p)))**(-1./p) 233:  ! inf-norm 234:  else if (p==0) then 235:  mag = 0.0 236:  do row = 1,m 237:  mag = max(mag,sum(abs(this%data(row,:)))) 238:  end do 239:  ! 1-norm 240:  else if (p==1) then 241:  mag = 0.0 242:  do col = 1,n 243:  mag = max(mag,sum(abs(this%data(:,col)))) 244:  end do 245:  ! p-norm 246:  else if (p>0) then 247:  vec = random(n) 248:  vec = vec/vec%norm(p) 249:  last_vec = zeros(n) 250:  mag = 0.0 251:  do iter = 1,max_iters 252:  last_vec = vec 253:  vec = this*last_vec 254:  vec = vec/vec%norm(p) 255:  if (vector_norm(vec-last_vec,p)<tol) exit 256:  end do 257:  mag = vector_norm(this*vec,p) 258:  end if 259:  end function matrix_norm 260:  261:  ! structured vectors/matrices 262:  function random_matrix(m,n) result(mat) 263:  integer :: m,n 264:  type(matrix) :: mat 265:  call mat%construct(m,n) 266:  call random_seed() 267:  call random_number(mat%data) 268:  end function random_matrix 269:  270:  function random_vector(n) result(vec) 271:  integer :: n 272:  type(vector) :: vec 273:  call vec%construct(n) 274:  call random_seed() 275:  call random_number(vec%data) 276:  end function random_vector 277:  278:  function zeros_vector(n) result(vec) 279:  integer :: n 280:  type(vector) :: vec 281:  call vec%construct(n) 282:  vec = 0.0 283:  end function zeros_vector 284:  285:  function zeros_matrix(m,n) result(mat) 286:  integer :: m,n 287:  type(matrix) :: mat 288:  call mat%construct(m,n) 289:  mat = 0.0 290:  end function zeros_matrix 291:  292:  function identity(m,n) result(mat) 293:  integer :: m,n,i 294:  type(matrix) :: mat 295:  call mat%construct(m,n) 296:  do i = 1,min(m,n) 297:  mat%data(i,i) = 1.0 298:  end do 299:  end function identity 300:  301: end module matrix_module 
 1: program main  2:  use matrix_module  3:  implicit none  4:   5:  type(vector) :: vec1, vec2  6:  type(matrix) :: mat1, mat2  7:  real :: x  8:  integer :: i  9:  10:  ! 0s, id, random 11:  mat1 = zeros(3,3) 12:  call mat1%destruct() 13:  mat1 = identity(3,3) 14:  mat2 = random(3,3) 15:  mat1 = mat1*mat1 16:  vec1 = zeros(3) 17:  call vec1%destruct() 18:  vec1 = random(3) 19:  vec2 = random(3) 20:  ! +,-,*,/,** 21:  mat1 = mat1+mat2 22:  vec1 = vec1+vec2 23:  mat1 = mat1-mat2 24:  vec1 = vec1-vec2 25:  vec1 = mat1*vec2 26:  mat1 = mat2*mat1 27:  mat1 = 2.0*mat1 28:  vec1 = 2.0*vec1 29:  mat1 = mat1/2.0 30:  vec1 = vec1/2.0 31:  mat2 = mat1**3 32:  ! norm 33:  x = norm(vec1,0) 34:  x = norm(vec1,1) 35:  x = norm(mat1,-1) 36:  x = norm(mat1,0) 37:  x = norm(mat1,1) 38:  x = norm(mat1,2) 39:  call vec1%destruct 40:  call vec2%destruct 41:  call mat1%destruct 42:  call mat2%destruct 43: end program main 
./main 
 2 14 1.03737926 3.62866984E-07 3 10 0.976910591 1.07453801E-07 4 9 1.84589541 2.06476543E-07 5 10 2.39860868 2.45756240E-07

:snippets:

:end: