[Next] [Previous] [Top]

4 Array Processing

4.1 Terminology and Specifications

REAL, DIMENSION(50) :: w
REAL, DIMENSION(5:54) :: x
REAL y(50)
REAL z(11:60)

Here, w, x, y and z are all arrays containing 50 elements.

The rank of an array is the number of dimensions. Thus, a scalar has rank 0, a vector has rank 1 and a matrix has rank 2.

The extent refers to a particular dimension, and is the number of elements in that dimension.

The shape of an array is a vector consisting of the extent of each dimension.

The size of an array is the total number of elements which make up the array. This may be zero.

Two arrays are said to be conformable if they have the same shape. All arrays are conformable with a scalar.

Take, for example the following arrays:

REAL, DIMENSION :: a(-3:4,7)
REAL, DIMENSION :: b(8,2:8)
REAL, DIMENSION :: d(8,1:8)
INTEGER :: c

The array a has

• rank 2

• extents 8 and 7

• shape (/8,7/)

• size 56

The general form of an array specification is as follows:

type [[,DIMENSION (extent-list)] [,attribute]... ::] entity list

This is simply a special case of the form of declaration given in "Specifications".

Here, type can be any intrinsic type or a derived type (so long as the derived type definition is accessible to the program unit declaring the array). DIMENSION is optional and defines default dimensions in the extent-list, these can alternatively by defined in the entity list.

The extent-list gives the array dimensions as:

• integer constants

• integer expressions using dummy arguments or constants

• ':' to show the array is allocatable or assumed shape

PARAMETER
PUBLIC
PRIVATE
POINTER
TARGET
ALLOCATABLE
DIMENSION(extent-list)
INTENT(inout)
OPTIONAL
SAVE
EXTERNAL
INTRINSIC

Finally, the entity list is a list of array names with optional dimensions and initial values.

The following examples show the form of the declaration of several kinds of arrays, some of which are new to Fortran 90 and will be met later in this chapter:

1. Initialisation of one-dimensional arrays containing 3 elements:

1. Declaration of automatic array logb. Here loga is a dummy array argument, and SIZE is an intrinsic function which returns a scalar default integer corresponding to the size of the array loga:

1. Declaration of dynamic (allocatable) arrays a and b. The dimensions would be defined in a subsequent ALLOCATE statement:

1. Declaration of assumed shape arrays a and b. The dimensions would be taken from the actual calling routine:

4.2 Whole Array Operations

In order for whole array operations to be performed, the arrays concerned must be conformable. Remember, that for two arrays to be conformable they must have the same shape, and any array is conformable with a scalar. Operations between two conformable arrays are carried out on an element by element basis, and all intrinsic operations are defined between two such arrays.

For example, if a and b are both 2x3 arrays

a = , b =

the result of addition is

a + b =

and of multiplication is

a x b =

If one of the operands is a scalar, then the scalar is broadcast into an array which is conformable with the other operand. Thus, the result of adding 5 to b is

b + 5 = + =

Such broadcasting of scalars is useful when initialising arrays and scaling arrays.

An important concept regarding array-valued assignment is that the right hand side evaluation is computed before any assignment takes place. This is of relevance when an array appears in both the left and right hand side of an assignment. If this were not the case, then elements in the right hand side array may be affected before the operation was complete.

The advantage of whole array processing can best be seen by comparing examples of Fortran 77 and Fortran 90 code:

1. Consider three one-dimensional arrays all of the same length. Assign all the elements of a to be zero, then perform the assignment a(i) = a(i)/3.1 + b(i)*SQRT(c(i)) for all i.

   Fortran 77 Solution

Fortran 90 Solution

REAL, DIMENSION(20) :: a, b, c
...
a=0
...
a=a/3.1+b*SQRT(c)

Note, the intrinsic function SQRT operates on each element of the array c.

1. Consider three two-dimensional arrays of the same shape. Multiply two of the arrays element by element and assign the result to the third array.

   Fortran 77 Solution

Fortran 90 Solution

REAL, DIMENSION (5, 5) :: a, b, c
...
c = a * b

