!  Tests and basic timing for programs in lineqn.fcm.
!
!	Nx  -  array dimension -- defined in lineqn.hp
!
!  Frequently used arguments
!
!	N   -  matrix size (i.e., problem size)
!	I   -  coordinates, axis 1
!	J   -  coordinates, axis 2
!
!  All of these are determined at run time.
!
!  File:  gtest.fcm
!						Ernest Sibert, September 1996

!  Construct random matrix with diagonal D + random

      function rmat(N, D, I, J)
      INCLUDE 'lineqn.hp'
      integer, dimension(Nx, Nx) :: I, J
      real,    dimension(Nx, Nx) :: rmat
      integer N
      real D

      integer r
      real RR(Nx,Nx), RV(Nx*Nx)

      rmat = 0.0
      call random_number(RV)
      RR = RESHAPE(RV, (/ Nx, Nx /) )
      rmat(1:N, 1:N) = RR(1:N, 1:N) - 0.5
      where ( I == J .and. I <= N)  rmat = rmat + D
      end function rmat

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

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

      integer ir, cx

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

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

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

      integer ir, cx

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

!  Elementary test of Gauss routines

      program Gausstest
      INCLUDE 'lineqn.hp'
      INCLUDE 'lineqn.hI'	! Interfaces

      INTERFACE
        function rmat(N, D, I, J)
        INCLUDE 'lineqn.hp'
        integer, dimension(Nx, Nx) :: I, J
        real, dimension(Nx, Nx) :: rmat
        integer N
        real D
      END INTERFACE

      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,    dimension(Nx,Nx) :: Xx, A0, A1, B0, B1, B2
      real,    dimension(Nx,Nx) :: C0, C1
      integer, dimension(Nx,Nx) :: I, J
      integer N, Pm(Nx), seed(3), up, k

      real D, err, mxa, mxx
      logical more
      character*50 buff

      print '(" Three integer seeds for RNG > "$)'
      read *, seed(1:3)
      call random_seed(SIZE=N);    print *, "Seed size = ", N
      call random_seed(PUT=seed);
      I = SPREAD( (/ (k, k=1,Nx) /), 2, Nx)
      J = SPREAD( (/ (k, k=1,Nx) /), 1, Nx)
      print *, "Specify dimension N <= 0 to quit.  D is diagonal addition."
      print *, "Maximum problem dimension is ", Nx
      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, I, J)
        A1 = A0
!        print *, "A0(1, 8, 1, 6) ="	! OK for large N
!        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
!  Test Gauss and Rowperm with inverse calculation
!       print *, "A1:";  call prmf(A1, 1, N, 1, N)
!       print *, "B1:";  call prmf(B1, 1, N, 1, N)
        call Gauss(A1, B1, Pm, N)
!       print *, "After Gauss -- A1:";  call prmf(A1, 1, N, 1, N)
!       print *, "B1:";  call prmf(B1, 1, N, 1, N)
	call Rowperm(Xx, B1, Pm, N)
	print *, "Permutation:"
        do, ip = 1,N,8
          up = min(ip + 7, N)
          print 22, Pm(ip:up)
 22       FORMAT ( 8(1X, I4:) )
        end do
!       print *, "Xx:";  call prmf(Xx, 1, N, 1, N)
        C0 = matmul(A0, Xx)	      ! Should be Ident
        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(A0) = ", E11.4, "   max(Xx) = ", E11.4, &
                "   max err = ", E11.4)
!  Test GaussX by recalculating inverse
        B2 = B0
        call GaussX(A1, B2, Pm, N)
	print *, "GaussX finished:"
	err = maxval(abs(B1 - B2))
	print 24, err
 24	FORMAT (" Max abs diff from Gauss: ", E11.4)
! Test GaussI by repeating inverse calculation
	C0  = GaussI(A0, N)
	print *, "GaussI finished:"
	err = maxval(abs(C0 - Xx))
	print 25, err
 25	FORMAT (" Max abs diff from Gauss: ", E11.4)
      end do
 99   continue			! For exit
      end program Gausstest
