! matmult: multiplies two matrices and prints the result to standard output

PROGRAM matmult
    IMPLICIT none

    INTEGER, PARAMETER :: MAXN = 100, MAXM = 100, MAXL = 100
    REAL, DIMENSION(MAXN,MAXM) :: mat
    REAL, DIMENSION(MAXM,MAXL) :: nat

    INTEGER :: n, m, l, i, j

    WRITE(*,'("Enter n,m,l:")',ADVANCE='NO')
    READ(*,*)  n,m,l

    IF(n > MAXN .OR. m > MAXM .OR. l > MAXL)  THEN
        WRITE(0,*)  "n, m, or l is too large"
        STOP
    END IF

    WRITE(*,'("Enter ", I2, "x", I2, " matrix:")')  n, m
    READ(*,*)  ((mat(i,j), j = 1,m), i = 1,n)
    WRITE(*,'("Enter ", I2, "x", I2, " matrix:")')  m, l
    READ(*,*)  ((nat(i,j), j = 1,l), i = 1,m)

    CALL domult(n, m, l, mat(1:n,1:m), nat(1:m,1:l))

END PROGRAM matmult

! domult: subroutine to do the multplication

SUBROUTINE domult(n, m, l, mat, nat)
    IMPLICIT none
    INTEGER, INTENT(IN)                 :: n, m, l
    REAL,    INTENT(IN), DIMENSION(n,m) :: mat
    REAL,    INTENT(IN), DIMENSION(m,l) :: nat

    REAL :: dot
    INTEGER :: i, j
    REAL, DIMENSION(n,l) :: mnat

    DO i = 1,n
        DO j = 1,l
            mnat(i,j) = dot(m, mat(i,:), nat(:,j))
        END DO
    END DO

    WRITE(*,'("The product of ",I2,"x",I2," matrix:")') n, m
    WRITE(*,*)  ((mat(i,j), j = 1,m), i = 1,n)
    WRITE(*,'("With ", I2, "x", I2, " matrix:")') m, l
    WRITE(*,*)  ((nat(i,j), j = 1,l), i = 1,m)
    WRITE(*,'("is the ", I2, "x", I2, " matrix:")'), n, l
    WRITE(*,*)  ((mnat(i,j), j = 1,n), i = 1,l)
END SUBROUTINE domult

! dot: function to compute a dot product

REAL FUNCTION dot(n, v, w) RESULT(dnv)
    IMPLICIT none
    INTEGER, INTENT(IN)               :: n
    REAL,    INTENT(IN), DIMENSION(n) :: v, w

    INTEGER :: i

    dnv = 0
    DO i = 1,n
        dnv = dnv + v(i) * w(i)
    END DO
END FUNCTION dot