1. Consider a three-dimensional array. Find the maximum value less than 1000 in
   this array.
   In Fortran 77 this requires triple DO loop and IF statements, whereas the Fortran 90 code is:

1. Find the average value greater than 3000 in an array.
   In Fortran 77 this requires DO loops and IF statements, whereas Fortran 90 code is:

Note in the last two examples the use of the following array intrinsic functions:

MAXVAL - returns the maximum array element value.
SUM - returns the sum of the array elements.
COUNT - returns the number of true array elements.

4.3 Elemental Intrinsic Procedures

The following are examples of elemental intrinsic procedures:

1. Find the square roots of all elements of an array, a. (Note that the SQRT function has already been seen in an example in "Whole Array Operations".)

1. Find the string length excluding trailing blanks for all elements of a character array ch.

4.4 WHERE Statement

A simple example is to avoid division by zero:

REAL, DIMENSION(5,5) : ra, rb
...
WHERE(rb>0.0) ra=ra/rb

The general form is

WHERE(logical-array-expression) array-variable=array-expression

The logical expression is evaluated, and all those elements of array expression which have value true are evaluated and assigned to array variable. The elements which have value false remain unchanged. Note that the logical array expression must have the same shape as the array variable.

It is also possible for one logical array expression to determine a number of array assignments. The form of this WHERE construct is:

WHERE (logical-array-expression)
array-assignment-statements
END WHERE

or

WHERE (logical-array-expression)
array-assignment-statements
ELSEWHERE
array-assignment-statements
END WHERE

In the latter form, the assignments after the ELSEWHERE statement are performed on those elements that have the value false for the logical array expression.

For example, the WHERE construct can be used to divide every element of the array ra by the corresponding element of the array ra avoiding division by zero, and assigning zero to those values of ra corresponding to zero values of rb.

REAL, DIMENSION(5,5) :: ra,rb
...
WHERE(rb>0.0)
ra=ra/rb
ELSEWHERE
ra=0.0
END WHERE

4.5 Array Sections

Array sections can be extracted using either:

• A simple subscript.

• A subscript triplet.

• A vector subscript.

4.5.1 Simple Subscripts

= ra(2,2)

4.5.2 Subscript Triplets

[lower bound]:[upper bound][:stride]

If either the lower bound or upper bound is omitted, then the bound of the array from which the array section is extracted is assumed, and if stride is omitted the default stride=1 is used.

The following examples show various array sections of an array using vector subscripts. The elements marked with x denote the array section. Let the defined array from which the array section is extracted by a 5x5 array.

= ra(2:2,2:2); Array element, shape (/1/)

= ra(3,3:5); Sub-row, shape(/3/)

= ra(:,3); Whole column; shape(/5/)

= ra(1::2,2:4); Stride 2, shape(/3,3/)

4.5.3 Vector Subscripts

An example of an integer expression of rank 1 is:

(/3,2,12,2,1/)

An example showing the use of a vector subscript iv is:

REAL, DIMENSION :: ra(6), rb(3)
INTEGER, DIMENSION (3) :: iv
iv = (/ 1, 3, 5 /) ! rank 1 integer expression
ra = (/ 1.2, 3.4, 3.0, 11.2, 1.0, 3.7 /)
rb = ra(iv) ! iv is the vector subscript
! = (/ ra(1), ra(3), ra(5) /)
! = (/ 1.2, 3.0, 1.0 /)

Note that the vector subscript can be on the left hand side of an expression:

iv = (/1, 3, 5/) ! vector subscript
ra(iv) = (/1.2, 3.4, 5.6/)
! = ra((/1, 3, 5/)) = (/1.2, 3.4, 5.6/)
! = ra(1:5:2) = (/1.2, 3.4, 5.6/)

It is also possible to use the same subscript more than once and hence using a vector subscript an array section that is bigger than the parent array can be constructed. Such a section is called a many-one array section. A many-one section cannot appear on the left hand side of an assignment statement or as an input item in a READ statement, as such uses would be ambiguous.

