!
! Sparse Matrix 
!   Charles O'Neill July 3, 2008
!    www.caselab.okstate.edu
!
! Copyright (c) 2008 Charles O'Neill
! 
! Permission is hereby granted, free of charge, to any person
! obtaining a copy of this software and associated documentation
! files (the "Software"), to deal in the Software without
! restriction, including without limitation the rights to use,
! copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following
! conditions:
! 
! The above copyright notice and this permission notice shall be
! included in all copies or substantial portions of the Software.
! 
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
! OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
! HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
! WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
! OTHER DEALINGS IN THE SOFTWARE.

module Sparse
  implicit none

  integer,parameter :: DP = 8

  type Matrix
    integer :: CurrentMax
    real(DP),allocatable :: M(:)
    integer,allocatable  :: i(:),j(:)
  end type

contains

  subroutine testSparse()
    type (Matrix) :: Mtx
    call setSize(Mtx,10)
    
    write(*,*)
    write(*,*) "Store"    
    call store(Mtx,1,1,1.0d0)
    call store(Mtx,1,2,1.0d0)
    call store(Mtx,1,3,1.0d0)
    call store(Mtx,1,3,1.0d0)
    call store(Mtx,1,10,1.0d0)
    call store(Mtx,14,10654,1.0d0)
    call store(Mtx,1992,103,1.0d0)
    call store(Mtx,12,10,1.0d0)
    call store(Mtx,12,10,1.0d0)
    
    write(*,*)
    write(*,*) "Matrix Raw"
    write(*,*) Mtx%M
    write(*,*) Mtx%i
    write(*,*) Mtx%j
    
    write(*,*)
    write(*,*) "Retrieve"
    write(*,*) Mtx%CurrentMax
    write(*,*) retrieve(Mtx,1,1)
    write(*,*) retrieve(Mtx,1,2)
    write(*,*) retrieve(Mtx,1,3)
    write(*,*) retrieve(Mtx,1,4)
    write(*,*) retrieve(Mtx,1,5)
    
    
    
  end subroutine

  ! Sets the maximum number of storage elements
  subroutine setSize(This,MaxSize)
    type (Matrix) :: This
    integer :: MaxSize  
    allocate(This%M(MaxSize))
    allocate(This%i(MaxSize))
    allocate(This%j(MaxSize))
    This%CurrentMax = 0
  end subroutine

  ! Store Value in i,j location in Sparse Matrix
  subroutine store(This,i,j,Value)
    type (Matrix) :: This
    integer :: i,j
    real(DP) :: Value
    integer :: index
      
    ! Unique Location
    index = search(This,i,j)
    
    if(index==0)then
      ! New Unique location
      This%CurrentMax = This%CurrentMax + 1
      This%M(This%CurrentMax) = Value
      This%i(This%CurrentMax) = i 
      This%j(This%CurrentMax) = j
    else
      ! Old non-unique Location
      This%M(index) = This%M(index) + Value
    endif
      
  end subroutine
  
  ! Retrieve Value in i,j location in Sparse Matrix
  function retrieve(This,i,j)
    type (Matrix) :: This
    integer :: i,j,Index
    real(DP) :: retrieve
    
    ! Search for location of i,j in data array
    Index = search(This,i,j)
    
    if(index==0)then
      ! Location is not in sparse indexing 
      retrieve = 0.0d0
    else
      ! Location is in sparse indexing
      retrieve = This%M(Index)
    endif
    
  end function
  
  ! Search for location of i,j in Sparse Matrix
  function search(This,i,j)  ! SLOW SLOW SLOW SLOW AND STUPID
    type (Matrix) :: This
    integer :: search
    integer :: i,j
    integer :: index
    
    ! Search for index
    do index=1,This%CurrentMax
      if( (This%i(index)==i) .AND. (This%j(index)==j))then
        search = index
        return
      endif
    enddo
    
    ! Nothing Found
    search = 0
    
  end function
  
end module