## 13. Extra Stuff

This lesson just contains some extra material(s) to be discussed in the final lecture. The discriptions here are quite brief.

### Generic subprograms

In Fortran 77 many of the intrinsic functions (sin, cos, abs) return the same type as they receive as an argument. These are known as generic functions. There was no way in Fortran 77 to do the same for user routines. Here is an example of Fortran 90 generic subroutines from Fortran 90 for the Fortran 77 Programmer. This example defines a subroutine, SWAP, which works for reals, integers and characters. Note that while this example only defines a genereic subroutine, generic functions are written in the same way.

```MODULE SWAPPER
INTERFACE SWAP
MODULE PROCEDURE SWAP_R, SWAP_I, SWAP_C
END INTERFACE
CONTAINS

SUBROUTINE SWAP_R(A, B)
IMPLICIT NONE
REAL, INTENT (INOUT)      :: A, B
REAL                      :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_R

SUBROUTINE SWAP_I(A, B)
IMPLICIT NONE
INTEGER, INTENT (INOUT)   :: A, B
INTEGER                   :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_I

SUBROUTINE SWAP_C(A, B)
IMPLICIT NONE
CHARACTER, INTENT (INOUT) :: A, B
CHARACTER                 :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_C
END MODULE SWAPPER
```

Here is a simple program to use the SWAP subroutine:

```PROGRAM SWAP_MAIN
USE SWAPPER
IMPLICIT NONE
INTEGER                    :: I, J, K, L
REAL                       :: A, B, X, Y
CHARACTER                  :: C, D, E, F

I = 1   ;  J = 2         ;     K = 100 ; L = 200
A = 7.1 ;  B = 10.9      ;     X = 11.1; Y = 17.0
C = 'a' ;  d = 'b'       ;     E = '1' ; F = '"'

WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F
CALL  SWAP (I, J)  ;  CALL SWAP (K, L)
CALL  SWAP (A, B)  ;  CALL SWAP (X, Y)
CALL  SWAP (C, D)  ;  CALL SWAP (E, F)
WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F
END
```

Here is (long) an example which shows both modules procedures and operator overloading. The purpose of this example is to define a matrix type so that the following things work as they would in Matlab:

```     type(matrix), dimension(n,n) :: A, B, C
type(matrix), dimension(n) :: x, b

A = B * C
b = A * x
x = b / A
```

A normal array of reals in Fortran 90 would do addition, subtraction and assignment correctly. However, it would do multiplication and division element-wise as explained in class. The module below defines the matrix type and the =, +, -, * and / operators for it. Note that the operators are defined in the same manner as other generic subprograms.

