!-*-f90-*-
!
!  API: Sparse matrices
!
!> \page "Comments on sparse matrix routines"
!> Please go to api/spmatrix.finc for the API documentation.

function fgsl_spmatrix_alloc(n1, n2)
  integer(fgsl_size_t), intent(in) :: n1, n2
  type(fgsl_spmatrix) :: fgsl_spmatrix_alloc
  fgsl_spmatrix_alloc%gsl_spmatrix = gsl_spmatrix_alloc(n1, n2)
end function fgsl_spmatrix_alloc
function fgsl_spmatrix_alloc_nzmax(n1, n2, nzmax, flags)
  integer(fgsl_size_t), intent(in) :: n1, n2, nzmax, flags
  type(fgsl_spmatrix) :: fgsl_spmatrix_alloc_nzmax
  fgsl_spmatrix_alloc_nzmax%gsl_spmatrix = &
  gsl_spmatrix_alloc_nzmax(n1, n2, nzmax, flags)
end function fgsl_spmatrix_alloc_nzmax
subroutine fgsl_spmatrix_size(m, n1, n2)
  type(fgsl_spmatrix), intent(in) :: m
  integer(fgsl_size_t), intent(inout) :: n1, n2
  call gsl_spmatrix_size(m%gsl_spmatrix, n1, n2)
end subroutine
subroutine fgsl_spmatrix_free(m)
  type(fgsl_spmatrix), intent(in) :: m
  call gsl_spmatrix_free(m%gsl_spmatrix)
end subroutine fgsl_spmatrix_free
function fgsl_spmatrix_realloc(nzmax, m)
  integer(fgsl_size_t), intent(in) :: nzmax
  type(fgsl_spmatrix), intent(inout) :: m
  integer(fgsl_int) :: fgsl_spmatrix_realloc
  fgsl_spmatrix_realloc = gsl_spmatrix_realloc(nzmax, m%gsl_spmatrix)
end function fgsl_spmatrix_realloc
function fgsl_spmatrix_set_zero(m)
  type(fgsl_spmatrix), intent(inout) :: m
  integer(fgsl_int) :: fgsl_spmatrix_set_zero
  fgsl_spmatrix_set_zero = gsl_spmatrix_set_zero(m%gsl_spmatrix)
end function fgsl_spmatrix_set_zero
function fgsl_spmatrix_nnz(m)
  type(fgsl_spmatrix), intent(in) :: m
  integer(fgsl_size_t) :: fgsl_spmatrix_nnz
  fgsl_spmatrix_nnz = gsl_spmatrix_nnz(m%gsl_spmatrix)
end function fgsl_spmatrix_nnz
function fgsl_spmatrix_memcpy(dest, src)
  type(fgsl_spmatrix), intent(inout) :: dest
  type(fgsl_spmatrix), intent(in) :: src
  integer(fgsl_int) :: fgsl_spmatrix_memcpy
  fgsl_spmatrix_memcpy = gsl_spmatrix_memcpy(dest%gsl_spmatrix, src%gsl_spmatrix)
end function fgsl_spmatrix_memcpy
function fgsl_spmatrix_get(m, i, j)
  type(fgsl_spmatrix), intent(in) :: m
  integer(fgsl_size_t), intent(in) :: i, j
  real(fgsl_double) :: fgsl_spmatrix_get
  fgsl_spmatrix_get = gsl_spmatrix_get(m%gsl_spmatrix, i ,j)
end function fgsl_spmatrix_get
function fgsl_spmatrix_set(m, i, j, x)
  type(fgsl_spmatrix), intent(in) :: m
  integer(fgsl_size_t), intent(in) :: i, j
  real(fgsl_double), intent(in) :: x
  integer(fgsl_int) :: fgsl_spmatrix_set
  fgsl_spmatrix_set = gsl_spmatrix_set(m%gsl_spmatrix, i ,j, x)
end function fgsl_spmatrix_set
function fgsl_spmatrix_compcol(T)
  type(fgsl_spmatrix), intent(in) :: T
  type(fgsl_spmatrix) :: fgsl_spmatrix_compcol
  fgsl_spmatrix_compcol%gsl_spmatrix = gsl_spmatrix_compcol(T%gsl_spmatrix)
end function fgsl_spmatrix_compcol
subroutine fgsl_spmatrix_cumsum(n, c)
  integer(fgsl_size_t), intent(in) :: n
  integer(fgsl_size_t), dimension(:), intent(inout), target, contiguous :: c
  if (size(c, dim=1) /= n+1) then
    call fgsl_error('c must have n+1 elements', __FILE__, __LINE__, fgsl_ebadlen)
    return
  endif
  call gsl_spmatrix_cumsum(n, c_loc(c))
