! test_sparse.f90   Test Sparse.f90
!
!      use sparse    compile module  sparse.f90
!      integer :: i, j
!      double precision :: val
!      integer, parameter :: nrow = 2
!      type(cell), pointer, dimension(nrow) :: A
!      call init(A,nrow)
!      call put(A, i, j, val)
!      call write_all(A)

program test_sparse
  use sparse
  implicit none
  integer :: i, j
  double precision :: val
  integer :: n = 4
  type(cell), dimension(:), pointer :: A

  ! for testing
  double precision, dimension(4) :: Y
  double precision, dimension(4) :: Yout
  double precision, dimension(4,5) :: M
  double precision, dimension(4) :: X

  print *, 'test_sparse.f90 running'
  call init(A,n)
  i = 1
  j = 1
  val = 1.0 
  call put(A, i, j, val)
  call put(A, 1, 3, 1.3d0)
  call put(A, 1, 5, 1.5d0)
  call put(A, 2, 2, 2.0d0)
  call put(A, 2, 4, 4.0d0)
  call put(A, 2, 5, 2.5d0) ! RHS
  call write_all(A)
  print *, 'change 1,3 to 3.3'
  call put(A, 1, 3, 3.3d0)
  call write_all(A)
  print *, 'change 1,1 to 1.1'
  call put(A, 1, 1, 1.1d0)
  call write_all(A)
  print *, 'insert 1,2,2.2  between'
  call put(A, 1, 2, 2.2d0)
  call write_all(A)
  print *, 'insert at head 2,1,2.1'
  call put(A, 2, 1, 2.1d0)
  call write_all(A)
  val = get(A,1,1)
  print *, 'get A(1,1) = ', val
  val = get(A,2,5)
  print *, 'get A(2,5) = ', val
  print *, 'fill in RHS'
  call put(A, 3, 3, 1.0d0)
  call put(A, 3, 5, 3.5d0)
  call put(A, 4, 5, 4.5d0)
  call put(A, 4, 4, 1.0d0)
  call write_all(A)
  call getrhs(A, Y)
  print *, 'get RHS  using getrhs'
  print "('RHS=',10F7.3)",(Y(i), i=1,n)
  call export(A, M)
  print *, 'export matrix to non sparse form'
  do i=1,n
    print "('row',I4,'=',10F7.3)",i, (M(i,j), j=1,n+1)
  end do

  do i=1,n
    Y(i) = -Y(i)
  end do
  call putrhs(A,Y)
  print *, 'check putrhs'
  call write_all(A)

  do i=1,n
    do j=1,n+1
      M(i,j) = -M(i,j)
    end do
  end do
  call import(A, M)
  print *, 'check import'
  call write_all(A)
  print *, 'check multiply'
  call multiply(A,Y,Yout)
  print "('Yout=',10F7.3)",(Yout(i), i=1,n)
  call add(A,1,1,98.0d0)
  print *, 'check add(A,1,1,98.0) ', get(A,1,1)
  print *, 'solve simultaneous equations'
  call export(A, M) ! save it
  call getrhs(A, Y) ! save it
  call simeq(A, X) ! A * X = Y RHS
  print "('simeq=',10F7.3)",(X(i), i=1,n)
  call import(A, M) ! must restore
  call multiply(A, X, Yout)
  print *, 'check simeq A * X = Yout against Y'
  do i=1,n
    print *, 'check Y,Yout,err',Y(i), Yout(i), abs(Y(i)-Yout(i))
  end do

  print "('simeq =',10F7.3)",(Y(i), i=1,n)


  print *,'test_sparse.f90 ending'
end program test_sparse