      PROGRAM MATRIX
      IMPLICIT NONE
      
      INTEGER DIM
      PARAMETER (DIM = 3) 
      REAL A(DIM, DIM)
      REAL B(DIM, DIM)
      REAL C(DIM, DIM)
      INTEGER I, J, K
      
      CALL READM(A, DIM)
      CALL READM(B, DIM)
       
      CALL MATRIXMUL(A, B, C, DIM)
       
      PRINT *, 'MATRIX A'
      CALL PRINTM(A, DIM)
      PRINT *, 'MATRIX B'
      CALL PRINTM(B, DIM)
      PRINT *, 'MATRIX C'
      CALL PRINTM(C, DIM)
      
      STOP
      END
      
      SUBROUTINE READM(MATRIX, DIM)
      INTEGER DIM
      REAL MATRIX(DIM, DIM)
      
      INTEGER I, J
      DO 10 I = 1, DIM
        READ *, (MATRIX(I, J), J = 1, DIM)
10    CONTINUE        
      RETURN
      END
      
      SUBROUTINE PRINTM(MATRIX, DIM)
      INTEGER DIM
      REAL MATRIX(DIM, DIM)
      
      INTEGER I, J
      DO 10 I = 1, DIM
        PRINT *, (MATRIX(I, J), J = 1, DIM)
10    CONTINUE        
      RETURN
      END
      
     
      SUBROUTINE MATRIXMUL(A, B, C, DIM)
      IMPLICIT NONE
      INTEGER DIM
      REAL A(DIM, DIM), B(DIM, DIM), C(DIM, DIM)
      
      INTEGER I, J, K
      REAL SUM
            
      DO 10 I = 1, DIM
        DO 20 J = 1, DIM
          SUM = 0.0
	  DO 30 K = 1, DIM
	    SUM = SUM + A(I, K) * B(K, J)
30        CONTINUE	    
          C(I, J) = SUM
20      CONTINUE
10    CONTINUE   

      RETURN
      END