! sparse.f90     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 matrix
!      type(cell), pointer, dimension(nrow) :: A
!      call init(A,nrow)
!      call put(A, i, j, val) ! j=nrow+1  is RHS
!      call write_all(A)

module sparse
  ! 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
  integer :: nrow
  
  public :: init
  public :: get
  public :: put
  public :: write_all
  public :: getrhs
  public :: putrhs
  public :: export
  public :: import
  public :: add
  public :: multiply
  public :: simeq
  private :: nrow

contains

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

    if(associated(A)) deallocate(A)
    allocate(A(n))
    nrow = size(A)
    allocate(cinit)
    cinit%j = -1
    cinit%val = 0.0d0
    nullify(cinit%next) 
    do i=1,nrow
      A(i) = cinit
    end do
  end subroutine init 

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

    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
    if(A(i)%j.eq.-1) then  ! first
      A(i)%j = j
      A(i)%val = val
      nullify(A(i)%next)
    else if(A(i)%j.eq.j) then  ! change first
      A(i)%val = val
    else if(j .lt. A(i)%j) then ! insert before first
      allocate(cinit) ! clone A(i)
      cinit%j = A(i)%j
      cinit%val = A(i)%val
      nullify(cinit%next)
      if(associated(A(i)%next)) cinit%next => A(i)%next
      A(i)%j = j
      A(i)%val = val
      A(i)%next => cinit
    else if(.not. associated(A(i)%next)) then ! add to first
      allocate(cinit)
      cinit%j = j
      cinit%val = val
      nullify(cinit%next)
      A(i)%next => cinit
    else ! find where to put
      prev => A(i)
      next => A(i)%next
      do while(associated(next))
        if(j.eq.next%j) then
          next%val = val ! just change value
          return
        else if(j .gt. next%j) then ! keep looking
           prev => next
           next => next%next
        else
          exit ! next is greater, insert here
        end if
      end do
      allocate(cinit)
      cinit%j = j
      cinit%val = val
      nullify(cinit%next)
      if(associated(next)) cinit%next => next ! may be null
      prev%next => cinit ! link previous to new cell
    end if
  end subroutine put

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

    nrow = size(A)
    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
    if(A(i)%j.eq.-1) then  ! first
      val = 0.0d0
      return
    else if(A(i)%j.eq.j) then  ! first
      val = A(i)%val
      return 
    else if(j .lt. A(i)%j) then ! before first
      val = 0.0d0
      return
    else if(.not. associated(A(i)%next)) then ! add to first
      val = 0.0d0
      return
    else ! find where to get
      next => A(i)%next
      do while(associated(next))
        if(j.eq.next%j) then
          val = next%val
          return
        else if(j .lt. next%j) then ! sparse zero
          val = 0.0d0
          return
        else
          prev => next
          next => next%next
        end if
      end do
    end if
  end function get

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

    nrow = size(A)
    do i=1,nrow
      j = A(i)%j
      if(j.eq.-1) cycle
      val = A(i)%val
      print "('(',I4,',',I4,')=',1PE15.6 )", i, j, val
      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(cell), dimension(:), pointer :: A
    double precision, dimension(:), intent(out) :: Y 
    integer :: i, j
    type(cell), pointer :: next

    nrow = size(A)
    if(size(Y).lt.nrow) then
      print *, 'sparse getrhs error, Y too small'
      return
    end if
    do i=1,nrow
      Y(i) = 0.0d0 ! some may not exists
    end do
    do i=1,nrow
      j = A(i)%j
      if(j .eq. -1) cycle
      if(j .eq.nrow+1) then
        Y(i) = A(i)%val
        cycle
      end if
      next => A(i)%next
      do while(associated(next))
        j = next%j
        if(j.eq.-1) exit
        if(j .eq.nrow+1) then
          Y(i) = next%val
          exit
        end if
        next => next%next
      end do ! while
    end do ! i
  end subroutine getrhs


  subroutine putrhs(A, Y)
    implicit none
    type(cell), dimension(:), pointer :: A
    double precision, dimension(:), intent(in) :: Y 
    integer :: i, j
    type(cell), pointer :: next

    nrow = size(A)
    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(cell), dimension(:), pointer :: A
    double precision, dimension(:,:), intent(out) :: M 
    integer :: i, j
    type(cell), pointer :: next
    double precision :: val

    nrow = size(A)
    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) = 0.0d0
      end do
    end do
    do i=1,nrow
      j = A(i)%j
      if(j.eq.-1) cycle
      val = A(i)%val
      M(i,j) = val
      next => A(i)%next
      do while(associated(next))
        j = next%j
        if(j.eq.-1) exit
        val = next%val
        M(i,j) = val
        next => next%next
      end do
    end do
  end subroutine export

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

    nrow = size(A)
    if(size(M).lt.nrow) then
      print *, 'sparse import error, M too small'
      return
    end if
    allocate(cinit)
    cinit%j = -1
    cinit%val = 0.0d0
    nullify(cinit%next) 
    do i=1,nrow ! start new
      A(i) = cinit
    end do
    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(cell), dimension(:), pointer :: A
    double precision, dimension(:), intent(in) :: Yin 
    double precision, dimension(:), intent(out) :: Yout 
    integer :: i, j

    nrow = size(A)
    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
      do j=1,nrow ! make more efficient
        Yout(i) = Yout(i) + get(A,i,j)*Yin(j)
      end do
    end do
  end subroutine multiply

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

    nrow = size(A)
    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
    implicit none
    type(cell), 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
    integer :: i, j, k, n, m

    n = size(A)
    if(size(X).lt.n) then
      print *, 'sparse simeq error, X too small'
      return
    end if
    ! quick and dirty version, needs to be rewritten
    m = n+1
    ! SET UP row  INTERCHANGE VECTORS 
    do k=1,n
      row(k) = k
    end do ! k

    ! BEGIN MAIN REDUCTION LOOP 
    do k=1,n

      ! FIND LARGEST ELEMENT FOR pivot 
      pivot = get(A,row(k),k) ! B(row(k),k)
      abs_pivot = abs(pivot)
      ipivot = k
      do i=k,n
        tmp = abs(get(A,row(i),k))
        if( tmp .gt. abs_pivot) then
          ipivot = i
          pivot = get(A,row(i),k)
          abs_pivot = tmp
        end if
      end do ! i

      ! HAVE pivot, INTERCHANGE row POINTERS 
      hold = row(k)
      row(k) = row(ipivot)
      row(ipivot) = hold

      ! CHECK FOR NEAR SINGULAR 
      if( abs_pivot < 1.0D-52 ) then
        do j=k+1,n+1
          call put(A,row(k),j,0.0D0)
        end do ! j
        print *, 'redundant row (singular) ', row(k)
      else

        ! REDUCE ABOUT pivot 
        do j=k+1,n+1
          ! B(row(k),j) = B(row(k),j) / B(row(k),k)
          call put(A,row(k),j,(get(A,row(k),j)/pivot))
        end do ! j

        ! INNER REDUCTION LOOP 
        do i=1,n
          if( i .ne. k) then
            tmp = get(A,row(i),k)
            do j=k+1,n+1
              call put(A,row(i),j,(get(A,row(i),j)-tmp*get(A,row(k),j)))
            end do ! j
          end if
        end do ! i
      end if 
      ! FINISHED INNER REDUCTION 
    end do ! k

    ! END OF MAIN REDUCTION LOOP 
    ! BUILD  X  FOR RETURN, UNSCRAMBLING rowS 
    do i=1,n
      X(i) = get(A,row(i),n+1)
    end do ! i
  end subroutine simeq 

  ! possible to simulate constructor  new cell(j, val, next)
  function make_cell(j, val, next) result(new_cell)
    ! usage:  new_cell_ptr => make_cell(j, val, next)
    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)
  end function make_cell        ! return pointer to constructed cell

end module sparse