168 6 Linear Algebra
One of the better features of Fortran is dynamic storage allocation. That is, the size of an
array can be changed during the execution of the program. To see how the dynamic allocation
works in Fortran, consider the following simple example where we set up a 4 × 4 unity matrix.
IMPLICITNONE
! The definition of the matrix, using dynamic allocation
REAL,ALLOCATABLE,DIMENSION(:,:) :: unity
! The size of the matrix
INTEGER:: n
! Here we set the dim n=4
n=4
! Allocate now place in memory for the matrix
ALLOCATE( unity(n,n) )
! all elements are set equal zero
unity=0.
! setup identity matrix
DOi=1,n
unity(i,i)=1.
ENDDO
DEALLOCATE( unity)
We always recommend to use the deallocation statement, since this frees space in memory.
If the matrix is transferred to a function from a calling program, one can transfer the dimen-
sionalitynof that matrix with the call. Another possibility is to determine the dimensionality
with theSIZEfunction. Writing a statement liken=SIZE(unity,DIM=1)gives the number of rows,
while using DIM=2 gives the number of columns. Note however that this involves an extra
call to a function. If speed matters, one should avoid such calls.
6.4 Linear Systems
In this section we outline some of the most used algorithms tosolve sets of linear equations.
These algorithms are based on Gaussian elimination [24,28]and will allow us to catch several
birds with a stone. We will show how to rewrite a matrixAin terms of an upper and a lower
triangular matrix, from which we easily can solve linear equation, compute the inverse ofA
and obtain the determinant. We start with Gaussian elimination, move to the more efficient
LU-algorithm, which forms the basis for many linear algebraapplications, and end the discus-
sion with special cases such as the Cholesky decomposition and linear system of equations
with a tridiagonal matrix.
We begin however with an example which demonstrates the importance of being able to
solve linear equations. Suppose we want to solve the following boundary value equation
−
d^2 u(x)
dx^2
=f(x,u(x)),withx∈(a,b)and with boundary conditionsu(a) =u(b) = 0. We assume thatfis a continuous
function in the domainx∈(a,b). Since, except the few cases where it is possible to find ana-
lytic solutions, we will seek approximate solutions, we choose to represent the approximation
to the second derivative from the previous chapter
f′′=fh−^2 f^0 +f−h
h^2
+O(h^2 ).