Parallel Computer Centre

A major new feature of Fortran 90 are the array processing capabilities. It is possible to work directly with a whole array or an array section without explicit

`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:

- Initialisation of one-dimensional arrays containing 3 elements:

`INTEGER, DIMENSION(3) :: ia=(/1,2,3/), ib=(/(i,i=1,3)/)`

- 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*:

`LOGICAL, DIMENSION(SIZE(loga)) :: logb`

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

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

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

`REAL, DIMENSION(:,:,:) :: a,b`

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:

- 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**

...

DO 10 i=1,20

a(i)=0

10 CONTINUE

...

DO 20 i=1,20

a(i)=a(i)/3.1 + b(i)*SQRT(c(i))

20 20 CONTINUE

**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.**

- 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**

...

DO 20 i = 1, 5

DO 10 j = 1, 5

c(j, i) = a(j, i) * b(j, i)

10 CONTINUE

20 CONTINUE

**Fortran 90 Solution**

` REAL, DIMENSION (5, 5) :: a, b, c`

...

c = a * b

- 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:

`REAL, DIMENSION(10,10,10) :: a`

amax=MAXVAL(a,MASK=(a<1000))

- 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:

`av=SUM(a,MASK=(a>3000))/COUNT(MASK=(a>3000))`

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.

The following are examples of elemental intrinsic procedures:

- 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".)

`b=SQRT(a)`

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

`length=LEN_TRIM(ch)`

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

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

Array sections can be extracted using either:

- A simple subscript.
- A subscript triplet.
- A vector subscript.

= ra(2,2)

`[`

*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/)

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/)

`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)

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)/)`

`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

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

`(/ 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*=

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"**

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 *n*x*n* array, but the subroutine requires a *n1*x*n1* 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*.

`[`

*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

*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:

`MAXVAL (score) ! = 90`

- Largest score for each student:

`MAXVAL (score, DIM = 2) `

! = (/ 90, 80, 66 /)

- Student with largest score:

`MAXLOC (MAXVAL (SCORE, DIM = 2))`

! = MAXLOC ((/ 90, 80, 66 /)) = (/ 1 /)

- Average score:

`average = SUM (score) / SIZE (score) ! = 62`

! average is an INTEGER variable

- 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:

`...`

INTEGER, ALLOCATABLE, DIMENSION (:) :: &

score_gt_average

...

ALLOCATE (score_gt_average(n_gt_average)

scores_gt_average = PACK (score, above)

! = (/ 85, 71, 66, 76, 90, 80 /)

- Did any student always score above the average?

`ANY (ALL (above, DIM = 2)) ! = .FALSE.`

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

`ANY (ALL (above, DIM = 1)) ! = .TRUE.`

`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`

- 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.

- 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*)

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

- 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*)

- 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*)

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

- 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*)

- 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*)

- 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?

- 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*)

[Next] [Previous] [Top]

All documents are the responsibility of, and copyright, © their authors and do not represent the views of The Parallel Computer Centre, nor of The Queen's University of Belfast.

Maintained by Alan Rea, email A.Rea@qub.ac.uk Generated with CERN WebMaker