[ previous ] [ Contents ] [ 1 ] [ 2 ] [ 3 ] [ 4 ] [ 5 ] [ 6 ] [ 7 ] [ 8 ] [ 9 ] [ 10 ] [ 11 ] [ 12 ] [ 13 ] [ next ]
Fortran 90 Lessons for Computational Chemistry
The main aims of this lesson are the following:
presenting storage ordering of multidimensional arrays.
presenting how to manipulate whole matrices or arrays sections in
Fortran.
matrix definition using the WHERE statement.
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 )
FORTRAN allows to define multidimensional arrays, being seven is
the max number of indices. The code excode_4_2.f90, Section 4.3.2 an array is
fully characterized making use of several inquiry type functions (see
excode_8_2.f90, Section
8.3.2).
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.
Assigning a particular value to the full array:
V1 = 0.5
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
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.
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
....
ELSEWHERE (mask_expr_2) [name]
....
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.
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: '
READ *, num_days
! Dynamic storage definition
ALLOCATE(t_worked(1:num_days))
!
PRINT *,' Worked hours per day in ', num_days, ' days.'
! I/O
READ *, t_worked
!
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
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
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
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
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
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
mailto:francisco.perez@dfaie.uhu.es