iv = (/1, 3, 1/)
ra(iv) = (/1.2, 3.4, 5.6/) ! not permitted
! = ra((/1, 3, 1/) = (/1.2, 3.4, 5.6/)

rb = ra(iv) ! permitted
! = ra(/1, 3, 1/) = (/1.2, 3.4, 1.2/)

4.6 Array Assignment

REAL, DIMENSION(5,5) :: ra,rb,rc
INTEGER :: id
.
.
.
! Shape(/5,5/) and scalar
ra = rb + rc*id

! Shape(/3,2/)
ra(3:5,3:4) = rb(1::2,3:5:2) + rc(1:3,1:2)

! Shape(/5/)
ra(:,1) = rb(:,1) + rb(:,2) + rc(:,3)

4.7 Recursion

For example, the code:

DO i=2,n
x(i) = x(i) + x(i-1)
END DO

is not the same as:

x(2:n)= x(2:n) + x(1:n-1)

In the first case, the assignment is:

x(i) = x(i) + x(i-1) + x(i-2) + ... + x(1)

whereas in the second case the assignment is:

x(i) = x(i) + x(i-1)

In order to achieve the recursive effect of the DO-loop, in Fortran 90 it would be appropriate to use the intrinsic function SUM. This function returns the sum of all the elements of its array argument. Thus the equivalent assignment is:

x(2:n) = (/(SUM(1:i), i=2,n)/)

4.8 Element Renumbering

REAL, DIMENSION (1:8) :: ra
REAL, DIMENSION (-3:4) :: rb
INTEGER, DIMENSION (1) :: locmax1, locmax2
REAL :: max1, max2
ra = (/ 1.2, 3.4, 5.4, 11.2, 1.0, 3.7, 1.0, 1.0/)
rb = ra
! To find location of max value

locmax1 = MAXLOC(ra) ! = (/ 4 /)
locmax2 = MAXLOC(rb) ! = (/ 4 /)

! To find maximum value from location

max1 = ra(locmax(1))
! OK because ra defined with 1 as lower bound

max2 = rb(LBOUND(rb) + locmax2(1) - 1)
! general form required with lower bound not equal to 1

4.9 Zero Sized Arrays

Zero sized arrays are useful for boundary operations. For example,

DO i=1,n
x(i)=b(i)/a(i,i)
b(i+1:n)=b(i+1:n) - a(i+1:n,1)*x(i)
! zero sized when i=n
END DO

4.10 Array Constructors

(/ array-constructor-value-list /)

For example,

REAL, DIMENSION(6) :: a
a=(/array-constructer-value-list/)

where, for example, array-constructer-value-list can be any of:

(/(i,i=1,5)/)
! = (/1,2,3,4,5/)

(/7,(i,i=1,4),9/)
! = (/7,1,2,3,4,9/)

(/1.0/REAL(i),i=1,4)/)
! = (/1.0/1.0,1.0/2.0,1.0/4.0/)

(/((i+j,i=1,3),j=1,2)/)
! = (/2,3,4,3,4,5/)

