Analysis of a FORTRAN Matrix Program - III

The code, as a reminder:

! program to calculate whether a graph is strongly connected
! it calculates this my using matrix multiplication on the
! adjacency matrix representing the graph
PROGRAM matrix

INTEGER :: length
PARAMETER (length = 4)
INTEGER :: a,b, c, m(length,length), copy(length,length), total, d, e, f

m(1, :) = (/ 1, 0, 1, 0/)
m(2, :) = (/ 1, 0, 0, 0/)
m(3, :) = (/ 1, 1, 0, 1/)
m(4, :) = (/ 0, 0, 1, 0/)

!initialize a copy of m so that we may calculate N^2
copy = m

DO a = 1, length
DO b = 1, length
total = 0
DO c = 1, length
total = m(a,c) * copy(c,b) + total
END DO
IF (total > 0) THEN
m(a,b) = 1
ELSE
m(a,b) = 0
END IF
END DO
END DO

DO f = 1, length
PRINT *, m(f,1:4)
END DO

DO e = 1, length
total = 0
DO d = 1, length
IF( m(d,e) .EQ. 1 .AND. m(e,d) .EQ. 1 ) total = total + 2
END DO
IF ( total == 8 ) THEN
PRINT *, 'Strongly Connected At: ', e
ELSE
PRINT *, 'Not Strongly Connected At: ', e
END IF
END DO

END PROGRAM matrix

Now, while the example I used in class had only 1s and 0s in the result, in reality we may have higher numbers. This is fine. A zero in position (i,j) means that there is no path from i to j, but a number greater than zero, such as 1 or 3 or 8, means that there is a path. So we need not check specifically for a nonzero and set the result to 1.

One we know this, we may use the built-in matmul function that I suggested you use, rather than doing the matrix multiplication by hand. Doing matmul it by hand is certainly still good programming exercise, but if we have this function built in, it may just be incredibly optimized. Further, it saves programmer work. Further, it makes the code smaller and thus more readable and less error-prone.

There was actually a bug in this code, which I should comment upon before proceeding:

DO a = 1, length
DO b = 1, length
    total = 0

        ! the following code calculates a single dotproduct
    DO c = 1, length
        total = m(a,c) * copy(c,b) + total
    END DO

        ! and here we assign the result to m
    IF (total > 0) THEN
        m(a,b) = 1
    ELSE
        m(a,b) = 0
    END IF
END DO
END DO

The problem is twofold. Firstly, we are only doing a single matrix multiplication in this entire code - calculating M^2 rather than M^4. Secondly, we are not even calculating M^2 correctly! This is because we are modifying M in the middle of the multiplication! M(1,1) is being set, and then that value is being used to compute M(1,2), M(1,3), M(1,4), and M(2,1), M(3,1), and M(4,1). We could solve this by making a third matrix, result, and assigning to that matrix, and then stating m = result at the end. But let us skip this entire step and just use the built in matmul function. We no longer need the copy matrix for this purpose, but I will leave it in in case we need it for some other purpose. I will change the code, commenting out the previous code so we can see what changed. I will still only calculate M^2. We will deal with M^4 and the associated problem later.

! program to calculate whether a graph is strongly connected
! it calculates this my using matrix multiplication on the
! adjacency matrix representing the graph
PROGRAM matrix

INTEGER :: length
PARAMETER (length = 4)
INTEGER :: a,b, c, m(length,length), copy(length,length), total, d, e, f

m(1, :) = (/ 1, 0, 1, 0/)
m(2, :) = (/ 1, 0, 0, 0/)
m(3, :) = (/ 1, 1, 0, 1/)
m(4, :) = (/ 0, 0, 1, 0/)

!initialize a copy of m so that we may calculate N^2
copy = m

!DO a = 1, length
! DO b = 1, length
!  total = 0
!
!  DO c = 1, length
!   total = m(a,c) * copy(c,b) + total
!  END DO
!
!  IF (total > 0) THEN
!   m(a,b) = 1
!  ELSE
!   m(a,b) = 0
!  END IF
! END DO
!END DO

! calculate M^2
m = matmul(m, m)

DO f = 1, length
PRINT *, m(f,1:4)
END DO

DO e = 1, length
total = 0
DO d = 1, length
IF( m(d,e) .EQ. 1 .AND. m(e,d) .EQ. 1 ) total = total + 2
END DO
IF ( total == 8 ) THEN
PRINT *, 'Strongly Connected At: ', e
ELSE
PRINT *, 'Not Strongly Connected At: ', e
END IF
END DO

END PROGRAM matrix