[ previous ] [ Contents ] [ 1 ] [ 2 ] [ 3 ] [ 4 ] [ 5 ] [ 6 ] [ 7 ] [ 8 ] [ 9 ] [ 10 ] [ 11 ] [ 12 ] [ 13 ] [ next ]

# `Fortran 90` Lessons for Computational Chemistry Chapter 4 - More on Arrays

## 4.1 Objectives

The main aims of this lesson are the following:

1. presenting storage ordering of multidimensional arrays.

1. presenting how to manipulate whole matrices or arrays sections in `Fortran`.

1. matrix definition using the WHERE statement.

## 4.2 Main items.

• Storage ordering

Multidimensional arrays are stored in memory by `Fortran` in such a way that the first subindex varies faster than the second, that varies faster than the third and so on and so forth. This is known as column major order.

For example, if we define a 4 x 2 matrix as

REAL , DIMENSION(1:4,1:2) :: A,

the A array has eight elements stored into memory as follows

```     A(1,1), A(2,1), A(3,1), A(4,1), A(1,2), A(2,2), A(3,2), A(4,2)
```

The A matrix initialization can be carried out in several ways. Assuming that each element should be initialized with a number equal to the index of the corresponding row, we could use two loops[5]

```     DO I_col = 1, 2
DO I_row = 1, 4
A(I_row, I_col) = I_row
ENDDO
ENDDO
```

An array constructor can also be of help, though the seemingly simple solution

```     A = (/ 1, 2, 3, 4, 1, 2, 3, 4 /)
```

does not work. The array constructors produce vectors and not matrices. The vector defined above is of dimension 8, but not a matrix 4x2. The vector and the array A have identical sizes, but are not conformal. The statement RESHAPE gives a possible solution. The sytax of this statement is

```     output_array = RESHAPE(array_1, array_2)
```

Where array_1 is a matrix that would be reshaped and array_2 is a vector with the dimensions of the new matrix output_array. The total number of elements of array_1 and output_array needs to be identical. In the previous example a correct array constructor is

```     A = RESHAPE( (/ 1, 2, 3, 4, 1, 2, 3, 4 /), (/ 4, 2 /) )
```

Another example can be found in code excode_4_3.f90, Section 4.3.3. The RESHAPE command can be used in the array declaration

```     INTEGER, DIMENSION(1:4,1:2) :: A = &
RESHAPE( (/ 1, 2, 3, 4, 1, 2, 3, 4 /), (/ 4, 2 /) )
```

The data ordering in storage is specially important in I/O operations. The command

```     PRINT*, A
```

will give as a result

```     A(1,1), A(2,1), A(3,1), A(4,1), A(1,2), A(2,2), A(3,2), A(4,2)
```

It is necessary to take this into account also when making use of the READ statement to fill with values the elements of a multidimensional array: READ(unit,*) A. The implicit DO statement allow to change the standard reading sequence

```     READ(unit,*) ( ( A(row,col), col = 1, 2 ), row = 1, 4 )
```
• The usage of whole matrices is a great advantage. Definig floating point vectors V1, V2, V3 y V4 as

```     REAL , DIMENSION(1:21) :: V1, V2
REAL , DIMENSION(-10:10) :: V3
REAL , DIMENSION(0:10) :: V4
```

The following tasks are simply performed using this `Fortran 90` feature.

1. Assigning a particular value to the full array:

```     V1 = 0.5
```
1. Equating matrices:

```     V1 = V2
```

Making each V1 element equal to the corresponding element of V2. This is only valid when both matrices are conformal. It is also valid

```     V3 = V2
```

but it is not valid

```     V1 = V4
```
1. All arithmetic operation for scalars can be also applied to conformal matrices, though they may not be the expected mathematical operations.

```     V1 = V2 + V3
V1 = V2*V3
```

In the first case V1 is the sum of two vectors, but in the second case each V1 element is the product of the corresponding V2 and V3 elements, which is not the scalar product. In the two-dimensional matrices case, if we define

```     REAL , DIMENSION(1:4,1:4) :: A, B, C
```

The following are valid statements in `Fortran 90`

```     A = A**0.5
C = A + B
C = A * B
```

The last case is not the matrix product but a matrix having each element as the result of the product of the corresponding A annd B elements.

1. A matrix can also be read without a DO loop, as in example excode_4_1.f90, Section 4.3.1, where also the intrinsic function SUM is presented.

• The definition of array slices is possible using the index syntax liminf:limsup:step

```     V1(1:10) = 0.5
B(1,1:4) = 100.0
```

In the first case the first ten elements of the V1 array take the value 0.5, while in the second elements in the first row of B take the value 100.0. See example excode_4_1.f90, Section 4.3.1.

