C  Data-parallel LU decomposition for smaller problems.
C  Note that load balancing is impossible in this case.
C  The program expects to transform an array of right-hand
C  sides at the same time it carries determines the LU
C  decomposition.
C
C  This version uses cannon.o for matrix multiplication
C  rather than the matmul intrinsic function.
C
C  File:  LUcan.f
C
C  Common information in /LUFE/, /LUPE/ is described in LU.h:
C
C    Nx     -  Array dimension
C    eps    -  Value too small for pivot
C    I, J   -  "physical" array coordinates

C  Transform A to LU decomposition while transforming
C  B to [L^-1][B] (matrix product).  `Row' is left
C  with logical row number in all entries.  P records
C  permuatation, so Row(P(i),:) == i .
C  Problem is taken to have N rows, with N <= Nx.
C  Uses partial pivoting.

      subroutine LUd(A, B, Row, P, N)
      INCLUDE 'LU.h'
      real, array(Nx,Nx) :: A, B
      integer, array(Nx,Nx) :: Row
      integer N, P(Nx)
CMPF  ONDPU A, B, Row
CMPF  ONFE P

      real, array(Nx, Nx) :: Lj, Uj
      integer r, rp, mxl(2)
      real amx

      Row = I		! Unpermuted ordering
      where ( Row > N ) Row = 0
      do, r = 1, N-1
        mxl = maxloc(abs(A), J == r .and. Row >= r)
        rp = mxl(1)
        P(r) = rp            ! Record permutation
        amx = A(rp, r)       ! Pivot value
        if ( abs(amx) < eps ) then
          print 21
 21       FORMAT (" Matrix appears to be singular.  Failed with")
          print 22, amx, r, rp
 22       FORMAT (" pivot ", E14.5, " at step ", I4, "row ", I4)
          go to 99        ! Exit
        end if
        where ( Row == r) Row = Row(rp, 1)
        Row(rp, :) = r    ! Swap rows implicitly
D        print *, "Row(1:N, 1:N) ="
D        call prmi(Row, 1, N, 1, N)
        where ( Row(:, r) > r ) A(:, r) = A(:, r)/amx
        Lj = spread(A(:, r), 2, Nx)
        Uj = spread(A(rp, :), 1, Nx)
        where ( Row(:, r+1:N) > r )  ! Update A's
     *     A(:, r+1:N)  
     *        = A(:, r+1:N) - Lj(:, r+1:N)*Uj(:, r+1:N)
        Uj = spread(B(rp, :), 1, Nx)
        where ( Row > r ) B = B - Lj*Uj   ! Update B's
D        print *, "A(1:N, 1:N) ="
D        call prmf(A, 1, N, 1, N)
D        print *, "B(1:N, 1:n) ="
D        call prmf(B, 1, N, 1, N)
      end do ! r = 1, N-1
      mxl = maxloc(Row)
      P(N) = mxl(1)
 99   continue    ! in case of failure
      end subroutine LUd

C  Carry out back-substitution to obtain X.  A, B, Row, P
C  as left by LUd.  X(j, 1:N) is the solution corresponding
C  to B(1:N, j)

      subroutine LUbk(X, A, B, Row, P, N)
      INCLUDE 'LU.h'
      real, array(Nx,Nx) :: X, A, B
      integer, array(Nx,Nx) :: Row
      integer N, P(Nx)
CMPF  ONDPU X, A, B, Row
CMPF  ONFE P

      real, array(Nx, Nx) :: Ujj
      integer r, rp

      Ujj = minval(A)             ! Ujj will have "diagonal"
      where ( Row == J ) Ujj = A  ! entries spread across cols.
      Ujj = spread(maxval(Ujj, 1), 1, Nx)
D      print *, "Ujj(1:8, 1:6) ="
D      call prmf(Ujj, 1, 8, 1, 6)
      X(:, N) = B(P(N), :)/Ujj(:, N)
      do, r = N-1, 1, -1
        X(:, r) =
     *    ( B(P(r), :)
     *     - SUM(X(:, r+1:N)*SPREAD(A(P(r), r+1:N), 1, NX), 2) )
     *       / Ujj(:, r)
      end do
      end subroutine LUbk

C  Construct random matrix with diagonal D + random

      function rmat(N, D)
      INCLUDE 'LU.h'
      real, array(Nx, Nx) :: rmat
      integer N
      real D
CMPF  mpl MPF_RANDOM, MPF_SRANDOM
      real, parameter :: Two31 = 2147483648.0
      integer r, RR(Nx)

      call mpf_srandom([1:Nx] + 118)
      rmat = 0.0
      do, r = 1, N
        call mpf_random(RR)
        rmat(r, 1:N) = RR(1:N)/Two31 - 0.5
      end do
      where ( I == J .and. I <= N)  rmat = rmat + D
      end function rmat

C  Print real A(ra:rb, ca:cb), at most 6 columns

      subroutine prmf(A, ra, rb, ca, cb)
      INCLUDE 'LU.h'
      real, array(Nx, Nx) :: A
      integer ra, rb, ca, cb
CMPF  ONDPU A
      integer ir

      if ( cb >= ca+6 ) cb = ca+5
      do, ir = ra, rb
       print 21, A(ir, ca:cb)
 21    FORMAT ( 6(1X, E11.4:) )
      end do
      end subroutine prmf

C  Print integer A(ra:rb, ca:cb), at most 8 columns

      subroutine prmi(A, ra, rb, ca, cb)
      INCLUDE 'LU.h'
      integer, array(Nx, Nx) :: A
      integer ra, rb, ca, cb
CMPF  ONDPU A
      integer ir

      if ( cb >= ca+8 ) cb = ca+7
      do, ir = ra, rb
       print 21, A(ir, ca:cb)
 21    FORMAT ( 8(1X, I4:) )
      end do
      end subroutine prmi

C  Elementary test of LU routines

      program LUtest
      INCLUDE 'LU.h'

      INTERFACE
        subroutine LUd(A, B, Row, P, N)
        INCLUDE 'LU.h'
        real, array(Nx,Nx) :: A, B
        integer, array(Nx,Nx) :: Row
        integer N, P(Nx)
CMPF    ONDPU A, B, Row
CMPF    ONFE P
      END INTERFACE

      INTERFACE
        subroutine LUbk(X, A, B, Row, P, N)
        INCLUDE 'LU.h'
        real, array(Nx,Nx) :: X, A, B
        integer, array(Nx,Nx) :: Row
        integer N, P(Nx)
CMPF    ONDPU X, A, B, Row
CMPF    ONFE P
      END INTERFACE
      INTERFACE
        function rmat(N, D)
        INCLUDE 'LU.h'
        real, array(Nx, Nx) :: rmat
        integer N
        real D
      END INTERFACE
      INTERFACE
        subroutine prmf(A, ra, rb, ca, cb)
        INCLUDE 'LU.h'
        real, array(Nx, Nx) :: A
        integer ra, rb, ca, cb
CMPF    ONDPU A
      END INTERFACE
      INTERFACE
        subroutine prmi(A, ra, rb, ca, cb)
        INCLUDE 'LU.h'
        integer, array(Nx, Nx) :: A
        integer ra, rb, ca, cb
CMPF    ONDPU A
      END INTERFACE
      INCLUDE 'cannon.h'

      real, array(Nx,Nx) :: Xx, A0, B0, A1, B1, C0, C1
      integer, array(Nx,Nx) :: Rw
      integer N, Pm(Nx), eltm, ip, mmtm
CMPF  ONDPU Xx, A0, B0, A1, B1, C0, C1, Rw
CMPF  ONFE  Pm
      real D, err, mxa, mxx
      logical more
      character*50 buff

      print *, "Maximum N is ", Nx
      I = spread([1:Nx], 2, Nx)
      J = spread([1:Nx], 1, Nx)
      print *, "Specify N <= 0 to quit."
      more = .TRUE.
      do while ( more )   ! Test specified value of N
        print '(" N, D > ", $)'
        read *, N, D
        if ( N <= 0 ) go to 99
        A0 = rmat(N, D)
        A1 = A0
D        print *, "A0(1, 8, 1, 6) ="
D        call prmf(A0, 1, 8, 1, 6)
        B0 = 0.0
        Xx = 0.0
        C0 = 0.0
        where ( I == J .and. I <= N) B0 = 1.0
        B1 = B0    ! Identity matrix
        call mpTimerStart()
        call LUd(A1, B1, Rw, Pm, N)
        call LUbk(Xx, A1, B1, Rw, Pm, N)
        eltm = mpTimerElapsed()
        print 21, eltm/1000.0
 21     FORMAT (" Finished in ", F8.3, " sec.")
        print *, "Pm ="
        do, ip = 1,N,8
          print 22, Pm(ip),Pm(ip+1),Pm(ip+2),Pm(ip+3),
     *            Pm(ip+4),Pm(ip+5),Pm(ip+6),Pm(ip+7)
 22       FORMAT ( 8(1X, I4:) )
        end do
D        print *, "A1(1:N, 1:N) ="
D        call prmf(A1, 1, N, 1, N)
D        print *, "B1(1:N, 1:N) ="
D        call prmf(B1, 1, N, 1, N)
D        print *, "Rw(1:N, 1:N) ="
D        call prmi(Rw, 1, N, 1, N)
D        print *, "Xx(1:N, 1:N) ="
D        call prmf(Xx, 1, N, 1, N)
        call MpTimerStart()
        call Cansqc(C0, A0, transpose(Xx), N)
        mmtm = mpTimerElapsed()
        print 31, mmtm
 31     FORMAT (" Cannon mmpy used ", I6, " msec.")
D        print *, "C0(1:N, 1:N) ="
D        call prmf(C0, 1, N, 1, N)
        C1 = C0
        where ( I == J .and. I <= N) C1 = C1 - 1.0
        mxa = maxval(abs(A0))
        mxx = maxval(abs(Xx))
        err = maxval(abs(C1))
        print 23, mxa, mxx, err
 23     FORMAT (" max(A) = ", E11.4, "   max(X) = ", E11.4,
     *          "   max err = ", E11.4)

      end do
 99   continue
      end program LUtest