(/a(i,2:4),a(1:5:2,i+3)/)
! = (/a(i,2),a(i,3),a(i,4),a(1,i+3),a(3,i+3),a(5,i+3)

It is possible to transfer a rank-one array of values to an array of a different shape using the RESHAPE function. The RESHAPE function has the form

RESHAPE(SOURCE,SHAPE,[,PAD][,ORDER])

where the argument SOURCE can be an array of any sort (in this case a rank-one array), and the elements of source are rearranged to form an array RESHAPE of shape SHAPE. If SOURCE has more elements than RESHAPE, then the unwanted elements will be ignored. If RESHAPE has more elements than SOURCE, then the argument PAD must be present. The argument PAD must be an array of the same type as SOURCE, and the elements of PAD are used in array element order, using the array repeatedly if necessary, to fill the missing elements of RESHAPE. Finally, the optional argument ORDER allows the elements of RESHAPE to be placed in an alternative order to array element order. The array ORDER must be the same size and shape as SHAPE, and contains the dimensions of RESHAPE in the order that they should be run through.

A simple example is:

REAL, DIMENSION(3,2) :: ra
ra=RESHAPE((/((i+j,i=1,3),j=1,2)/),SHAPE=(/3,2/))

which creates ra=

If the argument ORDER is included as follows

ra=RESHAPE((/(((j-1)*3+i,i=1,3),j=1,2)/),SHAPE= &
(/3,2/),ORDER(2,1))

then the result would be ra=

4.11 Allocatable Arrays

Allocatable arrays allow large chunks of memory to be used only when required and then be released. This produces a much more efficient use of memory than Fortran 77, which offered only static (fixed) memory allocation.

An allocatable array is declared in a type declaration statement with the attribute ALLOCATABLE. The rank of the array must also be specified in the declaration statement and this can be done by including the appropriate number of colons in the DIMENSION attribute. For example, a two dimensional array could be declared as:

REAL, DIMENSION(:,:), ALLOCATABLE :: a

This form of declaration statement does not allocate any memory space to the array. Space is dynamically allocated later in the program, when the array is required, using the ALLOCATE statement. The ALLOCATE specifies the bounds of the array and, as with any array allocation, the lower bound defaults to one if only the upper bound is specified. For example, the array declared above could be allocated with lower bound zero:

ALLOCATE (a(0:n))

The bounds may also be integer expressions, for example:

ALLOCATE (a(0:n+1))

The space allocated to the array with the ALLOCATE statement can later be released with the DEALLOCATE statement. The DEALLOCATE statement requires only the name of the array concerned and not the shape. For example:

DEALLOCATE (a)

Allocatable arrays make possible the frequent requirement to declare an array having a variable number of elements. For example, it may be necessary to read variables, say nsize1 and nsize2, and then declare an array to have nsize1 x nsize2 elements:

INTEGER n
REAL, DIMENSION(:,:), ALLOCATABLE :: ra
.
.
.
READ(*,*) nsize1,nsize2
ALLLOCATE (ra(nsize1,nsize2))
.
.
.
DEALLOCATE (ra)


Note that both ALLOCATE and DEALLOCATE statements can allocate/deallocate several arrays in one single statement.

An allocatable array is said to have an allocation status. When an array has been defined but not allocated the status is said to be unallocated or not currently allocated. When an array appears in an ALLOCATE statement then the array is said to be allocated, and once the array has been deallocated it is said to be not currently allocated. The DEALLOCATE statement can only be used on arrays which are currently allocated, and similarly, the ALLOCATE statement can only be used on arrays which are not currently allocated. Thus, ALLOCATE can only be used on a previously allocated array if it has been deallocated first.

It is possible to check whether or not an array is currently allocated using the intrinsic function ALLOCATED. This is a logical function with one argument, which must be the name of an allocatable array. Using this function, statements like the following are possible:

IF (ALLOCATED(a)) DEALLOCATE(a)

or

IF (.NOT.ALLOCATED(a)) ALLOCATE(a(5,20))

An allocatable array has a third allocation status, undefined. An array is said to be undefined if it is allocated within a procedure and the program returns to the calling program without deallocating it. Once an array is undefined, it can no longer be used. Hence it is good programming practice to deallocate all arrays that have been allocated. There are, however, two other ways around this problem. Firstly, an allocatable array can be declared with the SAVE attribute:

REAL, DIMENSION(:), ALLOCATABLE, SAVE :: a

This permits the allocatable array to remain allocated upon exit from the procedure and preserves the current values of the array elements. Secondly, the allocatable arrays could be put into modules, and in this case the arrays are preserved as long as the executing program unit uses the modules. The array can also be ALLOCATED and DEALLOCATED by any program unit using the module which the array was declared in.

Finally, there are three restrictions on the use of allocatable arrays:

• Allocatable arrays cannot be dummy arguments of a procedure and must, therefore, be allocated and deallocated in the same program unit

• The result of a function cannot be an allocatable array

• Allocatable arrays cannot be used in a derived type definition

It is frequently necessary, particularly in numerical analysis, to make use of temporary work arrays. In Fortran 77, this presented serious problems for providers of subroutine libraries, who had to resort to requiring the calling sequence to include the work arrays along with the genuine parameters. The parameter list was further lengthened by the need to pass information about the dimensions of the array. This example, where b is a work array, shows how an allocatable array can be used to overcome this problem.

Fortran 77 code:

SUBROUTINE invert(a,n,inv,b)
INTEGER n
REAL DIMENSION a(n,n), inv(n,n), b(n,2*n)
m = 2*n
.
.
.
END

Fortran 90 code:

SUBROUTINE invert(a,inverse_a)
REAL DIMENSION(:,:), INTENT(IN) :: a
REAL DIMENSION(:,:), INTENT(OUT) :: inverse_a
INTEGER :: n,m
REAL DIMENSION(:,:), ALLOCATABLE :: b
.
.
.
n = SIZE(a)
m = 2*n
ALLOCATE (b(n,m))
.
.
.
DEALLOCATE
END SUBROUTINE invert

Whenever this procedure is called, only two genuine parameters are needed as the work array, b, is declared, used, and discarded entirely within the procedure, and we can use the intrinsic function SIZE to remove the need to pass the parameter n. SIZE returns the scalar size of the array a. SIZE is described in more detail in "Automatic Arrays"

4.12 Automatic Arrays

Automatic arrays are automatically created (allocated) upon entry to the procedure in which they are declared, and automatically deallocated upon exit from the procedure. Thus the size of the automatic array can be different in different procedure calls.

The intrinsic function SIZE is often used when declaring automatic arrays. SIZE has the form:

SIZE(ARRAY [,DIM])

This returns the extent of ARRAY along dimension DIM, or returns the size of ARRAY if DIM is absent.

Note that an automatic array must not appear in a SAVE or NAMELIST statement, nor be initialised in a type declaration.

The following example shows the automatic arrays, work1 and work2 which take their size from the dummy arguments n and a:

SUBROUTINE sub(n,a)
INTEGER :: n
REAL, DIMENSION(n,n) :: a
REAL, DIMENSION(n,n) :: work1
REAL, DIMENSION(SIZE(a,1)) :: work2
END SUBROUTINE sub

The next example shows automatic array bounds dependent on a global variable defined in a module. Both use association and host association are shown:

MODULE auto_mod
INTEGER :: n
CONTAINS
SUBROUTINE sub
REAL, DIMENSION(n) :: w
WRITE (*, *) 'Bounds and size of a: ', &
LBOUND(w), UBOUND(w), SIZE(w)
END SUBROUTINE sub
END MODULE auto_mod

PROGRAM auto_arrays
! automatic arrays using modules instead of
! procedure dummy arguments
USE auto_mod
IMPLICIT NONE
n = 10
CALL sub
END PROGRAM auto_arrays

In the example below the power of dynamic arrays can be seen when passing only part of an array to a subroutine. Suppose the main program declares a nxn array, but the subroutine requires a n1xn1 section of this array a. In order to achieve this in Fortran 77, both parameters n and n1 must be passed as subroutine arguments:

PROGRAM array
PARAMETER (n=10)
REAL a(n,n)
REAL work(n,n)
...
READ(*,*) n1
CALL sub(a,n,n1,res,work)
...
END PROGRAM array


SUBROUTINE sub(a,n,n1,res,work)
REAL a(n,n1)
REAL work(n1,n1)
...
res=a(...)
...
END SUBROUTINE sub

Note that the work array is required only as local array, but for flexibility standard procedure in Fortran 77 is to also pass it as an argument.

Using dynamic arrays in Fortran 90, this can be achieved with much simplified argument passing:

PROGRAM array
REAL, ALLOCATABLE, DIMENSION(:,:) :: a
...
READ(*,*) n1
ALLOCATE(A(n1,n1))
CALL sub(a,n1,res)
DEALLOCATE(a)
...
END


SUBROUTINE sub(a,n1,res)
REAL, DIMENSION(n1,n1) :: a
REAL, DIMENSION(n1,n1) :: work
...
res=a(...)
...
END SUBROUTINE sub

Notice that using an allocatable array a, the array is exactly the size we require in the main program and so we can pass this easily to the subroutine. The work array, work, is an automatic array whose bounds depend on the dummy argument n1.

4.13 Assumed Shape Arrays

[lower_bound]:

where the lower bound defaults to 1 if omitted.

Assumed shape arrays make possible the passing of arrays between program units without having to pass the dimensions as arguments. However, if an external procedure has an assumed shape array as a dummy argument, then an interface block must be provided in the calling program unit.

For example, consider the following external subprogram with assumed shape arrays ra, rb and rc:

SUBROUTINE sub(ra,rb,rc)
REAL, DIMENSION (:,:) :: ra ! Shape (10, 10)
REAL, DIMENSION (:,:) :: rb ! Shape (5, 5)
! = REAL, DIMENSION (1:5,1:5) :: rb
REAL, DIMENSION (0:,2:) :: rc ! Shape (5, 5)
! = REAL, DIMENSION (0:4,2:6) :: rc
.
.
.
END SUBROUTINE sub

The calling program might include:

REAL, DIMENSION (0:9,10) :: ra ! Shape (10, 10)

INTERFACE
SUBROUTINE sub(ra,rb,rc)
REAL, DIMENSION (:,:) :: ra,rb
REAL, DIMENSION (0:,2:) :: rc
END SUBROUTINE sub
END INTERFACE
.
.
.
CALL SUB (ra,ra(0:4,2:6),ra(0:4,2:6))

The following example uses allocatable, automatic and assumed shape arrays, and shows another method of coding the final example in "Automatic Arrays":

PROGRAM array
REAL, ALLOCATABLE, DIMENSION(:,:) :: a

INTERFACE
SUBROUTINE sub(a,res)
REAL, DIMENSION (:, :) :: a
REAL, DIMENSION (SIZE(a, 1),SIZE(a, 2))
END SUBROUTINE sub
END INTERFACE
...
READ (*, *) n1
ALLOCATE (a(n1, n1)) ! allocatable array
CALL sub(a,res)
...
END PROGRAM array

SUBROUTINE sub(a,res)
REAL, DIMENSION (:, :) :: a ! assumed shape array
REAL, DIMENSION (SIZE(a, 1),SIZE(a, 2)) :: work
! automatic array
...
res = a(...)
...
END SUBROUTINE sub

4.14 Array Intrinsics

ALL(MASK[,DIM])
True if all elements true

ANY(MASK[,DIM])
True if any element true

COUNT(MASK[,DIM])
Number of true elements

MAXVAL(ARRAY[,DIM][,MASK])
Maximum element value

MINVAL(ARRAY[,DIM][,MASK])
Minimum element value

PRODUCT(ARRAY[,DIM][,MASK])
Product of array elements

SUM(ARRAY[,DIM][,MASK])
Sum of array elements

ALLOCATED(ARRAY)
True if array allocated

LBOUND(ARRAY[,DIM])
Lower bounds of array

SHAPE(SOURCE)
Shape of array (or scalar)

SIZE(ARRAY[,DIM])
Size of array

UBOUND(ARRAY[,DIM])
Upper bounds of array

MERGE(TSOURCE,FSOURCE,MASK)
Merge arrays subject to mask

PACK(ARRAY,MASK[,VECTOR])
Pack elements into vector subject to mask

SPREAD(SOURCE,DIM,NCOPIES)
Construct an array by duplicating an array section

UNPACK(VECTOR,MASK,FIELD)
Unpack elements of vector subject to mask

RESHAPE(SOURCE,SHAPE[,PAD][,ORDER])
Reshape array

MAXLOC(ARRAY[,MASK])
Location of maximum element

MINLOC(ARRAY[,MASK])
Location of minimum element

CSHIFT(ARRAY,SHIFT[,DIM])
Perform circular shift

EOSHIFT(ARRAY,SHIFT[,BOUNDARY][,DIM])
Perform end-off shift

TRANSPOSE(MATRIX)
Transpose matrix

DOT_PRODUCT(VECTOR_A,VECTOR_B)
Compute dot product

MATMUL(MATRIX_A,MATRIX_B)
Matrix multiplication

The following example shows the use of several intrinsic functions:

Three students take four exams. The results are stored in an INTEGER array:

score(1:3,1:4) =

• Largest score:

• Largest score for each student:

• Student with largest score:

• Average score:

• Number of scores above average:

above = score > average
! above(3, 4) is a LOGICAL array
! above =
n_gt_average = COUNT (above) ! = 6
! n_gt_average is an INTEGER variablE

• Pack all scores above the average:

• Did any student always score above the average?

• Did all students score above the average on any of the tests?

4.15 Array Example

PROGRAM conjugate_gradients

IMPLICIT NONE
INTEGER :: iters, its, n
LOGICAL :: converged
REAL :: tol, up, alpha, beta
REAL, ALLOCATABLE :: a(:,:),b(:),x(:),r(:),u(:),p(:),xnew(:)

READ (*,*) n, tol, its
ALLOCATE ( a(n,n),b(n),x(n),r(n),u(n),p(n),xnew(n) )

OPEN (10, FILE='data')
READ (10,*) a
READ (10,*) b

x = 1.0
r = b - MATMUL(a,x)
p = r

iters = 0

DO
iters = iters + 1
u = MATMUL(a, p)
up = DOT_PRODUCT(r, r)
alpha = up / DOT_PRODUCT(p, u)
xnew = x + p * alpha
r = r - u * alpha
beta = DOT_PRODUCT(r, r) / up
p = r + p * beta
converged = ( MAXVAL(ABS(xnew-x)) / &
MAXVAL(ABS(x)) < tol )
x = xnew
IF (converged .OR. iters == its) EXIT
END DO

WRITE (*,*) iters
WRITE (*,*) x

END PROGRAM conjugate_gradients

4.16 Exercises

1. Run the program matrix.f90 which declares a 2-dimensional integer array, with extents (n,n), where n is set to 9 in a PARAMETER statement.
   
   This program uses DO loops to assign elements of the array to have values r c, where r is the row number and c is the column number, e.g., a(3,2) = 32, a(5,7) = 57. It writes the resulting array to the file matrix.dat for later use.
   

2. From the array constructed in exercise 1, use array sections to write out:
   (a) the first row
   (b) the fifth column
   (c) every second element of each row and column, columnwise
   ( a(1,1), a(3,1),...)
   (d) every second element of each row and column, rowwise
   ( a1,1), a(1,3),...)
   (e) the 3 non-overlapping 3x3 sub-matrices in columns 4 to 6
   (section.f90)
   

