! sparse1.f90    develop  Sparse storage of a matrix
!                specifically designed for solving PDE
!                simultaneous equations.
!                one (1) based subscripting
!                column (nrow+1) is right hand side
!
!      use sparse    compile module  sparse.f90
!      integer :: i, j
!      double precision :: val
!      integer :: nrow = number_of_rows_in_sparse_matrix
!      type(cellptr), dimension(:), pointer :: A
!      call init(A,nrow)
!      call put(A, i, j, val)
!      call write_all(A)

module global
  ! 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
  
  public :: init
  public :: put
  public :: write_all
  private :: nrow

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 

  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

    nrow = size(A)
    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

  ! needed to simulate constructor  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

  ! 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
end module global

program sparse1
  use global
  implicit none
  integer :: i, j
  double precision :: val
  integer :: n = 99
  type(cellptr), dimension(:), pointer :: A;

  print *, 'sparse1 running'
  call init(A,n)
  i = 1
  j = 1
  val = 1.0 
  call put(A, i, j, val)
  call put(A, 1, 3, 3.0d0)
  call put(A, 1, 5, 5.0d0)
  call put(A, 2, 2, 2.0d0)
  call put(A, 2, 4, 4.0d0)
  call put(A, 2, 99, 99.0d0) ! 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,1.1'
  call put(A, 2, 1, 1.1d0)
  call write_all(A)

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