Arrays

Multidimensional Arrays

Matrices

• arrays can have more than one dimension

integer, dimension(10,10) :: a

• defines a 2-dimensional integer array a with 100 elements

• matrix with 10 rows and 10 columns

• element of array a is accessed by specifying two indices

• a(i,j) is the element of array a in the i^th row and j^th column, 1 <= i,j <= 10

Bounds on Indices

• separate lower and upper bounds may be specified for each dimension

integer, dimension(-10:10,0:19,51:100) :: a

• defines a three dimensional array a

• the first index varies from –10 to 10, second from 0 to 19 and third from 51 to 100

• total number of elements 21*20*50 = 21000

• Fortran allows arrays with up to 7 dimensions

Terminology

• rank is the number of dimensions

• size is the total number of elements

• bounds in a particular dimension specify the range of valid index values in that dimension

• extent in a particular dimension is the number of valid index values in that dimension

• shape is a one dimensional array of extents

• type of an array is the type of its elements

Examples

integer, dimension(-10:10,0:19,51:100) :: a

integer, dimension(0:20,-19:0,50) :: b

• rank of a = 3 and size = 21000

• bounds in dimension 1 are –10 (lower) and 10

• shape of a = (/21,20,50/) = shape of b

• arrays with same shape are said to conform

• conforming arrays can be used in expressions

• a = b is a valid assignment statement

Element Ordering

• matrix elements are conceptually ordered in column wise order

• need not be physically stored in that order

integer, dimension(2,2,2) :: a

_• a_1,1,1 a_2,1,1 a_1,2,1 a_2,2,1 a_1,1,2 a_2,1,2 a_1,2,2 a_2,2,2

• first index varies most frequently, then second..

• important for initialization and i/o of arrays

• index values used to directly access element

Reading Arrays

• arrays can be read in one element at a time

integer, dimension(10,10,10) :: a

do k = 1,10 ; do j = 1, 10 ; do i = 1,10

read *, a(i,j,k) ! separate read for each element

end do; end do; end do

• short notation for this using implied do

read *, (((a(i,j,k),i=1,10),j=1,10),k=1,10) ! or

read *, a

• all elements read from one line with one read

Array Assignment

• individual elements can be assigned values

• a conforming array may be assigned

– values assigned to corresponding variables

• scalar expression assigned to entire array

• operations on elements extended to array

• arithmetic operations applied to entire array

• operations for corresponding elements done in “parallel” ( matrix addition)

Array Sections

• array sections defined as for rank 1 arrays

• section defined separately for each dimension

• index may be specified for some dimensions

• rank of section <= rank of array

• sections may be used in array operations

integer, dimension(10,10,10) :: a

• a(3, :5, 2::2) is an array with shape (/5,5/)

a(3,1,2) a(3,2,2)….a(3,5,2) a(3,1,4)….a(3,5,10)

Examples of Sections

• assign colors to chessboard

integer, dimension(8,8) :: chessboard

chessboard( : :2, : :2) = 1 ! 1 for white

! white color for squares with two odd indices

chessboard(2: :2, 2: :2) = 1

! white color if both indices even

chessboard( : : 2, 2: :2) = 2 ! 2 for black

! black if first index is odd and second is even

chessboard(2: : 2, : :2) = 2

! black if first index is even and second odd

Array Initialization

• arrays may be initialized in declaration

• reshape function defines array constants

integer, dimension(2,2) :: a = &

reshape ( (/ 1, 2, 3, 4 /) , (/ 2, 2 /) )

! assigns a(1,1) = 1, a(2,1) = 2, a(1,2) = 3, a(2,2) = 4

• reshape function creates an array of desired shape from another array of same size

• required shape is second parameter

• default element ordering used

Allocatable Arrays

• many times required array size is not known to the programmer

• number of students in a course may vary

• be conservative and declare a “large” array

• leads to waste of memory locations

• if user needs a larger array, program would have to be recompiled

• Fortran 90 allows user to choose array size

• main new features of Fortran 90

Allocatable Arrays

• array may be given the allocatable attribute

integer,dimension( : ),allocatable :: a

integer,dimension( : , : ),allocatable :: b

• declares a to be a rank 1 integer array

• b is a rank 2 integer array

• only the rank and type of array is declared

• index bounds are not specified

• compiler does not allocate memory

• array does not actually “exist”

• allocatable arrays created during execution

• can slow down execution considerably

allocate(a(n), stat = i)