The most general syntax to define a slice is lowlimit:upplimit:step, the first slice element has index lowlimit, the last one is less than or equal to upplimit and step is the index variable increment. The default value of step is step=1. Examples:

```     V1(:)       ! the whole vector
V1(3:10)    ! elements V1(3), V1(4), ... ,   V1(10)
V1(3:10:1)  !   ""      ""     ""             ""
V1(3:10:2)  !   ""     V1(3), V1(5), ... , V1(9)
V1(m:n)     ! elements V1(m), V1(m+1), ... , V1(n)
V1(9:4:-2)  !   ""     V1(9), V1(7), V1(5)
V1(m:n:-k)  ! elements V1(m), V1(m-k), ... , V1(n)
V1(::2)     !   ""     V1(1), V1(3), ... , V1(21)
V1(m:m)     ! 1 x 1 array
V1(m)       ! Scalar
```
• The assignment of values to an array can be done making use of a logic mask, with the WHERE statement. The use of the mask allows to select those array elements that should undergo the initialization. If, e.g., we need to compute the square root of the elements of a floating point array called data_mat and store them in the array sq_data_mat, we can skip the use of loops and conditionals as in the following code

```     DO j_col = 1, dim_2
DO i_row = 1, dim_1
IF ( data_mat(i_row, j_col) >= 0.0 ) THEN
sq_data_mat(i_row, j_col) = SQRT( data_mat(i_row, j_col) )
ELSE
sq_data_mat(i_row, j_col) = -99999.0
ENDIF
ENDDO
ENDDO
```

The WHERE statement greatly simplifies this task. The statement syntax is

```     [name:] WHERE (mask_expr_1)
....
Array assignment block 1
....
....
Array assignment block 2
....
ELSEWHERE
....
Array assignment block 3
....
ENDWHERE [name]
```

where mask_expr_1 and mask_expr_2 are boolean arrays conformal with the array being assigned. The previous example is therefore simplified to

```     WHERE ( data_mat >= 0.0 )
sq_data_mat = SQRT( data_mat )
ELSEWHERE
sq_data_mat = -99999.0
ENDWHERE
```
• These aspects are treated in the different given examples. Example excode_4_3.f90, Section 4.3.3 shows how to initialize vectoras and matrices, in the last case making use of the RESHAPE statement. The example also introduces the `Fortran` intrinsics DOT_PRODUCT (scalar product) and MATMUL (matrices product).

Example excode_4_4.f90, Section 4.3.4 exemplifies the use WHERE in combination with a logical mask.

Example excode_4_5.f90, Section 4.3.5 stress the fact the the elimination of DO loops can sometimes bring surprising results about.

Example excode_4_6.f90, Section 4.3.6 shows how to use the RESHAPE statement in the definition of a matrix and how to use slicing in the defined array.

## 4.3 Example Codes.

### 4.3.1 excode_4_1.f90

```     PROGRAM ex_4_1
!
! VARIABLE DEFINITION
IMPLICIT NONE
REAL :: Total=0.0, Average=0.0
REAL , DIMENSION(:), ALLOCATABLE :: t_worked
! Correction Factor
REAL :: correction =1.05
INTEGER :: day, num_days
!
PRINT *,' Number of workdays: '
! Dynamic storage definition
ALLOCATE(t_worked(1:num_days))
!
PRINT *,' Worked hours per day in ', num_days, ' days.'
! I/O
!
t_worked(num_days-1:num_days) = correction*t_worked(num_days-1:num_days)
!
DO day=1,num_days
Total = Total + t_worked(day)
ENDDO
Average = Total / num_days
!
PRINT *,' Average daily hours of work in ',num_days, ' days : '
PRINT *, Average
!
END PROGRAM ex_4_1
```

### 4.3.2 excode_4_2.f90

```     PROGRAM ex_4_2
!
! Program to characterize an array making use of inquiry functions
!
IMPLICIT NONE
!
REAL, DIMENSION(:,:), ALLOCATABLE :: X_grid
INTEGER :: Ierr
!
!
ALLOCATE(X_grid(-20:20,0:50), STAT = Ierr)
IF (Ierr /= 0) THEN
STOP 'X_grid allocation failed'
ENDIF
!
WRITE(*, 100) SHAPE(X_grid)
100 FORMAT(1X, "Shape :        ", 7I7)
!
WRITE(*, 110) SIZE(X_grid)
110 FORMAT(1X, "Size :         ", I7)
!
WRITE(*, 120) LBOUND(X_grid)
120 FORMAT(1X, "Lower bounds : ", 7I6)
!
WRITE(*, 130) UBOUND(X_grid)
130 FORMAT(1X, "Upper bounds : ", 7I6)
!
DEALLOCATE(X_grid, STAT = Ierr)
IF (Ierr /= 0) THEN
STOP 'X_grid deallocation failed'
ENDIF
!
END PROGRAM EX_4_2
```