end subroutine fgsl_spmatrix_cumsum
function fgsl_spmatrix_scale(m , x)
  type(fgsl_spmatrix), intent(inout) :: m
  real(fgsl_double), intent(in) :: x
  integer(fgsl_int) :: fgsl_spmatrix_scale
  fgsl_spmatrix_scale = gsl_spmatrix_scale(m%gsl_spmatrix, x)
end function fgsl_spmatrix_scale
function fgsl_spmatrix_scale_columns(a , x)
  type(fgsl_spmatrix), intent(inout) :: a
  type(fgsl_vector), intent(in) :: x
  integer(fgsl_int) :: fgsl_spmatrix_scale_columns
  fgsl_spmatrix_scale_columns = gsl_spmatrix_scale_columns(a%gsl_spmatrix, x%gsl_vector)
end function fgsl_spmatrix_scale_columns
function fgsl_spmatrix_scale_rows(a , x)
  type(fgsl_spmatrix), intent(inout) :: a
  type(fgsl_vector), intent(in) :: x
  integer(fgsl_int) :: fgsl_spmatrix_scale_rows
  fgsl_spmatrix_scale_rows = gsl_spmatrix_scale_rows(a%gsl_spmatrix, x%gsl_vector)
end function fgsl_spmatrix_scale_rows
function fgsl_spmatrix_minmax(m , min_out, max_out)
  type(fgsl_spmatrix), intent(in) :: m
  real(fgsl_double), intent(out) :: min_out, max_out
  integer(fgsl_int) :: fgsl_spmatrix_minmax
  fgsl_spmatrix_minmax = gsl_spmatrix_minmax(m%gsl_spmatrix, min_out, max_out)
end function fgsl_spmatrix_minmax
function fgsl_spmatrix_min_index(m , imin, jmin)
  type(fgsl_spmatrix), intent(in) :: m
  real(fgsl_double), intent(out) :: imin, jmin
  integer(fgsl_int) :: fgsl_spmatrix_min_index
  fgsl_spmatrix_min_index = gsl_spmatrix_min_index(m%gsl_spmatrix, imin, jmin)
end function fgsl_spmatrix_min_index
function fgsl_spmatrix_csc(dest, src)
  type(fgsl_spmatrix), intent(inout) :: dest
  type(fgsl_spmatrix), intent(in) :: src
  integer(fgsl_int) :: fgsl_spmatrix_csc
  fgsl_spmatrix_csc = gsl_spmatrix_csc(dest%gsl_spmatrix, src%gsl_spmatrix)
end function fgsl_spmatrix_csc
function fgsl_spmatrix_csr(dest, src)
  type(fgsl_spmatrix), intent(inout) :: dest
  type(fgsl_spmatrix), intent(in) :: src
  integer(fgsl_int) :: fgsl_spmatrix_csr
  fgsl_spmatrix_csr = gsl_spmatrix_csr(dest%gsl_spmatrix, src%gsl_spmatrix)
end function fgsl_spmatrix_csr
function fgsl_spmatrix_compress(src, sptype)
  type(fgsl_spmatrix), intent(in) :: src
  integer(fgsl_int), intent(in) :: sptype
  type(fgsl_spmatrix) :: fgsl_spmatrix_compress
  fgsl_spmatrix_compress%gsl_spmatrix = gsl_spmatrix_compress(src%gsl_spmatrix, sptype)
end function fgsl_spmatrix_compress
function fgsl_spmatrix_add(c, a, b)
  type(fgsl_spmatrix), intent(inout) :: c
  type(fgsl_spmatrix), intent(in) :: a, b
  integer(fgsl_int) :: fgsl_spmatrix_add
  fgsl_spmatrix_add = gsl_spmatrix_add(c%gsl_spmatrix, &
  a%gsl_spmatrix, b%gsl_spmatrix)
end function fgsl_spmatrix_add
function fgsl_spmatrix_add_to_dense(a, b)
  type(fgsl_matrix), intent(inout) :: a
  type(fgsl_spmatrix), intent(in) ::  b
  integer(fgsl_int) :: fgsl_spmatrix_add_to_dense
  fgsl_spmatrix_add_to_dense = gsl_spmatrix_add_to_dense(a%gsl_matrix, b%gsl_spmatrix)
end function fgsl_spmatrix_add_to_dense
function fgsl_spmatrix_d2sp(S, A)
  type(fgsl_spmatrix), intent(inout) :: S
  type(fgsl_matrix), intent(in) :: A
  integer(fgsl_int) :: fgsl_spmatrix_d2sp
  fgsl_spmatrix_d2sp = gsl_spmatrix_d2sp(S%gsl_spmatrix, A%gsl_matrix)