! i and n must be integer variables, n must have a value

• creates a rank 1 array with size = value of n

• stat indicates success of allocation

• i is set to 0 if successful and /= 0 otherwise

• value of i must be checked before proceeding

allocate(b(-10:10,51:100), stat = i)

• bounds on indices may be specified for each dimension (may be variables or constants)

• array with allocatable attribute cannot be used unless it has been allocated

• an allocated array may be deallocated

deallocate(a, stat = i)

• maybe allocated again with different size

• always decallocate arrays after using them

Array Intrinsic Functions

• Fortran 90 has many array intrinsic functions

• elemental, inquiry, transformational functions

• elemental functions operate on individual elements separately

– sqrt(a), abs(a), sin(a),.. (a is a real array)

– result is an array of the same shape

– values of elements are changed

Array Intrinsic Functions

• inquiry functions give properties of the array

– size(a), shape(a), lbound(a), ubound(a)

– allocated(a)

• shape, lbound, ubound return a rank 1 array containing the shape, lower and upper bounds of array a respectively

• lbound(a,dim=1), size(a,dim=1) give lower bound and extent in dimension 1 respectively

• allocated(a) is .true. if a is allocated

Array Intrinsic Functions

• transformational functions work on entire array and may compute scalars or arrays of different shapes

– sum, dot_product, reshape, count, matmul, maxloc, minloc

• maxloc(a) gives location of maximum value in integer/real array a as a rank 1 array

• arrays computed by functions and array sections cannot be indexed, must be used in whole array operations

Linear Systems of Equations

• most commonly occurring problem, the original motivation for matrices

• solve a system of n equations in n variables

a₁₁x₁ + a₁₂x₂ + …… + a_1nx_n = b₁

a₁₂x₂ + a₂₂x₂ + …… + a_2nx_n = b₂

……….

a_n1x₁ + a_n2x₂ + …… + a_nnx_n = b_n

• a_i’s are all real numbers

Gaussian Elimination

• choose any equation and any variable with non zero coefficient ( called pivot) in the equation

• eliminate this variable from all other equations by adding a multiple of chosen equation

• repeat with remaining equations

• solve one equation in one variable

• back substitute value of variable to find others

• two steps, elimination and back substitution

Program

program linear

implicit none

real, dimension(:,:), allocatable :: a

integer :: n, i, j

real, dimension(:), allocatable :: x

print *, "enter number of variables"

read *, n

allocate(x(n), stat=i)

if ( i /= 0) then

print *, "allocation of x failed"

stop

endif

allocate(a(n,n+1), stat=i)

if ( i /= 0) then

print *, "allocation of array a failed"

stop

endif

do i = 1,n

print *, "enter coefficients of equation", i

read *, a(i, 1:n+1)

end do

do i = 1,n

a(i, i:n+1) = a(i,i:n+1)/a(i,i)

! make coefficient of x_i 1.0

do j = i+1,n

a(j,i:n+1) = a(j,i:n+1) - a(i,i:n+1)*a(j,i)

end do ! eliminate x_i

end do

do i = n, 1, -1

x(i) = a(i,n+1) - sum(a(i,i+1:n)*x(i+1:n))

end do ! sum of a null array is 0.0

print *, "solution is"

print *, x(1:n)

deallocate(a, x, stat=i)

if ( i /= 0) then

print *, "deallocation failed"

stop

endif

end program linear

Elimination

• a(1:n+1,1:n+1) contains the coefficients, x(1:n) the computed values of variables

• assumes a(i,i) /= 0.0 at every step

• true for many special types of matrices

• numerical problems if abs(a(i,i)) is too small

• elimination needs to be done only once for different right hand sides

• only back substitution needs to be done

^• reduces operations from n³ to n²

Pivoting

• swap rows and columns to make a(i,i) large

• partial pivoting - swap row i with row containing maxval(abs(a(i:n,i)))

– what if it is 0.0?

• complete pivoting- swap row and column to make a(i,i) = maxval(abs(a(i:n,i:n)))

– what if it is still 0.0?

– this changes order of the variables

• gives more accurate results but more time

LU Factorization

• a(i,i) = 1.0 after elimination, instead store value used to divide i^th row in elimination

• a(i,j) = 0.0 for j < i, instead store the multiple of i^th row subtracted from j^th row during elimination

• all required multipliers stored in array a

• matrix a is written as a product of an upper and lower triangular matrix

• a permutation used to keep track of pivoting