program mmult
c
c -- matrix multiplication
c -- figure 3.1
c
      integer out, nrow, ncol, maxr, maxc
      real x(9,3), y(9), a(3,3), g(3)
      common /inout/ out
c
      out = 6
      maxr = 9
      maxc = 3
      call input(x, y, nrow, ncol)
      call square(x, y, a, g, nrow, ncol, maxr, maxc)
      call output(x, y, a, g, nrow, ncol)
      stop
      end
      subroutine input(x, y, nrow, ncol)
c -- get values for nrow, ncol, and arrays x, y
c
      integer i, j, nrow, ncol
      real x(9,3), y(9)
c
      nrow = 5
      ncol = 3
      do 20 i = 1, nrow
       x(i,1) = 1.0
       do 10 j = 2, ncol
        x(i,j) = i * x(i, j-1)
10     continue
       y(i) = 2 * i
20    continue
      return
      end
      subroutine output(x, y, a, g, nrow, ncol)
c
c -- print out the answers
c
      integer out, nrow, ncol, i, j
      real x(9,3), y(9), a(3,3), g(3)
      common /inout/ out
c
      write(out, 101)
      do 10 i = 1, nrow
       write(out, 102) (x(i,j), j = 1, ncol), y(i)
10    continue
      write(out, 103)
      do 20 i = 1, ncol
       write(out, 102) (a(i,j), j = 1, ncol), g(i)
20    continue
      write(out, 104)
      return
101   format('          x                  y')
102   format(3f8.0, '  : ', f7.0)
103   format('          a                   g')
104   format(/)
      end
      subroutine square(x, y, a, g, nrow, ncol, n, m)
c
c -- matrix multiplication routine
c -- a = transpose x times x
c -- g = y times x
c
      integer nrow, ncol, i, k, l, n, m
      real x(n,m), y(n), a(m,m), g(m)
c
      do 40 k = 1, ncol
       do 20 l = 1, k
        a(k,l) = 0.0
        do 10 i = 1, nrow
         a(k,l) = a(k,l) + x(i,l) * x(i,k)
         if (k .ne. l) a(l,k) = a(k,l)
10      continue
20     continue
       g(k) = 0.0
       do 30 i = 1, nrow
        g(k) = g(k) + y(i) * x(i,k)
30     continue
40    continue
      return
      end