How to Find the Product of Two 3x3 Matrices

Fortran has an inbuilt function to calculate the matrix product of two matrices. This function is MATMUL( MatrixA,MatrixB). In this post, we will study the simple logic behind how this works. For simplicity, we choose two square matrices  of order 3, so that we do not have to bother about whether the multiplication is feasible. Let us see how this works.

• First, we read inputs for matrices A and B from their respective input files and display them at the terminal. We need to create an output matrix AB of the same dimension.
• We need to assign values to each element of AB. Therefore, we need two DO loops i=1 to n and j=1 to n. For non-square matrices of different orders, this is slightly different.
• Now we need one more loop for calculating elements of AB. Before the loop, we initialize the sum, s=0. Inside the loop, we will take the sum of each A(i,k)+B(k,j) for k=1 to 3. Therefore, we need another DO loop inside.
• After the END DO statement of the innermost loop, we will have calculated the sum, so we define AB(i,j)=s. Note that the s=0 statement used earlier automatically resets the sum for calculating the next element in a row.
• We then display the result. Additionally, we can also display the result using the intrinsic function, MATMUL.
• We can use a flag to check if both the results are the same. If we set the flag to 1 initially, and check each element of AB against an element of the result, C given by the intrinsic function. If the elements do not match we change  flag=0, exit and display that the results do not match.
• If everything goes well we retain the value of the flag=1 and display that the results match.

 Let us illustrate matrix multiplication using two  3x3 matrices.

Matrix Multiplication of A(3,3) and B(3,3)

Now, for the Fortran program,

  !To find the product of two matrices
program matrix_3x3
    implicit none
    integer i,j,k,flag
    real A(3,3),B(3,3),AB(3,3),C(3,3),s
  
    open (unit=1,file='matrix1.dat') !To read matrix A
    write (*,*) "matrix A"
    read (1,*) A
  
    do i=1,3 !To display matrix A
      write (*,*) (A(i,j),j=1,3)
    end do
    
    open (unit=1,file='matrix2.dat') !To read matrix B
    write (*,*) "matrix B"
    read (1,*) B
    
    do i=1,3 !To display matrix B
        write (*,*) (B(i,j),j=1,3)
    end do

    !To calculate the result
    c=0
    do i=1,3 !To change rows
        do j=1,3 !To change columns
            s=0 !To reset the sum for the next element in the row
              do k=1,3 !Notice why the number of columns of A have to be equal to the number of rows of B
                s=A(i,k)*B(k,j)+s
            end do
             AB(i,j)=s !To assign values to each element of the product
        end do
    end do
       
    write (*,*) "matrix AB (the result)" !To display the result
    do i=1,3
        write (*,*) (AB(i,j),j=1,3)
    end do
    
    C=matmul(A,B) !To perform matrix multiplication by intrinsic function. Optional
    write(*,*) "The product of matrices A and B by intrinsic function"
    do i=1,3
        write(*,*) (C(i,j),j=1,3)
    end do
      
    flag=1 !To check whether AB=C
    do i=1,3
          do j=1,3
            if (AB(i,j).ne.C(i,j)) then
                  flag=0
                  exit
              end if
        end do
    end do

    if (flag==1) then !To display if AB=C or not    
        write (*,*) "The results match"
        else
        write(*,*) "The results do not match. Recheck the program"
    end if
      

end program matrix_3x3

Sample Output for Matrix Multiplication of two 3x3 matrices