### 4.3.3 excode_4_3.f90

```     PROGRAM ex_4_3
!
! VARIABLES DEFINITION
IMPLICIT NONE
REAL, DIMENSION(1:5) :: VA = (/1.0,1.0,1.0,1.0,1.0/), PMAT
INTEGER I
INTEGER, DIMENSION(1:5) :: VB = (/(2*I,I=1,5)/)
REAL :: PE
REAL , DIMENSION(1:5,1:5) :: MC
REAL , DIMENSION(25) :: VC = &
(/ 0.0,0.0,0.0,0.0,1.0,0.5,2.0,3.2,0.0,0.0, &
0.0,0.0,0.0,0.0,11.0,0.5,2.3,3.2,0.0,0.0, &
1.0,3.0,-2.0,-2.0,-0.6 /)
! Scalar Product
PE = DOT_PRODUCT(VA,VB)
!
PRINT *, 'Scalar Product (VA,VB) = ', PE
!
! Product of matrices VAxMC
! RESHAPE VC to make it a  5 x 5 matrix
MC = RESHAPE(VC,(/5,5/))
PMAT = MATMUL(VA,MC)
!
PRINT *, ' VA x MC = ', PMAT(1:5)
!
END PROGRAM ex_4_3
```

### 4.3.4 excode_4_4.f90

```     PROGRAM ex_4_4
IMPLICIT NONE
REAL , DIMENSION(-180:180) :: Time=0
INTEGER :: Degree, Strip
REAL :: Value
CHARACTER (LEN=1), DIMENSION(-180:180) :: LEW=' '
!
DO Degree=-165,165,15
Value=Degree/15
DO Strip=-7,7
Time(Degree+Strip)=Value
ENDDO
ENDDO
!
DO Strip=0,7
Time(-180 + Strip) = -180/15
Time( 180 - Strip) = 180/15
ENDDO
!
DO Degree=-180,180
PRINT *,Degree,' ',Time(Degree), 12 + Time(Degree)
END DO
!
WHERE (Time > 0)
LEW='E'
ELSEWHERE (Time < 0)
LEW='W'
ENDWHERE
!
PRINT*, LEW
!
END PROGRAM ex_4_4
```

### 4.3.5 excode_4_5.f90

```     PROGRAM ex_4_5
!
! VARIABLE DEFINITION
IMPLICIT NONE
REAL, DIMENSION(1:7) :: VA = (/1.2,2.3,3.4,4.5,5.6,6.7,7.8/)
REAL, DIMENSION(1:7) :: VA1 = 0.0, VA2 = 0.0
INTEGER I
!
VA1 = VA
VA2 = VA
!
DO I = 2, 7
VA1(I) = VA1(I) + VA1(I-1)
ENDDO
!
VA2(2:7) = VA2(2:7) + VA2(1:6)
!
! Previous two operations with VA1 and VA2 seem that
! should provide the same result. Which is not the case.
PRINT*, VA1
PRINT*, VA2
!
! To obtain the same effect without an explicit DO loop we can do
! the following
VA2 = VA
VA2(2:7) = (/ (SUM(VA2(1:I)), I = 2,7) /)
!
PRINT*, VA1
PRINT*, VA2
END PROGRAM ex_4_5
```

### 4.3.6 excode_4_6.f90

```     PROGRAM ex_4_6
!
! DEFINITION OF VARIABLES
IMPLICIT NONE
INTEGER, DIMENSION(1:3,1:3) :: A = RESHAPE( (/ 1,2,3,4,5,6,7,8,9 /), (/ 3,3 /) )
!
!
!      1  4  7
!  A = 2  5  8
!      3  6  9
!
PRINT*, "Matrix Element", A(2,3)
PRINT*, "Submatrix", A(1:2,2:3)
PRINT*, "Submatrix", A(::2,::2)
PRINT*, "Matrix Column", A(:,3)
PRINT*, "Matrix Row", A(2,:)
PRINT*, "Full Matrix", A
PRINT*, "Transposed Matrix", TRANSPOSE(A)
END PROGRAM ex_4_6
```

[ previous ] [ Contents ] [ 1 ] [ 2 ] [ 3 ] [ 4 ] [ 5 ] [ 6 ] [ 7 ] [ 8 ] [ 9 ] [ 10 ] [ 11 ] [ 12 ] [ 13 ] [ next ]

`Fortran 90` Lessons for Computational Chemistry

0.0

Curro Pérez-Bernal `mailto:francisco.perez@dfaie.uhu.es`