program Matrix_Arithmetic

  implicit NONE

  double precision A(7,8), VECTOR(7), B(7,7)
  double precision AA(56)
  equivalence (A, AA), (AA(50), VECTOR(1))

  logical :: DONE = .false.
  integer :: I, J, K

  do i = 1, 7
    do j = 1, 7
      b(i,j) = 1.0 / ( i + j - 1 )
      a(i,j) = b(i,j)
    end do
  end do

  print "('HILBERT MATRIX'/)"
  print 2, ((A(I,J), J = 1,7), I = 1,7)
2 format(1X,7E25.17)
  print "(/'INVERSE'/)"

  do
    do k = 1, 7
      vector(1) = 1.0
      vector(2:7) = 0.0
      a(1,2:8) = a(1,2:8) / a(1,1)
      aa(1:55) = aa(2:56)

      do i = 1, 6
        aa(56) = a(i,1)
        do j = 1, 7
          a(i,j) = a(i,j+1) - aa(56)*a(7,j)
        end do
      end do
    end do

    print 2, ((A(I,J), J = 1,7), I = 1,7)

    if ( done ) exit
    done = .true.

    print "(/'MATRIX PRODUCT'/)"

    do k = 1, 7
      do i = 1, 7
        vector(i) = dot_product(a(i,1:7),b(:,k))
      end do
      print "(1x,7f25.17)", vector
    end do
    print "(/'TWICE INVERTED'/)"
  end do

end program Matrix_Arithmetic