3. Write a program which generates an 8x8 chequerboard, with 'B' and 'W' in alternate positions. Assume the first position is 'B'. (board.f90)
   

4. From the array constructed in exercise 1, use the WHERE construct to create an array containing all of the odd values and 0 elsewhere (use elemental function, MOD). (where.f90)
   

5. Declare a vector subscript, iv, with extent 5. From the array constructed in exercise 1 create a 9x5 array containing only the odd values. (vec_subs.f90)
   

6. Generate the array constructed in exercise 1 using a single array constructor. (reshape.f90)
   

7. Look at the Fortran 77 code sum2.f. Rewrite it using Fortran 90 with allocatable and assumed-shape arrays. (sum4.f90)
   
   Is there any instrinsic function which can simplify the same job? (sum5.f90)
   

8. Create an integer array whose size is allocated dynamically (read size from terminal). Assign odd and even values to the array (same as matrix.f90). Pass the array to a subroutine which uses an assumed shape argument and returns all odd values of the array and 0 elsewhere.
   (odd_val.f90)
   

9. Run the program spread1.f90. Modify it to create an real array with element values 1.0/REAL(i+j+1), where i is the row number and j is the column number. (spread2.f90)
   
   Can you find another way using Fortran 90 array?
   

10. Look at the program m_basis.f90. Modify it to select all values greater than 3000 and find the number of them, the maximum, the minimum and the average. (munro.f90)