C  Cannon's algorithm for square matrices
C  Requires enclosing system to define
C           Nx - actual array dimension
C
C  File:  cannon.f

C  This is set up for use with LU.f, but is easily
C  adapted for other contexts.

C  Product of N by N matrices C := A B, both stored
C  in normal configuration.  Only the upper left
C  N by N submatrix of C is altered.  Initial
C  skewing is performed with CSHIFT, thinking this
C  will be fast on the DECmpp.  All arithmetic is
C  with REAL

      subroutine Cansqc(C, A, B, N)
      INCLUDE 'LU.h'
      real, array(Nx, Nx) :: C, A, B
      integer N
CMPF  ONDPU C, A, B

      real, array(Nx, Nx) :: Ax, Bx
      integer ix

C   Skew rows of A, columns of B
      Ax = A
      Bx = B
      do, ix = 2, N
        Ax(ix:N, 1:N) = CSHIFT(Ax(ix:N, 1:N), 1, 2)
        Bx(1:N, ix:N) = CSHIFT(Bx(1:N, ix:N), 1, 1)
      end do
C   Multiply
      C(1:N, 1:N) = Ax(1:N, 1:N)*Bx(1:N, 1:N)
      do, ix = 2, N
        Ax(1:N, 1:N) = CSHIFT(Ax(1:N, 1:N),-1, 2)
        Bx(1:N, 1:N) = CSHIFT(Bx(1:N, 1:N),-1, 1)
        C(1:N, 1:N) = C(1:N, 1:N)
     *               + Ax(1:N, 1:N)*Bx(1:N, 1:N)
      end do
      end subroutine Cansqc