end function fgsl_spmatrix_d2sp
function fgsl_spmatrix_sp2d(A, S)
  type(fgsl_matrix), intent(inout) :: A
  type(fgsl_spmatrix), intent(in) :: S
  integer(fgsl_int) :: fgsl_spmatrix_sp2d
  fgsl_spmatrix_sp2d = gsl_spmatrix_sp2d(A%gsl_matrix, S%gsl_spmatrix)
end function fgsl_spmatrix_sp2d
function fgsl_spmatrix_equal(a, b)
  type(fgsl_spmatrix), intent(in) :: a, b
  integer(fgsl_int) :: fgsl_spmatrix_equal
  fgsl_spmatrix_equal = gsl_spmatrix_equal(a%gsl_spmatrix, b%gsl_spmatrix)
end function fgsl_spmatrix_equal
function fgsl_spmatrix_transpose_memcpy(dest, src)
  type(fgsl_spmatrix), intent(inout) :: dest
  type(fgsl_spmatrix), intent(in) :: src
  integer(fgsl_int) :: fgsl_spmatrix_transpose_memcpy
  fgsl_spmatrix_transpose_memcpy = gsl_spmatrix_transpose_memcpy(dest%gsl_spmatrix, src%gsl_spmatrix)
end function fgsl_spmatrix_transpose_memcpy
function fgsl_spmatrix_transpose(m)
  type(fgsl_spmatrix), intent(inout) :: m
  integer(fgsl_int) :: fgsl_spmatrix_transpose
  fgsl_spmatrix_transpose = gsl_spmatrix_transpose(m%gsl_spmatrix)
end function fgsl_spmatrix_transpose
integer(fgsl_int) function fgsl_spblas_dgemv(transa, alpha, a, x, beta, y)
  integer(fgsl_int), intent(in) :: transa
  real(fgsl_double), intent(in) :: alpha, beta
  type(fgsl_spmatrix), intent(in) :: a
  type(fgsl_vector), intent(in) :: x
  type(fgsl_vector), intent(inout) :: y
  fgsl_spblas_dgemv = gsl_spblas_dgemv(transa, alpha, a%gsl_spmatrix, x%gsl_vector, &
       beta, y%gsl_vector)
end function fgsl_spblas_dgemv
integer(fgsl_int) function fgsl_spblas_dgemm(alpha, a, b, c)
  real(fgsl_double), intent(in) :: alpha
  type(fgsl_spmatrix), intent(in) :: a, b
  type(fgsl_spmatrix), intent(inout) :: c
  fgsl_spblas_dgemm = gsl_spblas_dgemm(alpha, a%gsl_spmatrix, b%gsl_spmatrix, &
       c%gsl_spmatrix)
end function fgsl_spblas_dgemm
!
! I/O
!
integer(fgsl_int) function fgsl_spmatrix_fwrite(stream, m)
  type(fgsl_file) :: stream
  type(fgsl_spmatrix), intent(in) :: m
  fgsl_spmatrix_fwrite = gsl_spmatrix_fwrite(stream%gsl_file, m%gsl_spmatrix)
end function fgsl_spmatrix_fwrite
integer(fgsl_int) function fgsl_spmatrix_fread(stream, m)
  type(fgsl_file) :: stream
  type(fgsl_spmatrix), intent(inout) :: m
  fgsl_spmatrix_fread = gsl_spmatrix_fread(stream%gsl_file, m%gsl_spmatrix)
end function fgsl_spmatrix_fread
integer(fgsl_int) function fgsl_spmatrix_fprintf(stream, m, format)
  type(fgsl_file) :: stream
  type(fgsl_spmatrix), intent(in) :: m
  character(kind=fgsl_char, len=*), intent(in) :: format
  fgsl_spmatrix_fprintf = gsl_spmatrix_fprintf(stream%gsl_file, m%gsl_spmatrix, format)
end function fgsl_spmatrix_fprintf
type(fgsl_spmatrix) function fgsl_spmatrix_fscanf(stream)
  type(fgsl_file) :: stream
  fgsl_spmatrix_fscanf%gsl_spmatrix = gsl_spmatrix_fscanf(stream%gsl_file)
end function fgsl_spmatrix_fscanf
!
! helper function
!
subroutine fgsl_spmatrix_getfields(m, i, p, d)
  type(fgsl_spmatrix), intent(in) :: m
  integer(fgsl_int), pointer, intent(inout) :: i(:), p(:)
  real(fgsl_double), pointer, intent(inout) :: d(:)

  type(c_ptr) :: ip, pp, dp
  integer(fgsl_size_t) :: nz, psize

  call gsl_aux_spmatrix_getfields(m%gsl_spmatrix, ip, dp, pp, psize)
  nz = gsl_spmatrix_nnz(m%gsl_spmatrix) 
  call c_f_pointer(ip, i, [nz])
  call c_f_pointer(dp, d, [nz])
  call c_f_pointer(pp, p, [psize])
end subroutine fgsl_spmatrix_getfields