```!       matrix.f90
!
!       Module for matrix type and operations
!       by Paul H. Hargrove
!       May 13, 1996
!
module operator

type matrix
real elem
end type matrix

interface assignment(=)
module procedure matrix_from_real, matrix_from_matrix, &
vector_from_real, vector_from_vector
end interface

interface operator(+)
end interface

interface operator(-)
module procedure matrix_sub, vector_sub
end interface

interface operator(*)
module procedure matrix_mul, vector_mul, matrix_vector_mul
end interface

interface operator(/)
module procedure matrix_div, matrix_vector_div
end interface

contains

!
! ASSIGNMENT OPERATORS: X = Y
!

subroutine matrix_from_real(X, Y)
! copy a 2D array of reals to a 2D array of type matrix
real, intent(in), dimension(:,:) :: Y
type(matrix), intent(out), dimension(size(Y,1),size(Y,2)) :: X

X(:,:)%elem = Y(:,:)
end subroutine matrix_from_real

subroutine matrix_from_matrix(X, Y)
! copy a 2D array of type matrix
type(matrix), intent(in), dimension(:,:) :: Y
type(matrix), intent(out), dimension(size(Y,1),size(Y,2)) :: X

X(:,:)%elem = Y(:,:)%elem
end subroutine matrix_from_matrix

subroutine vector_from_real(X, Y)
! copy a 1D array of reals to a 1D array of type matrix
real, intent(in), dimension(:) :: Y
type(matrix), intent(out), dimension(size(Y,1)) :: X

X(:)%elem = Y(:)
end subroutine vector_from_real

subroutine vector_from_vector(X, Y)
! copy a 1D array of type matrix
type(matrix), intent(in), dimension(:) :: Y
type(matrix), intent(out), dimension(size(Y,1)) :: X

X(:)%elem = Y(:)%elem
end subroutine vector_from_vector

!
! ADDITION OPERATORS: X = Y + Z
!

! add 2D arrays of type matrix
type(matrix), intent(in), dimension(:,:) :: Y
type(matrix), intent(in), dimension(size(Y,1),size(Y,2)) :: Z
type(matrix), dimension(size(Y,1),size(Y,2)) :: X

X(:,:)%elem = Y(:,:)%elem + Z(:,:)%elem

! add 1D arrays of type matrix
type(matrix), intent(in), dimension(:) :: Y
type(matrix), intent(in), dimension(size(Y,1)) :: Z
type(matrix), dimension(size(Y,1)) :: X

X(:)%elem = Y(:)%elem + Z(:)%elem

!
! SUBTRACTION OPERATORS: X = Y - Z
!

function matrix_sub(Y, Z) result(X)
! subtract 2D arrays of type matrix
type(matrix), intent(in), dimension(:,:) :: Y
type(matrix), intent(in), dimension(size(Y,1),size(Y,2)) :: Z
type(matrix), dimension(size(Y,1),size(Y,2)) :: X

X(:,:)%elem = Y(:,:)%elem - Z(:,:)%elem
end function matrix_sub

function vector_sub(Y, Z) result(X)
! subtract 1D arrays of type matrix
type(matrix), intent(in), dimension(:) :: Y
type(matrix), intent(in), dimension(size(Y,1)) :: Z
type(matrix), dimension(size(Y,1)) :: X

X(:)%elem = Y(:)%elem - Z(:)%elem
end function vector_sub

!
! MULTIPLICATION OPERATORS: X = Y * Z
!

function matrix_mul(Y, Z) result(X)
! multiply 2D arrays of type matrix
! NOTE: NAG's F90 demo won't deal w/ "half constrained" dimensions
type(matrix), intent(in), dimension(:,:) :: Y
type(matrix), intent(in), dimension(:,:) :: Z
type(matrix), dimension(size(Y,1),size(Z,2)) :: X

X(:,:)%elem = MATMUL(Y(:,:)%elem, Z(:,:)%elem)
end function matrix_mul

function vector_mul(Y, Z) result(X)
! multiply 1D arrays of type matrix
type(matrix), intent(in), dimension(:) :: Y
type(matrix), intent(in), dimension(size(Y,1)) :: Z
real X

X = DOTPRODUCT(Y(:)%elem, Z(:)%elem)
end function vector_mul

function matrix_vector_mul(Y, Z) result(X)
! multiply 2D array times 1D array of type matrix
type(matrix), intent(in), dimension(:,:) :: Y
type(matrix), intent(in), dimension(size(Y,2)) :: Z
type(matrix), dimension(size(Y,1)) :: X

X(:)%elem = MATMUL(Y(:,:)%elem, Z(:)%elem)
end function matrix_vector_mul

!
! DIVISION OPERATORS: X = Y/Z = INV(Z) * Y
!

function matrix_div(Y, Z) result(X)
! "divide" 2D arrays of type matrix
type(matrix), intent(in), dimension(:,:) :: Y
type(matrix), intent(in), dimension(:,:) :: Z
type(matrix), dimension(size(Y,1),size(Y,2)) :: X
real, dimension(size(Z,1),size(Z,2)) :: W
integer i, j, k, n

! copy arguments so they aren't modified
W(:,:) = Z(:,:)%elem
X(:,:)%elem = Y(:,:)%elem

! perform Gauss elimination on augmented matrix (W|X)
n = size(Z,2)
do k = 1,n-1
do i=k+1,n
W(i,k) = W(i,k)/W(k,k)
X(i,:)%elem = X(i,:)%elem - W(i,k) * X(k,:)%elem
end do
do j=k+1,n
do i=k+1,n
W(i,j) = W(i,j) - W(i,k) * W(k,j)
end do
end do
end do

! perform back substitution on X
do k = n,1,-1
X(k,:)%elem = X(k,:)%elem / W(k,k)
do i=1,k-1
X(i,:)%elem = X(i,:)%elem - W(i,k) * X(k,:)%elem
end do
end do
end function matrix_div

function matrix_vector_div(Y, Z) result(X)
! "divide" 1D array by 2D array of type matrix
type(matrix), intent(in), dimension(:) :: Y
type(matrix), intent(in), dimension(:,:) :: Z
type(matrix), dimension(size(Y,1)) :: X
real, dimension(size(Z,1),size(Z,2)) :: W
integer i, j, k, n

! copy arguments so they aren't modified
W(:,:) = Z(:,:)%elem
X(:)%elem = Y(:)%elem

! perform Gauss elimination on augmented matrix (W|X)
n = size(Z,2)
do k = 1,n-1
do i=k+1,n
W(i,k) = W(i,k)/W(k,k)
X(i)%elem = X(i)%elem - W(i,k) * X(k)%elem
end do
do j=k+1,n
do i=k+1,n
W(i,j) = W(i,j) - W(i,k) * W(k,j)
end do
end do
end do

! perform back substitution on X
do k = n,1,-1
X(k)%elem = X(k)%elem / W(k,k)
do i=1,k-1
X(i)%elem = X(i)%elem - W(i,k) * X(k)%elem
end do
end do
end function matrix_vector_div

end module operator
```