Lecture 9 11
Arrays
In mathematics, linear algebra, etc. we often define variables with subscripts. The dot product of two vectors A and B is given by:
The dot product can also be written in matrix notation, as can the general product of a row and column matrix. Define [A] to be a 1xn row matrix and {B} to be an nx1 column matrix as shown below:
The product of A and B is defined by:
Consider the product of two rectangular matrices [A] and [B]. [A] is defined to have n rows and m columns:
The product of [A] times [B] is defined by
or
Another example: suppose we need to define the pressure (call it P) at various positions along the cord of an airfoil (call this X) and say that we need 100 data points). It would be useful to have these variables subscripted so that we had the data pairs: , , … , .
In order to do such computations in Fortran, we need to define what are called arrays, which allow the use of subscripts.
Consider the following Fortran statement:
REAL, DIMENSION (1:50) :: A
or REAL :: A(1:50)
or REAL :: A(50)
The first form of the statement defines as an array the real variable A which can have subscripts that range from 1 to 50. The second form is an alternate, acceptable Fortran statement saying the same thing. The third form shows an implicit from of the initial subscript - since none is given, it is assumed to be 1.
Suppose we need a two dimensional array (like the matrix multiplication example above) that is of size 20 x 30. We could use the statement:
REAL, DIMENSION (1:20,1:30) :: A
or REAL :: A (1:20,1:30)
or REAL :: A(20,30) ß assumes subscripts
start at 1
or REAL, DIMENSION (20,30) :: A " "
The matrix multiplication
or
would look something like the following using subscripted variables.
REAL::A(100,100),B(100,100),C(100,100),SUM
(assume arrays A,B and C, n, m, and p are defined)
DO I=1,N
DO J=1,P
SUM=0.0
DO K=1,M
SUM=SUM+A(I,K)*B(K,J)
END DO
C(I,J)=SUM
END DO
END DO
There are two forms of Fortran statements that define arrays:
· Compile-time arrays – memory allocated at compile time
· Run-time arrays – memory allocated at execution time
Compile-Time Arrays
The general form of the compile-time array specification is given by:
type, DIMENSION (l:u, l:u, .., l:u) :: list-of-array-names
or
type :: array-name (l:u, l:u, .., l:u), … , array-name (l:u, l:u, .., l:u)
l = lower range of subscript
u = upper range of subscript
· l and u must be integer constants, but may be negative or positive.
· If l=1, it may be omitted.
Some examples of compile-time arrays
REAL :: A(10), B(1:10) ß A and B have subscripts 1, 2, … 10
REAL, DIMENSION (10) :: A, B ß same as above
INTEGER :: CAT (-5,5) ß CAT (integer) has subscripts:
-5, -4, …, +5
REAL :: K(10,20), Q(10,20) ß K and Q are 10x20 arrays
REAL :: K2(6,6,100) ß K2 is real and is a 6x6x100 array
In the above, A & B are one-dimensional arrays; and can be thought of as row or column matrices.
The array K is a two-dimensional real array; and is referred to as a 10x20 matrix with 10 rows and 20 columns. Range of 1st subscript is 1 to 10 and range of 2nd subscript is 1 to 20.
Another form of compile-time array specification uses the parameter statement to define the array size. For example,
INTEGER, PARAMETER :: N1=10, N2=20
REAL , DIMENSION (N1:N2) :: K ß K is a 10x20 array
or
INTEGER, PARAMETER :: N1=10, N2=20
REAL :: K(N1:N2) ß K is a 10x20 array
What happens, when using compile-time arrays, if ---
· If you use subscripts outside the allocated range?
· You get array overflow. (more on this on next page)
· The computer may or may not give you a run-time error!!!!!!!!!
· If you have array overflow and are unlucky enough to receive no run-time errors; your results will most likely be wrong but you may not know it unless you check your results very carefully.
· If you don’t use all the available subscripts allocated for an array?
· You waist computer memory!
· If you allocate a large amount of arrays and/or array space which, together with your program statements, requires more computer memory (RAM or virtual RAM) then you have available?
· If you are lucky, the link-edit step of BUILD will give you an error and you will not be able to execute the program. If the compile/link-editor does not give you an error, you may still have a code that will not execute; or will execute incorrectly!!
Computer memory is laid out sequentially, from location 1 to location n (last storage location). How then does Fortran decide where an element of a multi-dimensional array should be stored? It uses a formula based on the information in the dimension statement where the array is defined. Suppose we have REAL :: A(4),C(3,3). The compiler will allocate memory for these two arrays as shown on the right (usually). Note that the two-dimensional array is stored in memory by columns.
If you had the following statements: A(2) = 3.2, B(2,3)= - 4.5. The value "3.2" would be stored in the second location from the start of A. The value - 4.5 would be stored in the eighth location from the start of B. To get the location number, Fortran uses the following formula: Suppose the dimension information is B(N,M) and we want to store something in B(I,J), then the memory location in the computer is = N(J-1)+I; relative to beginning of the B array.
What happens if -- ?
· What happens if you have the statement A(6)=16.5? The value 16.5 gets stored in sixth location from start of A; or the memory location allocated to B(2,1). ß NOT GOOD! AND YOU MAY NOT GET A RUN-TIME ERROR OR WARNING!
· What happens if you have the statement B(4,2)=22.5? The value 22.5 gets stored in the 3(2-1)+4=7th location from start of B; or the memory location allocated to B(1,3). ß NOT GOOD! AND YOU MAY NOT GET A RUN-TIME ERROR OR WARNING!
· What happens if you have the statement B(1,4)=18.5? The value 22.5 get stored in the 3(4-1)+1=10th location from the start of B. The location is outside the dimensioned spaced allocated for B, but the number WILL BE STORED THERE! It may overwrite either a memory location reserved for another array, a variable, or perhaps even a Fortran instruction!! ß NOT GOOD! AND YOU MAY NOT GET A RUN-TIME ERROR OR WARNING!
Run-Time or Allocatable Arrays
Some of the problems noted above can be avoiding by allocating memory at run time by using allocatable or run-time arrays. The general form of the run-time (or allocatable) array specification requires using two statements:
type, DIMENSION (:,:,…,:), ALLOCATABLE :: list-of-array-names
and
ALLOCATE (list, STAT = status-variable)
where list is a list of array specifications of the from
array-name(l:u, l:u, …, l:u,) .
· The number of colons (say k colons) in DIMENSION indicates what dimension the array is (k-dimensional).
· The number (k) of l:u quantities specified in array-name must be consistent with that given in DIMENSION.
· Thus, (: , :) indicates a two-dimensional array.
· l & u may be constants or integer variables
Example:
REAL, DIMENSION (: , :) A
REAL, DIMENSION (: , :, :) B
ALLOCATE ( A(N,M) , B(10,10,10) , STAT=i_alloc_status )
The DIMENSION statement indicates that A is a 2-D array and B is a 3-D array. The ALLOCATE statement allocates memory for array A as NxM and B as 10x10x10. The actual size of A depends on the values of the variables N and M – they can be assigned, or read as input.
The STAT=status-variable is optional; if used the status-variable is an integer variable whose value will be set to zero if allocation is successful or set to an error value if allocation is not successful.
The run-time form of allocating arrays is useful because the program allocates space during execution (and only as much as you need). Thus, one does not have to recompile a working program when larger (or smaller) arrays are required because of different input data (or program execution requirements).
Example: Suppose you want an array A to have subscripts 1 through N but N depends on the particular problem being solved (i.e., depends on input data). We can use the following:
INTEGER :: N, istat
REAL, DIMENSION (:), ALLOCATABLE :: A
READ *, N
ALLOCATE (A(N), STAT=istat)
IF (istat .ne. 0) PRINT *, “STOP, not enough memory”
(rest of program which uses the array A)
· If N is input as 20, this allocates 20 storage locations for A (subscripts 1 to 20).
· If N is input as 10000, this allocates 10000 storage locations for A.
· The STAT parameter is used to check if allocation was successful (STAT variable set to zero by computer it enough memory was available). The parameter is optional.
· The ALLOCATE statement is an executable statement. If sufficient memory is not available, the program must be stopped. If you don’t stop execution, unpredictable results can be expected!!!!!
Example: Suppose we had N (x-y) data points and that we wanted to do the following: 1) numerically integrate the data to find the area under the curve: ; or 2) determine a curve fit y=y(x); or 3) plot the data; etc. For example, to read the data, and numerically integrate from X(1) to X(N) using the trapezoidal rule, we could use the following:
INTEGER :: N, istat
REAL ::Area, DeltaX
REAL, DIMENSION (:), ALLOCATABLE :: X, Y
READ *, N
ALLOCATE (X(N), Y(N), STAT=istat)
IF (istat .ne. 0) PRINT *, “STOP, not enough memory”
DO I = 1, N
READ *, X(I), Y(I)
END DO
! numerically integrate with trapezoidal rule
Area = 0.0
DO I = 1,N – 1
DeltaX = X(I+1) – X(I)
Area = Area +( (Y(I+1) + Y(I))/2.0 ) * DeltaX
END DO
(rest of program)
We could have done this with one two-dimensional array by storing the x values in the 1st column and the y values in the 2nd column.
INTEGER :: N, istat
REAL ::Area, DeltaX
REAL, DIMENSION (: , :), ALLOCATABLE :: XYdata
READ *, N
ALLOCATE (XYdata(N,2), STAT=istat)
IF (istat /= 0) PRINT *, “STOP, not enough memory”
DO I = 1, N
READ *, XYdata(I,1), XYdata(I,2) ! x in first column
END DO
! numerically integrate with trapezoidal rule
Area = 0.0
DO I = 1,N – 1
DeltaX = XYdata(I+1,1) – XYdata(I,1)
Area = Area +( (XYdata(I+1,2) + XYdata(I,2))/2.0 ) * DeltaX
END DO
(rest of program)