!  Data-parallel Gaussian elimination  and associated
!  routines for smaller problems, using N^2 "virtual
!  processors".  Note that load balancing is impossible
!  in this case.  The programs expect right-hand sides
!  and solutions to be two-dimensional arrays.  In
!  general, argument arrays are understood to be
!  Nx by Nx, with the coefficient matrix being N by N,
!  where N <= Nx.  If N < Nx, the first N rows (and
!  columns) of A contain the important data.
!
!  This version prints at each step of the process.  It is intended
!  only for very small examples, up to 8x8.
!
!  Source file:  lineqd.fcm     Header with array dimension:  lineqn.hp
!
!    Nx  -  Actual array dimension - defined in lineqn.hp
!
!  Usual argument conventions:
!
!    N   -  Problem size, N <= Nx
!    A   -  Input matrix, becomes "history"
!    B   -  Right hand side, becomes solution
!    P   -  Vector giving permutation, using N of Nx entries

!  A "canonical example" to clarify the relationships among the routines
!  Gauss, GaussX, and Rowperm is given below, following Rowperm.


! ------ subroutine Gauss(A, B, P, N) --------------

!  Transform A to "history" while transforming B to
!  [A^-1].[B] (matrix product), though with rows
!  permuted.  P records the permutation used, so row
!  P(r) is the r-th pivot row.  Problem is taken to have
!  N rows, with N <= Nx, but all columns of B are
!  transformed.  Uses partial pivoting without moving
!  data.  Prints error message to standard output and
!  exits with partial results if any pivot value has
!  |pivot| < 1.0E-20 .

      subroutine Gauss(A, B, P, N)
      INCLUDE 'lineqn.hp'
      integer N, P(Nx)
      real, dimension(Nx,Nx) :: A, B

      INTERFACE
        subroutine prmf(A, ra, rb, ca, cb)
        INCLUDE 'lineqn.hp'
        real, dimension(Nx, Nx) :: A
        integer ra, rb, ca, cb
      END INTERFACE

      INTERFACE
        subroutine prmi(A, ra, rb, ca, cb)
        INCLUDE 'lineqn.hp'
        integer, dimension(Nx, Nx) :: A
        integer ra, rb, ca, cb
      END INTERFACE

      real, parameter :: eps = 1.0E-20
      integer, dimension(Nx,Nx) :: I, J, Row
      real, dimension(Nx, Nx) :: Lj, Uj
      integer r, rp, mxl(2)
      real amx

      I = spread([1:Nx], 2, Nx)	! Coordinates
      J = spread([1:Nx], 1, Nx)
      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 pivot ")
          print 22, amx, r, rp
 22       FORMAT (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
	A(rp, r+1:N) = A(rp, r+1:N)/amx
	B(rp, :) = B(rp, :)/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
        print *, "Iteration ", r, ",  rp = ", rp
        print *, "Row:";  call prmi(Row, 1, N, 1, N)
        print *, "A:";    call prmf(A,   1, N, 1, N)
        print *, "B:";    call prmf(B,   1, N, 1, N)
      end do ! r = 1, N-1
      mxl = maxloc(Row)
      rp  = mxl(1)		! Last pivot
      P(N) = rp
      amx = A(rp, N)
      if ( abs(amx) < eps ) then
        print 21
        print 22, amx, r, rp
        go to 99        ! Exit
      end if
      B(rp, :) = B(rp, :)/amx
      Lj = spread(A( :, N), 2, Nx)
      Uj = spread(B(rp, :), 1, Nx)
      where ( Row <> N ) B = B - Lj*Uj	! Update B's
      print *, "Iteration ", N, ",  rp = ", rp
      print *, "A:";    call prmf(A,   1, N, 1, N)
      print *, "B:";    call prmf(B,   1, N, 1, N)      
 99   continue    ! in case of failure
      where ( I > N ) B = 0.0		! For safety
      end subroutine Gauss


! ------ subroutine GaussX(A, B, P, N) ---------------------------

!  A should be a "history" produced by Gauss from an original matrix A0.
!  Transforms B to [A0^-1].[B] (leaving rows out of order) using history
!  of elimination left by Gauss in A, P.

      subroutine GaussX(A, B, P, N)
      INCLUDE 'lineqn.hp'
      integer N, P(Nx)
      real, dimension(Nx,Nx) :: A, B

      INTERFACE
        subroutine prmf(A, ra, rb, ca, cb)
        INCLUDE 'lineqn.hp'
        real, dimension(Nx, Nx) :: A
        integer ra, rb, ca, cb
      END INTERFACE

      integer, dimension(Nx,Nx) :: I
      real,    dimension(Nx,Nx) :: Lj, Uj
      integer r, rp
      real amx

      I = spread([1:Nx], 2, Nx)	! Coordinates
      print *, "A (history):";  call prmf(A, 1, N, 1, N)
      do, r = 1, N
	rp = P(r)		! Pivot row (physical)
	amx = A(rp, r)		! Pivot value
	B(rp, :) = B(rp, :)/amx
        Lj = spread(A( :, r), 2, Nx)
        Uj = spread(B(rp, :), 1, Nx)
        where ( I <> rp ) B = B - Lj*Uj ! Update B's
        print *, "Iteration ", r, ",  rp = ", rp
        print *, "B:";  call prmf(B, 1, N, 1, N)
      end do ! r = 1, N-1
      where ( I > N ) B = 0.0		! For safety
      end subroutine GaussX


! ------ subroutine Rowperm(A, B, P N) ------------

!  Form A by permuting rows of B according to P, for
!  effectively N by N matrices A, B.

      subroutine Rowperm(A, B, P, N)
      INCLUDE 'lineqn.hp'
      integer N, P(Nx)
      real, dimension(Nx,Nx) :: A, B

      integer r

      A = B(P, :)
      end subroutine Rowperm


! ------ A canonical example for Gauss, GaussX, Rowperm -----------------
!
!  The example uses a number of distinct arrays so that the relationships
!  among them can be kept clear.  Most uses of these routines will need
!  many fewer arrays.  We assume no pivot failures.
!
!	real, dimension(Nx,Nx) :: A0, Ah, B1, B2, Xp, X1, X2
!	integer, dimension(Nx) :: P
!
!	!{ 1 <= N <= Nx, A0, B1, B2 have initial data }
!
!	Ah = A0				! Coefficients
!	Xp = B1				! Right-hand side (RHS)
!	call Gauss(Ah, Xp, P, N)
!
!	!{ Xp solves A0.X = B1, except rows are out of order, }
!	!{ P gives correct order for rows, Ah has history.    }
!
!	call Rowperm(X1, Xp, P, N)
!
!	!{ A0(1:N,1:N).X1(1:N,1:Nx) = B1(1:N,1:Nx) except for }
!	!{ rounding errors, Ah, P retain history.             }
!
!	Xp = B2
!	call GaussX(Ah, Xp, P, N)
!
!	!{ Xp solves A0.X= B2, except rows are out of order.  }
!
!	call Rowperm(X2, Xp, P, N)
!
!	!{ A0(1:N,1:N).X2(1:N,1:Nx) = B2(1:N,1;Nx) except for }
!	!{ rounding errors, Ah, P still retain history.       }


! ------ function GaussI(X, N) --------------------------------------

!  Returns explicit inverse of N by N matrix X computed
!  using Gauss.  X is not altered, history is lost.

      function GaussI(X, N)
      INCLUDE 'lineqn.hp'
      integer N
      real, dimension(Nx,Nx) :: X, GaussI

      INTERFACE
	subroutine Gauss(A, B, P, N)
        INCLUDE 'lineqn.hp'
	integer N, P(Nx)
	real, dimension(Nx,Nx) :: A, B
      END INTERFACE

      INTERFACE
	subroutine Rowperm(A, B, P, N)
        INCLUDE 'lineqn.hp'
	integer N, P(Nx)
	real, dimension(Nx,Nx) :: A, B
      END INTERFACE

      integer, dimension(Nx,Nx) :: I, J
      real,    dimension(Nx,Nx) :: A, B

      integer P(Nx)

      I = spread([1:Nx], 2, Nx)		! Coordinates
      J = spread([1:Nx], 1, Nx)
      A = X			! So input isn't altered
      B = 0.0			! Make identity
      where ( I == J .and. I <= N) B = 1.0
      call Gauss(A, B, P, N)
      call Rowperm(A, B, P, N)
      GaussI = A
      end function GaussI
