↜ Back to index Introduction to Numerical Analysis 2
Part b–Lecture 1
Goal of this quarter
During this quarter we will learn methods to solve Ax = b that are faster than the Gaussian elimination (and even the LU decomposition) for special matrices. The main goal is to implement the conjugate gradient method (共役勾配法) for symmetric positive definite matrices. The basic code will turn out to be easier than the Gaussian elimination.
But first we learn how to define and use functions in Fortran. It will be more convenient to write the code using functions.
Functions in Fortran
In mathematics, a function (関数) like \sin or \operatorname{abs} is a rule x \mapsto \sin x that given one or more inputs produces an output. We have to specify the domain and the range, that is, what type of inputs are allowed and what type of output is to be expected: \sin: \mathbb R \to \mathbb R.
A function in Fortran captures a similar idea. In the function definition, we specify the name of the function, the types of inputs and the output, and provide the code that produces the output based on the inputs. The inputs are also called function arguments. We then can use the function in other parts of the code to easily get the output for the given inputs.
The main advantage of using functions is that we do not need to remember how exactly the function is implemented to use it.
To illustrate the syntax, let us define a function add_two
that simply adds 2 to its argument and returns the new value, and a function foo
that multiplies its two arguments increased by two:
implicit none
print *, add_two(2.0)
print *, add_two(10.0)
print *, foo(1.0, 3.0)
contains
function add_two(x)
real x
real add_two
= x + 2
add_two end
function foo(x, y)
real x, y, foo
= add_two(x) * add_two(y)
foo end
end
Before the end
command that concludes the program we included a new section contains
where one can define functions that are then accessible from the main program (the part before contains
).
A function definition starts with the keyword function
and ends with the keyword end
. function
must be followed by the name of the function and the list of arguments in parentheses, separated by commas. The lines that follow then declare all the variables used in the function code, including the types of the arguments (x
in the above code) and the output (add_two
, the same as the name of the function). The rest of the code before the end of the function is the code used to produce the output value.
We can call functions from other functions.
Here is a more complicated example of a function that computes a factorial.
implicit none
integer i
do i = 1, 5
print *, 'Factorial of', i, 'is', fact(i)
end do
contains
function fact(n)
integer n, fact, i
= 1
fact do i = 1, n
= fact * i
fact end do
end
end
This time we are using a local variable i
within the function. Note one very important fact: since i
is declared inside the function fact
, it is a different variable from the variable i
in the main program. In particular, the do loop inside fact
will not change the value of i
in the main program and the produced result is correct.
Note that the return variable fact
can be used like any other variable. The value that is returned is the value of the variable when the function ends.
Exercise 1. Try removing i
from the line integer n, fact, i
and rerun the program.
Exercise 2. Implement a function maximum
that takes two real numbers and returns their maximum. Implement another function positive_part
that takes one real number x as an input and returns \max(x, 0) using the function maximum
.
Using arrays as function arguments
Functions in Fortran can accept arrays as arguments and even return arrays.
The following function takes a one dimensional array of any size (indicated by (:)
in the declaration of x
) and returns the array of the differences of the consecutive elements.
implicit none
real a(3), b(5)
= [3.0, 2.5, 5.0]
a = [1.0, 3.0, 2.0, 2.5, 5.0]
b
print *, diffs(a)
print *, diffs(b)
contains
function diffs(x)
real x(:)
real diffs(size(x) - 1)
integer i
do i = 1, size(diffs)
= x(i+1) - x(i)
diffs(i) end do
end
end
Here is a more complicated function that takes a matrix as an input and returns the array of the sums of each row of the matrix.
implicit none
real a(3,2)
! initialize a by setting its rows
1, :) = [1., -1.]
a(2, :) = [2., 3.]
a(3, :) = [6., 7.]
a(
print *, row_sum(a)
contains
function row_sum(a)
real a(:,:) ! matrix of any size
real row_sum(size(a, 1)) ! size equal to the number of rows of a
integer i, j
! iterate of rows
do i = 1, size(a, 1)
= 0
row_sum(i) ! iterate over columns
do j = 1, size(a, 2)
= row_sum(i) + a(i, j)
row_sum(i) end do
end do
end
end
The code takes as input a two dimension array of any size indicated as (:,:)
, and returns a one-dimensional array with the length equal to the number of rows in a
.
Built-in function size
You can read more about the built-in function size
here. In short, size(a)
returns the number of all the elements of array a
, while size(a, 1)
returns the number of rows and size(a, 2)
the number of columns in the matrix represented by a
.
real a(5, 3)
print *, size(a, 1)
print *, size(a, 2)
print *, size(a)
end
prints
5
3
15
Exercise 3.
Implement a function inner
that takes two one-dimensional arrays x
and y
of the same size representing two vectors x, y and returns their inner product x \cdot y = (x, y) := x_1 y_1 + \cdots + x_n y_n.
Given arrays
real a(3), b(3)
= [3.0, 2.5, 5.0]
a = [1.0, -3.0, 2.0]
b
print *, inner(a, b)
we get
5.50000000
Note. Fortran 95 already has a built-in function dot_product
that does the same thing.
Exercise 4.
Implement a function mv_mul
that takes a two-dimensional array a
representing a matrix A and a one-dimensional array x
representing a column vector x as input and return the one-dimensional array that contains the component of the product Ax.
Hint. Recall that (Ax)_i = a_{i1} x_1 + a_{i2} x_2 + \cdots + a_{in} x_n. This is the inner product of row i of matrix A and vector x. You can therefore use the function inner
from the previous exercise if you want.
With the following arrays:
real a(3,2), x(2)
1, :) = [1., -1.]
a(2, :) = [2., 3.]
a(3, :) = [6., 7.]
a(
= [3., -1.]
x
print *, mv_mul(a, x)
we get
4.00000000 3.00000000 11.0000000
Exercise 5.
Implement a function norm
that that takes a one-dimensional array x
and returns the Euclidean norm \|x\| := \sqrt{x \cdot x} of the represented vector x.
Hint. You can use the function inner
implemented before.
For array
real x(3)
= [1., 2., -3.]
x
print *, norm(x)
we get
3.74165750
Note. Fortran 2008 has a built-in function norm2
that does the same thing.
More on Fortran functions
Recursive functions
In programming, functions are often called recursively (再帰呼び出し): a function calls itself with different argument values.
An important example is the definition of a factorial (階乗): n! = \begin{cases} 1 & n = 0,\\ n \cdot (n-1)! & n > 0. \end{cases}
In Fortran, we had to declare the function as recursive
. In most other languages, all functions can be called recursively. Fortran is very old…
implicit none
print *, fact(5)
print *, fact(12)
contains
! 👇👇👇👇👇
recursive function fact(n) result(r)
integer n, r
if (n == 0) then
= 1
r return ! return from the function
endif
= n * fact(n - 1)
r end
end
We also had to specify the name for a return value using result(r)
.
return
command is useful if you want to return from a function early.
Exercise 6. Write a function fib(n)
for computing the Fibonacci number F_n using the recursive formula
F_n :=
\begin{cases}
0, & n = 0,\\
1, & n = 1,\\
F_{n-2} + F_{n-1}, & n \geq 2.
\end{cases}
How many times is fib
called when computing fib(5)
? How many times for fib(10)
?
Exercise 7. Write a function gcd(n, k)
that computes the greatest common divisor (最大公約数) of two positive integers n and k recursively using the Euclidean algorithm (ユークリッドの互除法):
gcd(n, 0)
is n- if k > 0, then
gcd(n, k)
is the same asgcd(k, mod(n, k))
.
mod(n, k)
is the Fortran’s built-in function for the remainder (剰余).
Extra exercises
Exercise 8. Write a program that attempts to find the eigenvector of a given symmetric matrix A corresponding to the largest eigenvalue using the following algorithm.
Set x^{(0)} to some nonzero vector.
Iterate for i = 1, 2, …
Compute x^{(i)} = \frac{Ax^{(i-1)}}{\|Ax^{(i-1)}\|} Here the program should report an error if \|Ax^{(i-1)}\| is too small, say 10^{-6}.
Print the value of x^{(i)} to track progress.
Stop the iteration if \|x^{(i)} - x^{(i-1)}\| < \varepsilon.
Use \varepsilon = 10^{-6}.
Hints.
Use the functions
mv_mul
andnorm
that you implemented in the previous exercises.In Fortran, you can divide all the elements of an array x by a real number y by simply writing x / y. No need to use
do
loops. Similarly, to compute the element-wise difference of two arrays v and w of the same size you can use v - w.In numerical algorithms we usually restrict the maximum number of iterations in case the algorithm does not converge. You can use for example
do i = 1, 100
to limit the number of iterations to 100.
Sample run of the program with
real a(3,3), x(3)
1,:) = [1.0, 2.0, 3.0]
a(2,:) = [2.0, 1.0, 4.0]
a(3,:) = [3.0, 4.0, 1.0]
a(
= [1., 0., 0.] x
is
0.267261237 0.534522474 0.801783681
0.549971938 0.628539383 0.549971938
0.490831792 0.557763398 0.669316173
0.511315823 0.596535206 0.618629038
0.503328860 0.578738987 0.641655087
0.506864071 0.586905181 0.631372452
0.505297840 0.583228707 0.636017621
0.506004035 0.584889233 0.633927822
0.505686522 0.584141433 0.634870052
0.505829632 0.584478498 0.634445608
0.505765200 0.584326684 0.634636879
0.505794227 0.584395111 0.634550691
0.505781174 0.584364235 0.634589553
0.505787075 0.584378183 0.634572029
0.505784392 0.584371865 0.634579957
0.505785584 0.584374666 0.634576380
0.505785048 0.584373474 0.634577930
0.505785286 0.584374011 0.634577274