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

module sparse1
  ! define the type "cell" to be a element in sparse array
  type cell
    integer :: j = -1                ! column of sparse matrix
    double precision :: val = 0.0d0  ! value of entry
    type (cell), pointer :: next     ! link to next column entry
  end type cell
  type cellptr
    type (cell), pointer :: next
  end type cellptr
  integer :: nrow                    ! set by init

  public :: init
  public :: put
  public :: get
  public :: write_all
  public :: getrhs
  public :: putrhs
  public :: export
  public :: import
  public :: add
  public :: multiply
  public :: simeq
  private :: nrow
  private :: find
  private :: reduce

contains

  subroutine init(A,n)
    implicit none
    type(cellptr), dimension(:), pointer:: A 
    integer, intent(in) :: n
    integer :: i
    type(cell), pointer :: cinit

    if(associated(A)) deallocate(A)
    allocate(A(n))
    nrow = size(A)
    do i=1,nrow
      nullify(A(i)%next)
    end do
  end subroutine init 

  ! needed to simulate constructor  new cell(j, val, next)
  function make_cell(j, val, next) result(new_cell)
    implicit none
    integer, intent(in) :: j
    double precision, intent(in) :: val
    type(cell), pointer :: next ! no intent(in)
    type(cell), pointer :: new_cell ! returns pointer

    allocate(new_cell)
    new_cell%j = j 
    new_cell%val = val
    nullify(new_cell%next)
    if(associated(next)) new_cell%next => next
  end function make_cell        ! return pointer to constructed cell

  function find(next, j) result(found)
    implicit none
    integer, intent(in) :: j
    type (cell), pointer :: next ! no intent(in)
    type (cell), pointer :: found ! returns pointer
    type (cell), pointer :: p

    nullify(found)
    p => next
    do while(associated(p))
      if(p%j .eq. j) then
        found => p
        return
      end if
      if(p%j .gt. j) return
      p => p%next
    end do
  end function find        ! return pointer to found j cell

  subroutine put(A, i, j, val)
    implicit none
    type(cellptr), dimension(:), pointer:: A 
    integer, intent(in) :: i, j
    double precision, intent(in) :: val
    type(cell), pointer :: cinit, prev, next

    if(i.lt.1 .or. i.gt.nrow .or. j.lt.1 .or. j.gt.nrow+1) then
      print *, 'sparse put error i=',i,', j=',j
      return
    end if
    nullify(next)
    if(associated(A(i)%next)) next => A(i)%next
    if(.not. associated(next)) then ! put first
      A(i)%next => make_cell(j,val,next)
    else ! find where to put
      prev => A(i)%next
      do while(associated(next))
        if(j.eq.next%j) then
          next%val = val ! just change value
          return
        else if(j .lt. next%j) then ! insert before
          cinit => make_cell(next%j,next%val,next%next) ! clone
          next%j = j
          next%val = val
          next%next => cinit
          return
        else if(j .gt. next%j) then ! keep looking
           prev => next
           next => next%next
        end if
      end do
      nullify(next)
      prev%next => make_cell(j,val,next) ! link at end
    end if
  end subroutine put

  function get(A, i, j) result(val)
    implicit none
    type(cellptr), dimension(:), pointer:: A 
    integer, intent(in) :: i, j
    double precision :: val
    type(cell), pointer :: cinit, prev, next

    if(i.lt.1 .or. i.gt.nrow .or. j.lt.1 .or. j.gt.nrow+1) then
      print *, 'sparse get error i=',i,', j=',j
      return
    end if
    nullify(next)
    if(associated(A(i)%next)) next => A(i)%next
    val = 0.0d0
    if(.not. associated(next)) return ! return zero
    next => find(next, j)
    if(associated(next)) val = next%val
  end function get

  ! subroutine to write_all sparse matrix
  subroutine write_all(A)
    implicit none
    type (cellptr), dimension(:), pointer :: A
    integer :: i, j
    double precision :: val
    type (cell), pointer :: next

    nrow = size(A)
    do i=1,nrow
      next => A(i)%next
      do while(associated(next))
        j = next%j
        if(j.eq.-1) exit
        val = next%val
        print "('(',I4,',',I4,')=',1PE15.6 )", i, j, val
        next => next%next
      end do
    end do
  end subroutine write_all

  subroutine getrhs(A, Y)
    implicit none
    type(cellptr), dimension(:), pointer :: A
    double precision, dimension(:), intent(out) :: Y 
    integer :: i

    if(size(Y).lt.nrow) then
      print *, 'sparse getrhs error Y too small'
      return
    end if
    do i=1,nrow
      Y(i) = get(A,i,nrow+1)
    end do
  end subroutine getrhs

  subroutine putrhs(A, Y)
    implicit none
    type(cellptr), dimension(:), pointer :: A
    double precision, dimension(:), intent(out) :: Y 
    integer :: i

    if(size(Y).lt.nrow) then
      print *, 'sparse putrhs error Y too small'
      return
    end if
    do i=1,nrow
      call put(A,i,nrow+1,Y(i))
    end do
  end subroutine putrhs

  subroutine export(A, M)
    implicit none
    type(cellptr), dimension(:), pointer :: A
    double precision, dimension(:,:), intent(out) :: M 
    integer :: i, j

    if(size(M).lt.nrow) then
      print *, 'sparse export error M too small'
      return
    end if
    do i=1,nrow
      do j=1,nrow+1
        M(i,j) = get(A,i,j)
      end do ! j
    end do ! i
  end subroutine export

  subroutine import(A, M)
    implicit none
    type(cellptr), dimension(:), pointer :: A
    double precision, dimension(:,:), intent(in) :: M 
    double precision :: ztol = 1.0e-20 ! could be a passed argument
    integer :: i, j

    if(size(M).lt.nrow) then
      print *, 'sparse import error M too small'
      return
    end if
    call init(A,nrow)
    do i=1,nrow
      do j=1,nrow+1
        if(abs(M(i,j)).gt.ztol) call put(A,i,j,M(i,j))
      end do ! j
    end do ! i
  end subroutine import

  subroutine multiply(A, Yin, Yout)
    implicit none
    type(cellptr), dimension(:), pointer :: A
    double precision, dimension(:), intent(in) :: Yin 
    double precision, dimension(:), intent(out) :: Yout 
    integer :: i, j
    type (cell), pointer :: next

    if(size(Yin).lt.nrow .or. size(Yout).lt.nrow) then
      print *, 'sparse multiply error Yin or Yout too small'
      return
    end if
    do i=1,nrow
      Yout(i) = 0.0d0
      next => A(i)%next
      do while(associated(next))
        if(next%j.le.nrow) Yout(i) = Yout(i) + next%val*Yin(next%j)
        next => next%next
      end do
    end do
  end subroutine multiply

  subroutine add(A, i, j, val)
    implicit none
    type(cellptr), dimension(:), pointer:: A 
    integer, intent(in) :: i, j
    double precision, intent(in) :: val

    if(i.lt.1 .or. i.gt.nrow .or. j.lt.1 .or. j.gt.nrow+1) then
      print *, 'sparse add error i=',i,', j=',j
      return
    end if
    call put(A,i,j,val+get(A,i,j)) ! make more efficient
  end subroutine add

  subroutine simeq(A, X) ! solve A * X = Y RHS
    !      PURPOSE : SOLVE THE LINEAR SYSTEM OF EQUATIONS WITH REAL
    !                COEFFICIENTS   (A) * |X| = |Y|
    !
    !      INPUT  : THE NUMBER OF EQUATIONS  n
    !               THE REAL MATRIX  A   should be A(i)(j) but A(i,j) 
    !               THE REAL VECTOR  Y as RHS of a
    !      OUTPUT : THE REAL VECTOR  X
    !
    !      METHOD : GAUSS-JORDAN ELIMINATION USING MAXIMUM ELEMENT
    !               FOR PIVOT.
    !
    !      USAGE  :     call simeq(A,X)
    !
    !    WRITTEN BY : JON SQUIRE , 28 MAY 1983
    !    ORIGONAL DEC 1959 for IBM 650, TRANSLATED TO OTHER LANGUAGES
    !    e.g. FORTRAN converted to Ada converted to C to f90 to sparse
    implicit none
    type(cellptr), dimension(:), pointer :: A
    double precision, dimension(:), intent(out) :: X 
    integer, dimension(size(A)) :: row             ! row INTERCHANGE INDICIES 
    integer :: hold , ipivot                 ! pivot INDICIES 
    double precision :: pivot                ! pivot ELEMENT VALUE 
    double precision :: abs_pivot, tmp, kval, kjval
    integer :: i, j, k
    type (cell), pointer :: next, nexti, nextk

    if(size(X).lt.nrow) then
      print *, 'sparse simeq error X too small'
      return
    end if
    ! set up row interchange vectors
    do k=1,nrow
      row(k) = k
    end do
    ! begin main reduction loop
    do k=1,nrow
      ! find largest element for pivot
      pivot = get(A,row(k),k)
      abs_pivot = abs(pivot)
      ipivot = k
      do i=k+1,nrow
        tmp = get(A,row(i),k)
        if(abs(tmp).gt.abs_pivot) then
          ipivot = i
          pivot = tmp
          abs_pivot = abs(pivot)
        end if
      end do
      ! have pivot, interchange row indicies
      hold = row(k)
      row(k) = row(ipivot)
      row(ipivot) = hold
      ! check for near singular
      if(abs_pivot.lt.1.0D-10) then
        print *, 'singular at k=',k,' redundant row ',row(k)
      else
        ! divide by pivot
        print *, 'row(k)=',row(k),', pivot=',pivot
        next => A(row(k))%next
        do while(associated(next))
          if(next%j.gt.k) then
            next%val = next%val/pivot
            print *, 'simeq divide k=',k,', row(k)=',row(k),', j=',next%j
          end if
          next => next%next
        end do 
        do i=1,nrow ! reduce all except row(k)
          if(i.ne.k) then
            nexti => A(row(i))%next
            do while(associated(nexti))
              if(nexti%j.gt.k) then 
                kval = 0.0d0
                exit
              else if(nexti%j.eq.k) then
                kval = nexti%val
                nexti = nexti%next
                exit
              else
                nexti = nexti%next
              end if
            end do
            print *, 'kval=',kval
            if(kval.eq.0.0d0) cycle
            nextk => A(row(k))%next
            do while(associated(nexti))
              do while(associated(nextk))
                kjval = 0.0d0
                if(nextk%j.eq.nexti%j) then
                  kjval = nextk%val
                  nextk => nextk%next
                  exit
                else if(nextk%j.gt.nexti%j) then
                  kjval = 0.0d0
                  exit
                else
                  nextk => nextk%next
                end if
              end do
              nexti%val = nexti%val - kval*kjval
              print *, 'simeq reduce k=',k,', j=',nexti%j,' kjval=',kjval 
              nexti => nexti%next  
            end do ! while
          end if
        end do ! i
      end if ! singular
    end do ! k
    ! build X for return
    do i=1,nrow
      X(i) = get(A,row(i),nrow+1)
    end do
  end subroutine simeq 

end module sparse1
 
program test_sparse1
  use sparse1
  implicit none
  integer :: i, j
  double precision :: val
  integer :: n = 4
  type(cellptr), 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_sparse1.f90 running'
  call init(A,n)
  print *, 'expect error'
  val = get(A,0,1)
  print *, 'expect error'
  val = get(A,1,0)
  print *, 'expect error'
  val = get(A,5,1)
  print *, 'expect error'
  val = get(A,1,6)

  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_sparse1.f90 ending'
end program test_sparse1